(591e) Using Genome-Scale Metabolic Models and Machine Learning to Design Condition-Specific Combination Therapies

Background: Antibiotic resistance is a rapidly growing problem worldwide, especially as the supply of novel antibiotics dwindles. A promising solution is combination therapy, but the current trial-and-error approach for finding optimal drug combinations is time consuming and expensive. In addition, the growth environment of the pathogen greatly influences the potency of antibiotics^1,2, and the duration of antibiotic treatment impacts the rise of both drug resistance and collateral sensitivity to other antibiotics³. To address these challenges and expedite drug discovery, here we present a novel approach that can evaluate drug combination therapies using multiple factors that play a critical role in impacting drug resistance. Our approach uses flux data inferred from genome-scale metabolic models (GEMs) in combination with machine learning (ML) to predict condition-specific combination therapy outcomes. This approach overcomes the limitations of existing methods that do not account for diverse pathogen growth conditions, time intervals between treatments, and lack the flexibility to simultaneously use diverse data types^4–6.

Methods: We developed a data-driven approach called Condition-specific Antibiotic Regimen Design in silico using GEMs (CARDiG), which involves a two-step modeling process. First, high-throughput omics data and metabolite composition of the extracellular environment serve as GEM inputs to determine flux profiles in response to drug treatment and growth in defined media, respectively (Figure 1A). Second, GEM-derived flux profiles and drug interaction data serve as inputs to train a ML model that predicts interaction outcomes for novel drug combinations (Figure 1B).

We determined metabolic flux profiles in response to drug treatment and condition-specific growth by constraining the E. coli GEM iJO1366⁷ and the Mycobacterium tuberculosis (M. tb) GEM iEK1011⁸. Specifically, we used chemogenomic data for E. coli⁹ and transcriptomic data for M. tb⁶ to infer differential single-gene knockout fitness and differential gene expression, respectively. This yielded individual sets of differentially regulated genes that were imposed as GEM constraints to determine flux profiles for single drug treatments. To determine growth media flux profiles, we referred to the availability of metabolites present within a media condition to constrain the GEMs.

Prior to ML model training, we processed drug and media flux profiles to determine joint profiles for all combinations of interest. These profiles were comprised of three pieces of information: (a) the combined effect of all treatments (i.e. sigma scores), (b) the unique effect of individual treatments (i.e. delta scores), and (c) the overall metabolic entropy (i.e. entropy scores) as defined by Zhu et al.¹⁰ for drug and media flux profiles. We trained the ML algorithm - Random Forests, to associate these three feature patterns to drug combination outcomes. Using the trained ML model, we predicted outcomes for new drug combinations based on their feature information alone. We then compared our predictions against experimental data by calculating the Spearman correlation and the area under the receiver operating curve (AUROC).

Results: We benchmarked CARDiG against previous approaches by directly comparing our prediction accuracy against those reported in literature and evaluated their performance for five different cases:

1. pairwise drug interaction outcomes for E. coli⁴
2. three-way drug interaction outcomes for E. coli⁵
3. pairwise drug interaction outcomes for E. coli cultured in M9 glycerol media⁵
4. pairwise and three-way drug interaction outcomes for M. tb⁶
5. interaction outcomes for pairwise to five-way TB regimens used in clinical trials¹¹

Of note, the first, second, and fourth cases tested the model’s ability to predict novel combinations involving drugs with unseen mechanisms of action. The third case assessed whether the model could predict drug interaction outcomes in a new growth environment, while the fifth case ascertained if predicted outcomes corresponded with clinical efficacy. Table 1 provides a summarized comparison of CARDiG results against those reported in literature. Overall, we found that our approach retained high accuracies in predicting combination therapy outcomes for a diverse set of test cases based on both E. coli and M. tb data. This is notable given that our approach uses data on metabolic genes alone, suggesting that metabolism plays an important role in determining drug interactions.

