Tuberculosis (TB) remains notoriously difficult to treat and, until the COVID-19 pandemic, was the leading cause of death by a single infectious agent, Mycobacterium tuberculosis (Mtb) (WHO, 2021). TB requires lengthy combination therapy to prevent the acquisition of heritable drug resistance [addressed by others in this collection, including (Bhagwat et al., 2022; Jones et al., 2022; Liebenberg et al., 2022)] and to effectively target inherent heterogeneity in infection that results in drug tolerance to prevent treatment failure and subsequent disease relapse (Fox et al., 1999; Kerantzas and Jacobs, 2017). Genotypic resistance is defined by the heritable ability to grow in the presence of high concentrations of antibiotics beyond the minimum inhibitory concentration (MIC) (Balaban et al., 2019). Antibiotic tolerance is defined as the non-heritable ability of bacteria to survive transient exposure to drugs at concentrations that would otherwise be lethal (Brauner et al., 2016; Balaban et al., 2019). The current regimen for drug-sensitive Mtb requires a 2-month intensive phase of rifampicin, isoniazid, pyrazinamide, and ethambutol followed by a four or seven-month continuation phase of rifampicin and isoniazid. The phase III clinical trial (“Study 31”) for a four-month regimen with rifapentine, moxifloxacin, isoniazid, and pyrazinamide recently demonstrated non-inferiority to the six-month standard of care and is now recommended for some patients with drug-sensitive TB (Dorman et al., 2021; Carr et al., 2022). This progress notwithstanding, shorter and more effective therapies are urgently needed for both drug-sensitive and particularly for drug-resistant TB.

A hallmark of TB pathogenesis is the formation of different lesion types that have varied structures and provide different microenvironmental conditions to resident Mtb ( Figure 1 ). Mtb readily adapts to these different niches to withstand environmental stressors, often by slowing or halting replication and metabolic activity. These adaptations enable Mtb to tolerate drug treatment (Lenaerts et al., 2015; Sarathy and Dartois, 2020). Because the bacteria occupy separate niches with different environments, multiple states of drug tolerance exist together, influencing potential drug activity (pharmacodynamics) (Cadena et al., 2017). Furthermore, the physical structures of lesions vary, and drug penetration and accumulation in each lesion are dependent on the lesion structure and chemical properties of each drug (pharmacokinetics) (Prideaux et al., 2015; Sarathy et al., 2016) ( Figure 1 ). Therefore, drug combinations with varied pharmacokinetic and pharmacodynamic properties are necessary to target the bacteria in all their locations.

More than twenty-five new TB drugs are at various stages across the developmental pipeline (newtbdrugs.org). There are 2,300 possible three-way combinations for twenty-five drugs, far too many to test in vivo. Single-drug responses may not be indicative of drug combination activity and potential to reduce relapse, and drugs may not behave in combination as they do alone due to drug interactions such as synergies and antagonisms. For example, pyrazinamide, a crucial member of the current standard of care regimen, is inactive against Mtb in many in vitro growth conditions and as a monotherapy in vivo (Steenken and Wolinsky, 1954; Cho et al., 2007; Lanoix et al., 2016). However, pyrazinamide contributes significant sterilizing activity when used for the first two months of treatment in combination with the other drugs in the standard of care (Fox et al., 1999). On the other hand, pyrazinamide unexpectedly exhibited antagonism when combined with rifampicin in a 14-day PET/CT clinical study (Xie et al., 2021). This finding was consistent with the prior evidence that the combination of rifampicin + pyrazinamide + isoniazid for the continuation phase of treatment displayed slightly higher relapse rates than rifampicin + isoniazid without pyrazinamide (East African- British medical research councils, 1973). Because combination therapy is necessary to treat TB and single-drug responses are not informative of combination activity or interactions, we need to be able to systematically screen and optimize combination therapies. However, it is logistically impossible to screen the entire drug combination space in animal infection models. Therefore, it is critical to be able to perform systematic combination screens using in vitro and computational tools that translate to animal and clinical outcomes. There is an urgent need for in vitro and computational tools that capture the complex pharmacodynamics (due to drug tolerance) and pharmacokinetics (due to lesion structure) in infection to design drug combinations that successfully target the bacteria that are hardest to access and kill.

