Estimating tuberculosis drug resistance amplification rates in high-burden settings

Background Antimicrobial resistance develops following the accrual of mutations in the bacterial genome, and may variably impact organism fitness and hence, transmission risk. Classical representation of tuberculosis (TB) dynamics using a single or two strain (DS/MDR-TB) model typically does not capture elements of this important aspect of TB epidemiology. To understand and estimate the likelihood of resistance spreading in high drug-resistant TB incidence settings, we used epidemiological data to develop a mathematical model of Mycobacterium tuberculosis (Mtb) transmission. Methods A four-strain (drug-susceptible (DS), isoniazid mono-resistant (INH-R), rifampicin mono-resistant (RIF-R) and multidrug-resistant (MDR)) compartmental deterministic Mtb transmission model was developed to explore the progression from DS- to MDR-TB in The Philippines and Viet Nam. The models were calibrated using data from national tuberculosis prevalence (NTP) surveys and drug resistance surveys (DRS). An adaptive Metropolis algorithm was used to estimate the risks of drug resistance amplification among unsuccessfully treated individuals. Results The estimated proportion of INH-R amplification among failing treatments was 0.84 (95% CI 0.79–0.89) for The Philippines and 0.77 (95% CI 0.71–0.84) for Viet Nam. The proportion of RIF-R amplification among failing treatments was 0.05 (95% CI 0.04–0.07) for The Philippines and 0.011 (95% CI 0.010–0.012) for Viet Nam. Conclusion The risk of resistance amplification due to treatment failure for INH was dramatically higher than RIF. We observed RIF-R strains were more likely to be transmitted than acquired through amplification, while both mechanisms of acquisition were important contributors in the case of INH-R. These findings highlight the complexity of drug resistance dynamics in high-incidence settings, and emphasize the importance of prioritizing testing algorithms which allow for early detection of INH-R. Supplementary Information The online version contains supplementary material available at 10.1186/s12879-022-07067-1.


Background
Despite being both a preventable and curable disease, more than 10 million people develop tuberculosis (TB) each year, with 1.4 million deaths in 2019 [1]. Although 63 million lives have been saved through improvements in programmatic TB management this century, the increase in drug-resistant (DR-TB) cases is increasingly concerning [1]. Multidrug-resistant TB [MDR-TB; defined as resistance to both first-line drugs isoniazid (INH) and rifampicin (RIF)] is a particular barrier to TB control efforts [2]. In 2019, 465,000 people were diagnosed with MDR-TB [1]. MDR-TB can be acquired by transmission (primary resistance) or develop in vivo through inadequate or incomplete treatment (secondary resistance), and the relative contribution of these mechanisms is likely to vary by context [3]. In all settings, Open Access *Correspondence: justin.denholm@mh.org.au 3 Victorian Tuberculosis Program and Department of Microbiology and Immunology, Doherty Institute of Infection and Immunity, University of Melbourne, 792 Elizabeth Street, Melbourne, VIC 3000, Australia Full list of author information is available at the end of the article though, careful optimization of both clinical and public health management of MDR-TB is required to ensure good outcomes.
Mathematical modeling is increasingly used to support programmatic optimization for TB [4][5][6]. Accounting appropriately for DR-TB in mathematical models of disease is critical, as it differs considerably from drugsensitive TB (DS-TB) in both epidemiological parameters and relevant outcomes. Some variation in disease characteristics is relatively well-understood, including the prolonged treatment duration [7], adverse event rates [8] and diagnostic pathway performance [9,10]. However, considerable uncertainties persist regarding important characteristics of MDR-TB, including pathogen's fitness, transmissibility and risk of resistance amplification related to treatment [11,12]. Attempts to better characterize these features of MDR-TB have been challenging, in part due to the diversity of gene mutations which may confer resistance, many of which have limited clinical and epidemiological outcome data to inform model parameterization. Computational biological approaches have recently been used to bridge this gap, providing tools to estimate the fitness and resistance impact of novel TB mutations [13][14][15].
Modeling also offers an opportunity to quantify amplification and transmission of drug-resistant TB, by fitting dynamic models to observed data. We therefore aimed to incorporate epidemiological data into an empirically calibrated model, in order to explore parameter estimation for drug resistance amplification and transmission associated with both INH and RIF.

