Mechanism-based prediction of drug synergy via network controllability analysis of therapeutic pathways in intractable diseases.

Introduction

Intractable disease conditions, such as cancers, are typically caused by the dysregulation of multiple biomolecules (e.g., proteins) and are difficult to treat with single-drug therapies that can regulate only a limited subset of targets. These diseases are commonly multifactorial, involving complex disruptions across multiple pathways 1 ; thus, single-drug therapies are often insufficient to achieve durable responses. To overcome this limitation, multidrug therapies (e.g., doublet or triplet combinations) that exhibit synergistic effects have emerged as a promising strategy 2 , 3 for simultaneously targeting multiple disease mechanisms. 4 , 5 For example, surgery is not suitable for some patients with colorectal cancer (CRC) because of the tumor location or the risk of colon perforation and intestinal obstruction. 6 In such cases, multidrug therapy is a potential alternative. 7 However, the identification of effective drug combinations remains challenging owing to the vast space of possible drug combinations and the absence of reliable strategies to identify combinations with clinical potential. As an example, there are 499,500 possible pair combinations of 1000 drugs. Thus, the experimental evaluation of the synergistic effects of all possible drug combinations across thousands of disease conditions is infeasible.

To address this bottleneck, computational approaches have been developed to prioritize the discovery of synergistic drug combinations. Existing methods include supervised learning using known drug pairs, 8 , 9 , 10 , 11 and unsupervised or network-based models that do not rely on prior knowledge. 12 , 13 The supervised learning approach has become mainstream in this field, and it is effective for identifying drug combinations in which the chemical structures and mechanisms of action are similar to those of known drug combinations. However, as the supervised approach depends heavily on prior knowledge of such combinations; therefore, it is unlikely to discover combinations in which the chemical structures and mechanisms of action differ from those of known drug combinations. In contrast, the unsupervised learning approach, which does not rely on prior information about existing combinations, is more desirable in clinical practice. The use of transcriptome data to predict drugs that can be combined with first-line treatments for a specific disease has been proposed. 12 However, the prediction is constrained because one of the drugs in the pair is fixed. The use of topological relationships between drug target modules in biomolecular networks was also proposed to predict drug combinations. 13 However, this method is not applicable to diseases with few known therapeutic target molecules and often focuses only on the largest network modules, limiting the ability to address complex disease conditions involving diverse molecular mechanisms. There is still substantial room for improving the accuracy and applicability of drug combination prediction, particularly for complex disease conditions. Moreover, previous approaches have rarely distinguished between mechanistically distinct types of synergy.

To address this unmet need, we introduce the mechanistic classification of drug synergy into two paradigms—“boost” and “complement”—to distinguish synergy modes based on their molecular targeting patterns. Boost synergy describes the amplified therapeutic effect of multiple drugs acting on the same therapeutic target molecules. There is a bias toward boost-type combination in clinical practice; most drug combinations currently used involve drugs approved for the same disease category. For example, combinations of antineoplastic drugs are commonly used to treat patients with cancer. However, boost-type drug combinations result in clinical challenges such as the risk of treatment discontinuation due to cumulative side effects and unexpected drug-drug interactions. 14 Conversely, in complement synergy, multiple drugs are used to regulate different therapeutic target molecules and enhance therapeutic effects in a complementary manner. Complement-type drug combinations for cancers achieve high therapeutic efficacy with minimal antineoplastic doses by leveraging complementary interactions with non-antineoplastic drugs. However, the mechanisms underlying complement-type synergy are complex and often unpredictable and the lack of intrinsic anticancer effects in non-antineoplastic drugs makes achieving synergy particularly difficult. Thus, the design of complement-type synergy regimens poses significant scientific and technical challenges; indeed, their development remains largely reliant on trial-and-error approaches. However, their potential to reduce adverse effects and expand treatment options—particularly for elderly patients and those with comorbidities—is central to their clinical importance. In the context of precision medicine, demand is increasing for predictive models that account for disease-specific regulatory complexity. These challenges highlight an unmet need for predictive models with the ability to mechanistically distinguish and identify clinically actionable drug synergies.

In this research, we developed a mechanism-aware computational framework that categorizes drug synergy into boost-type and complement-type paradigms; here, we present a network-based method called SYNERGIE (drug SYNERGy inference engine) to predict context-specific and clinically useful drug combinations for a variety of disease conditions. The algorithm incorporates network controllability analysis, which models how disease-specific perturbations propagate through molecular systems to estimate potential therapeutic pathways and target molecules, extending its applicability to intractable disease conditions without clearly characterized therapeutic target molecules. We validated the identified complement-type drug triplets for CRC using in vitro experiments and cell type-resolved multiomics analyses of human cohort data at bulk and single-cell levels, demonstrating its utility for mechanism-driven drug combination therapy design. This dual-mechanism framework is distinguished from conventional predictors by accounting for mechanistic diversity and enabling its application, even in clinically challenging contexts. Notably, the predicted triplets enhanced first-line therapies, suggesting that SYNERGIE can enhance existing treatment regimens to support patient-contextualized precision medicine and increase the available therapeutic options for patients with limited alternatives.

Results

Overview of SYNERGIE for predicting synergy-aware drug combinations

SYNERGIE enables the systematic prediction of synergy-aware drug combinations across diverse diseases, including those without known therapeutic target molecules. Although the overlap between disease therapeutic target molecules and drug target molecules is critical for the prediction of effective drug combinations, many diseases—particularly intractable ones—lack known therapeutic target molecules. Moreover, as most known therapeutic target molecules are associated with first-line therapies, there is a need to identify alternative target molecules for the development of effective combination strategies. To address these challenges, SYNERGIE estimates therapeutic pathways and potential target molecules to facilitate the prediction of synergistic drug combinations.

To identify upstream therapeutic pathways regulating disease-specific expression changes, we first extracted differentially expressed gene (DEG)-enriched pathways (DEGPs) from disease transcriptome data using weighted gene co-expression network analysis (WGCNA) 15 followed by subsequent pathway enrichment ( Figure 1 A). The DEGPs were mapped onto a directed pathway-pathway association network to model inter-pathway regulatory relationships. The upstream pathways predicted to exert strong regulatory control over downstream DEGPs and mediate therapeutic effects through drug-induced modulation (e.g., inhibition or activation) were defined as therapeutic pathways based on network topology (see STAR Methods ). Figure 1 Overview of SYNERGIE: Prediction of synergy-aware drug combinations across diverse diseases by inferring therapeutic pathways and target molecules using network controllability analysis (A) Therapeutic pathway estimation from patient-specific transcriptome data. Patient-specific transcriptomic data were used to identify differentially expressed genes (DEGs) and weighted gene co-expression network analysis (WGCNA), and enrichment analysis were used to define DEG-enriched pathways (DEGPs). These DEGPs were mapped onto a directed pathway-pathway association network to infer upstream pathways likely to regulate disease-specific expression programs. The upstream pathways that were expected to strongly regulate the downstream DEGPs were subsequently identified as “therapeutic pathways” based on the topology of the directed pathway network. (B) Therapeutic target molecule estimation using network controllability analysis. Network controllability analysis was applied to each therapeutic pathway within a protein-protein functional interaction (PPI) network to model how perturbations propagate to estimate essential nodes exerting control over disease-relevant signaling. (C) Disease-specific drug combination prediction via the Bayesian optimization-based integration of drug-disease and drug-drug interactions across multiomics layers, which incorporates both network-based proximity and perturbation-based expression effects. Four factors are integrated to predict disease-specific synergistic drug combinations (doublets): drug-disease target commonality, drug-disease DEG inversion, drug-drug target random walk, and drug-drug DEG random walk. Three different types of synergistic drug combinations are predicted separately based on the overlapping target molecules of two drugs: (1) boost type, (2) complement type, and (3) indistinct type. (D) Validation of the predicted drug doublet and triplet combinations both in vitro and in human samples. Synergistic effects were validated in vitro using cell survival assays, and human-level validation of the cell type-specific mechanisms of the synergistic effects was achieved using both bulk and single-cell RNA-seq data. Created with BioRender.com.

To identify therapeutic target molecules for each disease, we performed the network controllability analysis for each therapeutic pathway. This analysis modeled how disease-specific perturbations propagate through molecular systems and exert control over biological signal flow using a protein-protein functional interaction (PPI) network. From this analysis, we then identified essential genes—indicated by dominating nodes capable of exerting dynamic control over the therapeutic pathways—and considered them as therapeutic target molecules ( Figure 1 B; see STAR Methods ).

Then, we conceptually redefined drug combination strategies into two different types of synergy paradigms based on the therapeutic target molecule relationships between drugs ( Figure 1 C): (1) “boost type,” in which the drugs regulate the same therapeutic target molecules to amplify therapeutic efficacy, and (2) “complement type,” in which the drugs regulate different therapeutic target molecules to enhance efficacy in a complementary manner. We distinctively predicted the associated drug combinations by integrating drug-disease and drug-drug interactions across multiomics layers of a PPI network using a Bayesian optimization framework, to reflect biological response flow (see STAR Methods ).

Finally, we achieved the comprehensive prediction of drug doublets and triplets for CRC. By building upon existing treatment regimens to ensure clinical feasibility, we constructed a triplet consisting of a first-line drug and a predicted drug doublet, and evaluated its therapeutic relevance using in vitro experiments and human-level validation from disease cohort data at both the bulk and single-cell level ( Figure 1 D).

SYNERGIE reveals upstream therapeutic pathways that govern disease-specific gene expression programs across a range of diseases