Animal models have been a critical part of the preclinical development of drug combinations. Multiple models are used in the drug development pipeline to balance throughput, cost, and information gained. Heterogeneity across animal models and even within species complicates comparison across models. Mouse models have been used for decades; the most commonly used models until recently were the C57BL/6 and BALB/c models. These mice are genetically resistant to TB infection and make only cellular granulomas, thus presenting a major limitation in their ability to recapitulate human disease. The C3HeB/FeJ model has recently grown in popularity and use; this model is genetically susceptible to TB and makes both cellular and caseous granulomas, more similar to human disease (Driver et al., 2012). Larger animals (i.e., the New Zealand White rabbit, marmosets, cynomolgus and rhesus macaques) recapitulate human disease more closely than smaller animals but are logistically prohibitive for comprehensive early screens (Lenaerts et al., 2015). Given the massive combination landscape and logistical limitation of screening every combination in animal models, many in vitro models have been developed to mimic the different stressors encountered in lesions. The utility of these models in predicting animal and clinical outcomes is not always clear. Standard minimum inhibitory concentration assays using glucose-based media optimized for Mtb growth were the gold standard for a long time. However, it is well-established that glucose is not the primary carbon source for Mtb in vivo, and Pethe et al. demonstrated the importance of making measurements of drug response in more in vivo-like conditions when they identified a class of drugs with glycerol-specific activity against Mtb that was inactive in vivo (Pethe et al., 2010). This study also highlighted the importance of validating the predictive power of in vitro measurements of in vivo outcomes.

In this review, we will discuss how heterogeneous lesion microenvironments give rise to non-replication and subsequent drug tolerance, why this tolerance requires combination therapy for effective TB treatment, and the implications for screening assay design. We review the current state of in vitro and computational tools for drug screens and speculate on areas of improvement as they relate to drug tolerance and non-replication. New tools have been developed to allow for systematic study of higher-order drug combinations, including the DiaMOND methodology to reduce the number of measurements required in the higher-order combination screening space, a mathematical “dose” model to predict higher-order combination activity from drug pairs, and chemo-genomic and regulatory network models to predict drug interactions from single-drug transcriptomics (Peterson et al., 2016; Zimmer et al., 2016; Cokol et al., 2017; Ma et al., 2019). Beyond drug interactions, additional tools have been developed to consider single-drug and drug combination potency (Wooten et al., 2021; Li et al., 2022). Computational models have advanced to incorporate complex pharmacokinetics and pharmacodynamics to predict drug response outcomes informed by animal data and clinical data (Pienaar et al., 2017; Savic et al., 2017; Lyons, 2022). With all these developments, we have the opportunity to consider how to best utilize all these tools in the context of drug tolerance for the rational development of drug combinations that will be most effective in treating TB.

To identify drugs and combinations that will best target the bacteria that are most tolerant to treatment, it is critical to understand the underlying causes of the drug tolerance that is most difficult to treat, how we currently model environments that induce drug tolerance, and to fill the knowledge and technical gaps that will improve our current models. Heterogeneous lesions are the hallmark of TB infection, and these lesions present a variety of microenvironments to which Mtb adapts. TB granulomas are classically categorized as fibrotic, cellular/non-necrotic, caseous/necrotizing, or cavitating, each with distinct characteristics (Lenaerts et al., 2015). However, the reality of infection is more of a spectrum than discrete classes, and the microenvironmental components are not homogenous even within one lesion type (Barry et al., 2009). Mtb adaptation to these environments results in varied drug tolerance that requires combination therapy to treat effectively. A critical effect of many of the environmental components that the bacteria experience is slowing or halting replication. Non-replication is associated with extreme drug tolerance (Sarathy et al., 2018). In this section, we detail how different environments slow and halt Mtb growth, their influence on drug tolerance, and how they are modeled.

