Metabolomics-based phenotypic screens for evaluation of drug synergy via direct-infusion mass spectrometry.

Drugs used in combination can synergize to increase efficacy, decrease toxicity, and prevent drug resistance. While conventional high-throughput screens that rely on univariate data are incredibly valuable to identify promising drug candidates, phenotypic screening methodologies could be beneficial to provide deep insight into the molecular response of drug combination with a likelihood of improved clinical outcomes. We developed a high-content metabolomics drug screening platform using stable isotope-tracer direct-infusion mass spectrometry that informs an algorithm to determine synergy from multivariate phenomics data. Using a cancer drug library, we validated the drug screening, integrating isotope-enriched metabolomics data and computational data mining, on a panel of prostate cell lines and verified the synergy between CB-839 and docetaxel both in vitro (three-dimensional model) and in vivo. The proposed unbiased metabolomics screening platform can be used to rapidly generate phenotype-informed datasets and quantify synergy for combinatorial drug discovery.

Introduction

Combination therapy has become a cornerstone to combat complex diseases such as cancer, metabolic, and autoimmune disorders with increased efficacy, potency, and decreased occurrence of drug resistance (Ascierto and Marincola, 2011; Mokhtari et al., 2017). In the United States alone, there are currently more than 10,000 ongoing clinical trials investigating combinatorial therapies (Author Anonymous, 2017), and even more preclinical studies (Liu et al., 2019). However, the success rate of combination therapy development and application remains low, in part due to poor reproducibility of the preclinical findings in clinical trials and ill-informed pharmacodynamic data in the early stages of clinical studies (Day and Siu, 2016). Therefore, new strategies that are more robust and informative are needed to correctly identify the most effective combinations for rapid and successful translation into the clinic.

Typical preclinical drug combination studies rely on the application of dose-response methods for synergy quantification, including Loewe (1953), Bliss independence (Bliss, 1939), or Chou-Talalay (Chou, 2010) models, that utilize single biochemical readouts, such as in vitro cell viability assays (Folkesson et al., 2020; Tomska et al., 2018; Wood et al., 2017). While the Bliss independence and Loewe models allow for the assessment of hundreds to thousands of drug combinations and can easily be applied to multiple cell lines simultaneously, they can often produce inconsistent results (Meyer et al., 2020), and thus more recent attempts have been made to unify the drug synergy principles (Wooten et al., 2021). The Chou-Talalay model has been considered the gold standard for evaluating synergism, but due to its requirement of a large number of data points to assess a wide range of concentrations between only two drugs, it is not applicable to high-throughput platforms (Amzallag et al., 2019). Another limitation of these methods is represented by the relatively high number of top hits that are often generated when screening large drug libraries. Even with seemingly promising results derived from these models, it is difficult to assert with confidence whether the selected synergistic drug combinations are likely to be clinically relevant because these methods do not provide insight into the molecular response of the combinatorial treatment. Therefore, the application of data-driven methods with a high-dimensional input (Celebi et al., 2019) might prove valuable to further refine the search and identify the most promising synergistic combinations.

Phenotypic methods that evaluate drug effects using more elaborate multivariate readouts, such as metabolomics studies, have gained a resurgence of interest in recent years (Aulner et al., 2019) and have been increasingly used to investigate treatment response over the traditional single readout methods (Zampieri et al., 2017, 2018; Dubuis et al., 2018; Guijas et al., 2018). As an example, we developed a 96-well format NMR-based high-content metabolomics screen to elucidate the actions of 56 kinase inhibitors in mammalian cells (Tiziani et al., 2011). Creek et al. utilized untargeted metabolomics analysis to test 90 antimalarial compounds (Creek et al., 2016) and captured significant metabolic modulations of in vitro cultured parasites. A subsequent antimalarial drug discovery study by Murithi et al. screened a library of 33 compounds and annotated 195 metabolites from parasite extracts to elucidate their modes of action (Murithi et al., 2020). Along with these efforts to reveal the mode of actions of single drugs, an experimental-computational approach (Campos and Zampieri, 2019) was recently developed to predict the efficacy of a drug combination using high-throughput metabolomics by integrating the metabolic profiles induced by single drugs and previously generated metabolic profiles of gene-knockout strains (Murithi et al., 2020). Previous methods were only based on single drug responses, which have shown narrowness in their potential applications due to the limitations in biological knowledge and existing databases (Wu et al., 2021). This is due to the complexity of cellular regulatory networks, and unique and unpredictable combinatorial drug responses that diverge from the summated single responses (Diaz et al., 2019; Menden et al., 2019). Despite the extensively proposed computational approaches for drug synergy prediction (Madani Tonekaboni et al., 2018), no algorithm is available to calculate the extent of synergy resulting from both individual and combinatorial responses in multivariate datasets. Similar to proteomic and transcriptomic data (Yang et al., 2020; Diaz et al., 2020), phenotypic multivariate metabolomics studies can provide deep insight into the drug effects of synergistic combination with an increased likelihood of clinical efficacy (Guijas et al., 2018).

In this study, we developed a metabolomics-based phenotypic screen in a three-dimensional (3D) spheroid prostate cancer model to evaluate synergy in drug combinations. To evaluate synergy between two drug treatments, based on the responses to both the individual and the combined drug treatments, we implemented a two-step approach integrating: (i) a primary screening that uses a conventional univariate cell viability readout to identify top hits, with (ii) a secondary phenotypic screen that uses stable isotope-tracer direct-infusion mass spectrometry (SIT-DIMS), to further restrict the selection of the most promising drug combinations based on a multivariate readout. Our approach is generalizable to any -omics datasets; however, in consideration of the ever-increasing number of drugs being developed that interfere with cancer cell metabolic processes (Alkan et al., 2018; Bruntz et al., 2017; Cluntun et al., 2017; Gregory et al., 2019; Lee et al., 2020; Li et al., 2020; Matre et al., 2016; Molina et al., 2018; Pavlova et al., 2018a; Vanhove et al., 2019; Voelkel-Johnson et al., 2018; White et al., 2017), metabolomics datasets can be particularly advantageous to evaluate drug combinations that include such drugs. While several analytical techniques are available and have been used for large-scale metabolomics studies, most require either long acquisition times (at least 20 min per sample) or can only detect a limited number of targeted metabolites of interest (Ganna et al., 2014, 2015; Cao et al., 2020). DIMS is well suited to perform high-throughput metabolic analyses by virtue of its relatively short acquisition time and sufficient representation of the overall metabolic fingerprint of samples (Habchi et al., 2016; Dubuis et al., 2018).

In the context of the secondary screen, to evaluate synergistic drug combinations based on the phenomics datasets, we developed a new principal component analysis (PCA)-based Euclidean distance synergy quantification (PEDS) algorithm able to quantify drug synergy and identify the most promising drug combinations from multi-array datasets. Thus, the secondary screen was used to validate and further select the synergistic hits generated from the primary screen through the application of PEDS. Importantly, an additional advantage of our methodology is that the metabolic data can be further mined to provide additional metabolic information on the single and combinatorial treatments’ behavior which can be used to investigate the molecular drug action. Here, we report the results of a drug library screen including the primary screen based on cell viability; the resulting top-hits combinations from the primary screen were then further evaluated using SIT-DIMS and PEDS. In addition, a recently developed data mining method (Matsuda et al., 2020) that allows the in-depth analysis of stable isotope-tracer enrichment data (from the secondary screen) was used to investigate the metabolic reprogramming induced by drug combinations compared to the individual drug treatments. The single top-hit synergistic drug combination that emerged from the two-step screen was then further validated both in vitro (by the Bliss independence model, Loewe additivity model, and additional biochemical assays) and in vivo, as well as with an in-depth metabolic flux analysis (MFA) to fully elucidate the metabolic mechanism underlying the synergistic treatment response.

