Metabolomics-Driven Exploration of the Chemical Drug Space to Predict Combination Antimicrobial Therapies.

Alternative to the conventional search for single-target, single-compound treatments, combination therapies can open entirely new opportunities to fight antibiotic resistance. However, combinatorial complexity prohibits experimental testing of drug combinations on a large scale, and methods to rationally design combination therapies are lagging behind. Here, we developed a combined experimental-computational approach to predict drug-drug interactions using high-throughput metabolomics. The approach was tested on 1,279 pharmacologically diverse drugs applied to the gram-negative bacterium Escherichia coli. Combining our metabolic profiling of drug response with previously generated metabolic and chemogenomic profiles of 3,807 single-gene deletion strains revealed an unexpectedly large space of inhibited gene functions and enabled rational design of drug combinations. This approach is applicable to other therapeutic areas and can unveil unprecedented insights into drug tolerance, side effects, and repurposing. The compendium of drug-associated metabolome profiles is available at https://zampierigroup.shinyapps.io/EcoPrestMet, providing a valuable resource for the microbiological and pharmacological communities.

Introduction

Despite modern phenotypic screening technologies to find growth-inhibitory compounds at high throughput (Brown and Wright, 2016, Zheng et al., 2013), it is still extremely hard to discover radically new antibacterial classes that truly enlarge the arsenal of antimicrobial compounds in the fight against resistant pathogens (Katz and Baltz, 2016, Payne et al., 2007). The high rate of rediscovery of known drug modes of action (MoAs) might be a consequence of the type of synthetic compound libraries screened (Demain, 2014, Helfrich et al., 2018, Shen, 2015) and possibly reflects a limited space of essential cellular processes that can be inhibited by a single molecule (Albert et al., 2000). One alternative solution to the conventional search for single-target compound treatments is to expand the search for novel antibacterials to combination therapies (Wang et al., 2016, Gonzales et al., 2015, Stokes et al., 2017). However, due to the explosion in the combinatorial complexity, combination discovery remains a major challenge. Furthermore, classical antimicrobial combination screens focus on discovering synergistic pairs of already known antibiotics (Stokes et al., 2017), thereby ignoring that nonantibiotic molecules can synergistically interact and become potent antimicrobial alternatives (Ejim et al., 2011, Andersson et al., 2017, Brochado et al., 2018).

To go beyond combinations of currently available antibiotics, we exploit high-throughput metabolomics (Zampieri, 2018) to monitor the metabolic response of Escherichia coli to a library of 1,279 chemical compounds (Prestwick Library), most of which are human-targeted drugs that have little if no antimicrobial activity (Maier et al., 2018). By combining the newly generated drug metabolome profiles with previously published compendia of metabolic (Fuhrer et al., 2017) and fitness (Nichols et al., 2011) profiles in E. coli gene-knockout mutants, we make de novo predictions of drug MoAs and systematically predict epistatic drug interactions. We show that high-throughput metabolic profiling of bacterial response to small molecules can expand the search for new antimicrobial treatments to compounds with no growth-inhibitory activity in vitro. The proposed methodology is applicable to other biological systems and therapeutic fields and can open new opportunities in drug discovery and aid drug target deconvolution, assessment of side effects, and prediction of drug-drug interactions.

Results

Multidimensional Metabolic Profiling of Drug Response

Phenotypic drug screens only select lead compounds that can kill bacteria or inhibit their growth (Zheng et al., 2013). However, nonlethal compounds with dispensable targets could become invaluable tools in combination therapies (Brochado et al., 2018). To systematically identify compounds that interfere with nonessential functions in gram-negative bacteria, we exposed exponentially growing E. coli cultures to a library of 1,279 chemically diverse compounds (i.e., Prestwick Chemical Library). This library includes US Food and Drug Administration (FDA)-approved drugs for diverse therapeutic purposes, ranging from treatment of infectious diseases to cancer and cardiovascular pathologies (Figure 1A). Only 11% of the compounds are antibiotics, while the majority are human-targeted drugs. Individual compounds were administered at a single concentration of 100 μM in 96 deep-well plate cultivations, and the metabolome response was monitored by flow injection analysis in a time of flight mass spectrometer (FIA-TOFMS) 2 h after drug exposure (Zampieri et al., 2018) (Figure 1B). In parallel, the optical density of treated cultures was monitored up to 6 h after drug exposure (Figures 1B and S1). This workflow enabled rapid profiling of relative changes in the abundance of ∼39,000 ions, out of which 969 could be putatively annotated as deprotonated metabolites. In total, we monitored metabolic changes across 1,279 perturbed conditions and DMSO treatments as vehicle controls in three biological replicates.

Figure 1

Metabolic Profiling of the Drug Response

(A) Distribution (pie chart) of Prestwick chemical compounds across therapeutic classes.

(B) Illustration of the metabolic drug profiling workflow. Growth is monitored using a plate reader up to 6 h after treatment, while metabolomics samples are collected after 2 h of treatment and analyzed by FIA-TOFMS (Fuhrer et al., 2011).

(C) Inner pie chart shows the distribution of compounds inhibitory activity. Outer pie charts illustrate the number of compounds with at least one (green) significant change (absolute Z score ≥3 and p value ≤ 1e−5) and more than 20 (blue) significant affected ions. The percentage of drugs exhibiting a metabolic phenotype is estimated on (1) annotated ions, (2) detected ions common to metabolome profiles of E. coli knockout strains (Fuhrer et al., 2017), and (3) totality of detected ions.

(D) For each class of therapeutic agents (Table S1), we report the distribution of growth rates relative to the untreated DMSO condition and number of responsive metabolites (absolute Z score ≥3 and p value ≤ 1e−5). For each therapeutic class, the tops and bottoms of each box are the 25th and 75th percentiles, respectively, while the red line in the middle of each box is the samples median. The lines extending above and below each box are the whiskers. Whiskers extend from the ends of the boxes delimited by the interquartile to the largest and smallest observations excluding outliers (red crosses). Outliers have values that are more than three scaled median absolute deviations.

To estimate drug-induced metabolic changes, raw mass spectrometry data were normalized by correcting for instrumental and systematic biases (Zampieri et al., 2018). To account for the confounding effect of different growth inhibitions across treatments, we employ a non-parametric smoothing function that for each metabolite normalizes relative changes in concentrations to corresponding changes in growth rate (Figure S1). Finally, a Z score normalization was applied on the growth-rate-corrected metabolic profiles before estimating average and SD over the three biological replicates (Table S1; see STAR Methods for full details).

Out of the 1,279 drugs, only 15% exhibited antimicrobial activity (i.e., inhibited growth more than 20% with respect to the untreated condition). As expected, most of these compounds are known antibiotics (Figures 1C, 1D, and S1). However, the vast majority (70%) induced a strong metabolic response causing significant changes in the abundance of at least one metabolite (absolute Z score ≥3 and p value ≤ 1e−5 Bonferroni-adjusted threshold) (Figure 1C). For each individual perturbation, on average, 10 metabolites underwent significant changes (absolute Z score ≥ 3 and p value ≤ 1e−5) (Figures 1D and S1), and most metabolites (92%) exhibited a significant response in at least one of the treatments (Figure S1). Altogether, we found an unexpectedly large spectrum of metabolic effects in E. coli and revealed drugs potentially targeting nonessential bacterial processes. Hence, we next evaluated which bacterial functions were affected by nonantibiotic drugs.

