The entire workflow of this study is depicted in Fig. 1. The drug responses of 625 chemical compounds tested on 199 bone marrow samples from patients with AML were obtained from the FIMM AML data set [11]. The bipartite network was constructed using this data set, as explained in the materials and methods section. A bipartite network can be projected onto two different types of unipartite networks containing nodes of only one type. The projection of the bipartite network, onto the “drug” node set is considered here, called the drug similarity network. The Louvain community detection approach was used to find drugs that behaved similarly in terms of drug response [14]. The results gave us two communities (clusters) of drugs denoted by C1 and C2 with network sizes of 155 and 141, respectively (Table S1).

Fig. 1: Schematic outline of the study.

Data pre-processing began after data collection, which was followed by full matrix extraction, weighted bipartite network reconstruction, and computational validation. After the selection of the best combinations, bone marrow and peripheral blood samples from both healthy individuals (n = 5) and AML patients (n = 16) were subjected to drug sensitivity assessment. For ATP-based viability assay the study design contains 8 drugs and 28 combinations in 384-well plates, each drug with 5 different concentrations and two replicates. The single cell sensitivity assay using the iQue® Screener PLUS flow cytometer was performed in 384-well plates to monitor drug effects on cell sub-types. The study design contains 5 drugs and 3 combinations, all with two replicates and five concentrations. For sapanisertib, the drug concentrations are 0.1, 1, 10, 100, and 1000 nM, and for all other drugs are 1, 10, 100, 1000, and 10,000 nM.

Comparing AML drug clusters: evaluating protein target pathways and chemical structure similarity

We used two independent computational methods to determine how distinct the two clusters are: the first identifies the significant difference between biological pathways of drug protein targets in each cluster, and the second evaluates the chemical structure similarity of drugs in each cluster. We constructed a drug-target network using the drug target commons (DTC) database [15], which is also a bipartite network in which each link connects drugs to their protein targets. Let T1 and T2 represent the set of protein targets of drugs in the cluster C1 and C2, respectively, and T represents the union of T1and T2. In this study |\(_\)|\(=921,_|=842,=1055.\)

Proteins with a score S (explained in the methods) greater than log (2) are considered to be preferentially targeted by drugs in cluster 1, denoted by PPT1. Similarly, PPT2 proteins have a score of less than \(\log (0.50)\). We performed GSEA (gene set enrichment analysis) on PPT1 and PPT2 proteins based on their associated scoring functions. As expected, the biological processes and signaling pathways affected by drugs in Clusters 1 and 2 are distinct. This difference enables us to inhibit two different signaling pathways using one combination. Drugs in cluster 1 (PPT1), such as LY3009120 (a pan-RAF inhibitor), predominantly target proteins associated with the RAF-MEK-ERK signaling pathway. This pathway plays a crucial role in cell proliferation and growth, indirectly influencing processes like cell-substrate adhesion and ion trans-membrane transport, which are enriched in our analysis [27]. In contrast, JAK1/2 inhibitors like ruxolitinib target JAK proteins, involved in cytokine signaling and immune responses, impacting pathways related to neuroactive ligand-receptor interactions and the regulation of actin cytoskeleton [28]. Drugs like birabresib, which target proteins in the bromodomain and extra-terminal (BET) family, have a role in gene regulation through chromatin binding, affecting gene expression and pathways related to chemical reactions and collagen metabolism [29]. Plicamycin, which binds to guanine-cytosine-rich regions of DNA, may influence gene expression and regulation, impacting pathways related to collagen metabolism and other DNA-dependent processes (Fig. 2A) [30]. On the other hand, proteins targeted by drugs in cluster 2 (PPT2) (silmitasertib, ulixertinib, sapanisertib, and teniposide) are in the p53 signaling pathway, cell cycle, apoptosis, and pancreatic, colorectal and chronic myeloid leukemia cancers and related to tumorigenesis and progression pathways, including human immunodeficiency virus 1 infection [31,32,33,34].

Fig. 2: Gene Enrichment Analysis for Proteins in Clusters 1 and 2.