Combination therapy’s superiority over monotherapy in bactericidal activity and preventing relapse is well-established; however, drug screens and animal and clinical studies typically focus on outcomes and rarely investigate and understand the relationship between this superiority and targeting drug tolerance. Early studies in the development of the current standard of care (isoniazid + rifampicin + pyrazinamide + ethambutol, or HRZE) found that the addition of rifampicin or pyrazinamide to the treatment regimen with streptomycin and isoniazid reduced disease relapse and including both allowed for shorter treatment duration (Fox et al., 1999). Though it was initially thought that this success was due to drug synergy (East African- British medical research councils, 1974), it was ultimately attributed to activity against semi-dormant bacteria (Dickinson and Mitchison, 1981). This was demonstrated using in vitro studies of isoniazid- or rifampicin-treated Mtb whose growth rate was stalled by decreasing incubation temperature or culture in acidic conditions. Short exposure to optimal growth conditions resulted in short recovery bursts, during which rifampicin had greater bactericidal activity than isoniazid, suggesting that rifampicin is better at killing these bacilli that are semi-dormant, that is, dormant much of the time with occasional metabolic bursts. Pyrazinamide has also been demonstrated specifically to target dormant (non-replicating) bacilli and loses its sterilizing activity when metabolic activity resumes (Hu et al., 2006). Therefore, using multiple drugs in combination may reduce the incidence of relapse by targeting multiple bacterial states in infection. In support of this, Mitchison and colleagues used a guinea pig model and showed that the addition of rifampicin to isoniazid or isoniazid + ethambutol did not increase the bactericidal activity of the combination but reduced relapse, suggesting that the success of the combination was not in additional killing, but in killing particular cells that would have otherwise caused relapse (Dickinson and Mitchison, 1976). Walter et al. proposed a ribosomal RNA synthesis ratio as a metric to distinguish sterilizing and non-sterilizing drugs and drug combinations and demonstrated its utility both in vitro and in the relapsing mouse model (Walter et al., 2021). This method could be used to provide molecular insight into how sterilizing combinations modulate cellular processes in subpopulations of drug-tolerant Mtb. Furthermore, the combination of the ribosomal RNA synthesis ratio with CFU measurement in the mouse model was recently found to be more informative of treatment-shortening potential than either metric alone (Dide-Agossou et al., 2022). A combination of pharmacodynamic markers could enhance our understanding of how combination therapy targets drug-tolerant cells both in vitro and in vivo.

Successful treatment of TB requires killing the bacteria that would otherwise remain to cause relapse. These cells may reside in different lesions and therefore exhibit different types of drug tolerance due to adaptation to different microenvironments. Combination therapy with multiple modes of action offers greater potential over monotherapy to kill multiple types of drug tolerance, and the best combinations will be those that successfully target different types of drug tolerance. Though this has not been demonstrated systematically, empirical evidence supports this concept: bedaquiline and pretomanid, which have performed remarkably well in both animal and clinical studies and have been approved in combination with linezolid for MDR-TB (NiX-TB) (Conradie et al., 2020), target both replicating and non-replicating bacteria, and have been shown to be particularly potent against non-replicating bacteria induced by different microenvironmental conditions (Cho et al., 2007; Koul et al., 2008; Tasneen et al., 2011; Gold et al., 2015). Their superior sterilizing performance could be attributed to their ability to target multiple types of drug tolerance. Recently, a four-month regimen with rifapentine and moxifloxacin was found to be non-inferior to the standard of care (“Study 31”) (Dorman et al., 2021). In preclinical studies, moxifloxacin exhibited potent early bactericidal activity in treatment (Nuermberger et al., 2004). Moxifloxacin’s potent activity against both replicating and non-replicating Mtb may contribute to this more rapid clearance. Taken together, the evidence suggests that combination therapy is required to reduce relapse, and the most effective regimens will include drugs that target multiple types of drug tolerance.

Heterogeneous lesions also result in varied drug access, absorption, and metabolism (i.e., pharmacokinetics) dependent on lesion structure, which drives the necessity of combination therapy to ensure that one or more of the drugs can access all the different locations and be effective in those locations ( Figure 1 ). MALDI mass spectrometry imaging (MSI) has been used to show how different drugs penetrate caseum in vivo (Sarathy et al., 2016). Pyrazinamide and isoniazid penetrate fairly evenly, but bedaquiline, which is highly lipophilic, binds intracellular lipids and caseum macromolecules and penetrates caseum very poorly (Sarathy et al., 2016; Greenwood et al., 2019). Clofazimine, which is highly effective against Mtb in BALB/c mice (which make exclusively cellular lesions), is relatively ineffective against Mtb in C3HeB/FeJ mice (which develop caseous necrotic lesions) (Irwin et al., 2014). Pharmacokinetic modeling demonstrated that clofazimine accumulates in cellular layers and does not diffuse into necrotic foci, explaining the lack of efficacy in the C3HeB/FeJ mice (Prideaux et al., 2015). Therefore, modeling lesion drug access to ensure that sufficient drug(s) reach the different types of drug tolerance in all their locations is important for the design of effective combinations.