Mapping Drug Effects on Microbial Metabolism

To identify metabolic genes that are potential direct drug targets or that indirectly mediate drug response, we systematically mapped enzymes to drug-induced metabolic changes using a genome-scale metabolic network of E. coli (Orth et al., 2011). We assessed whether largest metabolic changes are located in the proximity of an enzyme by using a weighted scoring function (i.e., locality score) (Zampieri et al., 2018) in which metabolic changes are weighted by their respective network distance to the tested enzyme (i.e., length of the shortest path) (Figures S1 and S2; Table S2). Enzymes associated with large metabolic changes in their proximity (p value ≤ 0.001) were tested for overrepresentation against KEGG (Kyoto Encyclopedia of Genes and Genomes) metabolic pathways. Overall, we found 70 different metabolic pathways enriched for local metabolic changes (q-value ≤ 0.001) in at least one treatment, hinting at the diversity of metabolic responses across the drug library.

While distinguishing direct from indirect metabolic changes remains challenging (Kwon et al., 2008), we systematically show that the nine drugs with non-promiscuous enzyme targets in bacteria (e.g., allopurinol; Figure S1) induce largest metabolic changes in the vicinity of known drug targets (Figures S2 and S3). This result suggests that local metabolic changes can reveal the most direct drug effects on metabolism. Surprisingly, in many cases, local metabolic changes were detected also for drugs with known molecular targets in humans that are absent in E. coli (Figure S1). Typical examples are lovastatin, a commonly used statin for lowering cholesterol, and docetaxel, a cancer chemotherapy agent. Lovastatin induced strong accumulation of intermediates in steroid degradation pathway (Figure S3), possibly hinting at the ability of E. coli to degrade the drug, similarly to other members of the human gut microbiome (Yoo et al., 2014). Docetaxel is an antineoplastic agent shown to hyper-stabilize microtubules (Dumontet and Jordan, 2010) and increase the anti-cancer activity of inhibitors of fatty acid synthase (Menendez et al., 2004). Although cytoplasmic tubules have been observed in bacteria, no homologies to eukaryotic microtubules have so far been established in E. coli (Pilhofer et al., 2011). Nevertheless, in E. coli, this anticancer compound induces a distinct metabolic response in the β-oxidation of fatty acids (Figures 2A and 2B), a pathway enabling fatty acids to be used as carbon source. Because the pathway is transcriptionally repressed in glucose minimal medium (Pauli et al., 1974), we hypothesized that metabolic changes report on the role of docetaxel as a direct or indirect transcriptional activator of pathway genes. To validate our hypothesis, we used a GFP reporter plasmid (Zaslaver et al., 2004) to monitor the promoter activity of the fadAB operon (Figure 2C). Consistent with our metabolome-based predictions, we observed a mild (27%) but significant (p value = 3.7e−05) increase in fadB promoter activity upon treatment of E. coli cells with 2 mM docetaxel. To test the physiological significance of detected changes in fadB promoter activity, we incubated wild-type E. coli cells in glucose M9 medium with and without 2 mM docetaxel for 9 h and then transferred growing bacteria to M9 medium with the long-chain fatty acid oleate as sole carbon source. In agreement with our metabolome-based prediction, we observed that cells preincubated with docetaxel exhibited faster growth rate when forced to switch to oleate consumption than E. coli cells preincubated in the absence of the drug (Figure 2D).

Figure 2

Locality Score Analysis

(A) A volcano plot of the 966 annotated metabolites in response to docetaxel.

(B) A volcano plot of the locality analysis results in docetaxel-treated cells.

(C) FadB promoter activity in M9 minimal medium with 2 mM docetaxel or DMSO control or in the presence of oleate as sole carbon source (∗∗∗p value < 1e−3). As a negative control, we monitored the effect of docetaxel on recA promoter activity using the same pUA66 plasmid (Zaslaver et al., 2004) (brown bars), mean ± SD across three biological replicates.

(D) Optical density measurements of an isogenic strain of E. coli growing in M9 medium with oleate as sole carbon source (0.2%). Cells previously grown for 9 h in glucose M9 medium with (black) and without docetaxel (green) are compared. Thick line and shaded region represent the average and SD estimated across three biological replicates. Inset histogram represents the doubling times estimated from the two growth curves (∗∗p value < 1e−2).

(E) Distribution of p values estimated from locality score analysis of drug induced metabolic changes in the proximity of enolase enzyme (Eno) for compounds in five major classes of antibiotics.

(F) Disk diffusion assays testing the susceptibility to ofloxacin (7 μg) in wild-type E. coli (WT) and two knockout mutants, ΔaspC (negative control) and ΔrhlB. Bar plot represents mean ± SD across four biological replicates.

(G) Distribution of ΔrhlB fitness scores (Nichols et al., 2011) in response to 16 DNA-damaging agents (green boxplot) and 66 chemical perturbations with diverse MoAs (gray boxplot).

Local metabolic changes can also uncover enzymes that mediate a tolerance response to drug treatment. To test this hypothesis, we systematically searched for enzymes that consistently associated to local metabolic changes in response to antibiotics. Here, we considered five major antibiotic classes: protein, RNA, folic acid, cell wall biosynthesis, and DNA replication inhibitors (Table S2). Enzymes significantly associated (p value ≤ 1e−5) to inhibitors of cell wall biosynthesis were enriched in nucleotide biosynthesis, while inhibitors of DNA replication associated to enzymes in pentose phosphate pathways and glycolysis (Table S2). Consistent with our hypothesis, we found an almost exclusive enrichment of β-lactam antibiotics among drugs inducing local changes in the proximity of ribose-phosphate diphosphokinase (Prs) (p value < 1e−16), an enzyme that mediates tolerance and facilitates emergence of resistance to ampicillin (Fridman et al., 2014). Interestingly, quinolone antibiotics exhibited most significant changes in proximity of the glycolytic enzyme enolase (Eno) (Figure 2E) (p value < 1e−21). This association possibly relates the function of enolase in regulating RNA degradosome and cell filamentation (Murashko and Lin-Chao, 2017) to the observation that E. coli mutants that are unable to filament can survive longer upon quinolone treatments (Piddock and Walters, 1992). To test a potential link of RNA degradosome activity with tolerance to quinolones we measured ofloxacin susceptibility of an E. coli mutant strain deleted for the nonessential degradosome component RhlB (Figure 2F). Consistent with our metabolome-based prediction, disk diffusion assays showed an increased susceptibility of ΔrhlB to ofloxacin. This result is in agreement with previous data (Nichols et al., 2011) in which ΔrhlB exhibited significantly reduced fitness to 16 different DNA-damaging agents (p value = 5.5e−4) (Figure 2G). Taken together, the characterization of metabolic changes induced by diverse pharmacological compounds unraveled the potential to selectively interfere with bacterial metabolic processes and opens up a large spectrum of potentially druggable cellular processes in gram-negative bacteria.

Metabolome-Based Grouping of Chemically Diverse Compounds

