secDrug: a pipeline to discover novel drug combinations to kill drug-resistant multiple myeloma cells using a greedy set cover algorithm and single-cell multi-omics

Multiple myeloma, the second-most common hematopoietic malignancy in the United States, still remains an incurable disease with dose-limiting toxicities and resistance to primary drugs like proteasome inhibitors (PIs) and Immunomodulatory drugs (IMiDs). We have created a computational pipeline that uses pharmacogenomics data-driven optimization-regularization/greedy algorithm to predict novel drugs (“secDrugs”) against drug-resistant myeloma. Next, we used single-cell RNA sequencing (scRNAseq) as a screening tool to predict top combination candidates based on the enrichment of target genes. For in vitro validation of secDrugs, we used a panel of human myeloma cell lines representing drug-sensitive, innate/refractory, and acquired/relapsed PI- and IMiD resistance. Next, we performed single-cell proteomics (CyTOF or Cytometry time of flight) in patient-derived bone marrow cells (ex vivo), genome-wide transcriptome analysis (bulk RNA sequencing), and functional assays like CRISPR-based gene editing to explore molecular pathways underlying secDrug efficacy and drug synergy. Finally, we developed a universally applicable R-software package for predicting novel secondary therapies in chemotherapy-resistant cancers that outputs a list of the top drug combination candidates with rank and confidence scores. Thus, using 17AAG (HSP90 inhibitor) + FK866 (NAMPT inhibitor) as proof of principle secDrugs, we established a novel pipeline to introduce several new therapeutic options for the management of PI and IMiD-resistant myeloma.