Pharmacokinetics are also important when considering drug interactions; if an antagonistic combination of drugs do not act in the same location, the antagonism may not be realized. Lesion-specific and caseum-specific pharmacokinetics are studied using a combination of MSI in vivo, ex vivo rabbit caseum, and in vitro caseum surrogate derived from foamy THP-1 macrophages (Sarathy et al., 2016). Lesion pharmacokinetics are also modeled using computational models, typically informed by in vivo data (Kjellsson et al., 2012; Prideaux et al., 2015). Most pharmacokinetic/dynamic (PK/PD) models do not consider drug interactions. Recently INDIGO-MTB, a computational tool to predict drug interactions using transcriptomics, was integrated with GranSim, a multi-scale model of tuberculosis granuloma formation. GranSim incorporates host immunity, Mtb growth dynamics, and drug PK/PD into one computational framework that describes interactions between these entities through space and time in granuloma evolution (Cicchese et al., 2021). Expanding this and other models to consider drug combinations, drug interactions, varied drug tolerance, and lesion-specific pharmacokinetics could provide invaluable information for the development and improvement of effective drug combinations.

Though many lesion microenvironmental components are thought to affect drug susceptibility in some capacity, it may be that we do not need to model every single one to capture the drug tolerance that causes disease relapse, i.e., the bacteria that are hardest to treat. Instead, we should identify which of these environmental niches and resultant bacterial states represent the bacteria that are most difficult to treat. Using a guinea pig model, Lenaerts and colleagues showed that the bacteria that remained after treatment were extracellular, primarily in the hypoxic acellular rim of caseous necrosis (Lenaerts et al., 2007). A marmoset study showed that the difference between sterilizing and non-sterilizing regimens was rapid clearance of cavitating lesions by the sterilizing regimen, indicating that drug combinations must target the cavitating lesion for a better outcome (Via et al., 2015). Given the complexity of the tuberculosis heterogeneity, computational models are a critical tool to link in vitro drug responses to such treatment outcomes in animals and the clinic. Recently Larkins-Ford et al. used a machine learning approach to predict relapse outcomes in the BALB/c mouse model from a suite of drug combination potency and interaction measures made in a variety of growth conditions to model the different environments in TB lesions (Larkins-Ford et al., 2021). This study demonstrated that in vitro measurements of response to drug combinations made in a subset of growth conditions (as a “sum of parts”) could be predictive of in vivo outcomes. Some sets of conditions were more predictive of in vivo response than others, thereby identifying validated sets of in vitro models as suitable for drug combination screening. We have yet to understand whether or not measurement in more complex growth conditions that combine these “parts” into one condition (i.e., a multi-stress model) will be more predictive than the sum of parts approach. Improving in vitro models to capture the environments that result in drug tolerance that enables the bacteria to survive beyond treatment and correlating response in these models to in vivo outcomes is a clear next step. Computational and mathematical modeling offer the advantage that they are inherently optimal for iterative learning and can leverage in vivo outcome data.

To get the best prediction from computational models, the input from in vitro models must be optimized and streamlined. Several tools have been developed to improve the quality and utility of metrics from in vitro models of drug combination activity (summarized in Table 1 ). Diagonal measurement of n-way drug interactions (DiaMOND) is a method to measure drug interactions using only a fraction of the drug-dose combination matrix (“checkerboard”), which is logistically prohibitive to measure for higher-order of systematic studies (Cokol et al., 2017). DiaMOND measures dose response curves, enabling the collection of additional metrics beyond drug interaction, including combination potency metrics. MuSyC is another framework that calculates drug interactions and distinguishes between different types of combination effects (whether we consider dose or efficacy to evaluate drug synergy) to overcome conflicting assumptions of the widely used Loewe Additivity and Bliss Independence principles and to account for drug interactions and combination potencies (Wooten et al., 2021). This widening of combination effects from only measuring traditional synergy may be important in our ability to develop predictive models of in vivo outcomes from in vitro data. Larkins-Ford et al. recently demonstrated that combination potency metrics were important to accurately predict treatment outcomes in mouse models (Larkins-Ford et al., 2021; Larkins-Ford et al., 2022). Additional modeling techniques utilize pairwise drug response measurements to predict high-order drug interactions (Zimmer et al., 2016; Katzir et al., 2019) and in vivo treatment outcomes (Larkins-Ford et al., 2022) as a path to reduce the number of measurements required to study the increasingly large drug combination space and still obtain informative metrics beyond traditional drug interactions. Other new computational tools incorporate pathway-specific effects of drug action to predict drug interactions (Peterson et al., 2016). INDIGO (inferring drug interactions using chemo-genomics and orthology), for example, uses transcriptomic data from single-drug responses to predict drug interactions (Chandrasekaran et al., 2016; Ma et al., 2019). These integrated molecular approaches may help us understand the pathways underlying drug tolerance and combination drug response.