Constructing the mathematical model and defining epidemiological parameters
We designed a deterministic compartmental model of Mtb transmission to capture five mutually exclusive health states with regards to TB infection and diseasesusceptible (S), early latent (L A ), late latent (L B ), infectious (I) and recovered (R, 16). The model included four TB strains: drug-susceptible (DS-TB, compartment subscript S), isoniazid mono-resistant (INH-R, compartment subscript H), rifampicin mono-resistant (RIF-R, compartment subscript R) and MDR-TB (compartment subscript M). It is to be noted that the strains are not phylogenetically related.
We assumed homogenous mixing in a closed population: All deaths are replaced as new births (rate π) entering the susceptible compartment. This includes both deaths due (1) N = S + L AS + L AH + L AR + L AM + L BS + L BH + L BR + L BM + I S + I H + I R + I M + R to TB disease (μ i ), as well as a universal population-wide death rate (μ).
When individuals in a population are infected with Mtb, they transition from the susceptible compartment (S) to the early latent compartment (L A ). The force of infection (λ) associated with each strain (Eq. 2) is defined as: where "x" indexes the drug resistance pattern-S, H, R or M.
β is the "effective contact rate" for DS-TB, defined as the product of the average number of contacts between two individuals per unit time and the probability of DS-TB transmission per contact. The relative transmissibility of the different strains is denoted r X and uses the DS-TB strain's transmissibility as reference ( r S = 1 ). In other words, r X represents the TB strains' relative fitness.
People entering the early latent compartment (L A ) can either progress (within two-years) directly to the active disease compartment (I) at rate ε, or transition to the late latent compartment (L B ) at rate κ. Progression from L B to the active disease state occurs at a much slower rate (ν), and is referred to as reactivation. Once individuals have entered the infectious compartment, one of the following six processes can occur: (1) the person may be correctly identified as having active DS-TB and commenced on treatment (rate τ), thence progressing towards cure and transitioning to the recovered (R) compartment; (2) person may be correctly identified to have DS-TB or DR-TB and commenced on treatment but experiences treatment failure without experiencing resistance amplification to other drugs and stay in the same infectious compartment; (3) spontaneous recovery (rate γ) with transition to the recovered compartment (R); (4) TB-related death (μ i ) (5) dying of natural causes or (6) the infecting strain could acquire resistance (α H and/or α R ) to isoniazid (INH-R), rifampicin (RIF-R) or MDR-TB and move to I H , I R and ultimately to I M compartments. To capture the progressive accrual of resistance with each transition, only one level of additional resistance not already present can be obtained during a disease episode. People who have spontaneously recovered from past TB or successfully completed treatment are both represented as a single compartment (R) on the assumption that prognosis is equivalent regardless of the infecting strain from which each person recovers. Once treatment is complete, the recovered person can transition back to L A through reinfection, represented as δ. We define δ (Eq. 3) as: where, RR r is the "relative risk of re-infection once recovered".
Latently infected people also have a risk of re-infection with the same or other strains represented as θ in the model; and the re-infecting strain would "override" the existing strain. We define θ x (Eq. 4) similarly to δ x as: where, RR i is the "relative risk of re-infection once latently infected". Figure 1 and Additional file 1: figure S1 which is a representation of the model, along with Eqs. 5-18, captures all the intermediary steps and parameters as explained in the paragraph above.
It is to be noted that the figure does not show individuals who are latently infected with a given strain will have the same strain if they develop active disease. An elaborated diagram in presented in the supplementary sheets where all the compartments modelled have been shown (Additional file 1: Fig. S1).
Ordinary differential equations used to define the fourstrain model Structure of four strain Mtb transmission model. The symbols S, L A , L B , I and R represent uninfected/susceptible, early latent, late latent, infected and recovered health states, respectively. The subscript "X" used in L A and L B compartments and other parameters, indexes the drug resistance patterns, with S, H, R and M representing susceptible, isoniazid mono-resistance, rifampicin mono-resistance and multidrug resistance respectively. The infectious compartment is elaborated in the figure to show the amplification flows of INH and RIF respectively, parameterized with α H and α R (red arrows). The green arrows represent infection/transmission flows, black arrows represent constant progression flows. Compartments stratified according to resistance profiles are shown in blue. (π-birth rate, λ X -force of infection, ε-rate of early progression, κ-rate of late progression, νreactivation rate, γ-spontaneous recovery rate, θ X -risk of re-infection once latently infected, μ-mortality rate, μ i -TB-specific mortality rate, τ X -treatment rate, δ X -risk of re-infection after recovery)