Sankey plot of enriched (A) KEGG signalling pathways and (B) GO biological processes related to target protein clusters PPT1 and PPT2. Each rectangle on the right side represents a pathway or biological process, and the size of each rectangle illustrates the degree of connectivity of each pathway. Each biological process or pathway is represented by a unique color. GO and KEGG pathway enrichment analysis on proteins that are merely targets by drugs in one cluster. G1 (G2) includes proteins that are targeted by at least three drugs in cluster 1 (cluster 2) (155 and 141 drugs). C Biological processes (BPs) of G1, (D) Biological processes (BPs) of G2, (E) KEGG pathway related to G1 proteins, and (F) KEGG pathway related to G2 proteins. The size of the node corresponds to number of genes, the x-axis is Fold Enrichment and the color of bars indicates the negative logarithm of Fold Enrichment.

We also performed ShinyGO [17] Gene Ontology and KEGG pathway enrichment analysis on proteins that are merely targeted by drugs in one cluster. For this purpose, two protein sets G1 and G2 were selected such that G1 includes proteins targeted by at least three drugs in cluster 1 and at most two drugs in cluster 2, and similarly, G2, consist of proteins that are mostly targeted by drugs in cluster 2 (at least three drugs in cluster 2 and at most two drugs in cluster 1). REVIGO was also used to summarize the enriched GO terms, and the results are shown in Fig. 2 and Tables S2 and S3. The cAMP signaling pathway, lipids and atherosclerosis, steroid hormone biosynthesis, and rhythmic processes and circadian rhythm are biological processes related to G1 proteins, which are mostly targeted by LY3009120, birabresib, plicamycin, and ruxolitinib. Cell cycle, cellular senescence, T-cell leukemia virus 1 infection and cell division, mitotic cell cycle, and protein phosphorylation processes are related to G2 proteins, mostly targeted by silmitasertib, ulixertinib, sapanisertib, and teniposide. Therefore, we demonstrate that the protein targets of drugs in each cluster are involved in distinct pathways and biological processes.

To do homogeneity analysis of chemical structure of drugs, the dice similarity test was used to show how structurally similar the drugs are in each cluster. This measurement compares the number of chemical features shared by a pair of compounds to the average size of the total number of features present. Pairwise similarities were calculated for chemical compounds chosen from two drug clusters for inter-cluster comparison. Drugs from different clusters are less similar than drugs from the same cluster, as shown in Fig. 3A. According to the box plot, the inter-cluster similarities are less than the intra-cluster similarities in both clusters. The results of the t-test imply that the mean of inter-cluster similarities is less than the mean of intra-cluster similarities in clusters 1 and 2 (p-value < 2.2e-16 for both t-test).

Fig. 3: Comparative analysis of dice similarity and drug efficacy-toxicity profiles in AML therapy.

A Box plot of dice similarity coefficient indices comparing intra-cluster 1 and 2 to inter-cluster compound pairs. P-value is generated using Wilcoxon signed-rank test, shown in red color. B The toxicity and efficacy of 296 drugs. Inset plot shows the relationship between toxicity and efficacy. Top five percent of drugs whose toxicity is less than the average of all drug toxicity and whose efficacy is greater than the average of all drug efficacy are in blue, and their name is shown in rectangle labels.

Combination selection: balancing toxicity and efficacy across clusters

As a result, we demonstrated that clusters are well-separated and that the protein targets of drugs in each cluster are involved in distinct pathways. In this novel combination strategy, we aim to select two drugs from distinct clusters while taking both toxicity and efficacy into account. The optimal combinations are those that have lower toxicity than the average toxicity and higher efficacy values than the average efficacy value for all drugs. For each drug, the average drug response of healthy and AML patient samples in the data set are considered as toxicity and efficacy, respectively. We assume that the ideal drugs have no inhibitory effect on healthy samples but significantly influence blast cells in AML patient samples. We chose the top 5% of drugs whose toxicity is less than the average of all drug toxicity and efficacy is greater than the average of all drug efficacy. Figure 3B depicts the link between toxicity and efficacy values of 296 drugs on 81 samples. The top four selected small molecules in each cluster are summarized in Table 1 and Table S4. Four chemical compounds from cluster 1 including birabresib, LY3009120, plicamycin, and ruxolitinib as well as four drugs from cluster 2 including sapanisertib, silmitasertib, teniposide, and ulixertinib were chosen for drug combination testing. According to our experimental design, the combination of drugs within a single cluster is known as negative group or intra-cluster, and the combination of drugs between clusters is considered as positive group or inter-cluster.