More complex, mechanistic computational models integrate host-pathogen interactions, pharmacodynamics, and pharmacokinetics (summarized in Table 1 ) (Ernest et al., 2021). GranSim is a multi-scale, agent-based model that captures the temporal and spatial dynamics of immune cell activity, bacterial growth, and bacterial killing by drugs in different lesion types, informed by pharmacokinetic and pharmacodynamic data from non-human primates and rabbits. These dynamics are captured at the molecular, scale, and whole lesion scales. As new in vivo data are acquired, the model may be updated. GranSim has demonstrated how changes in cytokines influence early infection and lesion formation (Fallahi-Sichani et al., 2011; Wong et al., 2020). GranSim can also be integrated with other computational frameworks, including INDIGO-MTB and a constraint-based model (CBM) that predicts Mtb metabolism and growth. GranSim-CBM is used to predict how environmental influence on Mtb metabolism and growth influences granuloma development and outcomes (Pienaar et al., 2016). The hollow fiber system model of tuberculosis is an in vitro drug development tool to optimize drug regimens and dose selection to maximize drug or combination efficacy and minimize the emergence of genotypic resistance. The in vitro system has a pharmacodynamic compartment that houses the bacteria and a pharmacokinetic compartment with semipermeable hollow fibers that allow drugs to diffuse to the pharmacodynamic compartment in a manner that mimics the appropriate concentration-time profile for the drugs. The results of the in vitro experiments are used in Monte Carlo simulations to predict outcomes in clinical populations (Gumbo et al., 2004; Gumbo et al., 2015). Another PK/PD model incorporates human clinical data to simulate the outcome of various multi-drug regimen scenarios (available at http://www.saviclab.org/systems-tb/) (Fors et al., 2020). Savic and colleagues have developed population PK/PD models that were used to determine the optimal dose of rifapentine with the most potential to shorten the duration of treatment based on outcomes from phase II clinical trials (Savic et al., 2017). This information was used to inform the dosing strategy for the phase III clinical trial for Study 31, which was found noninferior to the 6-month standard of care (Dorman et al., 2021). Incorporating validated in vitro model data that capture the drug tolerance in vivo in PK/PD models offers the potential to predict treatment outcomes at the granuloma level and to fine-tune and optimize regimens to target the bacteria that would otherwise withstand treatment either due to drug tolerance or poor drug access and cause relapse.

Though the field of computational models and mathematical frameworks for measuring combination drug response has seen tremendous progress and innovation, gaps remain. Tools like INDIGO-MTB and the dose model enable the prediction of higher-order interactions from single-drug transcriptomics or pairwise drug interactions and raise the question of the necessity of higher-order combination screens. Recently Larkins-Ford et al. demonstrated that synergy and potency metrics derived from pairwise combination measurements could be used to build successful higher-order combinations (Larkins-Ford et al., 2022). Streamlining predictive metrics of treatment outcomes (both animal and clinical) and the definition of “systematic” (is pairwise sufficient or is higher-order necessary)? is important to improve the in vitro data that will be used in modeling. There is also a need for a system to report failed regimens to improve computational models. Understanding why regimens failed will help inform model parameters to identify other combinations that might also fail.

The development of improved, streamlined in vitro models that inform computational models is not the end of the line. Additional considerations that may improve models include single-cell heterogeneity and strain-to-strain differences. For example, innate growth and metabolic differences among closely related bacilli create subpopulations that are differentially susceptible to drug treatment (Aldridge et al., 2012; Manina et al., 2015; Rego et al., 2017; Shee et al., 2022). Antibiotic-stressed, replicating mycobacteria display drug tolerance at the single-cell level; antibiotic-induced tolerance and tolerance in replicating Mtb warrant further study (Wakamoto et al., 2013; Zhu et al., 2018). It is unknown if and how much single-cell heterogeneity contributes to the subpopulations of bacteria in humans that survive beyond treatment to cause relapse. Once we have improved and validated in vitro models of non-replicating Mtb that survive beyond treatment, single-cell studies can be used to determine if certain cells exhibit identifiable transcriptomic or morphological markers that pre-dispose them to greater drug tolerance. We might also consider the effects of resistance mutations on drug combination activity. Schrader et al. showed that mutations that arise due to exposure to one antibiotic can cause multiform antibiotic resistance, i.e., different types of drug tolerance that may not be heritable (Schrader et al., 2021).

Effective targeting of drug-tolerant Mtb will also require that we understand how Mtb lineage contributes to different manifestations of tolerance. Different Mtb strains exhibit differences in virulence, adaptation to environmental stressors, and drug tolerance (De Groote et al., 2012; Tizzano et al., 2021). Using genome-wide association studies (GWAS), Hicks et al. identified SNPs in clinical strains of the Beijing lineage associated with propionate metabolism that conferred drug tolerance to certain drugs (Hicks et al., 2018). Recently Li et al. utilized CRISPRi technology to quantify the expression of Mtb genes and bacterial fitness in the presence of different drugs and discovered genetic mechanisms of intrinsic drug tolerance. Upon comparing this data against genomics of Mtb clinical isolates, strain-dependent mutations were discovered that conferred specific drug susceptibility (Li et al., 2022). These studies offer the possibility to tailor drug combinations to Mtb lineage based on genetic differences across lineages. Another recent study showed that the resR transcriptional regulator is a frequent target of positive selection for mutations and that strains with these mutations exhibit antibiotic “resilience,” which describes the potential for antibiotic tolerance and genetic resistance. The bacilli harboring these mutations exhibited more rapid post-antibiotic recovery than wild-type cells (Liu et al., 2022). This new finding demonstrates that we must better understand the various strategies Mtb uses to recover from drug treatment beyond standard definitions of non-replication and dormancy so that we can design combination therapies that target these subpopulations. To ensure that drug combinations adequately target tolerance across lineages and resilient subpopulations, multiple (clinical and laboratory) strains could be used in in vitro models for combination drug screens. We propose an iterative learning process ( Figure 2 ) to develop improved drug combinations where in vitro models of non-replicating-induced drug tolerance are used to screen drug combinations. Streamlined drug interaction and potency metrics would be used as input for computational models of complex PK/PD to identify better drug combinations to target drug-tolerant Mtb. Computational models should then identify potential populations that could cause relapse and then loop back to the in vitro models to redesign the combinations to better target those populations. Multiple strains could then be tested and used as additional input in the computational model so that the best regimens that target drug tolerance across lesions and strains are carried forward to test in animals and the clinic.

TB requires lengthy multidrug treatment due to populations of drug-tolerant bacteria that arise from the heterogeneity of infection. Complex lesion structure and microenvironments induce multiple types of non-replicating bacterial populations with differential drug susceptibility and variable drug access. Combination therapy is required to target these populations in all their locations, and better tools for in vitro screens are urgently needed. The ability to screen combinations systematically will enable us to consider the best combinations that will target drug-tolerant cells in a data-driven manner. Incorporating these data in lesion-scale PK/PD models that simulate which bacteria evade kill will also help inform which niches in the lesions induce the drug tolerance that results in relapse. Understanding which bacterial subpopulations remain will inform better choices of compounds to target those types of tolerance. Those drugs can then be included in systematic combination screens to design optimized combinations that target the types of drug-tolerant bacteria that evade kill. To achieve this goal, we must fill the knowledge gaps delineated in this review about the drug-tolerant non-replicating populations, streamline measurements from combinations screens, and optimize computational models of PK/PD to account for drug combinations and multiple types of drug tolerance.

Writing and editing: TG and BA. All authors contributed to the article and approved the submitted version.