Results

The overall workflow, implemented within the drug synergy screening platform and the subsequent validation of the top hits, consists of 4 main steps (Figure 1). The 2-step screen begins with a primary single-readout high-throughput screen of a drug library across a panel of cancer cell lines; the potential synergistic drug combinations are identified based on the Bliss independence model (Bliss, 1939). The subset of synergistic drugs from the primary screen are then re-assessed using the SIT-DIMS secondary screen performed in a 96-well format using a 3D spheroid cell model system. The newly developed PEDS algorithm is then used to quickly and cheaply further evaluate the synergistic drug combinations identified by the primary screen and further restrict the search. The selected candidate combination can then be prioritized for more robust analyses and mechanistic studies.

Figure 1

Schematic of the metabolomic-based phenotypic screening platform for drug synergy discovery

The workflow to evaluate synergies between agents in a drug library consists of four main steps. The primary screening uses an ATP bioluminescence assay and selects top-hit drug candidates based on the Bliss independence model. The secondary stable isotope-tracer direct-infusion mass spectrometry (SIT-DIMS) metabolomics-based phenotypic screening refines the synergy discovery based on multi-readout DIMS metabolomics data analyzed using our PEDS algorithm. The PEDS-selected top candidates can then be validated in vitro and in vivo, and the metabolic synergy mechanism can be further investigated by metabolic flux analysis (MFA) and computational data mining.

To validate the synergy and optimize the relative drug concentrations of the combination, a serial drug dosing modulation was performed. Both the Bliss independence and Loewe additivity models at serial drug doses were used, followed by the apoptosis assay. An in vivo study was also performed to confirm that the drug combination exerts a synergistic effect on HMVP2 PCa allograft tumors (Saha et al., 2016, 2017). Finally, the top candidate combination identified by the PEDS algorithm and synergy validation was further explored using an in-depth MFA and data mining method to fully characterize the metabolic changes induced by the drug combination and clarify the synergy mechanism.

Primary high-throughput, single-readout drug screen: identification of synergies between CB-839 and a drug library

Several recent studies highlighted the upregulation of glutamine utilization and dependence in multiple blood (Matre et al., 2016; Gregory et al., 2019) and solid cancers (Vanhove et al., 2019; Cluntun et al., 2017; Lee et al., 2020), including prostate cancer (White et al., 2017; Lodi et al., 2017), as a possible target of treatment. Therefore, we chose to screen candidate synergistic combinations between the glutaminase inhibitor CB-839 and a drug library of 292 compounds. A panel of four prostate cancer (PCa) cell lines (PC3, DU145, 22Rv1, and HMVP2), and three normal prostate (NP) cell lines (RWPE-1, PWR-1E, and NMVP) were screened.

To evaluate an appropriate treatment dose, the cell lines were initially treated for 24 h with the library compounds, administered as single agents, at multiple concentrations (100, 250, 500 nM, and 1μM). Cell viability was assessed using an ATP bioluminescence assay to determine drug sensitivity. The selection of the treatment dose for the library compounds was based on the following criteria: the majority of the drugs exert a (i) moderate (between 10%–50%) reduction in cell viability, and (ii) significant response compared to the control. For all the PCa cell lines, the 250 nM treatment dose had the highest number of drugs that met the criteria (Table S1). Treatment with CB-839 at the lowest tested dose (100 nM) reduced the cell viability of PC3, DU145, 22Rv1, and HMVP2 by 21.11% ± 5.02%, 43.24% ± 4.88%, 11.37% ± 10.58%, and 66.95% ± 3.96%, respectively, closely matching the criteria outlined above, and was selected for the subsequent combinatorial drug screening. The NP cell lines treated with the individual agents at the 100 and 250 nM doses were marginally more resistant than the PCa cell lines (DU145 and 22Rv1) screened at the same drug library concentrations (Table S2); similarly, the murine NP NMVP line was more resistant than HMVP2 and the human RWPE-1 line was more resistant than all the human PCa cells (PWR-1E was more resistant than DU-145 but not the other PCa cell lines) to 100 nM CB-839. In summary, based on the PCa and NP cell viability response to single treatments, we chose to use 250 nM as the drug library dose in combination with 100 nM CB-839 to perform the primary screen and identify new synergistic drug combination candidates.

Cell viability inhibition (based on the ATP bioluminescence assay) of both single and combination treatments was used to calculate the Bliss index using bootstrap resampling. Out of the 292 drug combinations tested across a panel of four prostate cancer cell lines, 147, 96, 82, and 118 candidate drug pairs were potentially synergistic (2.5th percentile of Bliss index>0) with CB-839 in 22Rv1, PC3, HMVP2, and DU145, respectively. Of those, 20 of the drugs (marked in red in Figure 2A) were potentially synergistic in all four of the PCa cell lines. These were thereby selected for the subsequent secondary screening using DIMS high-content technology (Southam et al., 2016). To test the selectivity of the drug combinations on cancer cells over normal cells, we also confirmed that NP cells were only marginally affected by the combinatorial treatments compared to PCa cells (Figure 2B: murine lines, HMVP2 compared to NMVP; Figure 2C: human lines, DU145 compared to RWPE-1; Figure S1).

Figure 2

The primary screening identified twenty drugs potentially synergistic with CB-839 in four PCa cell lines

(A–C) (A) Bliss ndex value (2.5th percentile, based on bootstrapping resampling of relative cell viability data from the primary screening) for a library of cancer drugs combined with CB-839 and administered to 4 PCa cell lines. The 2.5th percentile of Bliss index value higher than zero indicates that the drug combination is significantly synergistic. Twenty drugs (marked in red on the left) have Bliss index values greater than zero across all four prostate cancer cell lines. Relative cell viability comparison in PCa versus NP cells following individual (blue or grey squares) or combined (with CB-839; red or grey circles) treatment in (B) murine (NMVP vs HMVP2) and (C) human (RWPE1 vs DU145) prostate cells. The 20 top-hit drugs (marked in red in (A)) are shown as yellow squares (individual) or green circles (combined treatment), while the other library drugs are shown in gray.

Interestingly, eight of the 20 top-hit candidates that synergized with CB-839 in PCa cells are currently used in standard of care chemotherapy regimens or clinical trials for several cancer types. However, out of the five drugs included in the 292-compound drug library that are currently approved for treating prostate cancer, only one, docetaxel, was synergistic with CB-839 in all four PCa cell lines (Figure S2). From the 20 selected drugs, we also found drugs sharing similar primary targets: three HSP90 inhibitors (BIIB021, AUY922, and 17-DMAG), three mTOR inhibitors (temsirolimus, AZD8055, and rapamycin), and three microtubule inhibitors (ABT-751, nocodazole, and docetaxel).

Secondary multiple readout SIT-DIMS screen: mining of drug-induced metabolic modulations

To validate the 20 top-hit synergistic drug combinations identified from the primary screen as well as to investigate the molecular/metabolic modulation they induce, we utilized an automated chip-based nanoelectrospray ionization (nESI) source connected to a high-resolution Orbitrap mass spectrometer capable of acquiring both positive and negative ionization modes data in approximately 1.5 min per sample (1500 combined positive and negative ionization mode acquisitions or about 750 samples per day). Samples were prepared in 96-well format (13 min per plate) to ensure the demand of large throughput and fast quenching for halting metabolism (Lu et al., 2019). Moreover, the high sensitivity of the SIT-DIMS method allows acquiring data on cell samples from a single well of a 96-well plate. We assessed the high sensitivity and reproducibility of the DIMS-based screening method using 0.1, 0.5, and 1 ppm internal standards (ISs) run sequentially in both positive and negative ionization modes with tMS2 (as detailed in the Methods section). All the 14 ISs were detected at 0.1 ppm, and all coefficient of variations (CVs) of peak intensities at low, medium, and high concentrations were lower than 15% (Figure S3). No carryover effects were observed when comparing the two blanks run before and after the sequence.