Table 1 The selected chemical compounds from two clusters of drugs in the drug similarity network.Enhanced efficacy and reduced toxicity in inter-cluster drug combinations on AML patient samples revealed by cell viability drug screening

In the testing of all 16 inter-cluster and 12 intra-cluster combinations at five different concentrations, the cell viability of 16 samples from AML patients and 5 samples from healthy donors were monitored. Patient samples with blast percentage more than 49% were chosen for testing with the CTG assay (Table S5). The average inhibition across dosages on 16 patient samples is regarded as efficacy, whereas the average inhibition across dosages on healthy samples is regarded as toxicity. The drug combinations with rectangular labels have higher efficacy and lower toxicity than the median. The proportion test (p-value = 0.006) revealed that the percentages of inter-cluster drug combinations with high efficacy (efficacy higher than the third quantile of efficacy values) and low toxicity (toxicities lower than the first quantile of toxicities) are significantly more than random choices.

The synergy and combination ratio (CR) of drug combinations on AML and healthy samples was then calculated using synergy scoring functions HSA [21], Bliss [22], Loewe [23], and ZIP [24] (Figs. 4 and S3). The same analysis was done on synergy scoring values, and it was discovered that inter-cluster drug combinations differ significantly from random choices (P-values shown in Fig. 4A-F). The drug combinations shown with rectangular labels have the highest synergy on AML patient samples, and the lowest synergy on healthy samples. Table 2 summarizes all six plots and the significant drug combinations according to different measures are highlighted by green (inter-cluster), yellow and purple (intra-clusters). Following CTG analysis, consensus across synergy scoring functions led to the selection of the five best drug combinations out of 28 to quantify blast-specific drug responses with flow cytometry. Additionally, we used one of the most extensive databases, the Probes & medications portal (PDP) dataset [35], to extract the protein targets of these selected drugs. Table S6 provides a summary of the hypergeometric test findings, which show that there is no discernible overlap between the protein targets of these drugs whether taken separately or in combination. The need for future work arises to assess relevant biomarkers of on-target activity for each single and combination approach.

Fig. 4: Drug combinations’ synergy scores on 16 AML samples and 5 healthy samples.

The X-axis depicts the synergy in AML samples and the Y-axis represents the synergy in healthy samples. The median inhibition on AML and healthy samples is shown by dashed lines in red and blue, respectively. There are three groupings: clusters 1, 2, and intercluster, and the color of each dot indicates each of these groups. The p values presented in each panel are associated with the proportion test, comparing the inter-cluster combination with the random selection of drugs. The average of inhibition of drug combinations on dosages (A) and several synergy scores were depicted in separate panels using synergy scoring functions ZIP (B), HSA (C), Bliss (D), Loewe (E), and combination ratio (CR) of drug combinations on AML and healthy samples (F).

Table 2 Selected drug combinations sorted by synergy scoring functions.Cell subtype viability analysis highlights low toxicity of selected combinations

Using the CTG assay, we measure the general BM-MNC sensitivity, whereas with flow cytometry analysis we measure the number of live cells among different cell populations. Following 72-hour treatment with the 5 selected combinations on 3 different samples, viability of different cell subtypes of interest was measured by flow cytometry. Sample selection was based on the inclusion of three biological replicates for each combination, considering available cell numbers to enhance statistical power and result reliability. For each sample, there is a specific plate layout which can be found in Supplementary Fig. S1. We used six cell surface markers (CD14, CD15, CD45, CD38, CD117, and CD34; Table S7) to identify the major leukocyte populations present in the AML BM-MNCs: monoblasts, myelocytes, leukemic blasts, leukemic stem cells, and myeloid progenitor cells (Fig. S2).

In the studied samples, the average of blasts out of CD45 positive leukocytes, was 70% in DMSO, while on average 36% ± 16% of the blasts were killed by the combinations (Table S8). Based on the results, the percentage of dead cells for all five combinations in lymphocytes is considerably lower than 25% (Fig. 5). More importantly selected combinations have lower synergistic effect on lymphocytes compared to the blast population, demonstrating the lower toxicity of combinations (Figs. 6 and S3). The combination of JAK1/2 inhibitor (ruxolitinib) with either ERK or CSNK2A1 inhibitor had the highest efficacy and lowest toxicity, demonstrating the important role of these targets in AML. Numerous studies show the significance of the JAK/STAT signaling system in determining how hematopoietic cells react to various cytokines and growth factors [36, 37]. Recently there has been increased interest in different drug combinations with ruxolitinib [38,39,40,41] and as our results show the combinations of this drug, by having the lowest toxicity, seem to be promising for AML treatment.