Next, we aimed to group compounds inducing similar metabolic changes. The comparative analysis of metabolome responses can be a powerful method to find clues about drug MoAs (Allman et al., 2016, Zampieri et al., 2018) and systematically explore modes of metabolic interference of nonantibiotic compounds (Zampieri, 2018). To account for the possibility that drugs with similar MoAs have different permeability or target affinity, and hence induce different magnitude of metabolic changes after the same time of drug exposure, we developed an alternative metric to estimate pairwise drug similarity. The herein defined similarity metric uses an iterative thresholding scheme in which for a given Z score threshold (thr), the metabolic effect of drug A is described by a ternary vector of increased (+1), decreased (−1), and unchanged (0) metabolites . For each pair of drugs A and B, we exhaustively searched for the combination of thresholds (thr and thr) that maximizes the dot product between ternary drug vectors, and , divided by the total number of metabolites that were found to be changed in at least one of the two drugs (Figure 3A). The significance of the percentage of metabolites that exhibit consistent changes between the two drugs is estimated using a hypergeometric test. To benchmark our approach, we tested the ability of different similarity metrics to group known antibiotics with similar MoAs. Remarkably, the herein proposed metric outperforms correlation (i.e., Spearman correlation) or dependency (i.e., mutual information) measures of similarity (area under the receiver operating characteristic [ROC] curve >0.8; Figures 3B and S3). The iterative thresholding approach was hence used in an affinity propagation algorithm (Frey and Dueck, 2007) to group drugs in 164 clusters exhibiting distinct metabolic profiles (Table S1; Figure S3). Notably, clustering of metabolome profiles could group antibiotics with similar MoAs regardless of whether they target proteins with metabolic or other functions and in spite of different growth inhibitory activities (Figure S3).

Figure 3

Metabolome-Based Similarity Analysis

(A) Schematic representation of the computational framework used to group drugs eliciting similar metabolic responses. For each pair of drugs A and B, the combination of thresholds on the respective drug-induced metabolic changes (i.e., Z scores) that corresponds to the largest similarity score is selected. Subsequently, an affinity propagation algorithm is used to group drugs on the basis of their metabolome pairwise similarity.

(B) Area under the receiver-operating characteristic curves (AUCs) measuring the ability to discriminate antibiotics with similar MoAs using Spearman correlation (SP), mutual information (MI), and the herein-presented iterative score (IS). Only antibiotics that belong to cell wall, folic acid, DNA replication, protein, and RNA synthesis inhibitors are considered here.

(C) Similarity heatmap of the cluster containing ascorbic acid. Elements on the diagonal indicate the corresponding drug therapeutic class reported in Figure 1D.

(D) Differential responsive metabolites characteristic of the drug cluster of (C) (ANOVA p value ≤ 1e−5). Green and blue boxplots illustrate the distribution of relative metabolic changes among drugs in the cluster and all remaining treatments, respectively.

(E) A volcano plot of the 966 annotated metabolites in response to ascorbic acid.

(F) Uric acid conversion rate was normalized to the activity measured with H2O vehicle only. All tested compounds were dissolved in DMSO (1 mM kanamycin, ascorbic acid, and allopurinol) (∗∗p value ≤ 0.01). Values represent mean ± SD across three replicates.

Most of the identified clusters consist of nonantibiotic drugs, potentially unraveling classes of molecules impinging on cellular processes different from those affected by classical antibiotics. Ascorbic acid (vitamin C) belongs to such a group of compounds and recently gained traction for its potential versatile role as anti-cancer (Yun et al., 2015) and anti-tuberculosis (Vilchèze et al., 2013) agent. Ascorbic acid is an essential nutrient for some mammals (Hodges et al., 1969) and a reducing agent that at high dosage can increase the production of ferrous iron from ferric ions and drive the overproduction of reactive oxygen species (Chen et al., 2008). Our metabolome-based similarity analysis found ascorbic acid to group with largely diverse therapeutic agents (Figure 3C). Because most of these compounds have citrate as an excipient (Figure S1), we deduced that the common metabolome response (Figure 3D) mainly reflects similar iron chelating properties and downstream regulatory effects between citrate and ascorbic acid (Vilchèze et al., 2013, Frawley et al., 2013, Imlay, 2013). Consistent with this hypothesis, common metabolic changes that drive drugs to cluster together (Figure 3D) mainly localized in the tricarboxylic acid cycle (TCA) cycle and pathways allosterically regulated by citrate, such as degradation of galactarate and glyoxylate (Blumenthal and Jepson, 1966, Seo et al., 2014, Imlay, 2013). Different from other compounds in the cluster, however, ascorbic acid induces a marked increase in the levels of hypoxanthine (Figure 3E). To test the hypothesis that ascorbic acid also acts as a xanthine oxidase inhibitor, we performed in vitro enzyme assays measuring the activity of E. coli xanthine oxidases XdhA-B-C with and without ascorbic acid. Indeed, we found ascorbic acid to directly inhibit xanthine oxidases (Figure 3F), potentially linking the toxicity of high dosages of ascorbic acid to its combined action as an oxidative stress agent (Chen et al., 2008, Vilchèze et al., 2013) and inhibitor of oxidases. These results further support the potential of nontargeted metabolomics in finding key characteristics of drug MoAs, bridging the gap between drug-induced phenotype and drug mechanism of action.

From Drug Metabolic Fingerprint to Drug-Gene Associations

Next, we compared our compendium of drug metabolome profiles with previously generated metabolic profiles in 3,807 individual E. coli gene knockouts (Fuhrer et al., 2017) in order to link drug metabolic phenotypes directly to drug targets or indirectly to genes mediating drug response (Subramanian et al., 2017). Because metabolome profiles of gene-knockout mutants were collected in complex glucose medium, we tested the comparability of drug-induced metabolic changes between minimal and complex media for nine antibiotics with four different MoAs (Zampieri et al., 2017). We found that all respective antibiotics are within the top 3% of Prestwick compounds with the highest similarity (Figure S4), suggesting that the characteristic metabolic response to inhibition of drug target is highly conserved across different conditions (e.g., nutrients). However, future studies with larger number of drugs and conditions are necessary to generalize this principle to compounds different from those tested here.

In light of these results, we employed the same metric used to compare the metabolome profiles between drugs and compared individual drug responses with metabolic changes in gene deletion mutants (Figure 4A). This comparison resulted in 10,227 (i.e., 0.16% of all possible pairwise combinations) significant similarities (p value ≤ 1e−5 and similarity score ≥0.15), reflecting a large space of drug-gene associations (Figures 4A and S5). Because antibiotic targets are essential genes, our approach could not be directly tested against known targets of antimicrobial compounds. However, gene knockouts that elicit significantly (p value ≤ 5e−5) similar metabolic changes to five major classes of antibiotics could be often related to the antibiotic MoA (Table S3). Emblematic examples are (1) the metabolic similarity between an E. coli mutant deficient in the DNA repair gene recA and quinolone antibiotics (Liu and Imlay, 2013) (p value = 9.72e−5), (2) deletion of the aminodeoxychorismate synthase (ΔpabA) and folate inhibitors (Chakraborty et al., 2013) (p value = 3.43e−5), or (3) the knockout mutant for sn-glycerol-3-phosphate dehydrogenase (ΔglpD) and cell wall biosynthesis inhibitors (Spoering et al., 2006) (p value = 1.74e−27). Unexpectedly, we observed that deletion of aspartate aminotransferase (aspC) features the most similar metabolic changes to inhibitors of bacterial DNA-dependent RNA polymerase (e.g., rifampicin) (Figures 4B and S5). Because resistance to these antibiotics is nearly always due to a point mutation in the β subunit of bacterial RNA polymerase (rpoB) (Campbell et al., 2001), and aspC has already been shown to be involved in the general stress response to antibiotics (Shan et al., 2015, Mathieu et al., 2016), we hypothesized AspC to play an indirect role in mediating the adaptive response to inhibitors of RNA synthesis. The significantly (p value ≤ 0.01) higher resistance of ΔaspC to rifampicin (Figures 4C and 4D) confirmed our hypothesis and suggested that E. coli may actively reduce AspC activity to endure exposure to this class of antibiotics, although the underlying mechanisms remain to be clarified.