Parameter estimation
Parameters can be categorized as universal, countryspecific and time-variant parameters, as presented in Table 1.

Universal parameters
From the literature we gathered information on disease-specific and epidemiological parameters to calibrate the Mtb transmission model. We considered these parameters to be universal to all TB settings and so assigned the same values for all strains and settings (Table 1A).
Once we defined the parameters in our model, we next reviewed literature for information on the prior distributions of uncertain parameters (Table 1B).

Defining time-variant model processes
To capture the rise of drug resistance over time, we allowed the case detection rate (CDR, a proportion, defined in Eq. 20) and the treatment success rate (TSR) to vary with time. TSR is the probability of a person being first tested and ultimately put on treatment to be cured, or simply put the probability of treatment success at presentation. This parameter was further varied by strain.
People diagnosed with active TB are commenced on treatment upon identification and move from the infectious compartments (I, I H , I R and I M ) to the recovered compartment (R). The transition from the infectious to the recovered compartment is represented using the where, "d" is calculated by solving the following CDR equation: which can be re-written as: (* "t" represents time in Eqs. 19, 20 and 21). Sigmoidal functions were used to model progressive increases for both the CDR and the TSR between 1950 and 2020. The final value of the TSR was set to the most recent TSR estimate reported by the WHO. In contrast, the final value of the CDR was varied during calibration. This allowed flexibility in simulating the historical dynamics of TB control in the countries considered.

Defining the amplification rate
Treatment for tuberculosis begins once individuals are detected with TB and the TB strain is correctly identified. Treatment then proceeds and may result in three possible outcomes: death, successful treatment or treatment failure. Treatment failure can further be associated with acquiring new resistance to one additional (19) drug that was not previously present in the infecting organism. INH and RIF are part of the standard regimen for the treatment of drug-susceptible strains. Gain in resistance to either INH or RIF is represented using amplification rates α H or α R respectively in the model. Mathematical representation of INH and RIF monoresistant amplification is shown in Eqs. 22 and 23. where, ρ H = Proportion of previously INH-susceptible individuals that acquire resistance on treatment failure, ρ R = Proportion of previously RIF-susceptible individuals that acquire resistance on treatment failure, and "t" stands for time.

Model calibration to prevalence and notification data Prevalence data
The model presented above was calibrated to countryspecific data. We fitted the models using TB prevalence estimates from national TB prevalence (NTP) surveys  Table 2.

Notification data
We used WHO-reported TB notifications as a calibration target for both models. For Viet Nam, in 2018, 102,171 cases were notified and for The Philippines 382,543 cases were notified and we calibrated to the per capita notification rates corresponding to these values.

Uncertainty analysis
An Adaptive Metropolis (AM) algorithm was used to estimate model parameters [30], including drug resistance amplification rates. An AM algorithm adapts continuously to the target distribution. It uses the history of the process to tune the effective proposal distribution suitably. The size and the spatial distribution of the proposal distribution is significantly affected due to adaptation. Moreover, AM is easy to implement and use, with no additional computational cost. As soon as the simulation process starts, AM algorithms start using the cumulating information. This gives the algorithm a major advantage because at an early stage of the simulation, the rapid start of adaptation enables the search to be more effective. This even diminishes the number of function evaluations needed [30].
The AM algorithm [30] was used to generate samples from the posterior distribution of the parameters from 25,000 iterations for each country. The primary estimates are reported as the posterior median value for all parameters of interest such as amplification proportions, CDR, relative fitness of each modelled strain and the relative risk of infection once recovered (δ). The intervals reported are obtained by calculating the 25th and 75th percentile of each parameter's posterior distribution. Programming was done in Python 3.7.3 and all code and associated data are publicly available on GitHub (github. com/malanchak/AuTuMN). Figure 2A-C show the model fits to reported INH-R, RIF-R and MDR levels for the high DR-TB settings, The Philippines and Viet Nam respectively. Table 3 shows the posterior distributions of all calibrated parameters.