Fig. 5: The cell viability assay (CTG) and response of different cell populations to 5 selected combinations using flow cytometry assay.

Response signifies the percentage of dead cells following a 72 h treatment. The number of cells in each well was counted and normalized by the min–max normalization method. For each combination, three different samples, distinguished by the color of points, were treated with three different doses (10, 100, and 1000 nM), which are illustrated by the different point shapes. The colors in each cell group facet corresponds to a specific drug combination.

Fig. 6: Characteristics of drug responses and correlation analysis in AML treatment: Comprehensive flow cytometry and CTG assessment.

A A heat map showing characteristics of single agent and combination responses measured by flow cytometry readout. Blast-specific response of single drugs is highlighted according to drug sensitivity score (DSS) values with dark blue corresponding to high DSS value and white to low DSS value. Blast-specific and lymphocyte-specific response combinations at 1000 nM are highlighted according to percentage of apoptotic/dead cells, with dark blue in blast and red for lymphocyte corresponding to high percentage and white to a low percentage of apoptotic/dead cells. The synergistic effect of the drug combination was assessed based on the HSA synergistic score in 1000 nM on blast cells shown in blue and lymphocytes shown in red. B The correlation between responses measured by CTG and flow cytometry on five single drugs ruxolitinib, silmitasertib, ulixertinib, LY3009120, and sapanisertib, and (C) five drug combinations sapanisertib-LY3009120, ulixertinib-ruxolitinib, silmitasertib-ruxolitinib, silmitasertib-LY3009120, and ulixertinib-LY300912.

Blast-specific drug responses in AML: Efficacy profiles of selected combinations

We were able to assess blast-specific drug combination responses and compare them to the other combinations within different samples. Among the five tested combinations, two combinations with ruxolitinib which targets JAK1/2 were among the most efficient combinations. The combination of ruxolitinib with ulixertinib, an ERK inhibitor, exhibits the strongest efficacy against blasts, according to the results. After treatment, the combination induced 47% ± 13% cell death in blasts (Fig. 5 and Table S8) with a more synergistic effect on the blast population compared to the lymphocyte population (Fig. 6A). We depicted the gating of 1000 nM concentration of each drug on sample AML_3 to better understand the impact of combination therapy vs. DMSO control and single drug treated samples in Fig. 7. The number of blast cells in the ruxolitinib and ulixertinib treated well was reduced to 37%, showing the largest reduction compared to all other treatments, as shown in Fig. 7A. The second combination of ruxolitinib and silmitasertib, a CSNK2A1 inhibitor, showed high efficacy on blasts. On average, this combination induced death to almost half ± 14% of the blast population but had less effect on lymphocytes (Fig. 5 and Table S8). Additionally, this combination had a substantially higher inhibition rate compared to each single drug and acted synergistically toward the blast population (Fig. 6A).

Fig. 7: Flow cytometry scatter plots showing the effects of drug combinations on cell populations, along with comparisons to DMSO and single drug treated samples.

A This figure illustrates the effects of ruxolitinib and ulixertinib combination and (B) LY3009120 and sapanisertib combination on blasts, monocytic cells (CD14+) and lymphocytes after 72 h drug treatment. Numbers represent the percentage of cell counts in each population in comparison with untreated control. The plot represents a concentration of 1000 nM on sample AML_3.

Given the importance of pan-RAF inhibition, we next examined LY3009120 in combination with three other drugs. The samples used for the combination of LY3009120 and sapanisertib (mTOR1/2 inhibitor), consist of 56% blast and the response for them is 40% ± 12% inhibition. To confirm that this combination is efficient, we analyzed the effect of LY3009120 and sapanisertib combination with single treated and DMSO-treated cells in AML_3. In the combination-treated sample, the blast cells were significantly reduced to 13% while in the individual drugs LY3009120 and sapanisertib reduced the blasts to 38% and 75%, respectively (Fig. 6B). These results indicate that this combination has substantially higher inhibition rate compared to each single drug and a greater synergistic effect on blasts than on lymphocytes (Fig. 6A). Ulixertinib (ERK inhibitor) is the second drug that was used in combination with pan-RAF inhibitor. Patient samples treated with this combination, on average, contained 60% blast cells and after treatment they are reduced to 28% ± 14%. Finally, we tested the combination of LY3009120 with silmitasertib, a CSNK2A1 inhibitor on three different samples. The average blast population for these three samples is 62% and the response was 21% ± 5%. Overall, as shown in Figs. 5 and 6A, all combinations have very little impact on the lymphocyte populations, demonstrating low toxicity, and significantly more impact on less differentiated malignant cells, demonstrating the efficacy of the combinations.