Figure 4

Prediction of Drug-Perturbed Genes

(A) Circos plot (Krzywinski et al., 2009) shows the distribution of gene knockout-drug associations found among genes in six major biological functions.

(B) Distribution of p values estimated from calculation of the metabolic similarity between deletion of aspartate aminotransferase (aspC) and metabolic changes induce by drugs in five major classes of antibiotics.

(C) Growth-rate inhibition in liquid culture of different concentrations of rifampicin and ofloxacin in wild-type E. coli (WT) and ΔaspC with respect to the untreated condition (mean ± SD across the three replicates).

(D) Disk diffusion assays testing the susceptibility to rifampicin (7 μg) in wild-type E. coli (WT) and ΔaspC. Bar plot represents mean ± SD across four biological replicates.

(E) Volcano plot of the 966 annotated metabolites in response to disulfiram.

(F) Volcano plot of the 3,807 gene similarity profiles with disulfiram treatment.

(G) Impact of disulfiram on B. subtilis Icd in vitro activity (mean ± SD across three replicates). The measured NADPH conversion rate was normalized to the activity measured with H2O vehicle only. All tested compounds (kanamycin, disulfiram, and phosphoenolpyruvate) were dissolved in DMSO to a 1 mM concentration. ∗p value ≤ 0.05; ∗∗p value ≤ 0.01.

Overall, most recurrent associations were found between drugs and genes involved in carbon metabolism and signaling (Figures 4A and S5), suggesting new targets to hamper adaptation to external stressors (Fridman et al., 2014, Zampieri et al., 2017, Allison et al., 2011). Moreover, the predicted network of drug-perturbed genes provides experimentally testable hypotheses for in-depth mechanistic studies of tolerance or resistance mechanisms and drug MoAs. This comparative approach is orthogonal to the previously described analysis of locality in drug-induced metabolic changes, which is purely based on stoichiometry, and does not reveal whether observed changes correspond to an increased or decreased enzyme activity. By comparing drug with gene-deletion metabolome profiles, we not only expand our analysis to non-metabolic genes but also more directly associate drug changes to inhibition of gene functions. Overall, 54 drugs out of the 416 with at least one significant gene-knockout association showed similar predictions in the locality and gene-knockout-based analysis. One interesting example is disulfiram, a human-targeted drug used to treat chronic alcoholism by inhibiting the human enzyme acetaldehyde dehydrogenase (Figures 4E and S2). By comparing disulfiram induced metabolic changes to metabolic profiles of gene-knockout mutants, we found the largest similarity with the deletion of isocitrate dehydogrenase (Δicd) (Figure 4F). Here, we tested the hypothesis that Icd is directly inhibited by disulfiram using the commercially available purified Bacillus subtilis Icd protein. In vitro enzyme assays revealed a mild but significant (p value ≤ 0.05) inhibitory activity of disulfiram that at 1 mM causes 25% inhibition of Icd, similar to the Icd allosteric inhibitor phosphoenolpyruvate (Figure 4G) (Ogawa et al., 2007). Altogether, we demonstrated that our metabolome-based strategy is capable of finding genes that directly or indirectly mediate the metabolic response to drugs, suggesting new alternative ways to interfere with drug stress response.

Predicting Drug-Drug Interactions from Drug Metabolome Profiles

Compounds for which their combined inhibitory effect is greater (e.g., synergistic) than the one predicted on the basis of their individual effects are attractive alternatives to monotherapies, because they can optimize efficacy, decrease toxicity, and lower the risk for resistance evolution (Stokes et al., 2017, Zhao et al., 2013). Equally important is drug antagonism, where one drug may block or reduce the effectiveness of complex multiple drug regimes in patients with multimorbidity. To predict drug epistatic interactions, we integrated metabolome-based prediction of drug-gene associations with previously published chemogenomic data (Nichols et al., 2011). In this dataset, the sensitivity of all viable E. coli knockout strains to 82 chemical perturbations was quantified by differences in colony size (i.e., fitness score). Because drug-induced gene function inactivation can render E. coli more or less vulnerable to the action of a second drug, we predict that Prestwick drugs exhibiting a significant metabolome similarity (similarity score ≥ 0.15 and p value ≤ 1e−5) to a gene deletion can epistatically interact with drugs for which the fitness of the corresponding gene knockout strain largely differs from the average (Nichols et al., 2011) (Figure S5). It is worth noting that increased antibiotic susceptibility (i.e., negative fitness score) or resistance (i.e., positive fitness score) could reveal genes important for survival, for which expression is not regulated in response to the stress, and genes that are actively induced or repressed upon drug perturbation (Palmer et al., 2018). Our framework cannot resolve the potential conflicts emerging from drugs that individually trigger opposite (i.e., inducing or repressing) gene-regulatory responses (Bollenbach and Kishony, 2011) and hence whether drug-drug interference will result in a synergistic or antagonist interaction.

To systematically validate our approach, we used different thresholds on the fitness score of knockout mutant strains to distinguish susceptible or resistant strains to a given treatment, and for each threshold predicted drug epistatic interactions by identifying drugs and knockout strains that exhibit similar metabolic profiles. Next, metabolome-based predictions were compared to experimentally determined interaction scores previously measured for 465 pairs of drugs (Brochado et al., 2018). For drug pairs predicted to epistatically interact, we computed the average of experimentally determined interaction scores (Nichols et al., 2011) and compared it to drug pairs selected at random. Despite not including deletions of genes encoding for antibiotic targets in our predictions and that most experimentally tested drugs were antibiotics applied in a complex medium (Brochado et al., 2018), we found that the larger the difference in fitness of knockout mutants, the stronger (p value < 0.05) the average epistatic interactions between drug pairs (Figure 5A).

Figure 5

Linking Metabolomics to Chemogenomic Drug Profiling

(A) Systematic comparison of metabolome-based prediction of epistasis and interaction scores experimentally measured in (Brochado et al., 2018). The average of absolute drug interaction scores is calculated for drug pairs predicted to interact via knockout mutant strains with an increasing fitness score (dashed black line). Significance is estimated using a permutation test. For increasing thresholds on the knockout fitness scores, the red shaded area represents the mean ± SD of 1,000 randomization of predicted interacting drug-pairs.

(B) The volcano plot represents the similarity between all Prestwick drugs and ΔrpoS metabolic profiles.

(C) Optical density measurements of wild-type E. coli in LB medium with (green) and without (black) 500 μM benzydamine.

(D) Optical density measurements of wild-type E. coli in M9 glucose minimal medium with (green) and without (black) 500 μM benzydamine. Thick line is the average across three biological replicates. Shaded region represents average ± SD estimated over three biological replicates.

(E) Sensitivity profile of ΔrpoS against compounds tested in the chemogenomic screen (Nichols et al., 2011). Error bars represent the standard deviation of fitness score measured for different drug concentrations in (Nichols et al., 2011).

(F) Response matrix of doxycycline and benzydamine at different doses. For each combination of drug doses maximum growth rate is measured.