Drug resistance amplification and transmission
We observed higher proportions of drug resistance amplification for INH compared to RIF for both the high DR-TB incidence settings we simulated (Fig. 3). The estimated risk of INH-R amplification when treatment fails was 0.84 (95% CI 0.79-0.89) for The Philippines and 0.77 (95% CI 0.71-0.84) for Viet Nam. The estimated risk of RIF-R amplification when treatment fails was 0.05 (95% CI 0.04-0.07) for The Philippines and 0.011 (95% CI 0.010-0.012) for Viet Nam. This meant approximately 84% and 77% of the people who failed treatment in The Philippines and Viet Nam respectively would end up with resistance to INH.
In The Philippines, the model estimates for amplification from DS to INH-R was 26 per 100,000 (95% CI 15-33) people, followed by 10 per 100,000 (95% CI 6-14) people then gain resistance to RIF and moving to the MDR compartment. Comparing this to acquiring RIF resistance first, we see 2 per 100,000 (95% CI 1-3) moving from DS to RIF-R, followed by only 0.05 (95% CI 0.02-0.08) people gaining resistance to INH to move to the MDR compartment. A similar observation was seen for Viet Nam, the model estimates for amplification from DS to INH-R shows 6 per 100,000 (95% CI 4-8) people followed by people 4 per 100,000 (95% CI 3-5) gaining resistance to RIF and moving to the MDR compartment. In case of DS to RIF-R transition the estimates were 0.08 per 100,000 (0.06-0.11) people, followed by 0.0007 per 100,000 (0.0004-0.001) people gaining resistance to INH to move to the MDR compartment.
Our study also provided information on estimates of CDR with high precision for both the settings, as inclusion of notification and prevalence of infection data for the analysis helped in constraining the parameter. The estimates obtained for The Philippines was 0.49 (95% CI 0.47-0.51) and for Viet Nam was 0.66 (95% CI 0.63-0.69).