To better recapitulate the in vivo tumor microenvironment, PCa cells were aggregated into 3D spheroids in 96-well plates and treated with the top-hit candidates either as single agents or in combination with CB-839 for 24 h. As proof of concept, one human (DU145) and one murine (HMVP2) PCa cell lines were selected for the secondary screening. During the treatment, PCa cells were maintained either in the regular (unlabeled) medium, or in media containing either 13C6-glucose or 13C5,15N2-glutamine to probe the incorporation of these nutrients into major metabolic pathways (Bruntz et al., 2017) using the high-resolution and high mass accuracy SIT-DIMS screen.

While this is not required to evaluate synergy between the drug combinations (which is described below), our methodology is based on the acquisition of metabolic data that can be used to further evaluate the molecular response of cells to the individual and combined treatments. Both the individual and combination treatments administered to prostate cancer cells induce notable modulations of the cellular metabolic profiles in both HMVP2 (Figure 3A) and DU145 (Figure 3B) cells. Hierarchical clustering of the unlabeled data (based on 69 and 108 identified metabolites in HMVP2 and DU145 cell lines, respectively) clearly segregates samples treated with CB-839 (either alone or in combination) from samples that did not receive CB-839. As anticipated, in response to glutaminase inhibition by CB-839, glutamine and glutamate levels were strongly modulated as well as other metabolites involved in glutaminolysis and the TCA cycle (highlighted in red in Figures 3A and 3B). In DU145 cells, the administration of CB-839 also induced increased levels of several precursors of sphingolipid metabolism (marked in cyan, Figure 3B) which have been known to play a key role in malignancy and chemoresistance (Voelkel-Johnson et al., 2018). Another interesting observation in the DU145 cells was the dramatic surge in the levels of several amino acids (marked in blue, Figure 3B) following treatment with AZD8055 (with or without CB-839) likely resulting from the restored amino acid balance by the mTOR (amino-acid-responsive serine/threonine kinase mechanistic target of rapamycin) inhibitor. Notably, most samples receiving CB-839-containing treatments resulted in the intracellular accumulation of asparagine (marked in green in Figures 3A and 3B). To clarify the origin of this specific CB-839-induced metabolic modulation, we also used the SIT-DIMS method (in spheroids in a 96-well format and with 13C4,15N2-asparagine as the tracer) and confirmed that, in line with previous reports (Pavlova et al., 2018a), the observed accumulation of asparagine is due to the upregulation of asparagine uptake (increased intracellular accumulation of fully labeled asparagine and undetected labeled aspartate) following the administration of CB-839-containing treatments (Figure S4), likely to maintain proper rates of protein biosynthesis (Pavlova et al., 2018a).

Figure 3

Untargeted metabolomics profiling by high-content DIMS demonstrates individual drug metabolic modulations that are further affected by the combined administration with CB-839

(A and B) Relative metabolite levels measured following individual or combined treatment in (A) HMVP2 and (B) DU145 PCa cells show that prominent metabolic modulations are dictated by the presence or absence of CB-839, as well as additional important metabolic modulations determined by the library drugs. Metabolites are clustered by Euclidean distance, and treatments are clustered by Spearman’s correlation.

In addition to the untargeted metabolomics results mentioned above, the SIT-DIMS-based secondary screening datasets acquired on the cell samples labeled with an isotope tracer provide a wealth of additional information on the drug mode-of-action. The number of identified metabolites (identified by both tMS2 and mass accuracy, and including, for some of the metabolites, multiple isotopically enriched forms) in samples prepared with either 13C6-glucose or 13C5, 15N2-glutamine labeling equaled 117 and 180, respectively, in DU145 and 82 and 138, respectively, for HMVP2 cells (Table S3).

The fractions of different isotopic forms (M0, M+1, …, M+6) of several representative metabolites, were depicted in Figure 4A, to illustrate the disparate metabolic responses elicited by the individual and combined treatments. In agreement with the unlabeled data, a general decrease in the degree of isotopic enrichment in glutaminolysis and TCA cycle metabolites is observed for all samples treated with the glutaminase inhibitor CB-839 (either alone or in combination). Interestingly, cells treated with CB-839 do not appear to resort to an increased incorporation of glucose carbons into glycolysis and TCA cycle intermediates (Figure 4A), thus the increased isotopic enrichment of glucose-6-phosphate (G6P) in most of the samples exposed to CB-839 likely results in an increased flux of glucose to the pentose phosphate pathway (PPP), for which further analysis is needed. Notably, the addition of CB-839 to treatment with the HSP-90 inhibitor AUY922 does not appear to alter an already much diminished incorporation of glucose into any of the depicted metabolites and related pathways. Moreover, the 13C5,15N2-glutamine-derived enrichment of the oncometabolite 2-hydroxyglutarate (2HG) was undetected in several samples receiving the combined treatment, but not the corresponding individual treatments (Figure 4A), possibly indicating the inhibition of the biosynthesis of this tumorigenesis promoter by some, but not all, of the combined treatments.

Figure 4

SIT-DIMS analysis combined with computational data mining reveals the metabolic modulations important for drug activity

(A) Fractional isotope enrichment forms of representative metabolites resulting from 13C5,15N2-glutamine or 13C6-glucose-traced SIT-DIMS analysis of HMVP2 PCa cells show distinct levels of incorporation of labeled substrates following drug individual and combined treatments. Isotopes were labeled in different colors, and white indicates that no isotopes were detected.

(B) The Cohen’s effect size (d) was calculated by comparing the data-mined Markov chains of metabolic fluxes’ ratios of different reactions over CS between selected (individual and combined) treatment groups and the control group in HMVP2 cells. d > 0.8: large difference, 0.5 < d < 0.8: medium difference, 0.2 < d < 0.5: small difference. Asp: aspartate; Glu: glutamate; Suc: succinate; 2HG: 2-hydroxyglutarate; Mal: malate; Cit: citrate; G6P: glucose 6-phosphate; Ala: alanine; Ctrl: control; fw: forward reaction; rv: reverse reaction; ex: excretion; in: ingestion; antero: cytosol to mitochondria; retro: mitochondria to cytosol; SubsGln: glutamine feeding, and more details about the reaction annotations can be found in Table S6.