In the following, we selected a conservative threshold of −5; i.e., only knockout strains with a strong increase in drug susceptibility are selected to predict drug epistatic interactions. This analysis resulted in 1,032 pairs of potentially interacting drug pairs (i.e., 0.8% of all possible pairwise combinations) via drug-mediated inactivation of 139 genes (Figure S6; Table S3). As a proof-of-concept, we selected five drug pairs that have not been previously tested and measured the dose-response matrix for different compound concentrations. Drug epistasis is defined as the deviation of measured growth inhibition from the product of the single-drug response curves using the Bliss approximation (Bliss, 1939). To test prediction robustness, we also used an alternative metric of epistasis, the fractional inhibitory concentration index (FICi) (Odds, 2003). An FICi <1 indicates synergy, while an FICi >1 identifies antagonistic drug combinations.

Strongest predictions of drug interaction mediated by regulatory genes involve the action of regulators that play key roles in in vivo cellular adaptation, such as the sigma factor RpoS and the ferric uptake regulator Fur (Griffiths, 1993, Dong and Schellhorn, 2010) (Table S3). These predicted interactions can potentially suggest promising drug candidates to enhance activity of known antibiotics in an in vivo setting (Dong and Schellhorn, 2009, Hryckowian and Welch, 2013). Here, we tested the RpoS-mediated interaction between the anti-inflammatory drug benzydamine and the protein synthesis inhibitor doxycycline. RpoS deletion induces metabolic changes similar to those produced by benzydamine (Figure 5B), which exhibits a metabolic profile similar also to that produced by deletions of crp and pflA, involved in the response to nutrient and oxygen availability, respectively (Figure S7). Phenotypically, we observed that cells treated with benzydamine exhibit a premature entry into stationary phase and that this effect is more prominent when cells are grown in rich medium (Figures 5C and 5D). This observation is consistent with the role of RpoS in catabolism of fermentative products, such as acetate, and with previous experimental evidence of a higher expression of rpoS during entry in stationary phase in complex media (Rahman et al., 2006). In parallel, ΔrpoS is significantly more sensitive to doxycycline (Figure 5E). Consistent with our metabolome-based predictions, the growth-inhibitory activity of different dose combinations of doxycycline and benzydamine revealed strong synergism (FICi = 0.2) (Figures 5F and S7).

Another important transcriptional regulator predicted to mediate drug epistatic interactions is the ferric uptake regulator Fur (Griffiths, 1993). As a transcriptional repressor, fur deletion causes an upregulation of genes involved in iron homeostasis, siderophore-mediated uptake, and enterobactin biosynthesis (McHugh et al., 2003). Δfur was found to be more sensitive to two cell-staining compounds, propidium iodide and acriflavine, and two oxidative stress agents, streptonigrin and pyocyanin (Figures 6A and 6B). The increased susceptibility of Δfur to fluorescent intercalating and oxidative stress agents possibly links the toxicity of compounds like propidium iodide to oxidative stress, which is known to be severely accentuated when iron levels are high (Touati et al., 1995, Imlay, 2013). Among Prestwick compounds exhibiting similar metabolic changes to Δfur, we found many metal chelators (e.g., deferoxamine, pentetic acid, and clioquinol) (Figure 6A). In particular, we observed the most significant similarity between Δfur and deferoxamine mesylate, an iron-chelating compound used in the clinic to treat acute iron poisoning. This finding is consistent with the expectation that lowering iron availability induces a compensatory response reducing Fur activity (i.e., increasing iron uptake). Even at high deferoxamine concentrations (>500 μM), E. coli can compensate for the presence of deferoxamine and scavenge enough iron to sustain normal growth (Figure 6C). However, while reduced Fur activity should aggravate propidium iodide toxicity, we found that deferoxamine antagonizes propidium iodide toxicity (FICi = 1.7) (Figure 6C). This finding suggests that the two compounds trigger opposite cellular responses regulating iron homeostasis. To verify this hypothesis, we monitored the promoter activity of the Fur-repressed entCEBAH operon (Zaslaver et al., 2004). This operon encodes for enterobactin biosynthetic enzymes required when iron is limiting. As expected, we found that deferoxamine induces an overexpression of the operon (Figure 6D). Vice versa, in response to propidium iodide, cells actively reduce the transcription of the operon (i.e., increase Fur activity) (Figure 6D). Therefore, consistent with our hypothesis, propidium iodide and deferoxamine trigger opposite regulatory effects. Hence, the next question is how this regulatory conflict is resolved upon perturbation of cells with both compounds. We experimentally show that propidium iodide at high concentrations exerts a stronger control than deferoxamine on entC promoter activity (Figures 6D and S7). Therefore, the combined action of the two compounds results in the simultaneous reduction of free available iron due to deferoxamine and decreased expression of iron uptake genes, ultimately reducing the toxicity of propidium iodide. Notably, we showed that a multiplicative model fits EntC promoter activities (Figure S7) better than an additive model (i.e., R2 of 0.9 and 0.6, respectively) (Figure 6E). This observation suggests that Fur activity is independently regulated by different signals triggered by the two compounds, possibly including changes in growth rates (Scott et al., 2010).

Figure 6

Conflictual Response to Drug Combination

(A) The volcano plot represents the similarity between all Prestwick drugs and Δfur metabolic profiles (Fuhrer et al., 2017).

(B) Sensitivity profile of Δfur against compounds tested in the chemogenomic screen (Nichols et al., 2011).

(C) Response matrix of deferoxamine mesylate and propidium iodide at different doses.

(D) EntC promoter activity measured in cells treated with different doses of deferoxamine and propidium iodide.

(E) Probability density function of residuals from fitting analysis of EntC promoter activity upon treatment of cells with different combinations of deferoxamine and propidium iodide (PEntC(drug x, drug y)) (Figure S7). Residuals of a multiplicative model (PEntC(drug x, drug y) − [ PEntC(drug x) • PEntC(drug y)]) are compared to those of an additive model (PEntC(drug x, drug y) − [ PEntC(drug x) + PEntC(drug y)]). Median and SD of the two distributions are reported on the top of the histogram.

Next, since enzymes identified to underlie the predictions of drug-drug interactions were significantly enriched in nucleotide metabolism (p value ≤ 0.05), we tested three drug pairs involving the mediating action of enzymes in purine, pyrimidine, and folate biosynthesis, namely PyrE, PurL, and NudB. We validated the predicted synergistic interactions of the β-lactam antibiotic cefoxitin combined with the folic acid biosynthesis inhibitor sulfamethizole (FICi = 0.42) and the cell-wall synthesis inhibitor meropenem (FICi = 0.44), mediated respectively by enzymes in purine and pyrimidine de novo biosynthesis pathways, purL and pyrE (Figures S6 and S7). Notably, these two enzymes are essential in glucose minimal medium. In contrast, NudB is a dispensable enzyme upstream of folic acid biosynthesis. Here, we show that zidovudine, a structural analog of thymidine and a potent inhibitor of HIV replication, exhibits similar metabolic patterns to ΔnudB (Figure 7A), and as predicted, a strong synergistic interaction with sulfamethizole (FICi = 0.1) (Figures 7B, 7C, and S7). To test whether deletion of nudB would be sufficient to abolish the synergistic interaction between the two drugs, we repeated the growth assay experiment with a nudB deletion strain (Figure S6). We found that in ΔnudB, zidovudine minimum inhibitory concentration (MIC) was reduced ∼800-fold with respect to wild-type E. coli (Figure S6). However, upon sulfamethizole treatment, ΔnudB exhibits only a 2-fold reduction in MIC, but the lag-phase was strongly prolonged with respect to wild type (Figure S6). Most important, the synergism between the two drugs is largely reduced in a ΔnudB genetic background (FIC = 0.84), suggesting that NudB plays a key role in mediating the interaction between zidovudine and sulfamethizole. Hence, our methodology not only predicts epistatic drug interactions but also provides insights on the mechanisms underlying the interaction.