Discussion
From this modeling study we were able to construct a model which successfully replicated epidemiological dynamics in two high burden TB settings, incorporating parameters drawn from microbiological fitness data. Using this model to explore the development of drug resistance in these contexts, we found that a much higher proportion of treatment failure resulted in amplification for INH-R rather than for RIF-R. This finding is consistent with observed higher rates of INH-R globally [31,32] and allows consideration of factors which might be mechanistically important for understanding and planning a programmatic response. For example, in a prevalence survey from 1975, the rate of INH-R in Canadians with TB following an initial course of therapy was 75.6% [33].
One factor likely to play a significant role in preferential INH-R amplification is current methods of DR-TB detection which prioritize RIF's resistance identification. According to WHO and many country guidelines, TB patients with strains found to be resistant to RIF need to start on a recommended MDR-TB treatment regimen. Longer MDR-TB regimens, and historical secondline therapy regimens, frequently have INH included in them, irrespective of resistance to INH being either undetermined or confirmed. Re-treatment regimens in particular have often incorporated prolonged durations of INH therapy-for example the category II regimen used in the Philippines comprised of 8 months of INH, RIF and ethambutol supplemented by streptomycin for the initial 2 months, and pyrazinamide for the initial 3 months (2SHRZE/HRZE/ 5HRE) [34], and older treatment regimens used in Viet Nam comprised of 8 months of INH and ethambutol supplemented by initial two months of streptomycin, pyrazinamide and rifampicin (2SHRZ/6HE, 28). These factors may be further amplified by use of isoniazid in the private sector and/or through community pharmacy settings, where worse guideline adherence and increased risk resistance development has been shown [35,36] but with poorer treatment outcomes compared to NTPs [37,38].
In addition to programmatic insights, our model provides novel information on parametrizing CDR. This is important, as this parameter cannot be measured directly yet plays a significant role in informing robust mathematical model of TB transmission. As with any mathematical representation our model has certain limitations. Our model was calibrated to TB prevalence and DR surveys to estimate the risk of resistance amplification. But the definition of a TB case may change between surveys, even within the same country. We have adopted a simplified model structure that does not capture factors such as age, comorbidities and other heterogeneity associated with TB epidemics. These factors may affect the risk of resistance amplification. Our model is primarily built for pulmonary tuberculosis and does not include extra-pulmonary TB data, as our primary focus was on transmission. For the same reason, this model has been parametrized from adult TB data given the limited TB transmission from young children to others. In our model we assumed the risk of INH-R amplification is the same starting from I S , as compared to starting from I R ; the same applies for RIF-R amplification. We even assumed the fitness cost of MDR-TB is independent of that of INH-R of RIF-R. Therefore, these limitations can potentially influence the estimated risk of resistance amplification.
Furthermore, it is difficult to compare the results generated by our model with existing studies because the findings are novel and not many researchers have tried  to estimate risk of monoresistance amplification of INH using a four strain TB transmission model. Two studies which reported similar findings were-(1) study conducted in the KwaZulu-Natal Province of South Africa for extremely drug resistant TB transmission, reported 84% of the participants failed treatment [39]. (2) A study from Kampala, Uganda found proportion of patients that acquired resistance were low, to which the author mentions that the study only reflects a single cycle of treatment in a heavily treated cohort of patients and therefore their results may underestimate the degree of drugresistance amplification [40]. With respect to estimates for transmission, the results are comparable to the study conducted by Kendall et al. where they report estimated percentage of MDR-TB resulting from transmission varied substantially with different countries notification data; for example, 48% (30-75%) in Bangladesh versus 99% (91-100%) in Uzbekistan [41].
Historically, diagnosis of MDR-TB has been reliant on culture-based phenotypic testing, which in highburden settings may be applied selectively, such as after treatment failure. As part of the global policy to control DR-TB, many high burden settings have pledged to deploy the molecular diagnostic assay Xpert MTB/RIF (detects resistance only in RIF), which is a nucleic acid amplification test that can be directly applied to sputum samples [42,43]. As the presence of RIF resistance is highly predictive of MDR-TB, these policies have led to significant improvements in the appropriate initiation of second-line therapy [44]. However, as our work highlights, these algorithms may also be associated with selecting for and further amplifying INH resistance. Alternative molecular tools, such as the line probe assay MTBDRplus [45], Abbott RealTime MTB RIF/INH (Abbott), FluoroType MTBDR [46], BD MAX MDR-TB assay [47], Roche cobas MTB [48] and whole genome sequencing can identify both RIF and INH resistance, and may offer the programmatic advantages of rapid MDR-TB diagnosis while avoiding this secondary effect [49]. Further research into the association between specific INH resistance mutations and differential risk of transmission will be helpful in better defining the public health impact of this effect [50].
An interesting observation was made between the estimates of the fitness cost of each resistant strain and their respective transmission dynamics for both the countries. We see that the RIF-R and MDR-TB have similar fitness in both countries indicating homogeneity of the population, while INH-R strains seem to have a considerable difference between the two countries. This is reflected in the transmission dynamics observed for the INH-R strains, where Viet Nam has a higher transmission rate compared to The Philippines. This could mean the INH-R strains in Viet Nam, although under significant drug pressure in the host, are more transmissible, perhaps due to strain-specific factors or local host/pathogen interaction. The same correlation cannot be made for RIF-R strains. They have similar fitness cost in both countries, but the transmission rate is higher in Viet Nam compared to The Philippines.
Comparing fitness cost of INH-R strains to the proportions of individuals who develop INH-R due to treatment failure (ρ H ) for both the high burden DR-TB setting, we see that the fitness cost is inversely related to ρ H . While INH-R strains in Viet Nam have a higher fitness cost than The Philippines, the number of individuals acquiring resistance due to treatment failure is lower in Viet Nam. The estimates for CDR were also higher in Viet Nam compared to The Philippines. As, ρ H and CDR is directly proportional to the amplification rate as defined in Eq. 22, we can say that the amplification rate α H is inversely proportional to the fitness cost of INH-R strains. Therefore, an increase in the mycobacterial fitness for INH-R can lead to a potential increase in transmission rate but a decrease in the amplification rate. In case of RIF-R strains, the fitness cost was similar for both the settings. The slight difference in the amplification rate can be explained by the higher CDR observed in Viet Nam. Factors such as age, behavior, and gender ratio of the overall population, could possibly be a reason to see this type of differences. But it is beyond the scope of this paper to predict the effects of these external influences on this model.

Conclusion
While rapid molecular diagnostics will continue to be important for programmatic adoption, it is also important to recognize that the principle of unrecognized resistance amplification demonstrated here can be repeated for any resistance not routinely addressed in diagnostic algorithms. It is therefore essential to incorporate genome sequencing into surveillance programs, to maximize the clinical and public health benefits [51]. With recent developments in next generation sequencing techniques, we have now have high-throughput diagnostic tools for the detection of DR-TB which are both