Given the wealth of metabolic information included in these datasets, we applied a recently published data mining method (Matsuda et al., 2020) to quantitatively pinpoint any significant changes of metabolic fluxes among the different individual and combined treatment conditions. The measured mass isotopomer distribution vectors (MDVs) (Wittmann and Heinzle, 1999) of a few representative metabolites are sufficient to perform this kind of isotopomer analysis. In line with our previously described results, the metabolic modulations related to the CB-839-induced glutaminase inhibition (regardless of whether CB-839 was administered alone or in combination) were the most prominent, as shown by the large Cohen’s effect size of the flux ratios of log([SubsGln (glutamine feeding)]/[CS]), which generally (though not always) had absolute values greater than 0.8, especially for combination treatments (Figures 4B and S5). Notably, while all of the library drugs had absolute d values of log([SubsGln (glutamine feeding)]/[CS]) generally close to zero, the combination treatments’ values showed strong drug-dependent modulations. For instance, docetaxel further enhanced the CB-839-driven change, while the combination of Panobinostat reduced the effect of the individual CB-839 administration, and Temsirolimus had virtually no effect (Figure 4B). Importantly, other notable, but milder, metabolic flux modulations (0.2 < absolute d value <0.5, representing smaller differences) are highlighted by this analysis for the different treatments, contributing to the disparate treatment efficacies and outcomes. Similar to d of log([SubsGln (glutamine feeding)]/[CS]) above, although to a more moderate extent and in the opposite direction, the administration of CB-839 affected d of log([PC (pyruvate carboxylase)]/[CS]) and the observed changes were affected differently by the library drugs (Figure 4B). Moreover, additional small changes that are not strictly CB-839-driven can also be observed (e.g. the decreased log([AKGDH (α-ketoglutarate dehydrogenase)]/[CS]) following temsirolimus + CB-839; Figure 4B). While it is beyond the scope of this report to look in detail into the molecular modulations induced by all the drug treatments, these results highlight the additional advantage of potentially using the metabolic data not only for evaluating synergy (as discussed in the section below) but also to investigate the metabolic reprogramming cells undergo following treatment that can potentially be further targeted.

Secondary multiple readout SIT-DIMS screen: validation of the top-hit drug combinations via PCA-based Euclidean distance synergy quantification (PEDS) analysis

Currently, there are no established statistical methods to quantify drug synergism from a multi-readout dataset. We developed a PCA-based algorithm to quantify synergy and prioritize the most promising synergistic combinations based on the phenotype-informed dataset generated from the secondary metabolomics screening. All the detectable features (dependent variables) obtained after filtering and normalization (details in the STAR Methods section) (Kirwan et al., 2014) of the dataset acquired on the unlabeled samples were simplified into independent variables (principal components, PCs) using PCA. Given that each PC is orthogonal to all the others, we applied the Euclidean distance in the n-dimensional space to measure the synergistic effect, using PC scores as coordinates and proportions of variances as weights to the axes (Equation 1). Drug-induced metabolic changes entailed discrepancies among groups in the multidimensional space. The two vectors between the control and each of the two individual treatment groups were used to predict the hypothetical additive drug combination effect (the subtrahend of Equation 1). Thus, synergy can be evaluated by comparing the actual (the minuend of Equation 1) and predicted weighted Euclidean distances between the control and combination groups. The nature of this approach is similar to the Response additivity model, where synergy was assessed, on univariate measurements, by comparing the observed combination effect to an expected effect (Foucquier and Guedj, 2015), while expanding it to a multivariate level.

In Equation 1, denotes the average PC scores from replicates in the different experimental treatment groups (ctrl: control; drug1 and drug2: individual drug treatments; comb: combined treatment with drugs 1 and 2), and is the proportion of variance for ith PC. Drug combinations with higher than zero are considered synergistic.

Before applying it to data from our metabolomic screen, we tested the PEDS algorithm using available omics datasets. To the best of our knowledge, there are no publicly available metabolomics datasets acquired on samples treated with synergistic drugs. However, we tested our algorithm using available transcriptomic data from drug combination studies (source: ArrayExpress (Athar et al., 2018)) with validated drug synergies (Paroni et al., 2012; Walsby et al., 2014; Girnun et al., 2007). PEDS values from the three studies were all higher than zero (Figure S6), thus confirming the synergy between the drugs and the value of the simple PEDS equation to evaluate treatment synergy. We therefore proceeded to apply PEDS to our metabolomics datasets acquired on unlabeled samples. The PEDS values were calculated for each of the 20 drug pairs and for both cell types (Figure 5A). The PEDS values for ten of the drugs (i.e. Afatinib, ABT-751, Vandetanib, 17-DMAG HCl, Obatoclax Mesylate, AZD8055, Daunorubicin HCl, Motesanib Diphosphate, Temsirolimus, and PHA-665752) were lower than zero for both cells lines, indicating that, based on the secondary multivariate analysis, their combination with CB-839 is not synergistic. For other drugs (including Rapamycin, AUY922, Idarubicin HCl, GSK461364, Panobinostat, TP-0903, BIIB021, and Nocodazole), the PEDS values were higher than zero, indicating synergy, but only for one of the two tested cell lines. Only the combinations of docetaxel and trichostatin A with CB-839 were synergistic (PEDS values higher than zero, Figure 5A) for both cell lines, based on the multivariate DIMS data. While trichostatin A resulted in higher PEDS values than docetaxel for DU145 cells, we chose to further investigate the combination of CB-839 and docetaxel in consideration of the fact that docetaxel is one of the most frequently prescribed chemotherapy compounds across lines of care for PCa (Wen et al., 2019). To visualize the magnitude of the metabolic differences emerging following the individual or the combined treatments, we performed a principal component analysis on the unlabeled DIMS data, and the scores plots (HMVP2: Figure 5B, DU145: Figure 5C) indicate clear separations among different treatment groups.

Figure 5

Application of PEDS to the analysis of DIMS datasets reveals the synergy of the combinations of docetaxel or trichostatin A with CB-839 in both HMVP2 and DU145 PCa cells

(A–C) (A) PEDS values were calculated from the DIMS-unlabeled metabolomic profiles for the 20 top-hit drug combinations administered to HMVP2 (green line) and DU145 (orange line) PCa cells. PEDS values greater than 0 indicate synergy. PCA score plots obtained from the untargeted DIMS metabolomics-based screening data for (B) HMVP2 and (C) DU145 PCa cells treated with docetaxel (DX) and/or CB-839 (CB) indicate clear group separation. Circles represent sample replicates in each group; stars represent the center of each group.

Validation of the synergy between docetaxel and CB-839 in a wide range of relative concentrations

To further validate the synergy between docetaxel and CB-839 in a range of different concentrations and relative drug ratios, we used two foundational and widely accepted methods to evaluate synergy, the Bliss independence and the Loewe additivity models. The synergy scores for the Bliss (Figures 6A and 6B) and the Loewe (Figures 6C and 6D) models were calculated using the measured relative cell viability in PCa cells grown as spheroids (3D culture) following treatment at multiple drug doses (HMVP2: 0, 5, 10, 25 nM CB-839 with 0, 2.5, 5, 10, 25 nM docetaxel; DU145: 0, 100, 250, 500, 1000 nM CB-839 with 0, 1, 2.5, 5 nM docetaxel). In all cases, the δ-scores had high excess response and confirmed the synergy between CB-839 and docetaxel in both cell lines due to drug interactions in ideal concentration windows. Similar outcomes were also observed when cells were cultured in a monolayer (2D, Figure S7).

Figure 6

The combination of docetaxel and CB-839 is synergistic in a wide range of concentrations in both HMVP2 and DU145 PCa cells and according to both the Bliss independence and the Loewe additivity models

(A–D) 3D δ-score synergy maps generated (using SynergyFinder) from ATP bioluminescence assay data based on the Bliss independence (A: HMVP2, B: DU145) and Loewe additivity models (C: HMVP2, D: DU145) on PCa cells cultured in 3D and following treatment with serial doses of docetaxel and CB-839 for 48 h. A δ-score higher than zero represents synergy.

To determine whether the drug combination induced synergistic apoptosis, an apoptosis assay for HMVP2 and DU145 cells in a 2D model was applied. Notably, cells treated with the drug combination exhibited significantly higher apoptotic rates compared to either single drugs or the control (HMVP2: Figure S7E, DU145: Figure S7F).

In vivo validation of the synergy between docetaxel and CB-839