Figure 7

Experimental Validation of Predicted Drug Epistatic Interactions

(A) The volcano plot represents the similarity between all Prestwick drugs and ΔnudB metabolic profiles.

(B) Sensitivity profile of ΔnudB against compounds tested in the chemogenomic screen.

(C) Response matrix of sulfamethizole and zidovudine at different doses.

Discussion

Finding effective multidrug combinations against bacterial infections represents a pressing challenge. However, screening all possible combinations in large compound libraries is not feasible, because the number of experiments grows exponentially with the number of drugs and doses. Moreover, despite the large efforts in integrating high-throughput data in computational models (Cokol et al., 2011, Ejim et al., 2011), predicting drug epistatic interactions remains challenging (Bansal et al., 2014, Chandrasekaran et al., 2016, Facchetti et al., 2012). By scaling with the typical size of compound libraries screened in large phenotypic drug-discovery campaigns, our methodology complements other molecular profiling technologies (Subramanian et al., 2017, Maier et al., 2018) and offers the possibility to expand the search for potentially novel antimicrobial treatments to compounds that alone exhibit poor growth inhibitory activity. By comparing the metabolic signatures associated to drug and genetic perturbations, we identified drug-perturbed genes that combined with chemogenomic data enable prediction of epistatic drug interactions. Hence, our approach can be extended to largely diverse compounds and potentially lead to the discovery of nonantibiotic compounds that, when combined, exhibit promising antimicrobial properties. A key strength of our approach is that during the very early steps in the drug discovery pipeline, it can identify the MoA and mode of drug-drug interactions of compounds that could be optimized by traditional chemistry to improve selectivity and reduce toxicity. Follow-up experiments confirmed new MoAs in E. coli for three compounds (i.e., docetaxel, ascorbic acid, and disulfiram), the key role of aspartate aminotransferase and enolase in mediating the response to inhibitors of RNA synthesis and DNA replication, respectively, and five novel epistatic interactions between drug pairs.

Recent advances in sequencing technology coupled with random transposon mutagenesis can expedite the phenotypic screening of large mutant libraries upon drug treatments (Wang et al., 2011). In parallel, expanding the inventory of reference metabolic profiles in genetically perturbed bacteria to a larger set of relevant conditions (e.g., starvation) and more naturalistic media composition can increase the predictive power of in vivo relevant drug-perturbed genes. Because here we compared drug with genetic perturbations under different nutritional environments, it is likely that we underestimated the number of gene-drug associations and that the space of potential synergistic drug interactions is larger than we estimated. Moreover, technologies for systematic genetic perturbations, such as RNAi or targeted genome editing, can facilitate generating large compendia of reference metabolome profiles and expand the comparative analysis also to knockdowns of essential genes. Finally, integrating our metabolome-based approach with transcriptional drug profiling (Subramanian et al., 2017) and gene regulatory models (Ortmayr et al., 2018) could help to resolve potential conflicts between compounds acting on the same cellular process and predict whether drugs would synergize with or antagonize each other.

By mapping the landscape of metabolic responses in E. coli to a representative sample of the current drug chemical space, we unraveled an unexpectedly broad impact of general medications on bacterial metabolism. Understanding interference of human-targeted drugs with bacterial metabolism can guide drug repurposing, inform on drug metabolism, and open new ways by which nonantibiotic compounds can create selective pressures altering the overall gut microbiome composition (Maier et al., 2018), ultimately affecting human health and therapeutic outcome. Overall, by identifying drug metabolic fingerprints on a large scale, we provide a proof-of-principle of how to enlarge the space of potential combinatorial mechanisms by which we can inhibit bacterial growth and propose an alternative solution to the antibiotic discovery problem. We envisage that the proposed approach will become a broadly useful tool in many other therapeutic areas, including anti-cancer discovery (Ortmayr et al., 2018).

STAR★Methods

Key Resources Table

Contact for Reagent and Resource Sharing

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Mattia Zampieri (zampieri@imsb.biol.ethz.ch).

Experimental Model and Subject Details

Strain and growth conditions

E. coli BW25113 was used throughout this study. Culture medium consisted of standard M9 minimal media with (per liter): 5g of glucose, 7.52 g Na2HPO4·2H2O, 3 g KH2PO4, 0.5 g NaCl, 2.5 g (NH4)2SO4, 14.7 mg CaCl2·2H2O, 246.5 mg MgSO4·7H2O, 16.2 mg FeCl3·6H2O, 180 μg ZnSO4·7H2O, 120 μg CuCl2·2H2O, 120 μg MnSO4·H2O, 180 μg CoCl2·6H2O and 1 mg thiamine·HCl. Cells were grown at 37°C.

Method Details

Metabolome profiling of drug response

E. coli overnight cultures growing on M9 minimal medium were diluted to an initial Optical Density (OD600) of 0.05. 700 μL cell cultures were distributed in 96 deep well plates and cells were grown at 37°C until exponential phase and an OD600of 0.4 before exposure to the drug treatment by the addition of 7ul of a drug solution or DMSO. Here we selected the Prestwick library, consisting of 1280 FDA approved drugs prepared in DMSO solution. It is worth noting that the final drug concentration of 100 μM is close to the estimated concentrations of human-targeted drugs in the small intestine and colon (Maier et al., 2018). Samples for metabolomics profiling were taken after 2 hours from drug exposure. 50 μl of whole cell broth were directly transferred to 120 μl extraction liquid solution containing 50% (v/v) methanol and 50% (v/v) acetonitrile at −20°C. The extraction was carried out by incubating the samples for 1 hour at −20°C. Samples were centrifuged for 5 minutes at 4000 RPM and 80 μl of the supernatant was transferred to 96 well storage plates and stored at −80°C. All treatments were measured in 3 biological replicates. DMSO control treatments were included in every plate as control. Because of technical issues during measuring, the drug Chicago sky blue 6B was excluded from further analysis. Bacterial culture growth was estimated by measuring the optical density at 600 nm (OD600) in a plate reader. OD600measurements were acquired immediately before drug exposure and after 60, 120, 240 and 360 minutes. Growth rates were calculated as a slope of a linear fits to the logarithmically transformed growth curves (Table S1).

Mass spectrometry measurementand processing

The analysis was performed on a platform consisting of an Agilent Series 1100 LC pump coupled to a Gerstel MPS2 autosampler and an Agilent 6550 Series Quadrupole Time of Flight mass spectrometer (Agilent, Santa Clara, CA) as described previously (Fuhrer et al., 2011). Mass spectra were recorded from m/z 50 to 1000 using the highest resolving power (4 GHz HiRes). All steps of mass spectrometry data processing and analysis were performed with MATLAB (The Mathworks, Natick). For each acquired spectra, ion-peaks were identified with the findpeaks function of the Signal toolbox. Identified peaks with less than 5000 ion counts were excluded from further analysis. Detected ions were matched to a list of metabolites based on the corresponding molar mass (Zampieri et al., 2018). For the full list of metabolites used for annotation see Table S1. The chemical formula of each metabolite was used to calculate the deprotonated monoisotopic molecular mass. Detected ions within a mass tolerance difference of less than 0.003 Da were associated to the nearest reference metabolites. In total 969 different reference masses could be detected (see Table S1).