Increased sensitivity of AML samples to combination therapies over single drugs, regardless of genetic mutations and prognosis categories

There is a significant correlation between CTG assay and blast specific results, indicating that reduction in cell number measured by CTG, is related to the malignant cell populations (Figs. 6B and 6C). The cell viability readout for a single drug is converted to a drug sensitivity score (DSS) which is a drug sensitivity metric based on area under the dose-response curve. A greater DSS indicates higher sensitivity [42]. Strikingly, by combining selected inter-cluster drugs, the blasts were targeted, and combinations showed a synergistic effect on this population (Fig. 6A). Considering the most prevalent mutations among AML patients [43, 44], we examined the existing mutations in selected samples to monitor the drug responses based on genetic changes (Fig. 6A). To evaluate the impact of the combinations on samples bearing genetic alterations, some mutations that are frequently found in AML patients were considered (Fig. 6A). Mutation to FLT3, a well-known driver gene in AML was represented in two samples. Other prevalent mutations occurred in NPM1, GATA2, DNMT3A, TET2, KMT2A, NRAS, SMC3, and SRSF2. The combinations induced a synergistic effect on the blast population, regardless of the genetic alterations. The European Leukemia Network (ELN) classifies patients into three prognosis categories: “favorable”, “intermediate”, or “adverse” [45]. AML patients are also classified using the French-American-British (FAB) classification [46], which is based on morphological features. Regardless of sample type, we observed a synergistic effect following treatment. Importantly, after therapy, we noticed a synergistic effect in all samples, indicating that these combinations are effective at combating the heterogeneity of AML. It has been demonstrated that drugs should target the less differentiated leukemic blasts to achieve the best response in patients [6]. Given these two observations—the presence of the most relevant mutations and the prevalence of blast cells in the samples— the combinations seem to be promising for treatment.

Efficacy and toxicity of the novel combinations compared to first-line treatment in AML

In the following analysis, we compared the proposed combinations in this study (ruxolitinib-ulixertinib and LY3009120-sapanisertib) with two FDA-approved combinations for AML (venetoclax-azacitidine and venetoclax-cytarabine), as well as the investigational combination of venetoclax-ruxolitinib. As illustrated in Fig. 8, venetoclax-ruxolitinib demonstrates the highest efficacy on both blast cells and lymphocytes compared to the other combinations. This dual efficacy profile is a noteworthy advantage; however, it comes at the cost of heightened toxicity, as indicated by our results.

Fig. 8: Flow cytometry assay of selected combinations compared to first-line AML combinations.

Response signifies the percentage of dead cells following 72 h treatment. The count of cells in each well was adjusted relative to the count in control wells featuring both positive (DMSO) and negative (BzCl) controls using the min–max normalization method. Each combination has been tested on different samples at 50 nM concentration for venetoclax and 1000 nM for the other drugs. Red asterisks define the average response for each combination and colored dots represent different samples. Each panel also represents six p-values resulting from the Wilcoxon signed-rank tests to compare statistically two proposed combinations with three other combinations (including two first-line treatments and one investigational combination) for AML.

Conversely, the novel combinations, ruxolitinib-ulixertinib and LY3009120-sapanisertib, showed comparable efficacy in targeting blast populations as the established combinations. Notably, there was no significant difference in terms of efficacy (p-values are shown in Fig. 8). However, these two combinations have a significant advantage in demonstrating lower toxicity compared to first-line combinations, particularly for lymphocytes. The effects of ruxolitinib-ulixertinib and LY3009120-sapanisertib on blast lymphocyte population were significantly lower than all other combinations except for venetoclax-cytarabine (p-value = 0.25) which is not significant but still lower. This reduction in toxicity suggests these combinations can offer effective treatment while minimizing side effects associated with current therapies.