To further validate the synergy between docetaxel and CB-839, the individual and combined treatments were tested in vivo in HMVP2 PCa allograft tumors. The combination of CB-839 and docetaxel produced a significant reduction of HMVP2 PCa tumor growth whereas neither compound produced a significant inhibition of tumor growth when administered alone at the doses used (Figure 7A). Notably, the combination of CB-839 and docetaxel produced a significant synergistic effect on tumor growth, as indicated by the Bliss index values (Table S4). There were no apparent signs of toxicity with the combination and there were no significant changes in body weight throughout the experiment in any of the treatment groups (Figure 7B).

Figure 7

Synergistic reduction in the growth of HMVP2 allograft tumors treated with the combination of CB-839 and docetaxel

(A) Growth curve of HMVP2 allograft tumors in syngeneic FVB/N male mice. Data indicate mean ± SEM of both flank tumors in mice with treatment of vehicle control (n = 14), CB-839 10 mg/kg 3x/week (n = 14), docetaxel 20 mg/kg 1x/week (n = 16), or CB-839 10 mg/kg 3x/week + docetaxel 20 mg/kg 1x/week (n = 16). Two-way, repeated-measure ANOVA, followed by Bonferroni’s multiple comparison test. ∗∗∗∗p < 0.0001 compared to control; $ and #, p < 0.05 compared to CB-839 or docetaxel alone, respectively.

(B) Average mice body weight in each treatment group over the treatment period.

In-depth MFA-based validation of the metabolic modulations detected by SIT-DIMS

The SIT-DIMS-based secondary metabolomics screening in the 96-well format represents a quick, inexpensive, and powerful approach that, especially when combined with data mining, provides important information underlying the drug molecular mechanism as well as identifies additional targetable metabolic pathways.

To validate the results of the secondary screen and investigate more in-depth the metabolic modulation of the combination of docetaxel and CB-839, metabolic profiling and flux analysis experiments were performed on PCa spheroids treated with either individual or combined docetaxel and CB-839 and grown in media containing either 1,2-13C2-glucose, to examine glucose-derived carbon incorporation in the glycolysis, PPP and related pathways, or 13C5,15N2-glutamine, to perform an in-depth flux analysis of glutamine and its downstream pathways. Intracellular extracts were analyzed using a high-resolution UHPLC-MS platform.

Incorporation of 1,2-13C2-glucose resulted in significantly increased enrichment of the glucose/fructose-6-phosphate M+2 fraction (Figure S8) in samples exposed to CB-839, thereby highlighting the increase in glucose uptake induced by CB-839 treatment (to a similar extent regardless of the presence or absence of docetaxel). However, the downstream glyceraldehyde 3-phosphate M+2 fraction levels remained unchanged regardless of what treatment the cells were exposed to. Concurrently, the fractions of erythrose-4-phosphate M+2 and sedoheptulose-7-phosphate M+4 increased significantly in cells exposed to treatments containing CB-839, indicating the carbon diversion from glycolysis to the non-oxidative PPP. In addition to these CB-839-driven modulations, the inhibition of the oxidative PPP (as indicated by the significantly decreased fractions of ribose-5-phosphate M+1) was observed in samples receiving docetaxel-containing treatments.

As expected, following the CB-839-induced glutaminase inhibition, the glutamine M+5 enriched fraction significantly increased. Moreover, the combinatorial treatment (when compared to the single-agent treatments) significantly decreased the enriched fractions of succinate, fumarate, malate, and citrate, suggesting that the energy production from the oxidative TCA cycle decreased (Figure 8). Even in the presence of important modulations in TCA cycle flux following glutaminase inhibition, no significant differences were observed in the relative contributions of oxidative vs reductive glutaminolysis (unchanged ratio of citrate and α-KG as well as the ratio between M+3 and M+4 forms of fumarate and malate and citrate M+5 and M+4).

Figure 8

Treatment with combined docetaxel and CB-839 in HMVP2 PCa cells results in inhibition of TCA cycle turnover as well as de novo nucleotide biosynthesis

Fractional isotopic enrichment (following labeling with 13C5,15N2-glutamine) of representative metabolites (carbon and nitrogen enrichment are shown for TCA cycle intermediates and nucleotides, respectively) in HMVP2 PCa cells spheroids following treatment with docetaxel (DX) and/or CB-839 (CB) for 24 h. Computational data mining highlighted several reactions with large Cohen’s effect sizes (d) indicated by the blue (d < −0.8) or red (d > 0.8) metabolic reaction arrows, representing a decreased/increased contribution from that reaction to citrate synthesis. Only reactions with large effect sizes are color-coded here. Glc: glucose; Pyr: pyruvate; PDH: pyruvate dehydrogenase; PC: pyruvate carboxylase; AcCoA: acetyl CoA; Lac: lactate; LacEx: extracellular lactate; Cit: citrate; α-KG: alpha-ketoglutarate; Suc: succinate; Fum: fumarate; Mal: malate; OAA: oxaloacetate; Glu: glutamate; Gln: glutamine; Asp: aspartate; AMP: adenosine monophosphate; ADP: adenosine diphosphate; ATP: adenosine triphosphate; NADH: reduced form of nicotinamide adenine dinucleotide; C: control; DX: docetaxel; CB: CB-839.

With the MDV of these TCA cycle metabolites and glutamate as inputs, the isotopomer data were mined using the same method used for the SIT-DIMS data and flux ratios of each reaction to CS and the Cohen’s effect sizes between control and docetaxel + CB-839 groups were calculated (Figure 8: blue lines: d < −0.8, red lines: d > 0.8). With these high-resolution data, more significant differences can be detected. The reduced contribution of glutamine and increased contribution from pyruvate by PC validated our SIT-DIMS data mining results. In addition, the flux ratios of log([MDHreverse (reverse reaction of malate dehydrogenase)]/[CS]) and log([FHreverse (reverse reaction of fumarate hydratase)]/[CS]) increased, suggesting a diversion of oxaloacetate from citrate biosynthesis. Moreover, the flux ratios of log([LDH (lactate dehydrogenase)]/[CS]) and log([Lac_ex (lactate excretion)]/[CS]) decreased in the docetaxel + CB-839 group which indicates the short supply of pyruvate for energy production (all results are shown in Table S5).

Interestingly, in addition to validating the results of the SIT-DIMS analysis, the in-depth MFA indicated that the nitrogen incorporation of nucleotides in the combination group significantly decreased compared to single-drug treatments, including metabolites containing adenosine (Figure 8, S9A), guanosine, inosine (Figure S9B), uridine, and cytidine (Figure S9C).

Discussion

The most widely utilized methods for assessing a drug pair’s performance have been vastly based on synergy metrics of univariate drug efficacy (Meyer et al., 2019). The complexity of cellular regulatory pathways and the mechanisms of action of drugs result in enriched biological modulations and suggest that new drug combinations can be better designed with the help of systems biology and omics-informed approaches. In this study, we proposed an unbiased metabolomics screening platform to generate phenotype-informed datasets and developed an algorithm, called PEDS, to quantify synergy and prioritize the most promising drug combinations.