In vitro enzyme assays

To test in vitro the activity of Xanthine Oxidase, we purchased the purified enzymes XdhA/B/C from Sigma (i.e., X2252) and performed the reaction in 150 uL volume following protocol instructions. ∼0.05U/ml of enzyme dissolved in 50mM Tris/HCl was incubated with 500 μM of Hypoxanthine in 0.1M Tris/HCl buffer at 37°C. A plate reader was used to continuously measure Optical Density at 290 nm every 13 s to measure the production rate of uric acid within the first 5 min. Rate of conversion were estimated using a linear regression scheme. Similarly, to test in vitro the activity of isocitrate dehydrogenase, we purchased the purified Bacillus enzyme Icd from Megaenzyme (i.e., E-ICDHBS) and performed the reaction in 150 uL volume following protocol instructions. ∼0.31U/ml of enzyme dissolved in 50mM Tris/HCl was incubated with 1 mM of D-/L-Isocitric acid, 7.1 mM MgCl2 and 1 mM of NADP+, in 0.1M Tris/HCl buffer at 37°C. A plate reader was used to continuously measure Optical Density (OD340) at 340 nm every 13 s to measure the production rate of NADPH in the first 2 min after reaction start. Rate of conversion were estimated using a linear regression scheme. Notably, Isocitrate dehydrogenase from Bacillus subtilis has vast sequence (70% calculated using the online tool PDBeFold) and structural similarities (87%) with Isocitrate dehydrogenase from E. coli (EcIDH). The active sites of the enzymes are almost completely conserved and they have similar steady-state kinetic parameters for the dehydrogenase reaction. E.coli icd reported binding and catalytic sites were not only conserved, but also in matching structural positions when compared to B. subtilis analogous residues.

Disk diffusion assay

Filter paper circles of 6 mm diameter were impregnated with 7 μl of antibiotic solution (2.8 mM for oxfloxacin and 1.21 mM for rifampicin) and left to dry at room temperature. 106–107 bacteria at exponential phase were plated on M9 glucose minimal medium agar plate. The plate was left to dry before applying the disk in the middle and incubating for 12 hours at 37°C. The amount of antibiotic in each disk shown in the figures is mentioned in the legend.

Measurement of promoter activity

We used GFP transcriptional reporter plasmids in which the promoter regions of fadAB and entCEBAH are fused to green fluorescent protein (Zaslaver et al., 2004). Promoter activity was measured in M9 minimal medium with glucose or oleate as sole carbon source using a plate reader recording GFP intensity and optical density. GFP levels were normalized by estimating their proportionality with optical density during exponential phase (i.e., at steady state). Instead of reporting dynamic changes of promoter activity throughout growth as in (Zaslaver et al., 2004), a linear model in which GFP signal is proportional to OD600 (i.e., GFP = α• OD600+ ϒ) is fitted on experimental measurements in exponential phase using least square regression analysis. To account for different offsets caused by differences in background fluorescence across conditions we introduce an offset parameter ϒ. Estimates of α report on the actual gene promoter activity at steady state.

Quantification and Statistical Analysis

Metabolomics data normalization

Data preprocessing and normalization consists of mainly three steps: (i) elimination of artifacts caused by drug-compounds with an identical mass to reference metabolites, (ii) detection and correction of systematic variations, (iii) removal of dependencies with growth rate inhibition. Procedures are described below:Where is the measured intensity of ion m upon treatment c, represents biomass (OD600) at the time of extraction, and represent quadratic and linear dependency between ion intensity and OD600. In our model the variance associated to plate effects across the large number of MS injections (∼4000) is captured by the parameter . By assuming that overall, each ion concentration within the cell is “directly” affected only by few drugs, equivalently to assume sparsity of D, the proportionality factors α, β and γ could be determined by multiple least square fitting analysis performed on all collected samples, for each ion individually.

To detect ions associated to external perturbing agents, each compound stock was diluted in water (to reach a concentration of 10 μM), injected in the mass spectrometer and profiled in parallel to the cell extracts. For each Prestwick drug, detected ions with a measured intensity greater than three times the standard deviation estimated across all pure drug spectra were removed from further analysis. The importance of this preliminary step is nicely illustrated by allopurinol. Allopurinol has a molecular mass identical to that of hypoxanthine, which is one of the substrates of the inhibited reactions (xanthine dehydrogenase). Because our mass spectrometry platform cannot distinguish between metabolites with an identical mass to drug compounds, drug related ions are removed from analysis of the corresponding metabolic response (see also Figure S1).

Raw data normalization was performed by using a multiple linear regression approach. Detected ion counts are considered to be a linear combination between: extracted biomass (measured by OD600 (O)), plate to plate variance (P), drug effect (D), and white noise , using the following model:

The function fitlm using ordinary least-squares, was used for the regression analysis. The relative difference between the model-derived expected intensity and the corresponding measured ion intensity under each condition was used to estimate the compound effect D:Variance in the measurements (ε) was estimated from the standard error between replicate fold changes: .