To demonstrate the power of using GEMs in predicting condition-specific combination therapy outcomes, we next applied CARDiG to predict pairwise drug interactions in multiple media conditions. For this task, we used experimental data of E. coli treated with four single drug treatments (AZTreonam, CEFoxitin, TETtracycline, TOBramycin) and two pairwise drug treatments (CEF + TET, CEF + TOB). Each treatment outcome was assessed for E. coli cultured in Biolog phenotype microarrays (PMs)¹² 1 and 5, which measured metabolic respiration in 96 carbon sources and 96 nutrient supplemented conditions, respectively (Figure 2A). Out of these 192 media conditions, 109 could be simulated based on the metabolites annotated in the E. coli GEM.

We constructed a ML model using the flux profiles for the four drug treatments as well as the 109 media conditions, and interaction outcomes for the four single drug treatments in 109 media combinations. We then evaluated our model by predicting outcomes for the pairwise treatments across all media combinations (Figure 2B). Overall, we found that model predictions significantly correlated with experimental outcomes (R = 0.79, p ~ 9E-48, Figure 2C). Given the nature of the test dataset, where drug-drug-media interactions fall into four broad groups defined by two distinct drug pairs and two Biolog array types, we further assessed group-based correlations. We found that model predictions still corresponded well with experimental data for all groups but one (CEF + TOB + PM05). Nevertheless, our model correctly predicted CEF + TOB treatment in PM05 conditions to be the most synergistic amongst the validation set.

In addition to predicting combination therapies where drug treatments are given simultaneously, we extended our approach to predict treatment efficacy for time- and order-dependent (i.e. sequential) interactions. For this task, we used data for E. coli evolved in single drug treatments for three timespans (10, 21, and 90 days) then subsequently treated with a second drug^13–15. To account for both time and order, we re-defined the delta scores for sequential joint profiles. Additionally, we added a time feature within joint profiles indicating the total treatment time for a combination. Our model achieved significant correlations in predicting drug interactions in unseen time points (R = 0.29, p ~ 1E-13) and using 10-fold cross-validation (R = 0.47, p ~ 5E-36).

Conclusions: In this study, we have predicted combination therapy outcomes based on simulated metabolic information from GEMs. We also showed the advantage of using a GEM-based approach for investigating drug interaction outcome in a myriad of growth conditions. Finally, we demonstrated how our proposed approach can be extended to predict outcomes for sequential interactions. Overall, CARDiG serves as a versatile tool that can be used to expedite the design of effective combination therapies using diverse data types. Considering that GEMs for hundreds of organisms are available, our proposed approach holds potential to predict efficacious combination therapies against a wide range of pathogens. With the additional flexibility in simulating various media conditions, our approach can also be leveraged to predict interaction outcomes in specific growth conditions such as the disease microenvironment. This could help expedite clinical translation by identifying combination therapy candidates that are effective in vivo. Finally, the ability to predict outcomes for sequential interactions can be a powerful tool to design cyclic drug regimens that mitigate resistance.

References:

1. Brown, S. A., Palmer, K. L. & Whiteley, M. Nature Reviews Microbiology vol. 6 657–666 (2008).
2. Zampieri, M. et al. Mol. Syst. Biol. 13, 917 (2017).
3. Nichol, D. et al. Nat. Commun. 10, 1–10 (2019).
4. Chandrasekaran, S. et al. Mol. Syst. Biol. 12, 872 (2016).
5. Cokol, M., Li, C. & Chandrasekaran, S. PLOS Comput. Biol. 14, e1006677 (2018).
6. Ma, S. et al. MBio 10, (2019).
7. Orth, J. D. et al. Mol. Syst. Biol. 7, 535 (2011).
8. Kavvas, E. S. et al. BMC Syst. Biol. 12, 25 (2018).
9. Nichols, R. J. et al. Cell 144, 143–56 (2011).
10. Zhu, Z. et al. Nat. Commun. 11, 1–15 (2020).
11. Bonnett, L. J., Ken-Dror, G., Koh, G. C. K. W. & Davies, G. R. Clin. Infect. Dis. 46, 46–54 (2017).
12. Bochner, B. R., Gadzinski, P. & Panomitros, E. Genome Res. 11, 1246–1255 (2001).
13. Imamovic, L. & Sommer, M. O. A. Sci. Transl. Med. 5, 204ra132-204ra132 (2013).
14. Oz, T. et al. Mol. Biol. Evol. 31, 2387–2401 (2014).
15. Horinouchi, T. et al. Sci. Rep. 7, 1–11 (2017).