We started by performing a comprehensive single-agent and combination drug screen against four PCa and three NP cell lines. This analysis yielded 20 candidate synergistic drug combinations. Next, we conducted a further in-depth exploration of the 20 selected drug combinations using a metabolic multi-readout, SIT-DIMS-based secondary screen. To improve the biological relevance of the screen, PCa cell spheroids were used as a model for the prostate tumors in the phenotypic secondary screening. While we acknowledge that this model still only represents a highly limited reconstruction of the native prostatic heterogeneity and complex in vivo architecture (Ellem et al., 2014), it represents a first step toward increasing the complexity of models used for drug screening purposes. To our knowledge, this is the first metabolomic drug screen of mammalian cell models for the discovery of synergies of drug combinations. The proposed workflow can be readily complemented by novel tools for high-throughput datasets (Lubbock et al., 2021). Moreover, we capitalized on the ability of the DIMS platform to reliably profile thousands of samples in a relatively short time, thereby meeting the requirements of high-throughput technology required for screening purposes and concurrently producing potentially more informative, multivariate phenomics data for synergy evaluation. In addition to the quantification of drug synergy, this phenotypic screening platform sheds light on the predictive biomarkers that indicate treatment efficacy. By examining the untargeted and isotope-tracer metabolic profiles from all the different treatments, we found pronounced dysregulated metabolites that could potentially be targeted to prevent the development of resistance and enhance drug potency. In parallel, the SIT-DIMS platform proved to be a powerful tool to answer controversial questions efficiently. For instance, the results of the secondary screening indicated increased intracellular accumulation of asparagine in CB-839-treated prostate cell lines. Alkan et al. (2018) reported that the mitochondrial aspartate exporter compensates for glutaminase inhibition in murine muscle and lung cancer cells. Biancur et al. (2017) found that asparagine synthetase was significantly increased after CB-839 treatment in a mouse pancreatic tumor model. Using 13C4,15N2-asparagine as the tracer, we were able to readily confirm, using the SIT-DIMS platform, that asparagine uptake was promoted upon glutaminase inhibition, in agreement with other reports (Pavlova et al., 2018b).

Given the lack of established methods to evaluate drug synergy based on the high content multi-readout data from both single and combinatorial conditions, we developed a new algorithm, PEDS, to quantify drug synergy from untargeted metabolomics profiling (or other -omics data), with the underlying rationale that drug combinations can achieve a greater magnitude of metabolic modulation than the expected sum of individual activities. PEDS showed robustness in confirming previously validated synergistic drug pairs with publicly available omics datasets (Paroni et al., 2012; Walsby et al., 2014; Girnun et al., 2007). The application of PEDS to our DIMS metabolomics datasets resulted in the selection of the combination of CB-839 and docetaxel as synergistic. The synergy was then further validated in vitro by conventional methods, as well as in vivo. With the aid of a new data mining method for isotopomer analysis, the metabolic reprogramming induced by the combined treatment was narrowed down to collaborative inhibition of glutamine and a complementary modulation of glucose metabolism, resulting in major disruptions of the TCA cycle, oxidative PPP, and de novo nucleotide biosynthesis.

Notably, docetaxel is among the top three treatments across all lines of therapy for patients with prostate cancer (Wen et al., 2019). Additionally, CB-839 is in multiple clinical trials. Thus, from the pharmacological perspective, the combination of CB-839 and docetaxel warrants further preclinical investigations.

In summary, we describe and validate a new metabolomics-based phenotypic screening platform which, in combination with our PEDS algorithm, represents a powerful tool for the discovery of drug synergism. This is a highly versatile and customizable methodology, as it is applicable to many other assays for primary screening or other omics datasets. Founded on the comprehensive profiling of metabolic modulations, this platform serves as a spur for the evolution of increasingly impactful and efficacious search methods for drug combinations.

Limitations of the study

Despite the advantages, PEDS has potential limitations when the effect of the drug combination (compared to control) is opposite to the theoretical sum of the effect of the individual agents (which is an uncommon, but not impossible circumstance). This is because the PEDS equation only takes into consideration the magnitude, but not the direction of the effect vectors in the PCA space. Future studies might lead to a more generalized PEDS equation able to correctly evaluate synergy under these uncommon conditions. Besides, there could be covariance between the generated orthogonal components that may be of importance for synergy, which the method would naturally not be able to capture.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for materials should be directed to and will be fulfilled by the lead contact, Dr. Stefano Tiziani (tiziani@utexas.edu).

Materials availability

This study did not generate new unique reagents.

Method details

Tissue culture

Prostate cancer cell lines, PC3, DU145, 22Rv1, and HMVP2 (established from ventral prostate tumors of one-year-old Hi-Myc mice (Saha et al., 2016)) were grown in RPMI 1640 supplemented with 10% FBS and 1% glutamine (Mahajan et al., 2017; Russell and Kingsley, 2003). NMVP (Saha et al., 2017), a nontumorigenic murine prostate epithelial cell line, were cultured in DMEM/F12 medium containing 10% FBS, BPE (28 μg/mL), EGF (10 ng/mL), insulin (8 μg/mL), gentamicin (80 μg/mL), and transferrin (5 μg/mL). Human normal prostate cell lines PWR-1E and RWPE-1 were cultured in Keratinocyte SFM. All cell lines were routinely checked to ensure that they were mycoplasma free and maintained in optimal growth at 37°C and 5% CO2.

Cell viability assay

Single agent drug screening assays were carried out on PC3, DU145, 22Rv1, HMVP2, NMVP, PWR-1E, and RWPE-1 in 384-well white square flat-bottom plates. Cells were seeded at 2,500 per well and allowed to adjust/adhere to the environment overnight, and then treated for 24 h using single agent, from a drug library (Cambridge cancer compound library [CCCL] and additional cherry-picked compounds amounting to a total of 292 drugs with a wide range of drug targets [Figure S10]) at doses of 100, 250, 500 nM, and 1 μM with two replicates. Combinatorial drug screening assays were tested on PC3, DU145, 22Rv1, HMVP2, NMVP, and RWPE-1 (2,500 cells per well in 384-well plates, 4 replicates) treated with the drug library at 250 nM with or without 100 nM CB-839 for 24 h.

The validation of the docetaxel + CB-839 synergy was performed in both 2D and 3D culture systems. For cells cultured in monolayer (2D), HMVP2 or DU145 cells were seeded in a 96-well white opaque plate at 2,000 cells/well and incubated overnight. PCa cell spheroids (3D) were generated seeding 20,000 cells/well in a 96-well cell-repellent plate until spheroid formation occurs. Cells were then treated with a matrix of 4 × 5 concentrations of the two drugs for 48 h in duplicates. Specifically, HMVP2 in 3D culture were treated by 0, 5, 10, 25 nM CB-839 with 0, 2.5, 5, 10, 25 nM docetaxel. DU145 in 3D culture were treated by 0, 100, 250, 500, 1000 nM CB-839 with 0, 1, 2.5, 5 nM docetaxel. HMVP2 in 2D culture were treated by 0, 25, 50, 100 nM CB-839 with 0, 1, 2.5, 5, 10 nM docetaxel. DU145 in 2D culture were treated by 0, 100, 250, 500, 1000 nM CB-839 with 0, 0.1, 0.25, 0.5 nM docetaxel. Drug concentrations were adjusted compared to the 24-h screening experiments to account for aspects including the number of cells seeded, growth rate of each cell type and previously tested responses to 48 h treatments with the single agents. In all cases, the final concentration of DMSO was kept at 0.1%.

After treatment, CellTiter-Glo 2.0 cell viability assay and CellTiter-Glo 3D cell viability assay was used, according to manufacturer’s instructions, to test cells in 2D and 3D culture, respectively. Relative cell viability (normalized to control) was calculated based on the intensity of ATP bioluminescence (proportional to the number of viable cells) (Lodi et al., 2017; Jansson et al., 2018; Karvonen et al., 2019; Crouch et al., 1993). Luminescence was recorded by Spark® Cyto (TECAN, Männedorf, Switzerland). Luminescence readings were normalized to the levels of DMSO-treated control samples to calculate cell viability in drug-treated samples.

To evaluate drug synergies between CB-839 and the drug library compounds, the 2.5th percentile of Bliss Index(Bliss, 1939) were calculated using bootstrapping by resampling the relative cell viabilities of the single and combination treatments from biological replicates. The drug library compounds with the 2.5th percentile of Bliss index values greater than zero (indicating synergy) were selected as drug candidates(Österroos et al., 2016). Hierarchical clustering with Spearman’s rank correlation was used to generate the heatmap of 2.5th percentile of Bliss Index values.