To correct for the dependency between relative changes in metabolite abundance and growth rate inhibition (Figure S1), a nonparametric function was adopted. Specifically, locally weighted smoothing regression (LOWESS) based normalization with a window size of 20% of total data was performed on each individual ion profile, using the MATLAB function “malowess.” Briefly, for a given growth rate the expected relative abundance of an ion is estimated as the average of the relative ion intensity for proximal growth rate inhibition. Subsequently expected ion changes in correspondence of a given growth inhibition are subtracted from original measurements: . Where corresponds to the abundance of metabolite m in conditions with growth inhibition g, and represents the LOWESS fit to the dependency between relative abundance and growth inhibition (red line in Figure S1D). A z-score normalization was applied to finally estimate average and standard deviation over biological replicates, yielding for each metabolite an estimate of D and ε, respectively (Table S1 and https://zampierigroup.shinyapps.io/EcoPrestMet)

Pairwise drug distance metric

To compare the metabolome response between two different drugs we developed an iterative thresholding algorithm. For each pair of drugs A and B, we selected two thresholds, thr and thr between 2 and 6. The metabolic effects of drug A and B are described by a ternary vector of increased (+1, Z-scores ≥ thr), decreased (−1, Z-scores ≤ -thr) and unchanged (0, -thr ≤ Z-scores ≤ thr) metabolites . Only annotated ions were used in this comparative analysis. For each pair of drugs A and B, we exhaustively searched for the combination of thresholds (thr and thr) that maximizes the dot product between ternary drug vectors, and , divided by the total number of metabolites that were found to be changed in at least one of the two drugs:The significance of the percentage of metabolites that exhibit consistent changes between the two drugs is estimated using a hypergeometric test:Where consistent changes in both conditions represent observed successes (s), the number of changes in one condition represents the maximum possible number of successes (K), the total number of metabolites represents the total population (N), and the number of changes in the second condition represents the number of events or draws (n). P values were calculated in MATLAB using the hygecdf function. In parallel, pairwise similarity between drug responses were estimated using mutual information, and Spearman correlation. The 3 different metrics were empirically compared by ROC curves performances determined from the set of antimicrobial compounds in the Prestwick library with known MoAs (Figure S3). Overall, we observed the best performance with the iterative similarity metric proposed in this study. The resulting matrix of drug pairwise similarities is then used to group drugs eliciting similar metabolic adaptive changes using an affinity propagation algorithm (Frey and Dueck, 2007) (Table S1).

Mutant drug association analysis

Drug-induced metabolic changes were compared to E. coli gene-knockout strains metabolic profiles previously collected in (Fuhrer et al., 2017). To this end, the same iterative similarity metric presented in this study was used. Here, we compared Z-scores estimated in 1996 ions detected at the same mass-to-charge ratio (m/z) in the compendium of drug metabolome profiles and the previously published large compendium of gene-deletion metabolome profiles. Hence, ions between our compound screens were matched to those detected in the knockout dataset, regardless of whether they could be annotated to a reference metabolite or not, with a mass difference threshold of 0.003 Da (Table S3).

Locality score

A genome-scale network model of E. coli K12 model iJO1366 (Orth et al., 2011) metabolism was used to determine the distance between each enzyme-metabolite pair. The resulting pairwise distance matrix (D) between metabolic enzymes and metabolites was estimated by means of the minimum number of reactions separating the two in a non-directional network. All highly connected metabolites were removed prior to calculation. Specifically, we removed all metabolites that participate in more than 20 reactions, as small changes in the selected threshold do not affect the overall distance between enzyme and metabolites (Figure S3 and Table S2). To assess whether largest metabolic changes were statistically more probable in the proximity of the deleted enzymes, we used a locality scoring function, in which all metabolic changes are weighted by the respective distance to a tested enzyme, as follows:where for each enzyme e, a weighted mean over corresponding Z (metabolite m Z-score corresponding to drug d) is computed. Weights are functions of the inverse of the squared distance between the enzyme i and the metabolite m. We perform a permutation test by randomly shuffling the distance matrix D (D) K times (K = 10000). For each gene we ensure to preserve original degree of connectivity in the stoichiometric network (Table S2).

Synergistic pair prediction validation

E. coli cultures were grown in glucose M9 medium with combinations of different concentrations of 2 drugs. For each combination of drugs concentrations, growth was monitored by OD600measurements every 5 minutes for 48 hours using a plate reader and maximum growth rate was estimated. Drug epistatic interaction is estimated as the deviation of growth rate from the BLISS independence model (Bliss, 1939) .For a given concentration of drug x and y, is the expected growth inhibitory effect of the combination of the two compounds, while and represent the growth inhibitions of individual compounds, measured as a percentage of growth rate with respect to the untreated condition. By measuring the dose response curves to the individual compounds and their combinations, the expected growth inhibition under each combination can be computed, and synergy can be calculated as the deviation of measured () from expected growth rates: is the epistasis score: positive values indicate antagonism (i.e., the drug combination yields a milder growth inhibitory effect than expected), while negative values indicate synergism (i.e., the drug combination yields a stronger growth inhibitory effect than expected).

Fractional Inhibitory Concentration index

As an alternative to the BLISS scoring formula for drug-drug interactions, we calculated the Fractional Inhibitory Concentration index (FICi) (Odds, 2003):Where MIC stands for Minimum Inhibitory Concentration, a,b represent two drugs and a.combined stands for the MIC of drug a when combined with drug b, and b.combined for the MIC of drug b when combined with drug a. An FICi < 1 indicates that paired combinations of drugs can exert inhibitory effects that are more than the sum of their effects alone (synergy), while an FICi > 1 identify combinations of drug for which their combined effect is less than the sum of their effects alone (antagonism) (Odds, 2003). To cope with the potential error in MIC estimation using a finite number of drug dilutions, we identify a set of boundaries around which the individual MICs should lay, and for each drug pair we calculated three FICi indexes based on three selection criteria:

Worst case scenario – This scenario penalizes the synergy estimation the most. We assume the individual drug MICs to be the highest concentration tested, for which we still observe bacterial growth. This score causes to be underestimated.

Best case scenario – This scenario assumes the individual drug MICs to be the lowest concentration point tested for which complete growth inhibition was observed.

Average case scenario – We assume the true value of to be the average between the estimates from the worst and best-case scenarios.

For sulfamethizole, the only compound for which we did not have an estimate for the best-case scenario the MIC was assumed to be always the highest concentration used (10μM), which hence corresponds to the worst case scenario.

Prediction of epistatic drug-pairs

The procedure consists of three main steps.

Association matrix between Prestwick drugs and knockout strains: each Prestwick drug is associated to gene knockouts exhibiting a significant similar metabolic profile (p value ≤ 1e-5 and similarity score ≥ 0.15). To this end, drug metabolic profiles generated here are compared to knockout strains metabolic profiles generated in (Fuhrer et al., 2017).

Association matrix between knockout strains and chemical compounds: each knockout strain is associated to drugs which exhibit a stronger or milder inhibitory activity (i.e., fitness score), based on the data reported in (Nichols et al., 2011).

Matrix of pairwise drug epistatic interactions: the dot product between the two association matrices is used to find pairs of epistatically interacting drugs.

Statistical parameters

Result definitions (e.g., point estimates and error bar explanations), as well as details on the number of replicates are reported in the figure legends.

Data and Software Availability

The data generated in this work are available as Supplementary Tables (Table S1). Metabolome data can be downloaded from https://www.ebi.ac.uk/biostudies with accession number S-BSST245. Results from data analysis can be found as Supplementary Tables (Tables S2 and S3).

Additional Resources

The results of our work are made available through an interactive public data portal: https://zampierigroup.shinyapps.io/EcoPrestMet.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Bacterial and Virus Strains
E. coli BW 25113	KEIO collection (Baba et al., 2006)	N/A
ΔrhlB, ΔaspC, ΔrpoS, ΔnudB E. coli	KEIO collection (Baba et al., 2006)	N/A
fadAB, recA and entCEBAH fluorescent transcriptional reporters	Zaslaver library (Zaslaver et al., 2004)	N/A
Chemicals, Peptides, and Recombinant Proteins
1,280 chemically diverse drugs	Prestwick chemical drug library	http://www.prestwickchemical.com
Ofloxacin	Sigma-Aldrich	O8757
Rifampicin	Sigma-Aldrich	R3501
Oleic acid	Sigma-Aldrich	O1008
Meropenem trihydrate	Sigma-Aldrich	M2574
Docetaxel	Sigma-Aldrich	01885
Ascorbic acid	Sigma-Aldrich	A5960
Disulfiram	Sigma-Aldrich	D2950000
Phosphoenolpyruvate	Sigma-Aldrich	P0564
Benzydamine	Sigma-Aldrich	BP610
Sulfamethizole	Sigma-Aldrich	S2050000
Critical Commercial Assays
Xanthine Oxidase (microbial)	Sigma-Aldrich	CAS # X2252
Isocitrate dehydrogenase (Bacillus subtilis)	Megazyme	CAS # E-ICDHBS
Deposited Data
Phenotypic screen of the KEIO library	http://ecoliwiki.net/tools/chemgen/ (Nichols et al., 2011)	N/A
Metabolic profiles of the KEIO library	https://www.ebi.ac.uk (Fuhrer et al., 2017)	S-BSST5
Metabolic profile of the Prestwick library	https://www.ebi.ac.uk (This paper)	S-BSST245
Software and Algorithms
Mutual Information	CLR package (Faith et al., 2007)	http://m3d.mssm.edu
Iterative Similarity Score	This paper	N/A
Prediction of interacting drug pairs	This paper	N/A
MATLAB	Mathworks	https://www.mathworks.com/products/matlab.html
Other