Using SYNERGIE, we identified therapeutic pathways that regulate disease-specific DEGPs in 14 diseases by analyzing the pathway-pathway association network (see STAR Methods and Table S1 ). The diseases studied were acute myeloid leukemia, breast cancer (BC), chronic myeloid leukemia (CML), CRC, melanoma, diffuse large B cell lymphoma (DLBCL), glioblastoma, multiple myeloma, non-small cell lung cancer, ovarian cancer, prostate cancer, renal cell carcinoma (RCC), T cell acute lymphoblastic leukemia, and type 2 diabetes mellitus. Figure 2 A presents the network modules containing both DEGPs and therapeutic pathways for CRC. The DEGPs included pathways such as “post-translational phosphorylation” and “MET signaling,” whereas the therapeutic pathways included “deubiquitination” and “PIP3-AKT signaling.” Using independent resources, we explored these therapeutic pathways. For example, the deubiquitination of lymphoid-specific helicase has been shown to suppress erastin-induced ferroptosis in CRC cells. 16 Additionally, PIP3-AKT signaling is a promising therapeutic target as it is crucial for CRC cell proliferation and migration. 17 Other pathways, such as Notch1 signaling, Toll-like signaling, and vitamin metabolism, are also strongly associated with CRC. 18 , 19 Together, these findings suggest that targeting the identified therapeutic pathways could offer a mechanism-based strategy for modulating CRC-specific gene expression programs. Figure 2 Estimation of therapeutic pathways and target molecules using network controllability analysis in conjunction with disease transcriptome data (A) The pathway-pathway association network in colorectal cancer (CRC) showing associations between therapeutic pathways and DEG-enriched pathways (DEGPs): pink nodes represent therapeutic pathways; blue nodes represent DEGPs; edges denote pathway-pathway associations. All therapeutic pathways and DEGPs for CRC and other diseases are presented in Table S1 and Figure S1 . (B) Chord diagram of the regulatory relationships between therapeutic pathways and DEGPs in CRC. Directed links represent pathway-level regulatory influence from therapeutic pathways to DEGPs, with link width scaled by random walk with restart (RWR) values. Each sector and color represents different pathway categories. The scale indicates the number of pathways for each pathway category. Results for other diseases are provided in Figure S2 . (C) Venn diagram comparing the causal genes, known target molecules, and estimated target molecules. The numbers represent the gene counts in each section. Source data are provided in Table S2 . Results for other diseases are presented in Figure S4 . (D) Partial visualization of the therapeutic pathways in the PPI network: pink circles denote the estimated target molecules; gray circles represent nonestimated target molecules; gray lines indicate PPIs. (E–G) Zoomed-in views of the Axon guidance pathway (E), GPCR downstream pathway (F), and MAPK1/MAPK3 signaling pathway (G), each mapped onto the PPI network. (H–J) Expression levels of the proteins EFNA5 (H), PDE4D (I), and MAP2K1 (J)—target molecules estimated for CRC—were assessed in five tissues and categorized as four expression levels (“not detected,” “low,” “medium,” and “high”) using The Human Protein Atlas ( https://www.proteinatlas.org/ ). Further details can be found in Figure S5 .

Next, we examined the regulatory relationships between therapeutic pathways and DEGPs based on their pathway categories in CRC ( Figure 2 B). Among therapeutic pathways, “Signal transduction” was found to regulate multiple DEGPs, including “programmed cell death.” Conversely, among DEGPs, pathways such as “DNA repair” and “cell-cell communication” were regulated by multiple therapeutic pathways. These findings suggest that targeting the identified therapeutic pathways could simultaneously modulate multiple DEGPs. Further details for CRC and other diseases are provided in the Figures S1 and S2 . In additionally, the identified therapeutic pathways were validated by analyzing the network topology between the therapeutic pathways and DEGPs ( Figure S3 ).

Network controllability-derived therapeutic target molecules reveal additional intervention opportunities across diseases

To uncover the therapeutic target molecules regulating each therapeutic pathway, we applied network controllability analysis using the PPI network ( Table S2 ). For CRC, we examined the overlap among the estimated therapeutic target molecules, the known therapeutic target molecules, and causal genes ( Figure 2 C), and identified six proteins, including EGFR and KIT, that were present in all three categories. Ten proteins, including MTOR and MAPK1, were both estimated and known therapeutic target molecules, and 23 proteins, including SRC and SMAD3, were both estimated therapeutic target molecules and causal genes. Notably, SRC has been recently identified as a potential inhibitory target molecule for CRC. 20 Estimated therapeutic target molecules that did not overlap with known therapeutic target molecules included EFNA5, PDE4D, and MAP2K1, which are potential inhibitory target molecules for CRC. 21 , 22 , 23 EFNA5, identified as a putative therapeutic target molecule in the Axon guidance pathway ( Figures 2 D and 2E), may be an effective treatment for CRC by interacting with causal genes such as SRC . 20 PDE4D, a candidate target molecule in the GPCR downstream pathway ( Figures 2 D and 2F), may support CRC treatment by interacting with proteins such as APP, which enhances proliferation via the phosphorylation of ERK pathway proteins. 24 MAP2K1, an estimated target molecule of the MAPK1/MAPK3 signaling pathway ( Figures 2 D and 2G), may participate in CRC treatment through interactions with causal genes including EGFR, KIT, RAF1, and RET. Importantly, these estimated therapeutic target molecules were expressed at medium to high protein expression levels in the colon and rectum ( Figures 2 H–2J). These findings suggest that the estimated therapeutic target molecules encompass both established and previously unrecognized candidates for therapeutic intervention. We also validated the therapeutic potential of the estimated target molecules for CRC and other diseases by examining their functional roles within the PPI network ( Figures S4–S9 and S17 ).

SYNERGIE achieves accurate prediction of mechanistically different drug synergies across diseases

We evaluated the predictive performance of SYNERGIE for two types of synergistic doublets across various diseases: (1) the boost type and (2) the complement type. To facilitate comparisons with previous methods that did not differentiate between these two synergy types, we also assessed the performance of SYNERGIE for (3) the “indistinct type,” which includes synergistic drug doublets combining both boost and complement types without distinction. We compared SYNERGIE with two baseline methods: the topology-based method 13 and the opposite-direction method. 12 SYNERGIE was tested in three variations: using causal genes (SYNERGIE-CG), known therapeutic target molecules (SYNERGIE-KT), and estimated therapeutic target molecules (SYNERGIE-ET). The performance was evaluated using gold standard data for 14 diseases (see STAR Methods ) and measured using three accuracy metrics: area under the receiver operating characteristic curve (AUC), area under the precision-recall curve (AUPR), and Boltzmann-enhanced discrimination of the AUC (BED AUC). 25

We first evaluated the performance of boost-type doublets across 14 diseases. No significant differences were observed among the three variations of SYNERGIE ( Figure 3 A). For BED AUC, SYNERGIE-ET was the most accurate for CRC and BC ( Figure 3 B). Although SYNERGIE-KT was not applicable to DLBCL owing to the lack of known therapeutic target molecules, SYNERGIE-ET still exhibited higher accuracy. These findings suggest that SYNERGIE-ET is effective for predicting boost-type doublets and remains applicable to diseases without known therapeutic target molecules. Figure 3 Performance evaluation of the prediction of three types of synergistic drug combinations in 14 diseases (A) Comparison of the proposed methods for predicting the boost-type drug combinations. The proposed methods were those utilizing estimated therapeutic target molecules (SYNERGIE-ET), known therapeutic target molecules (SYNERGIE-KT), and causal genes (SYNERGIE-CG). Boxplots display the AUC score distributions for diseases. The horizontal axis displays the proposed methods, and the vertical axis indicates the AUC scores. The colors correspond to the proposed methods. The dotted horizontal lines in the top panels show the AUC = 0.5. (B) Bar graph representation of the boost-type drug combinations for each disease. The horizontal axis represents diseases, and the vertical axis indicates the AUC, AUPR, and BED AUC scores. Other formats follow the same format as (A). (C) Comparison of the proposed methods for predicting the complement-type drug combinations. This follows the same format as (A). (D) Bar graph representation of the complement-type drug combinations for each disease. This follows the same format as (B). (E) Comparison of the proposed and baseline methods for predicting the indistinct-type drug combinations. The baseline methods were the topology-based method and the opposite-direction method. Statistical significance was assessed using Wilcoxon signed-rank test with Benjamini–Hochberg correction for multiple comparisons. Asterisks denote statistical significance: ∗ p < 0.05; ∗∗ p < 0.01; ∗∗∗ p < 0.001. The p values between the SYNERGIE-ET and topology-based method, or opposite direction method are p ∗ = 1.2 × 10 −4 and p ∗ = 5.3 × 10 −3 . Other formats follow the same format as (A). (F) Bar graph representation of the indistinct-type drug combinations for each disease. This follows the same format as (B).

Next, we evaluated the performance of complement-type doublets across 14 diseases. SYNERGIE-ET was particularly effective for CRC, RCC, and DLBCL, where known therapeutic target molecules were unavailable ( Figure 3 D). These findings indicate that SYNERGIE is highly effective for the prediction of complement-type doublets.

Finally, we compared the performance of indistinct-type doublets across 14 diseases between the proposed and baseline methods. SYNERGIE-ET outperformed the topology-based method ( p ∗ <1.2 × 10 −4 ) and the opposite-direction method ( p ∗ <5.3 × 10 −3 ), where p ∗ indicates the adjusted p -value ( Figure 3 E). Additionally, SYNERGIE-KT and SYNERGIE-CG outperformed the topology-based method (SYNERGIE-KT: p ∗ <9.6 × 10 −3 ; SYNERGIE-CG: p ∗ <2.4 × 10 −4 ) as well as the opposite-direction method (SYNERGIE-KT: p ∗ <1.2 × 10 −4 ; SYNERGIE-CG: p ∗ <3.4 × 10 −3 ). For AUPR and BED AUC, SYNERGIE-ET achieved the highest performance for CRC, BC, and DLBCL, even in the absence of known therapeutic target molecules ( Figure 3 F). These results demonstrate that SYNERGIE can robustly predict synergistic drug combinations, even when synergy types are not explicitly distinguished and known therapeutic targets are unavailable.

SYNERGIE identifies clinically relevant drug combinations across a range of diseases

We systematically predicted drug combinations of three mechanistic synergy types across 14 diseases. Using SYNERGIE-ET, we explored 500,500 possible drug pairs involving 1001 orally administered drugs ( Tables S5 , S6 , and S7 ).

For CRC, we focused on predicted complement-type drug combinations and cross-referenced clinical trial data from ClinicalTrials.gov ( Figure 4 A; Table S3 A). All of the predicted drugs except chlorotrianisene, ponatinib, and lucitanib have been tested for therapeutic or preventive effects against CRC in clinical trials ( Tables S3 B and S3C). Aspirin (Asp) has undergone multiple phase 1–3 clinical trials for its preventive effects in patients with CRC. Similarly, ibrutinib (Ibru), in combination with pembrolizumab, was tested in phase 1/2 clinical trials for metastatic CRC, showing limited antineoplastic effects but good tolerability. 26 Figure 4 In vitro experimental validation of the complement-type doublets and triplets that were predicted via SYNERGIE for CRC (A) The top 18 predicted complement-type drug combinations for CRC using the SYNERGIE-ET. The horizontal axis represents the drug combinations, whereas the vertical axis corresponds to the five types of prediction scores (see STAR Methods ). Blue and green indicate the prediction scores for drug A and drug B, respectively, with the color intensity reflecting the magnitude of the scores. Details can be found in Table S3 . (B–I) Comparison of the survival ratios of CRC cells among Asp, an antineoplastic drug, and their combination. Gray represents single drugs (singlets), and blue represents the combination (doublet). Data are presented as mean values from independent experiments. Standard deviation values are provided in the Table S3 . (J–O) Comparison of the survival ratios of CRC cells among PB, an antineoplastic drug, and their combination. These panels follow the same format as (B–I). (P) Heatmap of the dose-response synergy scores of Asp doublets. The vertical axis represents the doses of Asp, and the horizontal axis shows drugs predicted to be Asp partners. The color scale reflects the IA scores, with positive values indicating synergistic effects and negative values indicating nonsynergistic effects. Each box is annotated with the corresponding IA score. Data are presented as mean values from independent experiments. Standard deviation values are provided in the Table S3 . (Q) Comparison of the dose-response survival ratios between Oxa, a first-line drug for CRC, alone (singlet) and the Oxa-Asp doublet. In the Oxa-Asp doublet, the dose of Asp was fixed at 1 mM while the dose of Oxa was varied. The horizontal and vertical axes represent the doses of Oxa (μM) and the mean survival ratios of CRC cells, respectively. Error bars represent standard deviations (SD). The gray and blue points correspond to Oxa and the Oxa-Asp combination, respectively. The color bar reflects the IA scores for each drug combination, with values indicating the IA scores. (R) Comparison of CRC cell survival ratios among singlets, doublets, and the triplet combination: Oxa, Asp, and Ibru. The horizontal axis represents the drugs, and the vertical axis represents the mean survival ratios. Error bars indicate SD. The color bar reflects the IA scores for each drug combination, with values indicating the IA scores. Disease abbreviations: CRC, colorectal cancer. Drug abbreviations: Asp, aspirin; Chl, chlorotrianisene; Ced, cediranib; Das, dasatinib; Ibru, ibrutinib; Inf, infigratinib; Luc, lucitanib; Mir, mirdametinib; Oxa, oxaliplatin; PB, sodium phenylbutyrate; Pim, pimasertib; Pon, ponatinib; Sar, saracatinib; Tac, tacedinaline; Tan, tanespimycin; and Van, vandetanib.

We further validated the predicted drug combinations using independent resources not included in the learning process. For example, Asp use before the age of 70 years has been shown to reduce the risk of CRC. 27 Sodium phenylbutyrate (PB) exhibits antitumor effects on CRC, with enhanced cytotoxicity when combined with 5-fluorouracil and irinotecan. 28 Additionally, drugs such as pimasertib, saracatinib, ponatinib, tanespimycin, and cediranib have demonstrated therapeutic effects against CRC. 29 , 30 , 31 , 32 , 33 Further details on CRC and other diseases are provided in the Tables S3 and S4 . Overall, these results demonstrate SYNERGIE’s ability to uncover clinically plausible and mechanistically diverse drug combinations for a range of diseases, and provide validation using independent sources.

The Ibru-Asp doublet predicted by SYNERGIE synergizes with first-line CRC therapy and demonstrates clinical potential in in vitro validation

To validate the complement-type drug combinations predicted by SYNERGIE, we conducted in vitro cell survival assays using HT29 human CRC cells, focusing on 14 out of the 18 predicted pairs ( Table S3 ). Combinations involving chlorotrianisene were excluded from validation as the drug is now discontinued. HT29 cells were treated with these drug combinations for 72 h, and cell survival was measured using the water-soluble tetrazolium (WST) assay.

The survival ratios of singlets versus doublets predicted for CRC were compared. Figures 4 B–4I and 4J–4O present the survival ratios for the Asp-containing and PB-containing combinations, respectively. Asp alone slightly reduced the survival ratio, but its combination with Ibru, tacedinaline, or vandetanib further enhanced this reduction ( Figures 4 B–4I). Conversely, neither PB alone nor PB-containing combinations significantly affected survival ratios ( Figures 4 J–4O), potentially because HT29 cell were used, as PB is known to exhibit antitumor effects in CRC. 28 These findings suggest that Asp-containing complement-type combinations are promising for CRC treatment.

We then evaluated the significance of synergistic effects of the 14 drug combinations using the independent action (IA) model, 34 an indicator of drug synergy ( Table S3 D). Drug pairs with an IA score greater than 0 were considered synergistic. Among the combinations, Asp with Ibru, dasatinib, or vandetanib were clearly synergistic ( Figure 4 P). Although Asp is recognized for its preventive effects against CRC, 35 no evidence supports its effectiveness as a single-drug treatment. These results suggest that Asp, a non-antineoplastic drug, enhances the antitumor effects of antineoplastic drugs in multidrug therapies.

We also assessed the synergistic effects of combining Asp with Oxa, a first-line CRC drug ( Figure 4 Q). At doses lower than its IC50, synergistic effects were observed when Oxa was combined with 1 mM Asp. This indicates that Asp can potentiate the efficacy of low concentrations of Oxa, affording a potential reduction in Oxa-related side effects while enhancing its effectiveness against difficult-to-treat tumors.

Finally, we evaluated the triplet combination of Oxa with the Ibru-Asp doublet predicted for CRC ( Figure 4 R). Synergy was observed for all doublets involving these three drugs, and the triplet combination further reduced CRC cell survival ratios compared to the doublets, indicating strong synergistic effects ( Table S8 ). These findings suggest that the Oxa-Ibru-Asp triplet exhibits significant synergistic potential against CRC and holds considerable clinical promise, with the non-antineoplastic agent Asp enhancing the efficacy of standard therapies in specific combinations.

Mechanistic analysis reveals complementary and partial boost-like synergy of the Oxa-Ibru-Asp triplet in drug-disease and drug-drug contexts

We explored the synergistic mechanisms of the Oxa-Ibru-Asp triplet, validated using in vitro experiments, by analyzing the gene expression patterns in CRC and following induction of the triplet. Transcriptome data at 24 h and 9 months post-Oxa induction in HT29 cells were obtained from the Gene Expression Omnibus (GEO) database 36 (see STAR Methods ).

First, we examined the drug-disease associations between each drug in the triplet and CRC by comparing the gene expression patterns ( Figures 5 A–5C). For Ibru ( Figure 5 A), MYC , CDC25A , and IER3 , which are highly expressed in CRC, were downregulated in response to Ibru. MYC , overexpressed 7-fold in CRC owing to mutations in APC and β-catenin , as well as the activation of growth factors, signaling pathways, hormones, and inflammation, regulates 215 metabolic responses. 37 CDC25A , involved in the G2/M transition of the cell cycle, is linked to the poor prognosis of CRC. 38 IER3 plays a critical role in stress responses in colon epithelial cells. 39 For Asp ( Figure 5 B), CDC45 and CCNB1 , which are both upregulated in CRC, were downregulated in Asp. CDC45 , a target molecule of the KNK437 protein, suppresses CRC growth and metastasis when regulated. 40 CCNB1 , which is overexpressed via Chk1 activation in CRC, promotes tumor cell growth. 41 For Oxa ( Figure 5 C), CCNB1 and CDKN2A , which are upregulated in CRC, were downregulated by Oxa. CDKN2A , which is significantly overexpressed in CRC, is associated with poor prognosis. 42 These findings suggest that Ibru, Asp, and Oxa exert therapeutic effects by inversely regulating key genes involved in CRC progression and prognosis. Figure 5 Omics analyses of the synergistic effects of the predicted Ibru-Asp doublet combined with Oxa, a first-line drug in CRC (A–C) Drug-disease associations based on the transcriptome pattern comparisons (A) between Ibru and CRC, (B) between Asp and CRC, and (C) between Oxa and CRC. The horizontal axis represents the transcriptome signature score after drug treatment in HT29 cells, and the vertical axis represents the CRC transcriptome signature scores. The color intensity corresponds to the significance of the gene expression ratios in CRC. Gene expression changes were assessed using linear models with empirical Bayes moderation, with p values adjusted using the Benjamini–Hochberg method. Exact p values are provided in Table S14 . The pink line represents the regression line, and the light pink area around the line represents the 95% confidence interval (CI). (D–F) Drug-drug associations based on the transcriptome pattern comparisons (D) between Ibru and Asp, (E) between Ibru and Oxa, and (F) between Asp and Oxa. Horizontal and vertical axes represent the transcriptome signature scores after each drug treatment in HT29 cells. The color intensity corresponds to the significance of the gene expression ratios in CRC. Gene expression changes were assessed using linear models with empirical Bayes moderation, with p values adjusted using the Benjamini–Hochberg method. Exact p values are provided in Table S14 . The pink line represents the regression line, and the light pink area around the line represents the 95% confidence interval (CI). This panel follows the same format as (A–C). (G and H) Venn diagrams showing the distribution of (G) upregulated and (H) downregulated DEGs in CRC and the Oxa–Ibru–Asp triplet. (I and J) Heatmaps of all pathways significantly enriched for (I) upregulated and (J) downregulated DEGs in singlets, doublet, and triplet. The horizontal axis represents pathways, and the vertical axis represents drug combinations. The heatmap’s color intensity represents − log ⁡ ( p - value ) , with pathway categories shown in the color bar above the heatmap. The asterisks for each square reflect statistical significance (∗ p < 0.05; ∗∗ p < 0.01; ∗∗∗ p < 0.001). Enrichment analysis was performed by Fisher’s exact test. Corrections for multiple testing were applied, adjusting significance values based on the number of pathway terms. Exact p values are provided in Tables S15 and S16 . Further details can be found in Table S9 . (K and L) Synergistic mechanisms of the Oxa-Ibru-Asp triplet. (K) On the left, the diagram shows the regulatory flow from BCR signaling to calcium signaling, coordinated upstream by Ibru and downstream by the triplet. On the right, the diagram shows the regulatory flow from AMPK signaling to PIP3K/Akt/mTOR signaling, coordinated upstream by the triplet and downstream by Oxa. (L) Synergistic mechanism of the Oxa–Asp doublet. The diagram illustrates the regulatory flow from DNA synthesis to protein synthesis, coordinated upstream by Asp and downstream by Oxa. Created with BioRender.com . (M and N) Enrichment frequency of pathways important for CRC treatment with significantly (M) upregulated and (N) downregulated DEGs in each drug. Gray corresponds to singlets; red or blue correspond to the doublets or triplet. (O and P) Heatmaps of enriched pathways important for CRC treatment using (O) upregulated and (P) downregulated DEGs in the singlets, doublets, and triplet. The histograms on the horizontal and vertical axes represent the enrichment frequency of the pathways. The heatmap’s color intensity represents − log ⁡ ( p value ) , with pathway categories shown in the color bar above the heatmap. Enrichment analysis was performed by Fisher’s exact test. Corrections for multiple testing were applied, adjusting significance values based on the number of pathway terms. Exact p values are provided in Tables S17 and S18 . This panel follows the same format as (I–J). Drug abbreviations: Asp, aspirin; Ibru, ibrutinib; Oxa, oxaliplatin.

We then analyzed drug-drug associations by comparing gene expression patterns between drugs in the triplet ( Figures 5 D–5F). There were minimal or no significant correlations for all doublets involving the triplet (Ibru-Asp: cor = 0.093, p = 3.5 × 10 −4 ; Ibru-Oxa: cor = −0.026, p = 0.42; Asp-Oxa: cor = 0.014, p = 0.66). These low correlation coefficients suggest that the drugs act through different pathways, which is consistent with a complement-type mechanism. With Ibru and Asp ( Figure 5 D), ERBB3 was upregulated by Ibru but downregulated by Asp. Although ERBB3 is a therapeutic target molecule in metastatic CRC, 43 its upregulation by Ibru may be mitigated by Asp. With Ibru and Oxa ( Figure 5 E), Ibru downregulated MYC and CCND1 , both critical CRC targets, whereas Oxa had little effect. CCND1 is associated with the G1/S phase transition of cell cycle, and its upregulation is a marker of poor prognosis in CRC. 44 For Asp and Oxa ( Figure 5 F), CAPN1 , which is slightly downregulated by Oxa but upregulated by Asp, participates in apoptotic signaling and platinum resistance. 45 Using transcriptome data from 9-month post-Oxa treatment, we confirmed that Asp could alleviate CAPN1 downregulation and potentially address platinum resistance ( Figures S10 and S11 ). CDKN2A and CCNB1 , critical CRC target molecules, were downregulated by both Oxa and Asp. Together, these results show that all doublets involving the triplet achieve mechanistic synergy via the complementary regulation of CRC-related genes, with partial boost-like interactions, supporting their potential as a multi-target therapeutic strategy. We also confirmed the synergistic effects of the triplet in relation to transporters, mutations, and causal genes ( Figures S12–S14 ).

The Oxa-Ibru-Asp triplet broadly and coordinately regulates key pathways for CRC treatment beyond the reach of single-drug therapies

To elucidate the biological relevance of the triplet, we evaluated the associations between CRC and the triplet by comparing CRC-associated DEGs with drug-induced DEGs ( Figures 5 G and 5H). In terms of drug-disease associations, there was very little overlap of upregulated and downregulated DEGs between CRC and Ibru, or between CRC and Asp, suggesting that neither Ibru nor Asp exacerbates CRC. In terms of drug-drug associations, the minimal overlap between Ibru and Asp suggests that these drugs can complement each other. Furthermore, the small overlaps between Ibru and Oxa and between Asp and Oxa indicate that the triplet can exhibit boost-like effects in certain processes. In addition, we also confirmed that several DEGs were inversely expressed between CRC and each drug ( Figure S15 ); these results suggest that the Oxa-Ibru-Asp triplet can counteract CRC in a complementary manner.

We conducted pathway enrichment analyses of the DEGs of each drug to elucidate the mechanisms of action of the triplet ( Figures 5 I and 5J; Table S9 ). In the pathway analyses of upregulated DEGs, we focused on pathways that were not enriched by single drugs but instead on those that were enriched by the doublets or the triplet ( Figure 5 I). “Calcium reabsorption” and “vascular smooth muscle contraction” were not enriched by single drugs but were enriched by the triplet. For “calcium reabsorption,” the expression status of calcium-sensing receptors reduces the risk of CRC. 46 The calcium signaling pathway is a downstream process of BCR signaling, which is a target of Ibru, suggesting that Ibru regulates the upstream pathway, whereas the triplet regulates the downstream pathway ( Figure 5 K). Additionally, “vascular smooth muscle contraction” has been associated with CRC treatment. 47 These results suggest that these pathways may contribute to the mechanism of the synergistic effects of the triplet.

Next, we focused on pathways enriched in an Ibru- or Asp-specific manner in terms of the upregulated DEGs ( Figure 5 I). Ibru-specific enriched pathways included “inositol phosphate metabolism,” and “PPAR signaling,” both of which are associated with the oncogenic processes of CRC. 48 , 49 The Asp-specific enriched pathways included “ECM-receptor interaction,” and “valine, leucine and isoleucine,” which are also associated with the oncogenic processes of CRC. 50 , 51 These findings indicate that coordination between single-drug-specific and combination-specific pathways contributes to synergistic efficacy.

In the pathway analyses of downregulated DEGs, we focused on pathways that were not enriched by single drugs but rather instead by the doublets or triplet ( Figure 5 J). “ErbB signaling” and “MAPK signaling” were enriched specifically for the Ibru–Asp doublet, both of which are crucial in CRC treatment, 1 indicating their potential contribution to the synergistic effects. “Apoptosis,” “AMPK signaling,” and “pentose phosphate pathway (PPP)” were not enriched by either single drugs or the doublet but were enriched in the Oxa-Ibru-Asp triplet. Regarding “AMPK signaling,” inhibiting AMPK suppresses tumor growth and development by suppressing most anabolic processes 52 ; AMPK also inhibits PIP3K/Akt/mTOR, a key target for CRC therapy, 1 reducing protein synthesis. The PIP3K/Akt/mTOR signaling pathway was specifically enriched by Oxa ( Figures 5 I and 5J), indicating that the triplet regulates the upstream AMPK signaling, whereas Oxa regulates the downstream PIP3K/Akt/mTOR signaling ( Figure 5 K). Regarding “PPP,” its activation in CRC promotes the production of ribose 5-phosphate (R5P), which facilitates nucleotide synthesis. Metabolic pathways relating to “nucleotide,” “purine,” and “pyrimidine” were specifically enriched by Asp; DNA is a target of Oxa, and protein synthesis via mTOR is regulated by Oxa. Thus, the synthesis of purines, pyrimidines, and deoxyribonucleotides is regulated by Asp, ribonucleotide synthesis is regulated by the triplet, and DNA and protein synthesis are regulated by Oxa. This combination regulates the DNA replication and RNA synthesis pathways ( Figure 5 L). Consequently, results suggest that these pathways may contribute to the synergistic effects of the Oxa-Ibru-Asp triplet.

In addition, we analyzed pathways enriched in Ibru- or Asp-specific manners relating to the downregulated DEGs ( Figure 5 J). The Ibru-specific enriched pathways included “apelin signaling,” and “spliceosome,” both of which are involved in CRC growth and progression. 53 , 54 These results suggest that the coordination of singlet-specific, doublet-specific, and triplet-specific pathways enhances therapeutic efficacy.

Finally, we examined whether the enriched pathways were related to five processes critical for CRC treatment 1 ( Figures 5 M and 5N): (1) the WNT system, (2) the RAS system, (3) the mTOR system, (4) the TGFβ system, and (5) the p53 system. For the upregulated DEGs, these critical pathways were more frequently enriched by the triplet than by the single drugs ( Figure 5 M). For downregulated DEGs, these important pathways were most frequently enriched by the triplet, followed by the Ibru-Asp doublet, and then by the single drugs ( Figure 5 N). We also examined the distribution of the key enriched pathways ( Figures 5 O and 5P). Among upregulated DEGs, “PI3K-Akt signaling” and “cell cycle” were most significantly enriched by the triplet ( Figure 5 O). Among the downregulated DEGs, “MAPK signaling” and “ErbB signaling” were not enriched by the single drugs but were enriched by the Ibru-Asp doublet ( Figure 5 P). Additionally, “apoptosis” was not enriched by the single drugs or the doublet but was enriched by the triplet, whereas “cell cycle” and “p53 signaling” were most significantly enriched by the triplet. Collectively, these results demonstrate that the Oxa-Ibru-Asp triplet orchestrates the complementary regulation of essential CRC pathways, including apoptosis, cell cycle, MAPK, and p53 signaling, beyond what can be achieved with single agents or doublets.

The Ibru-Asp doublet exerts synergistic effects in CRC patient cohorts at both bulk and single-cell levels

To evaluate the clinical relevance of the Ibru-Asp doublet, we validated its effects using large cohorts of patients with CRC. The CRC cohorts included The Cancer Genome Atlas (TCGA) COAD dataset 55 (244 patients), which was used for predictions, and GEO datasets 56 ( GSE3958268 : 557 patients; GSE14333 : 290 patients) that were not part of the training process. Because the primary target of Oxa is DNA, this validation focused specifically on the Ibru-Asp doublet.

We compared the survival ratios of patients with upregulated and downregulated expression of the target molecules ( Figures 6 A–6L). We focused on PTGS1 ( COX-1 ) and PTGS2 ( COX-2 ), which are target molecules of Asp, and BTK , which is a target molecule of Ibru. Patients in whom PTGS1 or PTGS2 are downregulation have minimal differences in survival ratios compared with patients with upregulated PTGS1 or PTGS2 ( Figures 6 A–6F), suggesting that Asp alone has limited therapeutic effects in CRC. Similarly, patients with downregulated BTK had slightly higher survival ratios than those with upregulated BTK ( Figures 6 G and 6K), indicating that Ibru alone also provides limited therapeutic effects for CRC. Notably, patients with the simultaneous downregulation of PTGS1 , PTGS2 , and BTK had higher survival ratios compared with patients with upregulated expression of these target molecules ( Figures 6 J–6L). These findings suggest that the Ibru-Asp doublet may exert synergistic effects in CRC treatment. Figure 6 Validation of synergistic effects of the Ibru-Asp doublet using cohorts of patients with CRC at both bulk and single-cell levels (A) Comparison of survival ratios between patients with high and low expression of PTGS1 , a target molecule of Asp, using the TCGA COAD cohort (244 patients). The two groups were compared by two-sided Log rank test. The horizontal axis represents the survival time, and the vertical axis represents the survival ratios. The green and pink lines indicate the low and high PTGS1 expression groups, respectively. (B) As for (A) but for the GSE39582 cohort (557 patients). The two groups were compared by two-sided Log rank test. (C) As for (A) but for the GSE14333 cohort (290 patients). The two groups were compared by two-sided Log rank test. (D) As for (A) but for PTGS2 , a target molecule of Asp, in the TCGA cohort. The two groups were compared by two-sided Log rank test. (E) As for (D) but for the GSE39582 cohort. The two groups were compared by two-sided Log rank test. (F) As for (D) but for the GSE14333 cohort. The two groups were compared by two-sided Log rank test. (G) As for (A) but for BTK , a target molecule of Ibru, in the TCGA cohort. The two groups were compared by two-sided Log rank test. (H) As for (G) but for the GSE39582 cohort. The two groups were compared by two-sided Log rank test. (I) As for (G) but for the GSE14333 cohort. The two groups were compared by two-sided Log rank test. (J) As for (A) but for PTGS1 , PTGS2 , and BTK , the target molecules of Asp and Ibru, in the TCGA cohort. The two groups were compared by two-sided Log rank test. (K) As for (J) but for the GSE39582 cohort. The two groups were compared by two-sided Log rank test. (L) As for (J) but for the GSE14333 cohort. The two groups were compared by two-sided Log rank test. (M) Distribution of cell types in single-cell data from patients with CRC (371,223 cells; mismatch repair proficient [MMRp], 28 samples; mismatch repair deficient [MMRd], 34 samples). The scatterplot was generated following the application of uniform manifold approximation (UMAP) to the single-cell data. Colors indicate different cell types. The accompanying bar plots represent the total cell proportions for each donor group, with values indicating the proportion of each cell type. (N) Gene expression patterns of the drug target molecules in different cell types. The horizontal axis represents the target molecules, and the vertical axis represents the cell types. The color of the circle shows the mean gene expression levels, and the size of the circle indicates the fraction of cells expressing the gene. Square colors correspond to the cell types.

We further investigated the cell type-specific synergistic effects of the Ibru-Asp doublet using single-cell transcriptome data from patients with CRC (371,223 cells; mismatch repair proficient [MMRp], 28 samples; mismatch repair deficient [MMRd], 34 samples). 57 In both MMRp and MMRd samples, the proportions of B cells, epithelial cells, and plasma cells decreased, whereas the proportions of T/NK/ILC and myeloid cells increased ( Figure 6 M). Next, we analyzed the expression patterns of target molecules for Ibru and Asp ( Figure 6 N). PTGS1 , PTGS2 , and BTK were upregulated in mast cells, indicating that mast cells are regulated by the Ibru-Asp doublet. PTGS2 was also upregulated in myeloid cells. Furthermore, off-target molecules of Ibru, such as ERBB2 , IKBKB , and MAPK3 , were upregulated in epithelial and stromal cells. Together, these results indicate that Ibru and Asp act through both complementary and partial boost-like mechanisms by targeting distinct and common CRC-associated cell types. Further details are provided in Figure S16 .

Discussion

This study introduces a paradigm in drug combination prediction by categorizing and mechanistically modeling synergy into two distinct types—“boost” and “complement”. To implement this concept, we developed SYNERGIE, a computational method for predicting such synergy-aware drug combinations across a wide range of diseases, including intractable and multifactorial conditions. The SYNERGIE method integrates drug-disease and drug-drug interactions across multiple omics layers using a Bayesian optimization framework over a PPI network. By combining this with network controllability analysis to estimate disease-specific therapeutic pathways and target molecules, SYNERGIE effectively identified disease-context-specific synergistic combinations, even for diseases lacking known targets. Furthermore, SYNERGIE can expand beyond doublet synergies to triplet synergies by integrating predicted doublets with first-line drugs for a given disease, enhancing its clinical relevance.

Drug synergism can be categorized into boost-type and complement-type. In the treatment of cancers, boost-type combinations—consisting predominantly of antineoplastic drugs—are commonly used, partly owing to their established clinical use and well-characterized mechanisms. 8 These combinations often leverage independent mechanisms while exerting convergent synergistic effects on cancer cells, enhancing therapeutic outcomes, inhibiting proliferation, and reducing drug resistance. However, these combinations also pose significant clinical challenges that may lead to treatment discontinuation, including increased cumulative toxicity and unexpected drug-drug interactions. 14 These limitations underscore the need to explore alternative synergy paradigms that offer efficacy without compromising safety.

Conversely, complement-type drug combinations for cancers can achieve strong therapeutic effects with minimal doses of antineoplastic drugs that are enhanced through complementary interactions with non-antineoplastic drugs. These combinations are particularly beneficial for elderly patients and patients with pre-existing conditions who are more likely to experience adverse effects and require higher safety standards. Thus, complement-type combinations are anticipated to improve patient quality of life by balancing safety and therapeutic effects. Additionally, this approach may enable the reevaluation of antineoplastic drugs that are currently considered insufficiently effective or are restricted in usage, thereby expanding treatment options. Despite their potential, however, complement-type combinations have been historically overlooked owing to their mechanistic complexity and unpredictable efficacy. By distinguishing these synergy modes computationally, SYNERGIE reveals latent therapeutic opportunities that might otherwise be missed in conventional screens, thereby broadening the landscape of combination therapy.

Most previous computational methods for drug combination prediction do not differentiate between types of synergism, although SYNERGIE can distinguish between them. We validated our SYNERGIE-based predictions of complement-type drug combinations for CRC using in vitro experiments and human cohort data ( Figure 4 ) to demonstrate the synergistic effects of certain drug combinations. Notably, synergy was observed when combining the predicted drug pair with Oxa, a standard first-line treatment for CRC. This result emphasizes the practical utility of SYNERGIE: it does not aim to replace established treatments, but rather to enhance them through the addition of rationally selected complementary agents. This complementarity improves clinical feasibility by building upon existing therapeutic frameworks, which may facilitate exploratory use in clinical settings compared to entirely new regimens. Therefore, SYNERGIE represents a promising tool that clinicians can use to identify clinically actionable drug combinations that balance innovation with translational plausibility.

In the future, the stratified synergy strategies may reshape the landscape of combination therapy. Boost-type combinations are expected to remain the mainstay therapy owing to their broad applicability, but the risk of increased side effects may pose significant challenges for boost-type combinations, potentially hindering the feasibility of clinical trials. In contrast, complement-type combinations offer a promising avenue for patient-tailored therapies, particularly where safety and tolerability are critical. Complement-type combinations, with a lower risk of side effects, are less challenging to evaluate in clinical settings. Our findings suggest that the generation of complement-type synergies by incorporating non-antineoplastic drugs can broaden the therapeutic arsenal for intractable disease conditions and thus redefine combination therapy beyond conventional cytotoxic strategies.

Limitations of the study

In this study, we defined complement-type synergy as a drug combination that enhances therapeutic efficacy by regulating distinct target molecules in a complementary manner. We therefore evaluated synergistic effects using Bliss’s IA model, 34 which assumes independence between the mechanisms of action of the two drugs. It should be noted, however, that the IA model is not specifically designed to quantify complement-type synergy, and the development of synergy metrics tailored for such mechanisms would likely allow more accurate assessment in future studies.

Although many methods have been proposed for quantifying synergy in doublets, 34 , 58 , 59 , 60 , 61 triplet synergy remains difficult to evaluate because there are few established synergistic metrics available. We therefore applied the IA model for triplets because it provides simple calculation and clear interpretation, which is suitable under data-limited conditions. IA scores were computed from survival ratios obtained from single-drug and double-drug treatments, following the same procedure used for doublets (see Table S8 ). Since only oxaliplatin concentrations were varied in the triplet experiments, the CA model was not applied. We expect that future development of dedicated triplet synergy metrics will enable more precise evaluation.

In this study, we performed survival analysis primarily to evaluate the validity of the predicted doublet and triplet drug combinations. To do so, we included only those patient samples (donors) whose expression directions for all target genes within a given doublet or triplet were consistently aligned. This filtering step substantially reduced the number of eligible samples. To ensure robustness despite the small sample size, we performed survival analyses across three independent cohorts. As larger cohorts become available in the future, more accurate and statistically robust validation will be possible.

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Yoshihiro Yamanishi ( yamanishi@i.nagoya-u.ac.jp ).

Materials availability

This study did not generate new unique reagents.

Data and code availability

• All datasets analyzed in this study were obtained from publicly available repositories. Accession numbers and URLs are listed in the key resources table . • The source code for SYNERGIE is publicly available at GitHub repository [ https://github.com/YamanishiLab/SYNERGIE ]. • Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Acknowledgments

This work was supported by Grant-in-Aid for 10.13039/501100001691 JSPS fellows ( 22KJ2494 , and 25KJ0159 , S.N.). This study was also supported by 10.13039/501100001691 JSPS KAKENHI (grant nos. 20H05797 , 21H04915 , 21K18327 , and 24K21338 , Y.Y.). This work was also supported by the JSPS, KAKENHI Grant-in-Aid for Transformative Research Areas (B) ( 20H05799 , M. G.). This study was also supported by the 10.13039/100007428 Naito Foundation , an Intramural Research Grant, Grant-in-Aid for JSPS Fellows Restart Postdoctoral Fellow (RPD) ( 22J40019 , M.I.), and Grant-in-Aid for Scientific Research (C) ( 22K12265 , M.I.). Graphical abstract, Figures 1 , 5 K, and 5L were created with BioRender.com under a publication license.

Author contributions

S.N. and Y.Y. designed and supervised the study. S.N. carried out the methodological development and analyzed the data. M.G. and Y.K. carried out the experimental validation. S.N. and M.G. prepared the manuscript’s first draft and figures. M.I. and J.T. carried out the data preparation and provided and the methodological advice. K.I. and Y.Y. provided constructive suggestions for results and discussion in manuscript completion. All authors read and approved the final version of the manuscript.

Declaration of interests

The authors declare no competing interests.

STAR★Methods

Key resources table

REAGENT or RESOURCE SOURCE IDENTIFIER Chemicals, peptides, and recombinant proteins Pimasertib Selleck Chemicals Cat# S12345 Infigratinib Selleck Chemicals Cat# S2183 Aspirin Tokyo Chemical Industry Cat# A0001 Sodium phenylbutyrate Tokyo Chemical Industry Cat# O0511 Ibrutinib MedChemExpress Cat# HY-10997 Lucitanib MedChemExpress Cat# HY-18951 Saracatinib Cayman Chemical Cat# 16292 Tanespimycin Cayman Chemical Cat# 10008397 Cediranib Cayman Chemical Cat# 19841 Tacedinaline Cayman Chemical Cat# 10010336 Vandetanib Cayman Chemical Cat# 17061 Critical commercial assays Cell Counting Kit-8 DOJINDO Cat# CK04-11 Deposited data Disease patient transcriptomic data TCGA, UCSC Xena Platform https://portal.gdc.cancer.gov/ , https://xenabrowser.net/datapages/ Normal tissue gene expression data Genotype-Tissue Expression (GTEx) project https://gtexportal.org/home/ CRC bulk cohort GEO GSE39582 , GSE14333 CRC single-cell cohort GEO GSE178341 Oxaliplatin-induced transcriptome GEO GSE30011 , GSE4238 Protein–protein interaction data yeast-two-hybrid high-throughput datasets, protein complexes, kinase–substrate pairs, metabolic enzyme-coupled interactions, signaling interactions, innate immune response interactions, and 3D structurally resolved protein-protein interactions http://www.interactome-atlas.org/ , http://mips.helmholtz-muenchen.de/corum , https://www.phosphosite.org/homeAction.action , http://signalink.org/ , https://www.innatedb.com/ , https://commons.cri.uchicago.edu/instruct/ Drug induced transcriptome data L1000 https://clue.io/lincs#access Drug–protein interaction data ChEMBL, MATADOR, Drug Bank, Psychoactive Drug Screening Program Ki, KEGG DRUG, BindingDB, Therapeutic Target Database https://www.ebi.ac.uk/chembl/ , http://matador.embl.de , https://go.drugbank.com/ , https://pdspdb.unc.edu/pdspweb/ , https://www.genome.jp/kegg/drug/ , https://www.bindingdb.org/rwd/bind/index.jsp , https://idrblab.net/ttd/ Known drug combination data DrugCombDB http://drugcombdb.denglab.org/main Causal genes DisGeNet (v7.0), OMIM, ClinVar, NHGRI-EBI GWAS catalog, GWASdb, PheWAS catalog https://disgenet.com/ , https://omim.org/ , https://www.ncbi.nlm.nih.gov/clinvar/ , https://www.ebi.ac.uk/gwas/ , http://jjwanglab.org/gwasdb , https://phewascatalog.org/phewas/#home Experimental models: Cell lines HT-29 ATCC Cat# HTB-38 Software and algorithms Code This study https://github.com/YamanishiLab/SYNERGIE Python v3.10 Python Software Foundation (PyPI) https://www.python.org/downloads/ scikit-learn v1.6.0 (Python library) PyPI https://pypi.org/project/scikit-learn/ scipy 1.14.1 (Python library) PyPI https://pypi.org/project/scipy/ pandas 2.1.3 (Python library) PyPI https://pypi.org/project/pandas/ cloudpickle v.3.0.0 (Python library) PyPI https://pypi.org/project/cloudpickle/ numpy v.1.26.2 (Python library) PyPI https://pypi.org/project/numpy/ gpyopt v.1.2.6 (Python library) PyPI https://pypi.org/project/GPyOpt/ matplotlib v.3.8.2 (Python library) PyPI https://pypi.org/project/matplotlib/ NetworkX v.3.4.2 (Python library) Hagberg et al. https://networkx.org Pyrwr (Python library) Jinhong Jung https://github.com/jinhongjung/pyrwr UCSCXenaTools (R package) Wang and Liu https://cran.r-project.org/web/packages/UCSCXenaTools/vignettes/USCSXenaTools.html limma (R package) Ritchie et al. https://bioconductor.org/packages/release/bioc/html/limma.html survival (R package) Therneau https://cran.r-project.org/web/packages/survival/index.html

Experimental model and study participant details

Cell line and culture conditions

The HT-29 cell line was purchased directly from the American Type Culture Collection (ATCC) in 2021. The cells were authenticated by the supplier using STR profiling before shipment. Upon receipt, the cells were expanded and used at low passage numbers (less than ten passages after receipt) to maintain their original characteristics. The cells were tested for mycoplasma contamination using PCR-based detection and were confirmed to be negative throughout the experimental period.

The human colorectal cancer cell line HT-29 (HTB-38; ATCC) was used for in vitro validation experiments. Cells were maintained in Dulbecco’s Modified Eagle Medium (DMEM) with high glucose (4.5 g/L; NACALAI TESQUE, INC., Japan), supplemented with 10% fetal bovine serum (FBS; Funakoshi Co., Ltd., Japan). Cultures were incubated at 37°C in a humidified atmosphere with 5% CO 2 . Passages between 5 and 15 were used for all experiments.

No animals or human subjects were used in this study.

Method details

Disease- and patient-specific transcriptome data

In this study, we selected 14 diseases for which known drug combination data were available in DrugCombDB 62 and corresponding drug-induced gene expression profiles were available in LINCS database. 63 We obtained transcriptome data for patients with 14 diseases from TCGA 55 or GEO 36 ( Table S10 ). The diseases were acute myeloid leukemia, breast cancer (BC), chronic myeloid leukemia (CML), CRC, melanoma, diffuse large B-cell lymphoma (DLBCL), glioblastoma, multiple myeloma, non-small cell lung cancer, ovarian cancer, prostate cancer, renal cell carcinoma (RCC), T cell acute lymphoblastic leukemia, and type 2 diabetes mellitus. For TCGA, cohort data were downloaded from UCSC Xena Platform 64 using UCSCXenaTools R package, 65 including transcriptomic, clinical, and survival data. Donors were selected based on primary tissue site and histological type from the clinical annotations (see Table S11 ). Normal tissue gene expression data from Genotype-Tissue Expression (GTEx) project 66 were also obtained via UCSCXenaTools, 65 with donors selected based on primary tissue site (see Table S11 ).

Then, we constructed patient-specific and disease-specific gene expression signatures. Differential expression analysis was performed using limma R package, 67 which are referred to as “patient-specific gene expression signatures.” The patient-specific signatures were then averaged for each disease to generate “disease-specific gene expression signatures.”

Drug-induced transcriptome data

Gene expression profiles after chemical treatments were obtained from LINCS. 63 We selected cell lines corresponding to each disease (see Table S10 ) where available; otherwise, we averaged the gene expression profiles of all cell lines. In total, 1001 approved and orally administered drugs were selected. We imputed missing values using a tensor decomposition algorithm. 68 The tensor decomposition method was performed according to the default settings described in the reference paper. Specifically, we used the tensor-train weighted optimization (TT-WOPT) algorithm with TT-ranks set to {1, 30, 30, 30, 1}.

Drug–protein interaction data

Interactions between small compounds and their target proteins were collected from seven public databases (ChEMBL, 69 MATADOR, 70 Drug Bank, 71 Psychoactive Drug Screening Program Ki, 72 KEGG DRUG, 73 BindingDB, 74 and Therapeutic Target Database 75 ): we extracted 35,366 known compound–protein interactions involving 4956 compounds and 3062 target proteins. Additionally, from a computationally predicted chemical–protein interactome dataset, we extracted 182,021 estimated compound–protein interactions involving 7468 compounds and 4001 target proteins. 76

Drug combination data

Drug combinations for the 14 diseases were obtained from DrugCombDB, 62 a database in which drug combinations are based on high-throughput screening, manual curation, and the FDA Orange Book. Drug combinations were classified into three interaction types based on the overlap of their target molecules. Specifically, for each disease, all known drug combinations were initially regarded as indistinct type. Among them, if two drugs share at least one common target molecule, the combination was labeled as boost type, whereas if the two drugs have no common target molecules, the combination was labeled as complement type. Drugs for which drug-induced gene expression signatures were available were selected. Table S12 summarizes the number of known drug combinations for each disease.

Protein–protein functional interaction (PPI) undirected and directed networks

We constructed PPI undirected and directed networks. We combined different types of functional associations between proteins into a single one (yeast-two-hybrid high-throughput datasets, 77 protein complexes, 78 kinase–substrate pairs, 79 metabolic enzyme-coupled interactions, 80 signaling interactions, 81 innate immune response interactions, 82 and 3D structurally resolved protein-protein interactions 83 ). These PPIs yielded 219,192 interactions among 12,648 proteins. For the PPI-directed network, we obtained a directed graph using the KEGGgraph R package. 84 From this graph, 63,409 interactions between 6354 proteins were identified.

Pathway–pathway association network

A pathway–pathway association network was constructed by converting protein-level directed networks (PPI-directed network) to pathway-level directed networks accounting for the biological response flow. We used a general approach for measuring group-level centralities by averaging or summing the individual centralities of group members. 85

First, the edge betweenness centrality was calculated for all edges in the PPI-directed network. We then annotated the KEGG and Reactome pathways 86 for each node (protein), where the only third layer Reactome pathways in the pathway hierarchy were used to remove global pathways with many genes and small pathways with a few genes. The edge betweenness centrality of each protein pair was then summed for each pathway pair; this was defined as “edge weight” for each pathway pair. For the same pathway pair (e.g., pathway A–pathway A), the sum of the edge betweenness centralities in the pathway was defined as the “edge weight.” These operations revealed 52,246 interactions between 431 pathways.

Estimation of therapeutic pathways from differentially expressed gene (DEG)-enriched pathways (DEGPs)

We defined therapeutic pathways as those that regulate DEGPs and exert therapeutic effects through drug regulation (e.g., inhibition or activation) for each disease. Pathways upstream of multiple DEGPs are assumed to be promising therapeutic pathways. Therefore, we estimated controllable therapeutic pathways based on the topology of the pathway–pathway association network ( Figure 1 A).

First, a weighted gene co-expression network analysis (WGCNA), 15 which is used to identify co-expression network modules, was performed for each disease using the patient-specific gene expression signatures. For each identified disease module, pathway enrichment analyses were then performed, and significant DEGPs were used in subsequent procedures.

Therapeutic pathways were then estimated for each DEGP based on the network topology of the pathway–pathway association network. Let G Path be a pathway–pathway association directed network with node (pathway) set V = { V 1 , V 2 , ⋯ , V r } . The adjacency matrix of G Path is the r × r matrix E , and its ( k , l ) entry is denoted as follows: (Equation 1) e k l = e ( V k , V l ) = { 1 , if V k ∼ V l 0 , otherwise .

Let the DEGP subset be V DEGP = { V 1 DEGP , V 2 DEGP , ⋯ , V r DEGP DEGP } ( r DEGP < r ) . We estimated a pathway with the largest weight of the input edge to the l -th DEGP as a therapeutic pathway as follows: (Equation 2) V k Target ∶ = arg max V k ∈ V e ( V k , V l DEGP ) .

As the pathway–pathway association network contains self-loops, DEGPs themselves can be their therapeutic pathways. Finally, these operations yielded a therapeutic pathway subset V Target = { V 1 Target , V 2 Target , ⋯ , V r Target Target } ( r Target ≤ r DEGP ) .

Therapeutic target molecule estimation based on the network topology

Essential nodes for controlling network systems (i.e., dominating nodes) can be identified from a network controllability analysis based on network topology. 87 , 88 The minimum dominating set (MDS)-based controllability analysis 89 , 90 is a popular approach to estimate the MDS, the smallest set of dominating nodes, where nondominating nodes must be controlled by at least one dominating node adjacent to nondominating nodes. Therapeutic target molecules for each disease were assumed to correspond to the dominating nodes (proteins) for controlling therapeutic pathways in the PPI network. Here, the problem of estimating therapeutic target molecules is formulated in the framework of network controllability analysis based on MDS ( Figure 1 B).

Let G PPI be a PPI undirected network, and let G l PPI be a PPI undirected network of the l -th therapeutic pathway ( l = 1 , 2 , ⋯ , r Target ) with node (protein) set V l = { v 1 , v 2 , ⋯ , v q } . The adjacency matrix of the l -th network, G l PPI , is the q × q matrix H ( G l PPI ) , and its ( i , j ) entry is denoted as follows: (Equation 3) h i j = h ( v i , v j ) = { 1 , if v i ∼ v j 0 , otherwise .

Here, we solve the problem of estimating MDSs as a binary integer programming problem. Each node i is assigned a binary integer variable representing a dominant node. Let λ i ∈ { 0 , 1 } be the class label of the dominating set assigned to the i -th node, where λ i = 1 means that the i -th node belongs to the dominating set, and λ i = 0 means that the i -th node does not belong to the dominating set. To compute the MDS of the l -th therapeutic pathway, we consider the minimization of the following linear function: (Equation 4) Minimize ∑ i = 1 q λ i (Equation 5) subject to λ i + ∑ j = 1 q h i j λ j ≥ 1 (Equation 6) λ i ∈ { 0 , 1 } .

The estimated dominating nodes were defined as therapeutic target molecule candidates. To prioritize the therapeutic target molecule candidates, we filtered them based on their reachability to the DEGP corresponding to the l -th therapeutic pathway. Finally, these operations yielded the therapeutic target molecules, C l Target = { c 1 , c 2 , ⋯ , c q ′ } ( q ′ ≤ q ) for the l -th therapeutic pathway. Table S13 shows the number of estimated target molecules for each disease.

Prediction of different types of synergistic drug combinations considering target commonality, DEG inversion, and random walk (SYNERGIE-ET)

The treatment of a disease requires the consideration of drug-response pathways, which involve upstream drug–target interactions and downstream changes in gene expression. This is of particular importance for diseases with complex pathologies, in which therapeutic target molecules are often broadly distributed across PPI networks, necessitating comprehensive evaluation. To predict synergistic drug combinations, we considered five types of interactions encompassing both drug–disease and drug–drug interactions at the multiomics levels.

(I) Drug–disease target commonality score (target commonality)

Complex diseases can be treated by simultaneously regulating multiple therapeutic target molecules. Consequently, drug pairs that regulate many therapeutic target molecules are assumed to exhibit synergy.

Let L be an index set of therapeutic pathways for a disease and C l Target a set of disease therapeutic target molecules for the l -th therapeutic pathway ( l ∈ L ). Additionally, let A Target and B Target represent a set of target molecules for drugs A and B, respectively. To quantify the commonality of disease therapeutic target molecules for a drug pair, the similarity between the set of disease therapeutic target molecules ( C l Target ) and the sets of drug target molecules ( A Target or B Target ) is computed for each therapeutic pathway. To adjust for the number of drug target molecules, the Kendall tau rank correlation coefficient was used to calculate the commonality of the l -th therapeutic pathway as follows: (Equation 7) s l Target common ( A Target , B Target , C l Target ) ∶ = Kendall ( A Target , C l Target ) + Kendall ( B Target , C l Target ) 2 , where Kendall ( A Target , C l Target ) represents the Kendall tau rank correlation coefficient between the set of drug A target molecules A Target and the set of disease therapeutic target molecules C l Target on the l -th therapeutic pathway. Similarly, Kendall ( B Target , C l Target ) denotes the Kendall tau rank correlation coefficient between the set of drug B target molecules B Target and the set of disease therapeutic target molecules C l Target . Finally, we measured drug–disease commonalities at the protein level by calculating the overall “drug–disease target commonality score,” combining the Kendall tau rank correlation coefficients across all therapeutic pathways, as follows: (Equation 8) s Target common ( A Target , B Target , C Target ) = ∑ l ∈ L s l Target common r Target , where r Target represents the total number of therapeutic pathways.

(II) Drug–disease DEG inversion score (DEG inversion)

In diseased states, gene expression systems are disrupted. However, normal states can be restored if disease-specific gene expression patterns are counteracted by drug-induced gene expression patterns. Therefore, drug pairs that reverse a substantial number of disease-related DEGs may be able to exert synergistic effects.

Let the disease-specific gene expression signature be z l = ( z 1 , z 2 , ⋯ , z q ) T , the drug-A gene expression signature be x l A = ( x 1 A , x 2 A , ⋯ , x q A ) T , and the drug-B gene expression signature be x l B = ( x 1 B , x 2 B , ⋯ , x q B ) T in the l -th therapeutic pathway ( l ∈ L ). Drug–disease inversions for the l -th therapeutic pathway at the transcriptomic level are quantified as the “drug–disease DEG inversion score,” based on the number of DEGs that show an inverse direction between the disease and drug A or drug B. The score is zdefined as (Equation 9) s l DEG invers ( x l A , x l B , z l ) ∶ = 1 2 ( ∑ i q w i x i A ′ ∑ i q w i + ∑ i q w i x i B ′ ∑ i q w i ) with x i A ′ = { 1 if x i A · z i < 0 0 otherwise } and x i B ′ = { 1 if x i B · z i < 0 0 otherwise } . In this formulation, w i represents the weight for the i -th DEG, which is computed based on the shortest distance from the disease therapeutic target molecules C l Target = { c 1 Target , c 2 Target , ⋯ , c q ′ Target } on the PPI undirected network G l PPI ( l ∈ L ) to the disease DEGs C l DEG = { c 1 DEG , c 2 DEG , ⋯ , c q DEG } as follows: (Equation 10) w i = max c j ∈ C l Target exp ( − d ( c i DEG , c j Target ) ) , where d ( c i DEG , c j Target ) is the shortest path distance. The overall drug–disease DEG inversion score for drugs A and B is computed by combining the scores across all pathways as follows: (Equation 11) s DEG invers ( x A , x B , z ) = ∑ l ∈ L s l DEG invers r Target , where r Target represents the total number of therapeutic pathways.

(III) Drug–drug target random walk score (target random walk score)

Interactions between target molecules of drugs are a crucial determinant of synergistic effects. Here, these interactions are quantified using the random walk with restart (RWR) 91 algorithm, which models the propagation of biological responses to drugs across a PPI network. The drug–drug target random walk scores are calculated as follows: (Equation 12) s Target random = 1 | A Target | + | B Target | ∑ a ∈ A Target , b ∈ B Target max ( rwr ( a Target , b Target ) , rwr ( b Target , a Target ) ) , where A Target represents the set of drug A target molecules, B Target represents the set of drug B target molecules, | A Target | and | B Target | are the number of target molecules for drugs A and B, respectively. The terms r w r ( a Target , b Target ) and r w r ( b Target , a Target ) represent the RWR scores from a target molecule a Target of drug A ( a ∈ A Target ) to a target molecule b Target of drug B ( b ∈ B Target ) , and vice versa. RWR computations were performed using the Python pyrwr package [ https://github.com/jinhongjung/pyrwr ].

(IV) Drug–drug DEG random walk score (DEG random walk score)

Interactions between the DEGs induced by drugs are another key factor influencing synergistic effects. Here, these interactions are quantified using the RWR algorithm, which models the propagation of biological responses to drugs across the PPI network. The drug–drug DEG random walk scores are calculated as follows: (Equation 13) s DEG random = 1 | A DEG | + | B DEG | ∑ a ∈ A DEG , b ∈ B DEG max ( r w r ( a DEG , b DEG ) , r w r ( b DEG , a DEG ) ) , where A DEG represents the set of DEGs for drug A, and B DEG represents the set of DEGs for drug B. The terms | A DEG | and | B DEG | represent the number of DEGs for drugs A and B, respectively. The RWR scores, rwr ( a DEG , b DEG ) and rwr ( b DEG , a DEG ) , capture the interaction strength from a DEG a DEG of drug A ( a DEG ∈ A DEG ) to a DEG b DEG of drug B ( b DEG ∈ B DEG ) and vice versa.

(V) Synergy type-specific score (STS score)

1. Boost type . Boost-type synergy occurs when two drugs have common target molecules. Therefore, drug pairs with overlapping target molecules are considered to exhibit boost-type synergy. The STS score was calculated based on the overlap of target molecules between drugs as follows: (Equation 14) s STS = 1 − exp ( − α | A Target ∩ B Target | ) , where | A Target ∩ B Target | denotes the number of shared target molecules between drugs A and B, and α is a hyperparameter (default: α = 5 ). 2. Complement type . Complement-type synergy occurs when two drugs regulate different target molecules in a complementary manner. Thus, drug pairs with a minimal number of common target molecules are assumed to exhibit complement-type synergy. The STS score is calculated as (Equation 15) s STS = exp ( − α | A Target ∩ B Target | ) , where | A Target ∩ B Target | indicates the number of shared target molecules between drugs A and B, and α is a hyperparameter (default: α = 5 ). 3. Indistinct type . For this type, we merge the boost and complement synergy types. Drug pairs with or without shared target molecules are considered to exhibit synergistic effects. The STS score is calculated as follows: (Equation 16) s STS = 0 .

Unless otherwise noted, the calculations related to the networks were performed using the Python NetworkX package. 92

Integration of drug combination scores using bayesian optimization

We combine the target commonality score, DEG inversion score, target random walk score, DEG random walk score, and STS score for each disease; this integration is possible using Bayesian optimization to identify the most optimal drug combinations.

Assume there are m diseases and n candidate drug pairs. We aim to predict which drug pair would form a synergistic drug combination for the t -th disease ( t = 1 , 2 , ⋯ , m ) . For the u -th drug pair ( u = 1 , 2 , ⋯ , n ), we integrated the five types of prediction scores: s t , u Target common , s t , u DEG inverse , s t , u Target random , s t , u DEG random , and s t , u STS . The weighted sum of the prediction scores of the u -th drug pair for the t -th disease can be expressed as follows: (Equation 17) γ t · s t , u ( t = 1 , 2 , ⋯ , m ; u = 1 , 2 , ⋯ , n ) , where γ t = ( γ t Target common , γ t DEG inverse , γ t Target random , γ t DEG random , γ t STS ) T is the weight vector for the prediction scores for the t -th disease, and s t , u = ( s t , u Target common , s t , u DEG inverse , s t , u Target random , s t , u DEG random , s t , u STS ) T is the prediction score vector for the u -th drug pair for the t -th disease. The weight vector was optimized using Bayesian optimization with a radial basis function (RBF) kernel to maximize the prediction accuracy, measured using Boltzmann-enhanced discrimination of the AUC (BED AUC). 25

We assume that the function g ( γ t ) is drawn from a Gaussian process prior, with observations { γ t , κ , o κ } κ = 1 K , where γ t , κ ( t = 1 , 2 , ⋯ , m ; κ = 1, 2, …, K ) induces a multivariate Gaussian distribution on R K , and o κ represents the corresponding BED AUC value. The distribution is determined using a mean function and a positive definite covariance function. From the prior and the observations, we infer a posterior distribution over the function, referred to as the acquisition function. We employed the upper confidence bound (UCB) 93 as the acquisition function. Under the Gaussian process prior, the function a UCB ( γ t ; { γ t , κ , o κ } , τ ) is defined as follows: (Equation 18) a UCB ( γ t ; { γ t , κ , o κ } , τ ) = − μ ( γ t ; { γ t , κ , o κ } , τ ) + φ σ ( γ t ; { γ t , κ , o κ } , τ ) , where τ is a hyperparameter in the Gaussian process, φ is an exploration weight controlling the balance between exploration and exploitation (default: φ = 2 ), μ ( γ t ; { γ t , κ , o κ } , τ ) is a predictive mean function, and σ ( γ t ; { γ t , κ , o κ } , τ ) is a predictive variance function. We estimated γ t , the value that maximizes the BED AUC, and obtained the optimized weight vector γ t ∗ as follows: (Equation 19) γ t ∗ = argmax γ t g ( γ t )

The obtained weight vector was then used to compute the prediction score of the u -th candidate drug pair for the t -th disease as follows: (Equation 20) γ t ∗ · s t , u ( t = 1 , 2 , ⋯ , m ; u = 1 , 2 , ⋯ , n ) .

Drug pairs with high scores are predicted to have synergistic effects in treating the corresponding disease.

For triplet construction, each triplet was generated by pairing a predicted doublet with a first-line drug ( Figure 1 D).

Proposed drug combination prediction method with causal genes (SYNERGIE-CG)

In SYNERGIE-CG, we formulate SYNERGIE using known causal genes. We treated known causal genes as therapeutic target molecules and identified pathways with causal genes as therapeutic pathways. Then, based on the therapeutic pathways and causal genes, we predicted drug combinations using the methods described in Equations 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 , 15 , 16 , 17 , 18 , 19 , and 20 . Causal genes for diseases were acquired from six public databases (DisGeNet 94 (v7.0), OMIM, 95 ClinVar, 96 NHGRI-EBI GWAS catalog, 97 GWASdb, 98 and PheWAS catalog 99 ). Table S13 provides the number of causal genes for each disease.

Proposed drug combination prediction method with known therapeutic target molecules (SYNERGIE-KT)

In SYNERGIE-KT, we formulate SYNERGIE using known therapeutic target molecules. We identified pathways with known therapeutic target molecules as therapeutic pathways to predict drug combinations, as described in Equations 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 , 15 , 16 , 17 , 18 , 19 , and 20 . The therapeutic target molecule information was compiled by manually curating data from medical monographs 100 and the KEGG DISEASE database. 73 Table S13 provides the number of known target molecules for each disease.

Baseline drug combination prediction method with network topology (topology-based method)

The topology-based method 13 hypothesizes that drug–drug relationships are characterized by a network-based proximity measure between the target molecules of the two drugs. The network proximity between drugs A and B is calculated as described in Equation 2 in Cheng et al., 13 referred to as the “separation score.” If the separation score has a positive value, drugs A and B are considered to be separated, and the combination is predicted to have a synergistic effect for treating a disease.

Baseline drug combination prediction method with the transcriptome opposite direction (opposite-direction method)

The opposite-direction method 12 hypothesizes that drug combinations with synergistic effects are characterized by the directions of gene expression signatures between drugs and diseases. The relationship between drugs A and B was calculated using the drug-induced gene expression signatures as described in Equation 2 in Stathias et al. 12 The relationship among drug A, drug B, and a specific disease was calculated using drug-induced and disease-specific gene expression signatures, as presented in Equation 3 of Stathias et al. 12 Finally, the prediction score for drugs A and B was obtained by integrating the two relationships, as described in Equation 4 in Stathias et al. 12 Drug pairs with high prediction scores are predicted to have synergistic effects in the treatment of a disease.

Performance evaluation procedure

We used three accuracy measures: the area under the receiver operating characteristic (ROC) curve (AUC), the area under the precision–recall (PR) curve (AUPR), and the Boltzmann-enhanced discrimination AUC (BED AUC). 25 ROC curves and PR curves for the performance of classifiers over all possible cutoffs were generated by plotting the true positive rates (TPRs) against the false positive rates (FPRs) and precision against recall, respectively. AUC scores range from 0 to 1.0, where 1.0 indicates perfect inference (100% TPR, 0% FPR), and 0.5 indicates random inference. AUPR scores range from 0 to 1.0, with 1.0 indicating perfect inference (100% precision, 100% recall). The BED AUC scores range from 0 to 1.0, with 1.0 indicating perfect inference in the early retrieval performance.

Chemicals used for the cell survival assay

Pimasertib and infigratinib were purchased from Selleck Chemicals (Houston, TX, USA). Aspirin (acetylsalicylic acid) and sodium phenylbutyrate (sodium 4-phenylbutyrate) were purchased from Tokyo Chemical Industries Co., Ltd. (Tokyo, Japan). Saracatinib, tanespimycin, cediranib, tacedinaline, and vandetanib were purchased from Cayman Chemical (Ann Arbor, Michigan, USA). Lucitanib and ibrutinib were purchased from MedChemexpress (South Brunswick Township, NJ, USA). The Cell Counting Kit-8 (CCK-8) used for the WST-8 assay was purchased from DOJINDO (Tokyo, Japan).

Cell culture and reagents

The HT-29 human CRC cell line (HTB-38; ATCC) was obtained from the ATCC (Rockville, MD, USA). Cultures were used between passages 5 and 15. HT-29 cells were cultured in Dulbecco’s modified Eagle medium with high glucose (4.5 g/mL) (NACALAI TESQUE, INC., JAPAN) supplemented with 10% fetal bovine serum (Funakoshi Co., Ltd., JAPAN). All cells were incubated at 37°C in a humidified atmosphere containing 5% (v/v) CO 2 .

Cell survival assay

In vitro growth inhibition was evaluated according to the manufacturer’s standard protocol using CCK-8. Cells were seeded in 96-well plates at a density of 5,000 cells/well in a total volume of 100 μL and exposed to various drugs for 72 h at 37°C in a 5% CO 2 atmosphere. WST-8 was added, and after 3 h, the absorbance was measured at a wavelength of 450 nm (ref. 630 nm) using a microplate reader (Bio-Rad Laboratories, Inc., Hercules, CA). The results are expressed as percentages [i.e., as the ratio of the absorbance of treated cells to that of the control (drug untreated group, 100%)]. Percent survival was calculated using the following formula: percentage survival = (absorbance of treated wells − absorbance of blank wells)/(absorbance of untreated wells − absorbance of blank wells) × 100.

Statistical evaluation of drug synergy based on in vitro experiments

The Bliss’s Independent Action (IA) model 34 and Loewer’s additivity (CA) model 58 were used to evaluate the synergistic effects of drug combinations.

In the IA model, it is assumed that the effects of drugs are stochastic events, and thus, it is appropriate to apply the model when all drugs to be used at the same time are considered to have different modes of action. The effects of a drug combination were calculated as the joint probability of each effect. Drug A has v% effects and Drug B has w% effects at a given combination of concentrations. The total effect rate of the combination can be computed as C I m i x = 1 − ( 1 − v ) ( 1 − w ) under the IA model’s assumption. Let, in the experiment, the combination of Drug A and Drug B obtain z % effect and if the ( 1 − z ) − C I m i x > 0 , a given drug combination is considered to have a synergistic effect based on Bliss’s IA model. IA score = ( 1 − z ) − C I m i x is denoted in this study.

In the CA model, it is assumed that the toxic unit equals one under the non-interaction assumption between drugs. The concentrations weighted by the effect of each drug (toxic unit [TU]) are added together to yield the TU of a given drug combination. (Equation 21) Toxic Unit ( TU ) = C A E C u A + C B E C u B

For a given drug combination, C A and C B stand for the concentrations of drug A and drug B, respectively. E C u A and E C u B are the concentrations of the drugs causing u% effect by drug A and drug B, respectively. If TU = 1, the effect rate of the drug combination remains at u% under the additive assumption. Thus, if TU < 1, a given drug combination is assumed to have synergistic effect based on the CA model, when u% is selected as the results of the growth inhibition assay. TU is denoted as the CA score in this study.

Validation of predicted complement-type drug combinations for CRC using TCGA cohort

To stratify patients based on the gene-expression ratios of the target genes of drugs, we defined the donors with z -scores of a target gene >1 as “high-expression donors”, and donors with z -scores of a target gene < −1 as low-expression donors.” The survival curves were fitted with Cox models based on a tabulation of the number at risk and number of events for each unique death time using the “survival” R package. 101

Validation of predicted complement-type drug combinations for CRC using GEO bulk and single-cell cohorts

We validated predicted drug combinations for CRC using an independent resource, CRC bulk cohort (made up of 557 donors) and CRC single-cell cohort (made up of 28 MMRp donors and 34 MMRd donors) from GEO 56 under the accession code GSE39582 and GSE178341 , respectively. The single cell data were processed according to the original paper. 57

We compared the levels of gene expression for drug targets of predicted drug combinations between healthy and CRC donors. To stratify patients based on drug-target gene expression, we calculated the gene expression ratios between patient and healthy donors. We defined donors with a fold change of a target gene > β as “high-expression donors” and those with a fold change of a target gene < δ as “low-expression donors.” Here, β and δ were determined so that the number of donors was at least 3. The survival curves were fitted with Cox models based on a tabulation of the number at risk and number of events for each unique death time using the “survival” R package. 101

Quantification and statistical analysis

Bayesian optimization was employed to integrate five types of drug–disease and drug–drug interaction scores for predicting synergistic drug combinations. For each disease, the optimal weight vector was estimated to maximize prediction accuracy, measured using BED AUC. The Gaussian process prior and UCB acquisition function were used for optimization. Performance metrics included AUC, AUPR, and BED AUC. For in vitro experiments, absorbance values were analyzed using the CCK-8 assay and normalized to untreated controls. Drug synergy was assessed using the Bliss independence model. Random walk with restart (RWR) was implemented using the pyrwr package. Network analyses were conducted using NetworkX.

All statistical analyses were performed using Python (version 3.10) with custom scripts. For the performance evaluation of proposed methods, we calculated AUC scores for 14 diseases in the indistinct-type combination prediction ( Figures 3 E and 3F). The performance comparison was made between SYNERGIE-ET and Topology-based method or Opposite direction method, between SYNERGIE-KT and Topology-based method or Opposite direction method, and between SYNERGIE-CG and Topology-based method or Opposite direction method. The two groups were compared by one-sided Wilcoxon signed-rank test. We made corrections for multiple testing using the false discovery rate (FDR) by Benjamini–Homberg (BH) corrections, adjusting significance values for two tests per analysis stream. The test was considered nominally significant with p < 0.05 and FDR-corrected significant if FDR <0.05. Statistical significance in figures is indicated as follows: ∗ p < 0.05, ∗∗ p < 0.01, ∗∗∗ p < 0.001.

For the omics analyses, differential expression analysis was performed ( Figures 5 A–5F). This analysis was conducted using linear models with empirical Bayes moderation implemented in limma package. The test was considered nominally significant with p < 0.05 and FDR-corrected significant if FDR <0.05.

For the KEGG pathway analyses ( Figures 5 I, 5J, 5O, and 5P), we performed enrichment analyses of the DEGs in singlets, doublet, and triplet by Fisher’s exact test. Corrections for multiple testing were applied, adjusting significance values based on the number of pathway terms. The test was considered nominally significant with p < 0.05 and FDR-corrected significant if FDR <0.05. Statistical significance in figures is indicated as follows: ∗ p < 0.05, ∗∗ p < 0.01, ∗∗∗ p < 0.001.

For the validation of synergistic effects of the Ibru–Asp doublet using CRC cohorts ( Figures 6 A–6L), we performed Cox regressions for the survival curves. We compared the survival curves between donors with low expressed target genes and donors with high expressed target genes. The two groups were compared by two-sided Log rank test. The test was considered nominally significant with p < 0.05.

References

1 Muzny D.M. Bainbridge M.N. Chang K. Dinh H.H. Drummond J.A. Fowler G. Kovar C.L. Lewis L.R. Morgan M.B. Newsham I.F. Comprehensive molecular characterization of human colon and rectal cancer Nature 487 2012 330 337 10.1038/nature11252 22810696 PMC3401966

2 Sun X. Vilar S. Tatonetti N.P. High-throughput methods for combinatorial drug discovery Sci. Transl. Med. 5 2013 205rv1 10.1126/SCITRANSLMED.3006667 24089409

3 Jia J. Zhu F. Ma X. Cao Z. Cao Z.W. Li Y. Li Y.X. Chen Y.Z. Mechanisms of drug combinations: interaction and network perspectives Nat. Rev. Drug Discov. 8 2009 111 128 10.1038/NRD2683 19180105

4 Al-Lazikani B. Banerji U. Workman P. Combinatorial drug therapy for cancer in the post-genomic era Nat. Biotechnol. 30 2012 679 692 10.1038/NBT.2284 22781697

5 Lehár J. Krueger A.S. Avery W. Heilbut A.M. Johansen L.M. Price E.R. Rickles R.J. Short G.F. Staunton J.E. Jin X. Synergistic drug combinations tend to improve therapeutically relevant selectivity Nat. Biotechnol. 27 2009 659 666 10.1038/NBT.1549 19581876 PMC2708317

6 Chibaudel B. Tournigand C. André T. de Gramont A. Therapeutic strategy in unresectable metastatic colorectal cancer Ther. Adv. Med. Oncol. 4 2012 75 89 10.1177/1758834011431592 22423266 PMC3296081

7 Falcone A. Ricci S. Brunetti I. Pfanner E. Allegrini G. Barbara C. Crinò L. Benedetti G. Evangelista W. Fanchini L. Phase III trial of infusional fluorouracil, leucovorin, oxaliplatin, and irinotecan (FOLFOXIRI) compared with infusional fluorouracil, leucovorin, and irinotecan (FOLFIRI) as first-line treatment for metastatic colorectal cancer: the Gruppo Oncologico Nord Ovest J. Clin. Oncol. 25 2007 1670 1676 10.1200/JCO.2006.09.0928 17470860

8 Xu M. Zhao X. Wang J. Feng W. Wen N. Wang C. Wang J. Liu Y. Zhao L. DFFNDDS: prediction of synergistic drug combinations with dual feature fusion networks J. Cheminf. 15 2023 33 10.1186/S13321-023-00690-3/TABLES/9 PMC10022091 36927504

9 Ma J. Motsinger-Reif A. Prediction of synergistic drug combinations using PCA-initialized deep learning BioData Min. 14 2021 46 10.1186/S13040-021-00278-3/TABLES/9 34670583 PMC8527604

10 Preuer K. Lewis R.P.I. Hochreiter S. Bender A. Bulusu K.C. Klambauer G. DeepSynergy: predicting anti-cancer drug synergy with Deep Learning Bioinformatics 34 2018 1538 1546 10.1093/BIOINFORMATICS/BTX806 29253077 PMC5925774

11 Zhang H. Goedegebuure S.P. Ding L. DeNardo D. Fields R.C. Province M. Chen Y. Payne P. Li F. M3NetFlow: A multi-scale multi-hop graph AI model for integrative multi-omic data analysis iScience 28 2025 111920 10.1016/J.ISCI.2025.111920 PMC11872513 40034855

12 Stathias V. Jermakowicz A.M. Maloof M.E. Forlin M. Walters W. Suter R.K. Durante M.A. Williams S.L. Harbour J.W. Volmar C.H. Drug and disease signature integration identifies synergistic combinations in glioblastoma Nat. Commun. 9 2018 5315 10.1038/s41467-018-07659-z 30552330 PMC6294341

13 Cheng F. Kovács I.A. Barabási A.L. László A.-L. Lászlóbarabási L. Network-based prediction of drug combinations Nat. Commun. 10 2019 1197 10.1038/s41467-019-09186-x 30867426 PMC6416394

14 Yanagawa S. Karakuchi N. Mochizuki T. Kodama S. Takeshima Y. Sumimoto K. Drug-Induced Interstitial Pneumonia due to Application of FOLFOX as Adjuvant Chemotherapy after Rectal Cancer Surgery: A Case Report and Literature Review. Case Rep Oncol. Times 13 2020 768 773 10.1159/000507985 PMC7383214 32774274

15 Zhang B. Horvath S. A general framework for weighted gene co-expression network analysis Stat. Appl. Genet. Mol. Biol. 4 2005 17 10.2202/1544-6115.1128 16646834

16 Duan J. Huang D. Liu C. Lv Y. Zhang L. Chang F. Zeng X. Li L. Wang W. Shao G. USP11-mediated LSH deubiquitination inhibits ferroptosis in colorectal cancer through epigenetic activation of CYP24A1 Cell Death Dis. 14 2023 402 415 10.1038/s41419-023-05915-9 37414755 PMC10326026

17 Maharati A. Moghbeli M. PI3K/AKT signaling pathway as a critical regulator of epithelial-mesenchymal transition in colorectal tumor cells Cell Commun. Signal. 21 2023 201 215 10.1186/S12964-023-01225-X 37580737 PMC10424373

18 Tyagi A. Sharma A.K. Damodaran C. A Review on Notch Signaling and Colorectal Cancer Cells 9 2020 1549 10.3390/CELLS9061549 32630477 PMC7349609

19 Hu X. Fatima S. Chen M. Xu K. Huang C. Gong R.H. Su T. Wong H.L.X. Bian Z. Kwan H.Y. Toll-like receptor 4 is a master regulator for colorectal cancer growth under high-fat diet by programming cancer metabolism Cell Death Dis. 12 2021 791 10.1038/s41419-021-04076-x 34385421 PMC8360949

20 Chen J. Elfiky A. Han M. Chen C. Saif M.W. The Role of Src in Colon Cancer and Its Therapeutic Implications Clin. Colorectal Cancer 13 2014 5 13 10.1016/J.CLCC.2013.10.003 24361441

21 Li S. Hou X. Wu C. Han L. Li Q. Wang J. Luo S. MiR-645 promotes invasiveness, metastasis and tumor growth in colorectal cancer by targeting EFNA5 Biomed. Pharmacother. 125 2020 109889 10.1016/j.biopha.2020.109889 32036212

22 Kim D.U. Nam J. Cha M.D. Kim S.W. Inhibition of phosphodiesterase 4D decreases the malignant properties of DLD-1 colorectal cancer cells by repressing the AKT/mTOR/Myc signaling pathway Oncol. Lett. 17 2019 3589 3598 10.3892/ol.2019.9996 30867802 PMC6396110

23 Kopetz S. Grothey A. Yaeger R. Van Cutsem E. Desai J. Yoshino T. Wasan H. Ciardiello F. Loupakis F. Hong Y.S. Encorafenib, Binimetinib, and Cetuximab in BRAF V600E–Mutated Colorectal Cancer N. Engl. J. Med. 381 2019 1632 1643 10.1056/nejmoa1908075 31566309

24 Pandey P. Sliker B. Peters H.L. Tuli A. Herskovitz J. Smits K. Purohit A. Singh R.K. Dong J. Batra S.K. Amyloid precursor protein and amyloid precursor-like protein 2 in cancer Oncotarget 7 2016 19430 19444 10.18632/ONCOTARGET.7103 26840089 PMC4991393

25 Truchon J.F. Bayly C.I. Evaluating virtual screening methods: Good and bad metrics for the “early recognition” problem J. Chem. Inf. Model. 47 2007 488 508 10.1021/CI600426E/SUPPL_FILE/CI600426E.PDF 17288412

26 Kim D.W. Tan E. Zhou J.M. Schell M.J. Martinez M. Yu J. Carballido E. Mehta R. Strosberg J. Imanirad I. Kim R.D. A phase 1/2 trial of ibrutinib in combination with pembrolizumab in patients with mismatch repair proficient metastatic colorectal cancer Br. J. Cancer 124 2021 1803 1808 10.1038/S41416-021-01368-Z 33828254 PMC8144213

27 Keum N. Giovannucci E. Global burden of colorectal cancer: emerging trends, risk factors and prevention strategies Nat. Rev. Gastroenterol. Hepatol. 16 2019 713 732 10.1038/S41575-019-0189-8 31455888

28 Al-Keilani M. Darweesh R. Abstract 1188: Sodium phenylbutyrate has an antineoplastic effect and enhances the cytotoxicity of 5-fluorouracil and irinotecan in colorectal cancer cell lines Cancer Res. 77 2017 1188 10.1158/1538-7445.AM2017-1188

29 Queralt B. Cuyàs E. Bosch-Barrera J. Massaguer A. de Llorens R. Martin-Castillo B. Brunet J. Salazar R. Menendez J.A. Synthetic lethal interaction of cetuximab with MEK1/2 inhibition in NRAS-mutant metastatic colorectal cancer Oncotarget 7 2016 82185 82199 10.18632/ONCOTARGET.11985 27636997 PMC5347684

30 Arcaroli J.J. Touban B.M. Tan A.C. Varella-Garcia M. Powell R.W. Eckhardt S.G. Elvin P. Gao D. Messersmith W.A. Gene array and fluorescence in situ hybridization biomarkers of activity of saracatinib (AZD0530), a Src inhibitor, in a preclinical model of colorectal cancer Clin. Cancer Res. 16 2010 4165 4177 10.1158/1078-0432.CCR-10-0066 20682712 PMC3805460

31 Tan F.H. Putoczki T.L. Lou J. Hinde E. Hollande F. Giraud J. Stylli S.S. Paradiso L. Zhu H.J. Sieber O.M. Luwor R.B. Ponatinib Inhibits Multiple Signaling Pathways Involved in STAT3 Signaling and Attenuates Colorectal Tumor Growth Cancers (Basel) 10 2018 526 10.3390/CANCERS10120526 PMC6316865 30572654

32 Powers M.V. Valenti M. Miranda S. Maloney A. Eccles S.A. Thomas G. Clarke P.A. Workman P. Mode of cell death induced by the HSP90 inhibitor 17-AAG (tanespimycin) is dependent on the expression of pro-apoptotic BAX Oncotarget 4 2013 1963 1975 10.18632/ONCOTARGET.1419 24185264 PMC3875762

33 Chen E. Jonker D. Gauthier I. MacLean M. Wells J. Powers J. Seymour L. Phase I study of cediranib in combination with oxaliplatin and infusional 5-Fluorouracil in patients with advanced colorectal cancer Clin. Cancer Res. 15 2009 1481 1486 10.1158/1078-0432.CCR-08-0761 19228749

34 Bliss C.I. the Toxicity of Poisons Applied Jointly Ann. Appl. Biol. 26 1939 585 615 10.1111/j.1744-7348.1939.tb06990.x

35 An S. Gunathilake M. Lee J. Kim M. Oh J.H. Chang H.J. Sohn D.K. Shin A. Kim J. Relationship Between Aspirin Use and Site-Specific Colorectal Cancer Risk Among Individuals With Metabolic Comorbidity J. Kor. Med. Sci. 39 2024 e199 10.3346/jkms.2024.39.e199 PMC11231443 38978486

36 Barrett T. Troup D.B. Wilhite S.E. Ledoux P. Rudnev D. Evangelista C. Kim I.F. Soboleva A. Tomashevsky M. Edgar R. NCBI GEO: Mining tens of millions of expression profiles - Database and tools update Nucleic Acids Res. 35 2007 D760 D765 10.1093/nar/gkl887 17099226 PMC1669752

37 Satoh K. Yachida S. Sugimoto M. Oshima M. Nakagawa T. Akamoto S. Tabata S. Saitoh K. Kato K. Sato S. Global metabolic reprogramming of colorectal cancer occurs at adenoma stage and is induced by MYC Proc. Natl. Acad. Sci. USA 114 2017 E7697 E7706 10.1073/PNAS.1710366114 28847964 PMC5604037

38 Hernández S. Bessa X. Beà S. Hernández L. Nadal A. Mallofré C. Muntane J. Castells A. Fernández P.L. Cardesa A. Differential Expression of cdc25 Cell-Cycle–Activating Phosphatases in Human Colorectal Carcinoma Lab. Invest. 81 2001 465 473 10.1038/labinvest.3780254 11304565

39 Stachel I. Geismann C. Aden K. Deisinger F. Rosenstiel P. Schreiber S. Sebens S. Arlt A. Schäfer H. Modulation of nuclear factor E2-related factor-2 (Nrf2) activation by the stress response gene immediate early response-3 (IER3) in colonic epithelial cells: a novel mechanism of cellular adaption to inflammatory stress J. Biol. Chem. 289 2014 1917 1929 10.1074/JBC.M113.490920 24311782 PMC3900942

40 Yang S. Ren X. Liang Y. Yan Y. Zhou Y. Hu J. Wang Z. Song F. Wang F. Liao W. KNK437 restricts the growth and metastasis of colorectal cancer via targeting DNAJA1/CDC45 axis Oncogene 39 2020 249 261 10.1038/s41388-019-0978-0 31477839

41 Fang Y. Yu H. Liang X. Xu J. Cai X. Chk1-induced CCNB1 overexpression promotes cell proliferation and tumor growth in human colorectal cancer Cancer Biol. Ther. 15 2014 1268 1279 10.4161/CBT.29691 24971465 PMC4128869

42 Wang Q.Q. Zhou Y.C. Zhou Ge Y.J. Qin G. Yin T.F. Zhao D.Y. Tan C. Yao S.K. Comprehensive proteomic signature and identification of CDKN2A as a promising prognostic biomarker and therapeutic target of colorectal cancer World J. Clin. Cases 10 2022 7686 7697 10.12998/WJCC.V10.I22.7686 36158487 PMC9372836

43 Bent A.H. Maru D.M. Vauthey J.-N. Dasari A. Johnson B. Kee B.K. Parseghian C.M. Menter D. Overman M.J. Morris V.K. HER3 expression in metastatic colorectal cancer: Defining the clinicomolecular profile of an emerging target J. Clin. Oncol. 40 2022 3588 10.1200/JCO.2022.40.16_SUPPL.3588

44 Jun S.Y. Kim J. Yoon N. Maeng L.S. Byun J.H. Prognostic Potential of Cyclin D1 Expression in Colorectal Cancer J. Clin. Med. 12 2023 572 10.3390/JCM12020572 36675501 PMC9867305

45 Al-Bahlani S.M. Al-Bulushi K.H. Al-Alawi Z.M. Al-Abri N.Y. Al-Hadidi Z.R. Al-Rawahi S.S. Cisplatin Induces Apoptosis Through the Endoplasmic Reticulum-mediated, Calpain 1 Pathway in Triple-negative Breast Cancer Cells Clin. Breast Cancer 17 2017 e103 e112 10.1016/J.CLBC.2016.12.001 28089626

46 Yang W. Liu L. Masugi Y. Qian Z.R. Nishihara R. Keum N. Wu K. Smith-Warner S. Ma Y. Nowak J.A. Calcium intake and risk of colorectal cancer according to expression status of calcium-sensing receptor (CASR) Gut 67 2018 1475 1483 10.1136/GUTJNL-2017-314163 28676564 PMC5754263

47 Li W.W. Wang H.Y. Nie X. Liu Y.B. Han M. Li B.H. Human colorectal cancer cells induce vascular smooth muscle cell apoptosis in an exocrine manner Oncotarget 8 2017 62049 62056 10.18632/ONCOTARGET.18893 28977925 PMC5617485

48 Pandurangan A.K. Ismail S. Esa N.M. Munusamy M.A. Inositol-6 phosphate inhibits the mTOR pathway and induces autophagy-mediated death in HT-29 colon cancer cells Arch. Med. Sci. 14 2018 1281 1288 10.5114/AOMS.2018.76935 30393482 PMC6209706

49 Pazienza V. Vinciguerra M. Mazzoccoli G. PPARs Signaling and Cancer in the Gastrointestinal System PPAR Res. 2012 2012 560846 10.1155/2012/560846 PMC3458283 23028383

50 Kim M.S. Ha S.E. Wu M. Zogg H. Ronkon C.F. Lee M.Y. Ro S. Extracellular Matrix Biomarkers in Colorectal Cancer Int. J. Mol. Sci. 22 2021 9185 10.3390/IJMS22179185 PMC8430714 34502094

51 Jung M.K. Okekunle A.P. Lee J.E. Sung M.K. Lim Y.J. Role of Branched-chain Amino Acid Metabolism in Tumor Development and Progression J. Cancer Prev. 26 2021 237 243 10.15430/JCP.2021.26.4.237 35047449 PMC8749315

52 Guo B. Han X. Tkach D. Huang S.G. Zhang D. AMPK promotes the survival of colorectal cancer stem cells Animal Model. Exp. Med. 1 2018 134 142 10.1002/AME2.12016 30891558 PMC6388072

53 Chen T. Liu N. Xu G.M. Liu T.J. Liu Y. Zhou Y. Huo S.B. Zhang K. Apelin13/APJ promotes proliferation of colon carcinoma by activating Notch3 signaling pathway Oncotarget 8 2017 101697 101706 10.18632/ONCOTARGET.21904 29254197 PMC5731907

54 Adler A.S. McCleland M.L. Yee S. Yaylaoglu M. Hussain S. Cosino E. Quinones G. Modrusan Z. Seshagiri S. Torres E. An integrative analysis of colon cancer identifies an essential function for PRPF6 in tumor growth Genes Dev. 28 2014 1068 1084 10.1101/GAD.237206.113 24788092 PMC4035536

55 Cancer Genome Atlas Research Network Weinstein J.N. Collisson E.A. Mills G.B. Shaw K.R.M. Ozenberger B.A. Ellrott K. Shmulevich I. Sander C. Stuart J.M. The Cancer Genome Atlas Pan-Cancer Analysis Project Nat. Genet. 45 2013 1113 1120 10.1038/NG.2764 24071849 PMC3919969

56 Marisa L. de Reyniès A. Duval A. Selves J. Gaub M.P. Vescovo L. Etienne-Grimaldi M.C. Schiappa R. Guenot D. Ayadi M. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value PLoS Med. 10 2013 1001453 10.1371/JOURNAL.PMED.1001453 PMC3660251 23700391

57 Pelka K. Hofree M. Chen J.H. Sarkizova S. Pirl J.D. Jorgji V. Bejnood A. Dionne D. Ge W.H. Xu K.H. Spatially organized multicellular immune hubs in human colorectal cancer Cell 184 2021 4734 4752.e20 10.1016/J.CELL.2021.08.003/ATTACHMENT/6A06B6DB-8D91-4DCB-A53A-A74585AE0547/MMC7.XLSX 34450029 PMC8772395

58 Loewe S. Muischnek H. Effect of combinations: Mathematical basis of the problem Arch. f. Exp. Pathol. u. Pharmakol. 114 1926 313 326 10.1007/BF01952257

59 Berenbaum M.C. What is synergy? Pharmacol. Rev. 41 1989 93 141 10.1016/S0031-6997(25)00026-2 2692037

60 Wooten D.J. Meyer C.T. Lubbock A.L.R. Quaranta V. Lopez C.F. MuSyC is a consensus framework that unifies multi-drug synergy metrics for combinatorial drug discovery Nat. Commun. 12 2021 4607 10.1038/S41467-021-24789-Z PMC8322415 34326325

61 Ianevski A. Giri A.K. Aittokallio T. SynergyFinder 2.0: visual analytics of multi-drug combination synergies Nucleic Acids Res. 48 2020 W488 W493 10.1093/NAR/GKAA216 32246720 PMC7319457

62 Liu H. Zhang W. Zou B. Wang J. Deng Y. Deng L. DrugCombDB: a comprehensive database of drug combinations toward the discovery of combinatorial therapy Nucleic Acids Res. 48 2020 D871 D881 10.1093/NAR/GKZ1007 31665429 PMC7145671

63 Subramanian A. Narayan R. Corsello S.M. Peck D.D. Natoli T.E. Lu X. Gould J. Davis J.F. Tubelli A.A. Asiedu J.K. A Next Generation Connectivity Map: L1000 platform and the first 1,000,000 profiles Cell 171 2017 1437 1452.e17 10.1016/j.cell.2017.10.049 29195078 PMC5990023

64 Goldman M.J. Craft B. Hastie M. Repečka K. McDade F. Kamath A. Banerjee A. Luo Y. Rogers D. Brooks A.N. Zhu J. Haussler D. Visualizing and interpreting cancer genomics data via the Xena platform Nat. Biotechnol. 38 2020 675 678 10.1038/s41587-020-0546-8 32444850 PMC7386072

65 Wang S. Liu X. The UCSCXenaTools R package: a toolkit for accessing genomics data from UCSC Xena platform, from cancer multi-omics to single-cell RNA-seq J. Open Source Softw. 4 2019 1627 10.21105/JOSS.01627

66 Lonsdale J. Thomas J. Salvatore M. Phillips R. Lo E. Shad S. Hasz R. Walters G. Garcia F. Young N. The Genotype-Tissue Expression (GTEx) project Nat. Genet. 45 2013 580 585 10.1038/ng.2653 23715323 PMC4010069

67 Ritchie M.E. Phipson B. Wu D. Hu Y. Law C.W. Shi W. Smyth G.K. limma powers differential expression analyses for RNA-sequencing and microarray studies Nucleic Acids Res. 43 2015 e47 10.1093/NAR/GKV007 25605792 PMC4402510

68 Iwata M. Yuan L. Zhao Q. Tabei Y. Berenger F. Sawada R. Akiyoshi S. Hamano M. Yamanishi Y. Predicting drug-induced transcriptome responses of a wide range of human cell lines by a novel tensor-train decomposition algorithm Bioinformatics 35 2019 i191 i199 10.1093/bioinformatics/btz313 31510663 PMC6612872

69 Gaulton A. Bellis L.J. Bento A.P. Chambers J. Davies M. Hersey A. Light Y. McGlinchey S. Michalovich D. Al-Lazikani B. Overington J.P. ChEMBL: a large-scale bioactivity database for drug discovery Nucleic Acids Res. 40 2012 D1100 D1107 10.1093/NAR/GKR777 21948594 PMC3245175

70 Günther S. Kuhn M. Dunkel M. Campillos M. Senger C. Petsalaki E. Ahmed J. Urdiales E.G. Gewiess A. Jensen L.J. SuperTarget and Matador: resources for exploring drug-target relationships Nucleic Acids Res. 36 2008 D919 D922 10.1093/NAR/GKM862 17942422 PMC2238858

71 Knox C. Law V. Jewison T. Liu P. Ly S. Frolkis A. Pon A. Banco K. Mak C. Neveu V. DrugBank 3.0: a comprehensive resource for ‘Omics’ research on drugs Nucleic Acids Res. 39 2011 D1035 D1041 10.1093/NAR/GKQ1126 21059682 PMC3013709

72 Roth B.L. Lopez E. Patel S. Kroeze W.K. The Multiplicity of Serotonin Receptors Neuroscientist 6 2000 252 262 10.1177/107385840000600408

73 Kanehisa M. Goto S. Furumichi M. Tanabe M. Hirakawa M. KEGG for representation and analysis of molecular networks involving diseases and drugs Nucleic Acids Res. 38 2010 D355 D360 10.1093/nar/gkp896 19880382 PMC2808910

74 Liu T. Lin Y. Wen X. Jorissen R.N. Gilson M.K. BindingDB: a web-accessible database of experimentally determined protein–ligand binding affinities Nucleic Acids Res. 35 2007 D198 D201 10.1093/NAR/GKL999 17145705 PMC1751547

75 Qin C. Zhang C. Zhu F. Xu F. Chen S.Y. Zhang P. Li Y.H. Yang S.Y. Wei Y.Q. Tao L. Chen Y.Z. Therapeutic target database update 2014: a resource for targeted therapeutics Nucleic Acids Res. 42 2014 D1118 D1123 10.1093/NAR/GKT1129 24265219 PMC3964951

76 Sawada R. Iwata H. Mizutani S. Yamanishi Y. Target-Based Drug Repositioning Using Large-Scale Chemical-Protein Interactome Data J. Chem. Inf. Model. 55 2015 2717 2730 10.1021/ACS.JCIM.5B00330/SUPPL_FILE/CI5B00330_SI_001.PDF 26580494

77 Luck K. Kim D.K. Lambourne L. Spirohn K. Begg B.E. Bian W. Brignall R. Cafarelli T. Campos-Laborie F.J. Charloteaux B. A reference map of the human binary protein interactome Nature 580 2020 402 408 10.1038/s41586-020-2188-x 32296183 PMC7169983

78 Ruepp A. Brauner B. Dunger-Kaltenbach I. Frishman G. Montrone C. Stransky M. Waegele B. Schmidt T. Doudieu O.N. Stümpflen V. Mewes H.W. CORUM: the comprehensive resource of mammalian protein complexes Nucleic Acids Res. 36 2008 D646 D650 10.1093/NAR/GKM936 17965090 PMC2238909

79 Hornbeck P.V. Kornhauser J.M. Tkachev S. Zhang B. Skrzypek E. Murray B. Latham V. Sullivan M. PhosphoSitePlus: a comprehensive resource for investigating the structure and function of experimentally determined post-translational modifications in man and mouse Nucleic Acids Res. 40 2012 D261 D270 10.1093/NAR/GKR1122 22135298 PMC3245126

80 Shimizu Y. Hattori M. Goto S. Kanehisa M. Generalized reaction patterns for prediction of unknown enzymatic reactions Genome Inform. 20 2008 149 158 10.1142/9781848163003_0013 19425130

81 Fazekas D. Koltai M. Türei D. Módos D. Pálfy M. Dúl Z. Zsákai L. Szalay-Bekő M. Lenti K. Farkas I.J. SignaLink 2 - a signaling pathway resource with multi-layered regulatory networks BMC Syst. Biol. 7 2013 7 10.1186/1752-0509-7-7/TABLES/2 23331499 PMC3599410

82 Breuer K. Foroushani A.K. Laird M.R. Chen C. Sribnaia A. Lo R. Winsor G.L. Hancock R.E.W. Brinkman F.S.L. Lynn D.J. InnateDB: systems biology of innate immunity and beyond--recent updates and continuing curation Nucleic Acids Res. 41 2013 D1228 D1233 10.1093/nar/gks1147 23180781 PMC3531080

83 Meyer M.J. Das J. Wang X. Yu H. INstruct: a database of high-quality 3D structurally resolved protein interactome networks Bioinformatics 29 2013 1577 1579 10.1093/bioinformatics/btt181 23599502 PMC3673217

84 Zhang J.D. Wiemann S. KEGGgraph: a graph approach to KEGG PATHWAY in R and bioconductor Bioinformatics 25 2009 1470 1471 10.1093/BIOINFORMATICS/BTP167 19307239 PMC2682514

85 Everett M.G. Borgatti S.P. The centrality of groups and classes J. Math. Sociol. 23 1999 181 201 10.1080/0022250X.1999.9990219

86 Jassal B. Matthews L. Viteri G. Gong C. Lorente P. Fabregat A. Sidiropoulos K. Cook J. Gillespie M. Haw R. The reactome pathway knowledgebase Nucleic Acids Res. 48 2020 D498 D503 10.1093/NAR/GKZ1031 31691815 PMC7145712

87 Liu Y.Y. Slotine J.J. Barabási A.L. Controllability of complex networks Nature 473 2011 167 173 10.1038/nature10011 21562557

88 Liu Y.Y. Barabási A.L. Control principles of complex systems Rev. Mod. Phys. 88 2016 035006 10.1103/REVMODPHYS.88.035006/FIGURES/41/MEDIUM

89 Qiao M. Abbas J.J. Jindrich D.L. Jones D.O. Scolnic D.M. Foley R.J. Nacher J.C. Akutsu T. Analysis on controlling complex networks based on dominating sets J. Phys., Conf. Ser. 410 2013 012104 10.1088/1742-6596/410/1/012104

90 Nguyen M.H. Hà M.H. Nguyen D.N. Tran T.T. Solving the k-dominating set problem on very large-scale networks Comput. Soc. Networks 7 2020 1 15 10.1186/S40649-020-00078-5/TABLES/4

91 Pan J.Y. Yang H.J. Faloutsos C. Duygulu P. Automatic multimedia cross-modal correlation discovery KDD-2004 - Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining 2004 Association for Computing Machinery (ACM) 653 658 10.1145/1014052.1014135

92 Hagberg A. Schult D. Swart P. Exploring Network Structure, Dynamics, and Function Using NetworkX Proceedings of the 7th Python Science Conference 2008 11 16 10.25080/TCWV9851

93 AuerPeter Using confidence bounds for exploitation-exploration trade-offs J. Mach. Learn. Res. 3 2003 397 422 10.5555/944919.944941

94 Piñero J. Saüch J. Sanz F. Furlong L.I. The DisGeNET cytoscape app: Exploring and visualizing disease genomics data Comput. Struct. Biotechnol. J. 19 2021 2960 2967 10.1016/J.CSBJ.2021.05.015 34136095 PMC8163863

95 Amberger J.S. Bocchini C.A. Schiettecatte F. Scott A.F. Hamosh A. OMIM.org: Online Mendelian Inheritance in Man (OMIM®), an online catalog of human genes and genetic disorders Nucleic Acids Res. 43 2015 D789 D798 10.1093/NAR/GKU1205 25428349 PMC4383985

96 Landrum M.J. Lee J.M. Benson M. Brown G.R. Chao C. Chitipiralla S. Gu B. Hart J. Hoffman D. Jang W. ClinVar: improving access to variant interpretations and supporting evidence Nucleic Acids Res. 46 2018 D1062 D1067 10.1093/NAR/GKX1153 29165669 PMC5753237

97 Buniello A. Macarthur J.A.L. Cerezo M. Harris L.W. Hayhurst J. Malangone C. McMahon A. Morales J. Mountjoy E. Sollis E. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019 Nucleic Acids Res. 47 2019 D1005 D1012 10.1093/nar/gky1120 30445434 PMC6323933

98 Li M.J. Wang P. Liu X. Lim E.L. Wang Z. Yeager M. Wong M.P. Sham P.C. Chanock S.J. Wang J. GWASdb: a database for human genetic variants identified by genome-wide association studies Nucleic Acids Res. 40 2012 D1047 D1054 10.1093/NAR/GKR1182 22139925 PMC3245026

99 Bastarache L. Using Phecodes for Research with the Electronic Health Record: From PheWAS to PheRS Annu. Rev. Biomed. Data Sci. 4 2021 1 19 10.1146/ANNUREV-BIODATASCI-122320-112352 34465180 PMC9307256

100 Papadakis M.A. McPhee S.J. Rabow M.W. Current Medical Diagnosis and Treatment 2014 2014 McGraw Hill Medical

101 Therneau T.M. A Package for Survival Analysis in R. R Foundation for Statistical Computing 2023 https://CRAN.R-project.org/package=survival

Supplemental information

Document S1. Figures S1–S17, Tables S3, S4, and S8–S13 Table S1. DEGPs and therapeutic pathways estimated for diseases, related to Figure 2A and B DEGPs were identified from disease transcriptome data by weighted gene co-expression network analysis (WGCNA) 15 and pathway analysis. Therapeutic pathways regulating the DEGPs were estimated based on network topology of pathway–pathway interaction network. The “disease_id” and “disease_name” columns correspond to KEGG DISEASE IDs and disease names, respectively. The “pathway_type” column represents whether a pathway is DEGP or therapeutic pathway. The “pathway_id” and “pathway_name” columns correspond to REACTOME pathway IDs and pathway names, respectively. Table S2. Estimated therapeutic targets for various diseases by network controllability analysis, related to Figure 2C–J Gene set column represents the type of estimated target set. Table S5. Predicted boost-type drug combinations for various diseases, related to Figure 4A Top 20 drug combinations predicted for each disease by proposed method with estimated targets are listed. “Prediction rank” represents prediction rank based on prediction score. The “drug_id1” and “drug_id2” columns correspond to KEGG DRUG IDs of predicted drug combinations. The “drug_name1” and “drug_name2” columns correspond to drug names of predicted drug combinations. The “drug_efficacy1” and “drug_efficacy2” correspond to drug efficacy of drug combinations. The “score” column corresponds to prediction scores representing therapeutic druggability. Table S6. Predicted complement-type drug combinations for various diseases, related to Figure 4A Top 20 drug combinations predicted for each disease by proposed method with estimated targets are listed. “Prediction rank” represents prediction rank based on prediction score. The “drug_id1” and “drug_id2” columns correspond to KEGG DRUG IDs of predicted drug combinations. The “drug_name1” and “drug_name2” columns correspond to drug names of predicted drug combinations. The “drug_efficacy1” and “drug_efficacy2” correspond to drug efficacy of drug combinations. The “score” column corresponds to prediction scores representing therapeutic druggability. Table S7. Predicted indistinct-type drug combinations for various diseases, related to Figure 4A Top 20 drug combinations predicted for each disease by proposed method with estimated targets are listed. “Prediction rank” represents prediction rank based on prediction score. The “drug_id1” and “drug_id2” columns correspond to KEGG DRUG IDs of predicted drug combinations. The “drug_name1” and “drug_name2” columns correspond to drug names of predicted drug combinations. The “drug_efficacy1” and “drug_efficacy2” correspond to drug efficacy of drug combinations. The “score” column corresponds to prediction scores representing therapeutic druggability. Table S14. Differential expression analysis using CRC disease-specific transcriptome data, related to Figure 5A–D Differential expression analysis was performed between healthy controls and patients with CRC using CRC disease-specific transcriptome data. Table S15. All pathways significantly enriched for upregulated DEGs in singlets, doublet, and triplet, related to Figure 5I Exact −log(p-value) values obtained from pathway enrichment analysis of upregulated genes are provided. Table S16. All pathways significantly enriched for downregulated DEGs in singlets, doublet, and triplet, related to Figure 5J Exact −log(p-value) values obtained from pathway enrichment analysis of downregulated genes are provided. Table S17. Enriched pathways important for CRC treatment using upregulated DEGs in the singlets, doublets, and triplet, related to Figure 5O Exact −log(p-value) values obtained from pathway enrichment analysis of upregulated genes are provided. Table S18. Enriched pathways important for CRC treatment using upregulated DEGs in the singlets, doublets, and triplet, related to Figure 5P Exact −log(p-value) values obtained from pathway enrichment analysis of downregulated genes are provided.

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2026.115339 .