Cell viability data for synergy validation in a range of different drug combination doses were analyzed with SynergyFinder(Ianevski et al., 2020) to generate the 3D synergy maps according to the Bliss and Loewe methods.

Apoptosis assay

Cells were seeded in 10 cm dishes (DU145:1,000,000 cells/dish, HMVP2: 400,000 cells/dish due to their faster proliferation rate; in triplicates), incubated overnight, and then treated for 48 h. Cells were prepared with CF®488A Annexin V and PI apoptosis kit following the manufacturer’s instructions. Stained cells were analyzed using a flow cytometer (BD FACSAria II, San Jose, CA). Data were analyzed using FlowJo v10.6.1 Software. Quadrants were positioned based on the corresponding unstained samples. Student’s t test (two-tailed) was used to evaluate the significance of apoptotic cell fractions between control and treatment groups.

Stable isotope tracer direct infusion mass spectrometry metabolic profiling

For the SIT-DIMS experiments, cell spheroids were generated by adding magnetic nanoparticle-incorporated cells (HMVP2 or DU145, at a density of 100,000 cells/well) to cell-repellent microplates and leaving the plate on a magnet drive for 30 min to center the cells. The plate was then incubated overnight to complete spheroid formation and achieve strong cell-cell interactions. After the overnight incubation, spheroids were washed with 300 μL 1X PBS. A magnetic bead extractor (in a 96-well format) was used to move and wash spheroids simultaneously. The washed spheroids were then placed in either regular medium (unlabeled) or media containing either 13C6-glucose or 13C5,15N2-glutamine. Cells were treated immediately with the 20 drugs selected from the primary screen alone or in combination with CB-839 and 0.1% DMSO as the solvent control with five replicates in each group. The same procedure was followed to prepare the samples labeled with asparagine using medium containing 13C4,15N2-asparagine. Cells were treated immediately with CB-839 and docetaxel, either alone or in combination, 0.1% DMSO as the solvent control, and there were five replicates in each group.

After 24-h of treatment, the spheroids were washed with 300 μL ice cold 1X saline(Lu et al., 2019), and transferred to 96-well pre-washed (to eliminate chemical residues) filter plates with 200 μL ice cold 1X saline added to each well. The filter plates were centrifuged at 1,000 g for 3 min at 4°C to remove the saline while retaining the spheroids on the filters. Samples were extracted by adding 50 μL of ice cold 2:2:1 acetonitrile:methanol:water spiked with 14 internal standards to each well (extraction blanks were also included) and shaking the filter plates at 700 rpm for 2 min. The filter plates were centrifuged at 100 g for 1 min at 4°C and the extracts (filtrates) were collected in 96-well PCR plates. Samples were kept on ice as much as possible throughout the extraction. Small aliquots were removed from each sample, mixed and aliquoted to generate QC samples, with one QC solution generated by pooling all samples and 3 QC solutions each containing aliquots of unlabeled or 13C6-glucose or 13C5,15N2-glutamine labeled samples. When all samples, QCs and blanks were added, the plates were then heat-sealed for 2 s and stored in −80°C until data acquisition.

The DIMS-based data acquisition was performed as follows (Lu et al., 2019). The sample plates were thawed on ice and then centrifuged at 2,500 g for 10 min at 4°C immediately prior to the direct infusion. The analyses were performed on a Q Exactive Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Scientific, San Jose, CA) equipped with an automated chip-based nanoelectrospray ionization (nESI) source (Triversa NanoMate, Advion, Ithaca, NY). Nanoelectrospray conditions comprised 30 s acquisition time per ionization mode, 5 μL injection volume, 0.4 psi gas pressure, and 1.3 kV voltage controlled by ChipSoft software (version 8.3.3; Advion). MS was operated switching between positive and negative full-MS modes under the following parameters: capillary temperature, 250°C; microscans, 1; automatic gain control (AGC) target, 1e6; mass resolution, 70,000; and mass-to-charge ratio (m/z) range, 50–750. Calibration of the mass spectrometer for both ionization modes corresponding to the highest m/z range was achieved through the routine commercial calibration solutions provided by the manufacturer. In addition, customized calibrations were carried out at m/z 50–750 mass range using the following: for negative ionization mode, 87.00877 (pyruvic acid); 117.01624 (fumaric-d2 acid); 149.06471 (glutamic-d3 acid); 208.11399 (tryptophan-d5 indole); 265.14790 (sodium dodecyl sulfate) and 514.288441 (sodium taurocholate); for positive ionization mode, 74.09643 (n-butylamine), 138.06619 (caffeine fragment), 195.08765 (caffeine) and 524.26496 (MRFA).

Blanks were run before and after the entire acquisition. Sample sequence for injection was randomized, with one QC sample per corresponding labeling condition injected every six samples, and one total QC sample injected every twelve samples. Tandem mass spectrometry (tMS2) data were acquired on the unlabeled and differently labeled QC samples with 3 replicates for metabolite identification. Sample plates were kept at 4o C during the entire sequence.

In total for the SIT-DIMS metabolic screen, 6 blanks, 342 QC samples and 1260 cell samples were acquired. The raw files were converted to mzXML files by MSConvert(Kessner et al., 2008), and processed by custom written MATLAB scripts (available to academic laboratories upon request). Data from unlabeled, 13C6-glucose labeled, and 13C5, 15N2-glutamine labeled samples for the different cell types were processed separately. Metabolite assignment was achieved by comparing MS1 m/z of all samples to the Human Metabolome Database (Wishart et al., 2018) (filtered by detected and endogenous categories, 1387 different masses) and the tMS2 spectrum of QC samples to a library generated in-house (tMS2 spectrum for 369 metabolites including 14 IS, Table S6), within 5 ppm. According to previously published guidelines (Kirwan et al., 2014; Southam et al., 2016), the data processing workflow includes: (1) blank filtering: candidates are removed if their intensities in the sample are less than three times those of the extraction blanks; (2) QC filtering: assigned features with coefficient of variance (CV) higher than 25% in QCs are filtered out; (3) sample filtering: peaks are retained if they are present in >50% of samples; (4) replicate filtering: only peaks that are present in at least four out of five of the biological replicate analyses of the samples are retained; (5) Total spectral area normalization. Student’s t test (two-tailed) was used to evaluate the significance of metabolite intensities between different groups.

HMVP2 allograft tumor study

Male FBV/N mice, 6-7 weeks old, were purchased from Charles River Laboratories (Raleigh, NC) and allowed to acclimate for 1 week before the start of the tumor experiment. Mice were maintained on a semi-purified diet (AIN93M; 10kCal% fat; Custom Animal Diets, LLC) and water ad libitum. The animal protocol was approved by the Institutional Animal Care and Use Committee of the University of Texas at Austin. Mice were injected subcutaneously with 5 × 106 HMVP2 cells into both flanks. When tumors were palpable, mice were divided into 4 groups with approximately equal average tumor volumes and then treated intraperitoneally with i) vehicle control, ii) CB-839 alone (10 mg/kg, 3/week), iii) docetaxel alone (20 mg/kg, 1x/week) or iv) the combination of CB-839 (10 mg/kg, 3/week) and docetaxel (20 mg/kg, 1x/week). Tumor size was measured 1-2 times/week using a digital caliper and tumor volume was calculated by using the formula: 0.5236×D1×(D2)2, where D1 and D2 are the long and short diameter, respectively. Body weight of the mice was measured weekly. Experiments were continued until tumor sizes in the control group reached their maximum limit, as specified by the protocol.