INTRODUCTION
Multiple myeloma (MM) is the second-most common hematopoietic malignancy in the United States [1]. MM is an age-dependent plasma cell neoplasm characterized by clonal expansion of malignant antibody-producing post-germinal-center B cellderived plasma cells within the bone marrow with significant complexity and heterogeneity at the molecular level [1]. Proteasome inhibitors (PIs) are standard-of-care chemotherapeutic agents for myeloma that impede tumor metastasis and angiogenesis by accelerating unfolded protein response (UPR) or the ubiquitindependent proteolysis of important regulatory proteins involved in key physiological and pathophysiological cellular processes in cancer cells and by interfering with the NF-κB-enabled regulation of cell adhesion-mediated drug resistance [2][3][4][5][6]. Bortezomib (Bz/ Velcade) was the first PI to be approved by U.S. Food and Drug Administration (FDA) for clinical application in 2003 for the treatment of relapsed and refractory myeloma [1,7,8]. Other examples include second-generation PIs Carfilzomib (Cz/Kyprolis) and the oral medication Ixazomib (Ix /Ninlaro/MLN9708) [7][8][9]. PIs are effective anti-MM drugs when used alone or in combination with other anti-cancer agents like immunomodulatory drugs (IMiDs), alkylating agents, topoisomerase inhibitors, corticosteroids, and histone deacetylase inhibitors (HDACis) [1,3]. However, despite these and other recent improvements in therapies, myeloma still remains a difficult-to-cure disease with dose-limiting toxicities and drug resistance and a median survival rate of only around 7 years [7,10]. Not all patients respond equally well to treatment and those who do often develop resistance over the course of treatment. Drug resistance may therefore be categorized into (1) innate resistance already present in drug-naive patients who never respond to treatment, or (2) emerging/acquired resistance where a patient's tumor ultimately undergoes relapse or "acquires" the ability to resist therapy in the course of treatment despite good response to initial treatment [2,11]. Therefore, there is an urgent need to search for novel secondary therapeutic options where new agents may be combined with standard-of-care drugs to achieve synergistic effects for treating drug resistance in myeloma.
Deciphering key features within patients underlying tumor heterogeneity and personalized sensitivity to chemotherapy is essential to predict the efficacy of anti-cancer drugs and to prevent delay in the selection of more effective alternative strategies [2,[11][12][13][14]. However, assessing the survival endpoints in clinical applications requires the treatment of a large number of patients with these drugs that need to be measured in months to years. Therefore, developing prediction algorithms of response can be a long process. One alternative is to use collections of human cancer cell lines from patient tumors that represent a broad spectrum of the biological and genetic heterogeneity of cancer, commonly known as in vitro modeling of drug response. We have compiled a panel of >70 human myeloma cell lines (HMCLs) representing the broad spectrum of biological and genetic heterogeneity of myeloma patients [15].
In this study, we have developed a computational method called secDrug for discovering novel synergistic secondary drug combinations that may effectively reverse resistance as combination regimens and allow for reduced dosing and toxicity of FDAapproved myeloma drugs. Next, we introduced single-cell transcriptomics as a novel screening tool for prioritizing secDrug combinations based on the subclonal expression of the drug targets and observed that the 17AAG + FK866 combination is potentially highly efficacious.
Further, to validate our prediction results, we used our HMCL panel as in vitro model system representing inter-individual heterogeneity in drug-response/resistance to show that the top predicted secondary secDrugs are indeed effective against PI-and IMiD resistance as single agents or as a combination. Further, using 17-AAG (an HSP90 inhibitor) as the test secDrug, we added functional assays, next-generation RNA sequencing, CRISPR-based gene editing, and high-dimensional mass cytometry (CyTOF/cytometry time of flight) in primary bone marrow cells (PMCs; ex vivo model system) from myeloma patients to create a multi-pronged approach/pipeline to discover, validate and characterize novel drugs as potential secondary choices for circumventing resistance to primary drugs in myeloma and to generate better treatment outcomes. This also allowed the identification of differentially expressed (DE) genes and novel pathways associated with the successful drug combinations.

MATERIALS AND METHODS In silico prediction of secondary drugs
Design and development of the secDrug pipeline are non-trivial and mathematically involved (details provided in Supplementary Methods section). Briefly, we utilized we used the vast array of human cancer cell lines in the Genomics of Drug Sensitivity in Cancer (GDSC version GDSC1000) database and created a pharmacogenomics data-driven greedy algorithm-based set-covering computational optimization method followed by a regularization technique to seek all secondary drugs that could kill the maximum number of cell lines of the test disease (B-cell cancers) resistant to the test drug (PI) in a sequential manner ordered by the number of cell lines killed. A greedy algorithm constructs a solution to an optimization problem piece by piece through a sequence of choices to find the overall, or globally, optimal solution. The GDSC1000 database is the largest public collection of information on sensitivity to >250 drugs covering a wide range of targets and processes involved in cancer biology in >1000 human cancer cell lines [16,17] Drugs, reagents, antibodies, and kits Ixazomib (Ixa) was procured from Takeda (Takeda Pharmaceuticals Inc., Deerfield, IL, USA). All other drugs were purchased from Selleck Chemicals. Drugs were dissolved in dimethyl sulfoxide (DMSO) and stored at −20°C. Recombinant Human IL-6 was obtained from Peprotech. Cleaved caspase-3/8/9, HSP90, c-Myc, p65, and IRF4 antibodies were purchased from Cell signaling Technology (CST). Monoclonal Anti-β-Actin-Peroxidase antibody produced in mouse was purchased from Sigma-Aldrich (St Louis, MO). Goat anti-Mouse/Rabbit IgG (H + L) secondary antibody (HRP conjugated) was obtained from Thermo Fisher Scientific. DHE (Dihydroethidium) assay kit and JC-1 Mitochondrial Membrane Potential (MMP) assay kit were purchased from Abcam (Boston, MA). Caspase-Glo 3/7 Assay System and CellTiter-Glo 2.0 Assay were purchased from Promega (Madison, WI).

Human myeloma cell lines (HMCLs)
HMCLs generated through the immortalization of primary myeloma cells were used as in vitro model systems to screen top secDrugs against sensitive, innate resistant, and acquired (Parental/P vs. clonally-derived resistant/R pairs generated using dose escalation over a period of time) myeloma [15]. We have also generated in vitro drug response profile for the four PIs: Bz, Cz, Oprozomib (Opz), and Ix as single agents in all the HMCLs included in the panel. PI-sensitivity in these cell lines was highly correlated, which suggests that any of these four PIs could be used as surrogates [15]. Therefore, we used Ixazomib as the representative PI in this study. Further, we have used machine learning-based computational approaches to derive a gene expression signature predictive of baseline PIresponse in myeloma [15]. The creation of the ANBL6 N-ras (ANBL6/Ras) codon 61 activating mutant cell line has been described earlier [18]. The IMiD resistant cell line, MM1S LenR, was obtained as a gift from Dr. Keith Stewart, Mayo Clinic, AZ. All cell lines were authenticated at source and tested randomly at regular intervals at the AU Center for Pharmacogenomics and Single-Cell Omics (AUPharmGx) using GenePrint 24 System (Promega). All cell lines are mycoplasma negative. HMCLs were maintained in HMCL media supplemented with IL-6 [19].

Human primary myeloma cells (PMCs)
Bone marrow-derived CD138+ ve patient PMCs were obtained through Mayo Clinic, MN following written informed consent and used as ex vivo model systems. Prior IRB approval was obtained from the Mayo Clinic review board. Participants were identified by number, not by name.

Establishment of RPMI8226 Hsp90 CRISPR-knockout cell line
Chemically-modified synthetic single-guide RNA (sgRNA) were designed targeting the Hsp90AA1 gene and synthesized by Synthego (Synthego Corporation, Menlo Park, CA). The sgRNAs were required to meet strict offtarget requirements of at least 2 mismatches within an early exon and target a common exon present in the majority of annotated transcripts. The sgRNAs were complexed together with the spCas9 to form a ribonucleoprotein (RNP). The RNPs were then delivered to RPMI8226 cells via an optimized electroporation setting. The transfected cells were then recovered for 2 days before the edits created were evaluated. Positive control sgRNA (RELA) were transfected at the same time. The edited site was PCR-amplified, and Sanger sequencing was performed on the amplicons Sequencing data was then analyzed using Synthego's Inference of CRISPR Edits (ICE) software tool to determine the percentage of knockout (KO) sequences of the genetic target [20]. ICE identifies the editing frequency and the specific indels present in the pool. Additionally, ICE calculates the frequency of the desired KO, reported as the KO score. Finally, once minimum KO editing efficiency was confirmed, RPMI8226/ Hsp90KO cells were expanded and QC tested.

In vitro chemosensitivity assays and drug synergy analysis
Cells were treated with increasing concentrations of secDrugs and PIs (represented by Ixazomib) or IMiDs (represented by Lenalidomide) as single agents or in combination for 48 h, and cytotoxicity assays were performed using CellTiter-Glo ® Luminescent cell viability assay (Promega Madison, WI). Luminescence was recorded in a Neo2 Microplate Reader (Biotek), and halfmaximal inhibitory concentration (IC 50 ) values were determined using GraphPad Prism software by calculating the nonlinear regression using sigmoidal dose-response equation (variable slope). Drug synergy was calculated using Calcusyn software based on Chou-Talalay's combination index (CI) method and the isobologram algorithm (Biosoft, US) [21].

Apoptosis assays
Caspase-3/7 activity assay was performed on the HMCLs using Caspase-Glo 3/7 luminescent assay kit according to manufacturer's instructions H. Kumar et al.
(Promega Madison, WI) using Neo 2 Microplate Reader (Biotek). Cell death by apoptosis was also measured by immunoblotting analysis.

Determination of superoxide levels
Cells were incubated with 5 μM DHE (in RPMI) for 15 min in the dark at 37°C. Cells were then washed once with cell-based assay buffer, and red fluorescence was recorded by Synergy Neo2 multi-plate reader.

Measurement of mitochondrial membrane potential (MMP)
Cells were incubated with 5 μM JC-1 dye for 15 min in the dark at 37°C and washed twice in PBS, and then analyzed for red and green fluorescence by Synergy Neo2 multi-plate reader.

Mass cytometry (CyTOF)
Thirty-seven antibody targets directed against cell surface and intracellular markers were utilized as Immunphenotyping Panel for CyTOF analysis. The Antibody markers and respective metal conjugate are described in Supplementary Table S1. Panels were designed using the web-based panel designer software: Maxpar Panel Designer (www.fluidigm.com) for optimal signals, minimum background due to oxidation, isotopic purity, and sufficient sensitivity for each targeted marker. Prelabeled antibodies were purchased from Fluidigm. Purified antibodies from Biolegend and Santacruz Biotechnology were labeled using an X8 polymer MaxPAR antibody conjugation kit (Fluidigm) according to the manufacturer's instructions. CyTOF analysis was performed on PMCs treated with DMSO (vehicle/control), 0.2 µM 17AAG,1 µM 17AAG, and 5 µM 17AAG.

CyTOF data analysis
Cytobank software version 7.3.0 (Santa Clara, CA) was used for cleanup of cell debris, removal of doublets and dead cells. Cleaned fcs files were further gated and analyzed by Cytobank. Plasma cells were identified as CD19−, CD16−, CD3−, CD38+, and kappa OR lambda+ (based on each patient's kappa or lambda restriction from clinical flow data). If the plasma cells had diminished surface CD38 expression as a result of previous daratumumab exposure, CD229 was used as a positive selection marker. T-distributed stochastic neighbor embedding (t-SNE), viSNE, and FlowSom plots were generated to visualize the subpopulation architecture based on markers of interest. Relative marker intensities and cluster abundances per sample were visualized by a heatmap.

Single-cell RNA sequencing (scRNAseq)
Automated single-cell capture and cDNA synthesis were performed at 1500 tumor cells/sample using 10× Genomics Chromium platform that uses droplet-sequencing-based chemistry. Single-cell RNA sequencing was performed on Illumina HiSeq 2500 NGS platform (Paired-end. 2 × 125 bp, 100 cycles. v3 chemistry) at >50 million reads per sample. scRNAseq data analysis scRNAseq datasets were obtained as matrices in the Hierarchical Data Format (HDF5 or H5). A combination of Seurat and Partek Flow software packages was used to pre-process the data and perform single-cell transcriptomics analysis. Highly variable genes for clustering analysis were selected based on a graph-based clustering approach. The visualization of cell populations was performed by t-SNE.

Next-generation RNA sequencing (NGS)
HMCLs were plated at a density of 4 × 10 5 cells per mL, and 0.5 μM of 17-AAG was added as a single agent or in combination with 15 nM of Ixazomib. Baseline (untreated) and post-treatment (treated) cells were collected 24 h post-treatment. High-quality RNA was extracted using QIAshredder and RNeasy kit (Qiagen). RNA concentration and integrity were assessed using Nanodrop-8000 and Agilent 2100 Bioanalyzer and stored at −80°C. An RNA integrity number threshold of eight was applied, and RNA-seq libraries were constructed using Illumina TruSeq RNA Sample Preparation kit v2.
NGS Libraries were size-selected, and RNA sequencing (RNAseq) was performed on Illumina's NovaSeq platform using 150 bp paired-end protocol with a depth of >20million reads per sample.

RNAseq data analysis
Gene expression data were pre-processed, log 2 -transformed, and analyzed using a combination of command-line based analysis pipeline (DEseq2 and edgeR) and Partek Flow software to identify differential gene expression profiling (GEP) signatures. Genes with mean counts<10 were removed, and CPM (counts per million) data was used to perform differential expression testing to identify GEP signatures. Due to the small sample sizes, we used GSA to perform differential gene expression analysis between groups that applies limma, an empirical Bayesian method, to detect the DE genes (DEGs). Genes with mean fold-change > |1| and p < 0.05 were considered as the threshold for reporting significant differential gene expression. Heatmaps were generated using unsupervised hierarchical clustering (HC) analysis based on the top DEGs.

Pathway analysis
Ingenuity pathway analysis (IPA) software (QIAGEN) was used to identify the molecular pathways and upstream regulators predicted to be activated or inhibited in response to 17-AAG treatment (single-agent and combination with PIs) based on the list of significantly differentially regulated genes [22].

Western Blotting
HMCLs treated with 17-AAG alone, Ixa alone, or 17-AAG + Ixa combination were harvested, washed, and lysed using radioimmunoprecipitation assay (RIPA) lysis buffer containing 50 mM Tris-HCl, pH 7.5, 150 mM NaCl, 1% NP40, 5 mM EDTA, 1 mM DTT, phosphatase, and protease inhibitors cocktail (Sigma) and incubated on ice for 15 min. Samples were then centrifuged at 14,000 rpm at 4°C for 30 min. The supernatant was then aspirated and quantified using Pierce ™ BCA Protein Assay Kit (Thermo Scientific). Samples were solubilized in sodium dodecyl sulfatepolyacrylamide gel electrophoresis sample buffer, and equal amounts of protein were loaded per lane of 10% sodium dodecyl sulfatepolyacrylamide gels and transferred onto PVDF membranes (Millipore; Billerica, MA). Membranes were blocked in TBS with SuperBlock ™ blocking buffer (Thermo Fisher), incubated with primary antibodies and secondary antibodies in TBS with 0.2% Tween 20 and 2.5% bovine serum albumin. Immunoreactivity was detected by chemiluminescent HRP substrate (Bio-Rad), and the exposed image was captured using a ChemiDoc ™ MP Imaging System (Bio-Rad). Densitometry analysis was performed using Image J software.

Statistical analysis
All statistical analyses were performed using R (the project for statistical computing and graphics) and GraphPad Prism 9.0 software. All tests were two-sided, and p < 0.05 was considered statistically significant.

RESULTS
Identification of secondary drugs using secDrug algorithm Details on the design and development of the secDrug pipeline are provided in the Supplementary Methods section. Briefly, a novel modified greedy algorithm-based set-covering computational optimization-regularization pipeline was used to identify all secondary drugs that could kill the maximum number of cell lines in the GDSC1000 database belonging to the test disease (B-cell cancers including myeloma) and which are resistant to the PI/PI drug Bortezomib (Bz/Velcade; the primary anti-myeloma drug). A total of 1091 cells lines were present in the GDSC1000 database [16,17]. The following filtering criteria were applied to select computable B-cell lines: target cell-B-cell; cancer type-blood; tissue-blood; histology-lymphoid_neoplasm/haematopoietic_neoplasm; site-haematopoietic_and_lymphoid_tissue; no missing data). A total of 94 cell lines satisfied the above filtering criteria and were selected for further analysis. IC 50 values were processed, imputed, and categorized as S (PI-sensitive), R (PI-resistant), and N ("Neutral"/Intermediate PI IC 50 values) prior to analysis. We applied our computation algorithm to the GDSC1000 dataset and predicted the top secDrugs that can be best combined with PIs to achieve response in N and R lines. The predicted top secondary drug combinations in PI-resistant + PI-neutral B-cell cancers with a PI backbone are shown in Table 1. These include HSP90 inhibitor (17-AAG), Nicotinamide phosphoribosyltransferase or Nampt inhibitor (FK866), PIKfyve inhibitor (YM201636), Raf inhibitor (PLX-4720), Bcl2 inhibitor (Navitoclax), SB505124 (transforming growth factor-β type I receptor, ALK4, ALK7 inhibitor), S6K1- Percent coverage (cell lines predicted to be killed by each combination treatment regimen) of the cancer lines included in the prediction model is also provided.  Fig. 3A-C. We found that 17-AAG not only showed synergy with PIs and IMiDs, the combination of 17-AAG and FK866 also showed significant synergy, as depicted by the dose-response curves and CI values representing combination treatments. CI values were consistently less than 1, which indicates synergy [25].
In addition, FK866 also showed synergy with Ixazomib (Ixa + FK866 survival curves are shown in Fig. S1A). Cell survival curves representing other top secDrugs + Ixazomib combination in innate sensitive, innate resistant, and acquired resistant HMCLs are shown in Fig. S1B-F. These secDrugs also showed strikingly high synergy with Ixazomib, as depicted by the CI values. Further, based on dose reduction index (DRI) values, the IC 50 of Ixazomib in myeloma cell lines was predicted to be significantly reduced in the presence of these secDrugs. Figure S2 shows the relative decrease of the predicted effective IC 50 (nM concentration) of Ixazomib when used in combination with 17-AAG. The DRI values, calculated using CI theorem, demonstrated that 17-AAG improved the therapeutic index of PI and IMiD administration to the cells and decreased the amount of PI/IMiD required to achieve effective responses [25]. This points towards the possibility of reducing the dose and thereby toxicity of PIs when administered as 17AAG + PI combination. Druginduced apoptosis was confirmed in HMCLs using Caspase 3/7 activity assays (data not shown).

CyTOF analysis revealed 17-AAG-induced cell death of PMCs and key changes in myeloma-specific proteomic markers
High-dimensional mass cytometry or CyTOF analysis is a deep immunophenotyping method that combines flow cytometry and elemental mass spectrometry [26]. We performed CyTOF analysis on PMCs obtained from myeloma patients (n = 6) to assess 17-AAG-induced cell death through apoptosis as well as to evaluate changes in phenotypic and functional markers in MM cells at the single-cell/subclonal proteomics levels. As shown in Fig. 4A, CyTOF analysis following exposure to 17-AAG treatment revealed a distinct cluster of cells defined by elevated cleaved caspase levels in the primary samples, indicating treatment-induced cell death by apoptosis in the cells exposed to 17-AAG.
Myeloma cells are addicted to several proteins. Figure 4B shows results from our CyTOF-based differential expression (pre vs. post-17-AAG treatment) analysis, revealing shifts in IRF4 and pSTAT3 as well as myeloma cell survival markers like CD138 and phosphoproteins like pRB.
Western blotting analysis Treatment-induced protein expression of the phenotypic/ functional markers of myeloma (p65/NFκB, IRF4, and cMyc) and markers of apoptosis (including Cleaved Caspase-3, Cleaved Caspase-9, etc.) was confirmed using immunoblotting analysis in HMCLs. Figure S3A shows representative pre-vs. post-17-AAG treatment immunoblotting results on these myeloma cell survival and apoptotic pathways. Densitometry analysis results are provided in Fig. S3B. Our results show a substantial decrease in IRF4, p65, and cMyc following 17-AAG treatment and a concurrent increase in Cleaved Caspase-3, Cleaved Caspase-9 protein expression, which was also confirmed at the mRNA level using pre-vs. post-17AAG-treatment differential gene expression (RNAseq) analysis, along with several other myeloma protein/survival markers like STAT1, RELB, NFKBIA, NFKB2, and IKZF3.

17-AAG induces apoptosis via a mitochondrial-mediated pathway in myeloma
To investigate if 17-AAG imparts its cytotoxic effects in myeloma by generating reactive oxygen species (ROS), particularly superoxides and hydrogen peroxide (H 2 O 2 ), cellular superoxide anions were measured by using the fluorescent dye DHE (Abcam). MMP was assessed using JC-1 (Abcam). JC-1 is a cationic carbocyanine dye that accumulates in mitochondria. The dye exists as a monomer (green fluorescence) at low concentrations and changes color from green to red in energized mitochondria.
We observed induction of cellular superoxide anions (Fig. 5A) and intracellular ROS production (Fig. 5B) that causes mitochondrial membrane depolarization following 17-AAG treatment in myeloma cells representing sensitive and resistant disease.
17-AAG induced cell death was comparable with Hsp90 knockdown Next, we compared the effect of Ixa and Ixa+17-AAG combination therapy between wild-type and CRISPR-mediated HSP90AA1 gene knockdown cell lines. Dose-response curves in Fig. S4   toward an on-target effect of 17-AAG treatment leading to the anti-MM efficacy.
NRas mutant HMCL showed high sensitivity to 17-AAG treatment Finally, the myeloma cell lines ANBL6P, ANBL6VR, and ANBL6 NRas mutant were treated with single-agent 17-AAG, Ixa, and 17-AAG + Ixa combination. We have described earlier that these activating mutations of the ras oncogenes in ANBL6 (ANBL6 NRas) may lead to growth factor independence and suppression of apoptosis [18]. Notably, our ANBL6 NRas mutant cell line showed 20-times greater 17-AAG sensitivity (lower IC 50 ) compared to the ANBL6P or VR cell lines (Fig. S5).
Finally, IPA predicted upregulation of the microRNA let-7 (zscore 7.180; p-value 5.02e−12) as the top upstream regulator based on significantly DEGs (Fig. S8) Creation of secDrug software package Finally, we developed an R software package based on our secDrug pipeline for predicting novel secondary therapies in chemotherapy-resistant cancers. secDrug takes a query of any cancer type and any test/primary/standard-of-care drug and outputs a list of the top secondary drug combinations with a confidence score and biological pathway visualization. Thus, secDrug has potential application in clinical decision-making for discovering resistance-reversing cancer chemotherapy regimens. R codes for the package are available at https://github.com/Ujjal- Mukherjee/secDrug/tree/main/CombinationDrugMyeloma, and the datasets are available at GitHub repository.

DISCUSSION
Drug resistance is a major obstacle in achieving a complete and sustained therapeutic effect in cancer chemotherapy [2,11,14,27,28]. Chemo-resistance may also lead to overdosing and unwanted exposure to ineffective anti-tumor agents, thereby increasing the risk of negative side effects and the cost of drug development [29,30].
Further, we performed extensively in vitro, ex vivo, and functional validation in research models of refractory and resistant myeloma to validate 17-AAG + FK866, 17AAG + PI, and 17-AAG + IMID as combination treatment candidates that also served as a proof-of-principle for our secDrug pipeline. Overall, our validation results corroborated with our in silico prediction of secondary drugs based on secDrug analysis.
17-AAG/Tanespimycin has previously been shown to work against myeloma, in vitro, in vivo (animal studies) as well in clinical studies [31][32][33][34]. However, to our knowledge, this is the first study that specifically evaluates the use of 17-AAG combination therapy in relapsed and refractory myeloma models. Further, in our study, the ANBL6 Ras mutant cell line showed 20-times lower 17-AAG cytotoxicity compared to the ANBL6P/VR cell lines. An earlier study in metastatic malignant melanoma has shown that a patient harboring NRAS-activating mutation exhibited disease stabilization for 49months following administration of pharmacologically active doses of 17-AAG [35]. Mutations of NRAS have earlier been shown to be significantly associated with lower single-agent PIsensitivity and shorter time to progression in bortezomib-treated myeloma patients [36]. Thus, our study points towards a unique niche (NRas-mutant myeloma) where 17-AAG could be highly effective as a single agent as well as in combination with PIs and FK866, in addition to relapse and refractory myelomas. Moreover, we evaluated the molecular pathways involved in response to the top secondary drugs, which provided additional insights into the mechanism of action of 17-AAG as a secDrug.
Myeloma tumor cells have elevated intracellular NAD+ levels that support the high rate of energy metabolism for uncontrolled proliferation, tumor cell growth, and survival [37]. FK866 is a chemical inhibitor of Nampt (Nicotinamide phosphoribosyltransferase), a key enzyme in NAD+ metabolism [38]. Consequently, FK866 has been shown to reduce myeloma tumor growth in PIsensitive and PI-resistant myeloma through activation of autophagy and cell death in myeloma cells [39]. In this study, we showed that FK866 not only overcomes PI-resistance when used as a single-agent or as an Ixa combination, combining 17-AAG + FK866 is highly synergistic against our validation models of relapsed/refractory myeloma.
Our study introduces several novels secDrugs as potential synergistic partners of PIs that have never been studied as potent single-agent or combination therapy options in myeloma model systems, including KIN001-102 (A6730; Akt1/2 kinase inhibitor) and SB505124 (inhibitor of transforming growth factor-β type I receptor or ALK4, ALK7 that activates the SMAD2/3 pathway). These may serve as novel candidates for further studies on the pre-clinical and clinical validation in xenograft or mouse models of myeloma.
Although some of the other predicted secDrugs have earlier been shown to be effective against myeloma, very few studies have explored their efficacy as drug combinations with PIs/IMiDs in models of refractory/resistant myeloma. For example, PF-4708671 is a P70S6K1 isoform-specific inhibitor that has recently been shown to induce statistically significant apoptosis in HMCLs and PMCs in combination with several standard-of-care therapies [40]. NEDD8-activating enzyme/neddylation-inhibitor/MLN4924 has once earlier been shown effective against a subset of cell lines represented by cell surface expression of TNFR1 [41]. PLX4720 (a small-molecule, ATP-competitive inhibitor of Mutant BRAF kinase) was earlier shown to have a partial single-agent response in patients harboring subclonal BRAF mutations [42,43]. Our in silico predictions and single-agent cytotoxicity data thus builds a strong case to test these drugs as secDrug combination regimens in a wider panel of HMCLs representing refractory and clonally derived acquired resistant cell lines. Among the other secDrugs, Navitoclax is a high-affinity small-molecule BH3 mimetic that inhibits Bcl2 and Bcl-xl. Navitoclax has been shown to inhibit cell proliferation in myeloma leading to induction of apoptosis [44,45]. YM201636 is an inhibitor of PIKfyve, a mammalian protein involved in the regulation of crucial cellular functions, including nuclear signaling and autophagy. Few recent studies demonstrated the therapeutic efficacy of PIKfyve inhibitors in myeloma cell lines [46,47].
Overall, we present here a unique pipeline that introduces not only novel secDrugs but also provides additional niches for secondary drugs that are already under preclinical or clinical investigation in relapsed/refractory myeloma.
Our findings provide a strong case for combining the top predicted secDrugs with PIs and IMiDs to overcome resistance and thereby improve patient outcomes. This potentially introduces many more drugs as new and more effective therapeutic options for the management of resistant myeloma with a high probability of clinical success that promises to improve the quality of treatment, maximize drug efficacy, minimize toxicities and adverse drug reactions from over-dosing and decrease the rate of mortality in myeloma patients. A logical extension of this pipeline would be the development of model systems where a combination of more than two secDrugs can be effectively tested.
The integration of in silico modeling-based pipeline with singlecell technologies (scRNAseq and CyTOF analysis) introduces an innovative, evidence-based application in clinical decision-making that will minimize the number of test drugs required for discovering successful combination chemotherapy regimens against drug-resistant cancers.

DATA AVAILABILITY
All data used in the analysis are available on reasonable request from the corresponding author.

CODE AVAILABILITY
All code used in the analysis is available on reasonable request from the corresponding author.