High-resolution mass spectrometry-based metabolic flux analysis

HMVP2 cells were seeded in ultra-low attachment flasks at 500,000 cells/mL to generate spheroids for 24 h. Spheroids were then washed with PBS and added to 1,2-13C2-glucose labeled or 13C5,15N2-glutamine labeled or unlabeled fresh medium, and treated with 500 nM docetaxel and/or 100 nM CB-839 for 24 h with 3 replicates per group. Sample extraction was performed as previously reported (Sweeney et al., 2020). Briefly, metabolites were extracted by 1:1 methanol:water with 10 mM ammonium bicarbonate and equal parts of chloroform. Butylated hydroxytoluene (BHT) was added to the extraction buffer to preserve metabolites susceptible to oxidation. Both polar and apolar fractions were analyzed (Sweeney et al., 2020; Pandey et al., 2018) on a Q Exactive Hybrid Quadrupole-Orbitrap Mass Spectrometer equipped with an electrospray source, connected with a Vanquish UHPLC system (Thermo Scientific, Waltham, MA). In short, polar samples in 1:1 methanol:water were first analyzed using a SeQuant ZIC-HILIC 3.5 μm, 100 Å, 150 × 2.1 mm PEEK coated HPLC column (Millipore Sigma, Burlington, MA) with isobaric separation, and then dried and resuspended in water and analyzed using a Synergi 4 μm Hydro-RP 80 Å, 150 × 2 mm HPLC column (Phenomenex, Torrance, CA) with gradient separation. Apolar samples were dried and resuspended in 95:5 ethanol:6N HCl and then analyzed using a Synergi 4 μm Hydro-RP 80 Å, 150 × 2 mm HPLC column (Phenomenex, Torrance, CA) with isobaric separation. The raw files were processed using SIEVE 2.2.0 SP2 (Thermo Fisher Scientific) and an in-house script that operates in the MATLAB programming environment. Metabolite identifications were achieved by matching accurate masses and retention times to a library of standards generated in house (Lu et al., 2017). A pooled quality control was used to monitor instrument stability and a blank was used for background subtraction. Peaks were included in analysis if the coefficient of variation (CV) was less than 25% in the QC replicates and the signal to noise ratio was higher than 3. Student’s t test (two-tailed) was used to evaluate the significance of isotope fractions between different groups.

Quantification and statistical analysis

Data mining for isotopomer analysis

Data mining analysis was performed according to a recently published method (Matsuda et al., 2020). The metabolic model includes 53 reactions and 34 metabolites covering glycolysis, PPP, TCA cycle (oxidative and reductive), and glutaminolysis (Table S7). In line with previous studies and for simplicity, intracellular compartments were not considered. Glucose uptake was arbitrarily set as constant 100, while other fluxes were free with a lower and upper boundary. Isotope fractions (after natural abundance correction) with standard deviations served as the data mining input (SIT-DIMS data: Table S8, MFA data: Table S9). All analyses were performed using mfapy 0.5.8, where the metabolic model was used as a function to simulate the MDV for a given metabolite flux distribution, and the goal is to minimize the residual sum of square between the experimental and simulated fluxes by Metropolis-Hastings algorithm. For each condition, 8 rounds of 5,000,000 distributions were generated, where the initial 2,500,000 steps were discarded because of the burn-in process and 2,500 steps were sampled from the second half by every 1,000 steps. Thus, in total, a dataset of 20,000 steps was used to represent that condition. By calculating the flux ratios with a Log10 transform, metabolic reprogramming can be identified when the dataset distributions are distinct from each other between two conditions. Because of the large dataset, Cohen’s effect size was used to determine the difference between two populations (with absolute number at 0.2: small effect size, 0.5: medium effect size, and 0.8: large effect size).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Experimental models: Cell lines
Human: PC3	American type culture collection (ATCC; Manassas, VA)	CRL-1435
Human: DU145	ATCC (Manassas, VA)	HTB-81
Human: 22Rv1	ATCC (Manassas, VA)	CRL-2505
Human: LNCaP	ATCC (Manassas, VA)	CRL-1740
Human: RWPE-1	ATCC (Manassas, VA)	CRL-11609
Human: PWR-1E	ATCC (Manassas, VA)	CRL-11611
Mouse: HMVP2, NMVP	Laboratory of John DiGiovanni (Saha et al., 2016; Saha et al., 2017)	N/A
Chemicals, peptides, and recombinant proteins
All solvents and chemicals for mass spectrometry	Thermo Fisher Scientific (Waltham, MA)	LC/MS grade
Bovine pituitary extract (BPE)	VWR (Radnor, PA)	AAJ64417
Transferrin	Sigma-Aldrich (St. Louis, MO)	T8158
RPMI 1640	Thermo Fisher Scientific (Waltham, MA)	SH30096.01
DMEM/F12	Thermo Fisher Scientific (Waltham, MA)	SH30271.FS
Keratinocyte serum-free medium (SFM)	Thermo Fisher Scientific (Waltham, MA)	17005042
Fetal bovine serum (FBS)	Thermo Fisher Scientific (Waltham, MA)	SH3007103HI
Glutamine	Thermo Fisher Scientific (Waltham, MA)	SH30034.01
Epidermal growth factor (EGF)	Thermo Fisher Scientific (Waltham, MA)	PHG0311
Insulin	Thermo Fisher Scientific (Waltham, MA)	12585014
Gentamicin	Thermo Fisher Scientific (Waltham, MA)	15-710-064
Dimethyl sulfoxide (DMSO)	Thermo Fisher Scientific (Waltham, MA)	AC295522500
NaCl	Thermo Fisher Scientific (Waltham, MA)	S271-1
Phosphate buffered saline (PBS)	Thermo Fisher Scientific (Waltham, MA)	SH30028.02
Drug library	Selleckchem (Houston, TX)	L2300 & Cherry Picking
CB-839	Selleckchem (Houston, TX)	S7655
Docetaxel	Cayman Chemical (Ann Arbor, MI)	11637
CellTiter-Glo 2.0 cell viability assay	Promega (Madison, WI)	G9243
CellTiter-Glo 3D Cell Viability Assay	Promega (Madison, WI)	PAG9683
CF®488A Annexin V and PI apoptosis kit	Biotium (Fremont, CA)	30061
NanoShuttle	Greiner Bio-One	657846
1,2-13C2-glucose	Cambridge Isotope Laboratories (Tewksbury, MA)	CLM-504-1
13C4,15N2-asparagine	Cambridge Isotope Laboratories (Tewksbury, MA)	CNLM-3819-H-0.25
13C5,15N2-glutamine	Cambridge Isotope Laboratories (Tewksbury, MA)	CNLM-1275-H-PK
13C6-glucose	Cambridge Isotope Laboratories (Tewksbury, MA)	CLM-1396-10
Materials
384-well plates	Thermo Fisher Scientific (Waltham, MA)	264705
96-well plates	Thermo Fisher Scientific (Waltham, MA)	136101
96-well PCR plates	Thermo Fisher Scientific (Waltham, MA)	3482P
Heat sealing foil	Thermo Fisher Scientific (Waltham, MA)	AB-0745
Cell culture dishes	Thermo Fisher Scientific (Waltham, MA)	171099
Cell culture flasks	Thermo Fisher Scientific (Waltham, MA)	10-126-34
96-well cell-repellent plates	Greiner Bio-One	655976
96-well magnetic drive	Greiner Bio-One	655830
96-well magnetic bead extractor	V&P Scientific (San Diego, CA)	VP 407AM-N1
Ultra-low attachment flasks	VWR (Radnor, PA)	89089
96-well filter plates	Pall Corporation	8084