Inferences on the acquisition of multidrug resistance in \emph{Mycobacterium tuberculosis} using molecular epidemiological data

TL;DR

This study uses molecular epidemiological data and modeling to understand how multidrug resistance in tuberculosis develops and spreads, focusing on mutation rates, resistance acquisition pathways, and transmission versus evolution contributions.

Contribution

It introduces a novel model combining resistance evolution and VNTR loci tracking, applying approximate Bayesian computation to analyze incomplete data on MDR TB in Bolivia.

Findings

Double resistance can evolve directly from sensitive states within hosts.

Differences in mutation rates influence resistance patterns at the population level.

A significant proportion of MDR TB cases result from transmission rather than resistance acquisition.

Abstract

We investigate the rates of drug resistance acquisition in a natural population using molecular epidemiological data from Bolivia. First, we study the rate of direct acquisition of double resistance from the double sensitive state within patients and compare it to the rates of evolution to single resistance. In particular, we address whether or not double resistance can evolve directly from a double sensitive state within a given host. Second, we aim to understand whether the differences in mutation rates to rifampicin and isoniazid resistance translate to the epidemiological scale. Third, we estimate the proportion of MDR TB cases that are due to the transmission of MDR strains compared to acquisition of resistance through evolution. To address these problems we develop a model of TB transmission in which we track the evolution of resistance to two drugs and the evolution of VNTR loci. However, the available data is incomplete, in that it is recorded only {for a fraction of the population and} at a single point in time. The likelihood function induced by the proposed model is computationally prohibitive to evaluate and accordingly impractical to work with directly. We therefore approach statistical inference using approximate Bayesian computation techniques.

Tables

Table 1. Table 1: Molecular data set compiled from Monteserin et al. (2013). All isolates were classified according to their genotype and resistance profile. The symbol “a” represents 10 repeat units and “-” represents missing data. The entries in the four columns sum to the total number of isolates, 100.

Genotype	$\mathrm{𝟎}$	INH	RIF	MDR		Genotype	$\mathrm{𝟎}$	INH	RIF	MDR
253533233433427	4	0	0	0		243413342212437	1	0	0	0
253533233433327	3	0	0	0		233312442212437	1	0	0	0
253533233433527	11	1	0	0		233313441212437	1	0	0	0
253533233433525	3	0	0	0		233413442212248	1	0	0	0
143533233433527	2	1	0	0		233413442212249	1	0	0	0
253333244232232	1	0	0	0		233213442212349	1	0	0	0
25333324423-232	1	0	0	0		231413542212335	1	0	0	0
254333243232342	0	2	0	2		232433242212436	1	0	0	0
263532232423139	3	0	0	0		234413442212436	1	0	0	0
223413442212437	2	0	0	0		434433452212427	1	0	0	0
233413542212347	0	1	0	2		256432342122237	2	1	0	0
244333244232332	0	0	0	1		256433342123236	2	0	0	0
244333244232322	1	0	0	0		247432342122136	1	0	0	0
245333244242332	1	0	0	0		268432252122227	0	1	0	0
254333244232232	0	0	0	1		268632252122227	1	0	0	0
254333244232332	1	0	0	0		221313352122338	0	0	0	1
254333244242332	1	0	0	0		263532233423148	1	0	0	0
253333244242232	1	0	0	0		360332233423138	1	0	0	0
252333243232232	0	0	0	1		263513233523344	1	0	0	0
252333243232332	0	0	0	1		253523233433527	1	0	0	0
251333243242332	1	0	0	0		253533232433527	1	0	0	1
252333243262222	0	0	0	1		253533232433427	1	0	0	0
244233234222322	1	0	0	0		253523133433527	1	0	0	0
233373242232325	1	0	0	0		253533133433527	1	0	0	0
252343242232524	1	0	0	0		353533233433427	0	0	0	1
25233234423251a	1	0	0	0		253533233433837	1	0	0	0
35234234423251a	1	0	0	0		253533233433237	1	0	0	0
233413442212338	0	1	0	0		254533233433537	0	0	0	1
233413442212335	1	0	0	0		253533233433536	1	0	0	0
233413442212337	1	0	0	0		252533233433428	1	0	0	0
21341344221233a	1	0	0	0		253534233433325	1	0	0	0
213413442212327	0	0	0	1		243533232433737	1	0	0	0
233413442212437	3	0	0	0		242433433433436	1	0	0	0

Table 2. Table 2: Summary of linear algebra notation.

Symbol	Meaning
$𝐕$	$G\times L$ matrix describing the VNTR genotypes.
$𝐔$, $𝐓$, $𝐖$	$G\times 4$ matrices of untreated, treated and total cases respectively, with columns corresponding to resistance profiles.
${𝐔}_k$, ${𝐓}_k$, ${𝐖}_k$	$G\times 1$ column vector for resistance profile $k$ of untreated, treated and total cases respectively.
$U_k$, $T_k$, $W_k$	Total population sizes of untreated, treated and total with resistance profile $k$.
${𝐔}_{i,k}$, ${𝐓}_{i,k}$, ${𝐖}_{i,k}$	$(i,k)$ entries of the matrices $𝐔$, $𝐓$, $𝐖$: the number of cases in each category with genotype $i$ and resistance profile $k$.
${\mathrm{𝟏}}_i$	$i\times 1$ column vector whose entries are all 1.
$e_i$	Column vector whose entries are 0 except for 1 in position $i$. Dimension determined by context.

Table 3. Table 3: Summary of model parameters. The top set of parameters are given fixed values, whereas the bottom set of parameters are allocated prior distributions and estimated using ABC. Fixed values and priors are justified in Section  4.2 . Rates are in units of per capita per year, but the time unit is set to 1/12 year in simulations. ∗ Specifically, β 0 subscript 𝛽 0 \beta_{0} is assumed to follow a (shifted) Gamma prior defined as β 0 − 0.68 ∼ Gamma ​ ( shape = 2 , rate = 0.73 ) similar-to subscript 𝛽 0 0.68 Gamma formulae-sequence shape 2 rate 0.73 \beta_{0}-0.68\sim\text{Gamma}(\text{shape}=2,\text{rate}=0.73) . See Section  4.2 for further details.

Symbol	Meaning	Fixed
		value
$\delta$	rate of death and natural recovery	0.52
${\gamma }_0$	cure rate for resistance profile 0, when treated	0.5
${\gamma }_{\mathrm{INH}}$	cure rate for resistance profile INH, when treated	0.25
${\gamma }_{\mathrm{RIF}}$	cure rate for resistance profile RIF, when treated	0.25
${\gamma }_{\mathrm{MDR}}$	cure rate for resistance profile MDR, when treated	0.05
$N$	total susceptible population size in absence of disease	${10}^4$
$\tau$	treatment and detection rate	0.5
$c$	cost of resistance	0.1
Symbol	Meaning	Prior
${\beta }_0$	transmission rate for resistance profile $0$	Gamma∗
$\mu$	mutation rate of VNTR per locus per unit time	$U(0,1)$
${\rho }_{\mathrm{INH}}$	rate of acquisition of resistance to INH	$U(0,1)$
${\rho }_{\mathrm{RIF}}$	rate of acquisition of resistance to RIF	$U(0,1)$
${\rho }_{\mathrm{MDR}}$	rate of acquisition of resistance to INH and RIF	$U(0,1)$

Table 4. Table 4: ABC posterior means with lower and upper limits of the 95% HPD (highest posterior density) credibile intervals for each parameter of each fitted model.

	${\rho }_{\mathrm{INH}}$	${\rho }_{\mathrm{RIF}}$	${\rho }_{\mathrm{MDR}}$	$\mu$	${\beta }_0$
Full model
Posterior mean	$1.14\times {10}^{-3}$	$1.67\times {10}^{-4}$	$2.62\times {10}^{-4}$	$1.64\times {10}^{-3}$	2.85
CI lower limit	$3.40\times {10}^{-4}$	$3.82\times {10}^{-6}$	$3.93\times {10}^{-6}$	$1.11\times {10}^{-3}$	0.97
CI upper limit	$1.94\times {10}^{-3}$	$4.28\times {10}^{-4}$	$5.81\times {10}^{-4}$	$2.40\times {10}^{-3}$	5.33
Submodel l
Posterior mean	$1.60\times {10}^{-3}$	$6.37\times {10}^{-4}$	–	$1.59\times {10}^{-3}$	3.29
CI lower limit	$4.55\times {10}^{-4}$	$1.27\times {10}^{-4}$	–	$1.03\times {10}^{-3}$	1.20
CI upper limit	$2.49\times {10}^{-3}$	$1.24\times {10}^{-3}$	–	$2.19\times {10}^{-3}$	5.78
Submodel 2
Posterior mean	$3.46\times {10}^{-4}$	–	$1.56\times {10}^{-4}$	$1.70\times {10}^{-3}$	2.81
CI lower limit	$7.26\times {10}^{-5}$	–	$6.62\times {10}^{-7}$	$1.10\times {10}^{-3}$	0.86
CI upper limit	$6.90\times {10}^{-4}$	–	$3.76\times {10}^{-4}$	$2.54\times {10}^{-3}$	5.20

Table 5. Table 5: Contributions to MDR-TB from alternative sources. This table contains the posterior medians and means and lower and upper limits of the 95% HPD credibility intervals for the proportion of double resistance cases originating from each possible source.

Source	Median	Mean	95% Credible Interval
Transmission	0.9975	0.9655	(0.7826, 0.9999)
Conversion in one step	0.0023	0.0284	(0.0000, 0.1667)
Conversion in two steps	0.0000	0.0060	(0.0000, 0.0073)

Equations

$$S=N-1_G^{\top }\mathbf{W}{1}_4,$$

$$p=\frac{1}{N}{1}_G^{\mathrm{T}}\mathbf{W}{D}_{\beta }{1}_4$$

$$B\sim \text{Binomial}(S,\tilde{p}).$$

$$R\sim \mathrm{Binomial}(U,\delta )$$

$$C_k\sim \mathrm{Binomial}(T_k,\delta +{\gamma }_k),$$

$$D\sim \mathrm{Binomial}(U,\tau ).$$

$$A_{i,0\to \ast }$$

$$A_{i,\mathrm{INH}\to \mathrm{MDR}}$$

$$A_{i,\mathrm{RIF}\to \mathrm{MDR}}$$

$${\Delta }_{\rho }=\left({\Delta }_0\mid {\Delta }_{\mathrm{INH}}\mid {\Delta }_{\mathrm{RIF}}\mid {\Delta }_{\mathrm{MDR}}\right).$$

$${\Delta }_0$$

$${\Delta }_{\mathrm{INH}}$$

$${\Delta }_{\mathrm{RIF}}$$

$${\Delta }_{\mathrm{MDR}}$$

$$\frac{dU}{dt}$$

$$\hat{U}=N\left(1-\frac{\delta }{{\beta }_0}\right)\text{mbox}and\hat{S}=\frac{\delta N}{\beta }$$

$${\pi }_{ABC}(\theta |s_{obs})\propto \pi (\theta )\int K_h(\Vert s-s_{obs}\Vert )p(s|\theta )ds,$$

$$(1_g|\mathbf{X}{)}^{\top }(1_g|\mathbf{X}),$$

$$\Vert s-s_{obs}\Vert =\Vert S(\mathbf{X})-S({\mathbf{X}}_{obs})\Vert ={\left(\sum _{j=1}^{\dim (s)}{\left[S(\mathbf{X}{)}_j-S({\mathbf{X}}_{obs}{)}_j\right]}^{\frac{1}{2}}\right)}^2,$$

$$p_k=\frac{{\delta }_r+{\gamma }_k}{{\delta }_d+{\delta }_r+{\gamma }_k}.$$

$${\beta }_0-0.68\sim \text{Gamma}(\text{shape}=2,\text{rate}=0.73)$$

Full text

Inferences on the acquisition of multidrug resistance in Mycobacterium tuberculosis using molecular epidemiological data

Guilherme S. Rodrigues111School of Mathematics and Statistics, University of New South Wales, Australia.

222CAPES Foundation, Ministry of Education of Brazil, Brasília - DF 70040-020, Brazil, Andrew R. Francis333Centre for Research in Mathematics, Western Sydney University, Australia.,

Scott A. Sisson1 and Mark M. Tanaka444School of Biotechnology and Biological Sciences, and Evolution & Ecology Research Centre, University of New South Wales, Australia.

1 Introduction

Tuberculosis (TB) is a lung disease caused by the bacterium Mycobacterium tuberculosis which kills around 1.5 million people each year and remains a serious challenge for global public health [\citeauthoryearWHOWHO2015]. Antibiotic drugs for treating TB have been available since the mid 20th century, and currently implemented strategies for TB control rely on the efficacy of these drugs. Treatment of TB involves combination therapy – in which multiple drugs are administered together in part to improve killing efficacy. The “first-line” drugs used in combination to treat tuberculosis are rifampicin, isoniazid, pyrazinamide, ethambutol and streptomycin.

As with most other pathogens, resistance to antibiotic drugs has rapidly evolved in M. tuberculosis. Streptomycin was the first of the first-line drugs to be developed and deployed in 1943, but resistance was observed before the end of that decade [\citeauthoryearMitchisonMitchison1951, \citeauthoryearGillespieGillespie2002]. Of particular concern is the rise of bacterial strains resistant to multiple drugs, as cases caused by them are difficult to treat successfully. Multidrug resistance (MDR) is defined as resistance to both rifampicin and isoniazid. These are the two most effective drugs against tuberculosis (when the strain is not resistant). Currently, 3.3% of new TB cases are multi-drug resistant [\citeauthoryearWHOWHO2015]. The occurrence of MDR-TB strains that have additional resistance (called extensively drug resistant, XDR and totally drug resistant, TDR) are particularly problematic and have the potential to cause large outbreaks that are difficult to control \shortciteGandhi.06.Extensively. A better understanding of how multiple drug resistance evolves would aid efforts to contain resistance and control tuberculosis.

Genetic studies have established that many independent mutation events have led to resistance [\citeauthoryearRamaswamy and MusserRamaswamy and Musser1998]. Although this suggests that mutation of genes is an important source of resistance, model-based analysis of molecular data has revealed that among resistant cases, most are due to the transmission of already resistant bacteria \shortciteLuciani.09.epidemiological. It is therefore of interest to investigate whether or not this finding also holds for multi-drug resistant tuberculosis.

The rates at which resistance evolves against different drugs vary. For instance, isoniazid resistance is known to be acquired faster than rifampicin resistance \shortciteFord.13.Mycobacterium,Gillespie.02.Evolution,Nachega.03.Tuberculosis. The rates of mutation to resistance per cell generation are low in absolute value; for example for isoniazid the rate is around $3\times 10^{-8}$ and for rifampicin it is around $2\times 10^{-10}$ [\citeauthoryearDavidDavid1970, \citeauthoryearGillespieGillespie2002], although there is a high degree of variation across different lineages of M. tuberculosis \shortciteFord.13.Mycobacterium. One might therefore expect that double resistance of these drugs (MDR) evolves at an exceedingly low rate [\citeauthoryearNachega and ChaissonNachega and Chaisson2003]. However, MDR strains often occur at appreciable frequencies \shortciteZhao.12.National,Anderson.14.Transmission and a recent study has presented a theoretical model showing how double resistance can evolve rapidly within hosts \shortciteColijn.11.Spontaneous. It would be useful to establish whether such fast direct acquisition of double resistance can be detected in bacterial isolates from epidemiological studies.

To characterise patterns of TB transmission and drug resistance in a given geographic region, bacterial isolates from TB patients are often genotyped using molecular markers known as variable numbers of tandem repeats (VNTRs) which are repeated genetic sequences that exhibit variation across isolates. The source of this variation is mutation at the VNTR genetic loci which leads to the expansion or contraction of repeat numbers at those loci (Figure 1). A scheme for discriminating effectively among a set of isolates involves considering repeat numbers at multiple VNTR sites. This molecular typing scheme is called multi-locus VNTR analysis (MLVA); in the context of tuberculosis epidemiology it is often known as mycobacterial interspersed repetitive units-VNTR (MIRU-VNTR) \shortciteMazars.01.High-resolution,Supply.06.Proposal. Typing techniques such as MLVA have been useful for tracking particular strains and understanding how drug resistance evolves and disseminates at the epidemiological level \shortciteMonteserin.13.genotypes,Anderson.14.Transmission.

Here, we investigate the rates of drug resistance acquisition in a natural population using molecular epidemiological data from Bolivia \shortciteMonteserin.13.genotypes. First, we study the rate of direct acquisition of double resistance from the double sensitive state within patients and compare it to the rates of evolution to single resistance. In particular, we address whether or not double resistance can evolve directly from a double sensitive state within a given host. Second, we aim to understand whether the differences in mutation rates to rifampicin and isoniazid resistance translate to the epidemiological scale. Third, we estimate the proportion of MDR TB cases that are due to the transmission of MDR strains compared to acquisition of resistance through evolution. To address these problems we develop a model of TB transmission in which we track the evolution of resistance to two drugs and the evolution of VNTR loci. However, the available data (see Section 2) is incomplete, in that it is recorded only for a fraction of the population and at a single point in time. The likelihood function induced by the proposed model is computationally prohibitive to evaluate and accordingly impractical to work with directly. We therefore approach statistical inference using approximate Bayesian computation techniques.

2 Data

The data set we use is taken from a study of tuberculosis in Bolivia \shortciteMonteserin.13.genotypes. Bolivia has a population of 11 million people and a TB incidence of 120 per 100,000 per year. This rate is comparable to the global incidence of TB (133 per 100,000 per year) and to the rate in Peru, but is 3–6 times the TB incidence in neighbouring countries Brazil, Paraguay, Uruguay, Argentina and Chile [\citeauthoryearWHOWHO2015]. In the molecular epidemiological study, the investigators genotyped 100 isolates collected in 2010, which represented an estimated 1.1% of the cases in Bolivia at the time of the study \shortciteMonteserin.13.genotypes. Each isolate was tested for drug sensitivity to five drugs. Here, we focus on resistance against the two drugs isoniazid and rifampicin used to define multidrug resistance. Of the 100 isolates, 14 were found to be MDR, that is, resistant to both of these drugs, 78 were sensitive to both drugs and the remaining 8 were resistant to isoniazid but sensitive to rifampicin. No isolates were resistant to rifampicin while being sensitive to isoniazid.

In addition to these drug resistance profiles, each isolate was genotyped using 15 VNTR loci. For example, an isolate in the data set, which was resistant to isoniazid but sensitive to rifampicin, had the following 15 repeat numbers for its 15 VNTR loci: $143533233433527$, which together constitute its genotype. Variation in these genotypes occurs through a process of mutation in which repeat numbers increase or decrease (see Figure 1).

Let $g$ be the number of distinct genotypes present in a sample, and label the resistance profiles by ([math], INH, RIF, MDR), where [math] denotes sensitivity to both drugs, INH denotes resistance to isoniazid and sensitivity to rifampicin, RIF denotes resistance to rifampicin and sensitivity to isoniazid, and MDR denotes resistance to both drugs. The observed data $\mathbf{X}_{obs}$ are then a $g\times 4$ matrix of counts, such that each row gives the distribution of isolates across the four resistance profiles for a given genotype and each column gives the distribution of isolates across genotypes for a given resistance profile. The sum of entries in a particular row is the number of isolates with that genotype, while the sum of entries in a particular column is the number of isolates with that resistance profile. The data set also includes a $g\times 15$ matrix of repeat numbers from the VNTR genotyping.

The Bolivian data set is displayed in full in Table 1, which shows all $g=66$ distinct genotypes and classifies all 100 isolates according to genotype and resistance profile. The $\mathbf{X}_{obs}$ matrix is formed by combining the [math], INH, RIF and MDR columns.

3 Model

In this Section we introduce a model that incorporates both VNTR-based genotyping and drug resistance states. The dynamic variables of the model correspond to numbers of cases of untreated and treated tuberculosis, their resistance states and VNTR genotypes associated with these infections in the population. We will now briefly describe processes involved in the model, and provide further details in the following Subsections.

An untreated case of TB can become detected and treated, and treatment involves a combination of drugs including the two in question. Drug sensitive strains can acquire resistance under treatment with some probability and thereby change their resistance state. Treated and untreated cases can infect susceptible individuals and convert them to untreated cases. We disregard latent infections for simplicity (although latency is an important feature of the natural history of tuberculosis), and focus on active infections which are the larger source of new infections. Treated and untreated individuals can also recover or die. Treated individuals enjoy an additional probability of recovery that depends on the efficacy of the drugs, which in turn depends on the sensitivity or resistance of the infecting strain. Treated and untreated cases are also associated with a VNTR genotype, and this genotype evolves over time according to a stepwise mutation process for each locus. Figure 2 shows the broad structure of the model with respect to treatment and resistance states, while suppressing details of transmission, recovery, death and mutation of the VNTR loci.

At the end of the period of evolution, a simple random sample of 100 isolates is taken without replacement from the population, which matches the sample size of the Bolivian dataset. This provides a full description of the generative process for the observable data.

Let $G$ be the number of distinct genotypes in the population (the number of distinct genotypes in the sample is $g$) and $L$ be the number of VNTR loci used in the genotyping scheme. For the Bolivian dataset $L=15$. In the model, the variable $G$ is unknown and varies dynamically. We maintain three matrices which change through time: a $G\times L$ matrix, $\bf V$, which describes the VNTR genotypes; a $G\times 4$ matrix, $\bf U$, which describes the numbers of untreated cases of tuberculosis classified according to VNTR genotype and resistance state; and a $G\times 4$ matrix $\bf T$ which describes the numbers of treated cases of tuberculosis, again classified according to VNTR genotype and resistance state. It will be useful to define a $G\times 4$ matrix, $\bf W$, whose entries are the total numbers of both treated and untreated cases: ${\bf W}={\bf U}+{\bf T}$.

As it will be helpful to be able to pick out columns of these matrices, we adopt notation for the standard basis vectors of $\mathbb{R}^{n}$. Let $e_{i}$ denote the $i$-th basis (column) vector, so that $e_{i}=(0,\cdots,0,1,0\cdots,0)^{\top}$, with the 1 in the $i$-th position. This allows us, for instance, to write the columns of the matrix $\mathbf{T}$ corresponding to each resistance state as $\mathbf{T}_{0}=\mathbf{T}\,e_{1}$, $\mathbf{T}_{\mathrm{INH}}=\mathbf{T}\,e_{2}$, $\mathbf{T}_{\mathrm{RIF}}=\mathbf{T}\,e_{3}$ and $\mathbf{T}_{\mathrm{MDR}}=\mathbf{T}\,e_{4}$, with similar notation for other matrices (note that the dimension of the $e_{i}$ is left open but inferred from the matrix multiplication; in this case they are in $\mathbb{R}^{4}$).

Further, writing $\mathbf{1}_{i}$ for the column vector in $\mathbb{R}^{i}$ whose entries are all $1$, then the product $\mathbf{T}\,\mathbf{1}_{4}$ is a $G\times 1$ column vector whose entries are the numbers of individual cases of each VNTR genotype in the treated population, and the product $\mathbf{1}_{G}^{\top}\,\mathbf{T}\,\mathbf{1}_{4}$ is the sum of all the entries in $\mathbf{T}$ (the size of the treated population). Thus, we can write the size of the susceptible population, $S$, as

[TABLE]

where $\mathbf{1}_{G}^{\top}\,\mathbf{W}\,\mathbf{1}_{4}$ is the size of the infected population, and where $N$ is the total population size which remains constant. We treat $N$ as modelling the set of all individuals who come in contact with infectious cases and so we exclude individuals who either do not encounter infectious cases or are otherwise protected from infection. This variable therefore may be smaller than the actual population size.

The components of each vector $\mathbf{T}_{k}$ for $k=0$, INH, RIF or MDR, are integers, representing the number of individual cases for each genotype. In the schematic diagram of the model in Figure 2, we use $T_{k}=\mathbf{1}_{G}^{\top}\,\mathbf{T}_{k}$ to represent the total population number of treated individuals with resistance state $k$, with similar notation $U_{k}=\mathbf{1}_{G}^{\top}\,\mathbf{U}_{k}$ to represent the untreated populations. The matrix notation is gathered and shown in Table 2.

The arrows between populations in Figure 2 represent the directional rates of detection and treatment $\tau$ and acquisition of resistance to each drug or set of drugs, so that $\rho_{\mathrm{INH}}$ and $\rho_{\mathrm{RIF}}$ represent rates of acquisition of resistance to isoniazid and rifampicin respectively and $\rho_{\mathrm{MDR}}$ the rate of double acquisition.

In this model time is discrete, and during each time step the following events takes place in sequence.

Disease transmission giving rise to new cases; 2. 2.

Natural recovery, cure or death of cases; 3. 3.

Detection of cases which are then treated with drugs; 4. 4.

Conversion among resistant profiles in treated cases due to acquisition of resistance; and 5. 5.

Mutation of the genetic marker (multiple VNTR loci).

The remainder of this Section provides details of how each of these events are modelled. Readers wishing to focus on the statistical aspects of the ABC inference can skip these subsections and go directly to Section 4.

We regard the above process as a discrete-time stochastic model rather than a discrete-time approximation of a continuous-time stochastic process with rates approximating probabilities, although the latter interpretation becomes more appropriate as the time step length decreases. Here, rates of events will be treated as probabilities, which again is appropriate when time steps are short. The rate parameters are measured in years but we make time steps 1/12 of a year.

A summary of all model parameters, both fixed and to be estimated, and their meanings, are provided in Table 3.

3.1 New infections

In our model, new infections occur by mass action. The per capita rate at which a susceptible individual becomes infected by a case with resistance profile $k$ is given by $\beta_{k}/N$ times the number of infected cases in state $k$. The transmission parameters $\beta_{k}$ are scaled by $1/N$ for convenience since realistic values of $\beta_{k}/N$ are typically very small, and this ensures that the $\beta_{k}$ are on the natural “per person per unit time” scale.

The acquisition of resistance to antibiotics often comes at a cost to the fitness of the bacterium, and we implement this fitness cost by assuming the transmission rate is lower for cases that carry resistance. Specifically, we assume a cost of “$c$ per drug”, so that if $\beta_{0}$ is the transmission rate of sensitive cases, then cases resistant to one drug transmit at rate $\beta_{\mathrm{INH}}=\beta_{\mathrm{RIF}}=(1-c)\beta_{0}$ and cases resistant to two drugs at $\beta_{\mathrm{MDR}}=(1-c)^{2}\beta_{0}$. Here, $c=0.1$ is considered known and fixed based on previous analyses of molecular epidemiological data \shortciteLuciani.09.epidemiological.

We now construct an expression for the average transmission probability across the infected population. The matrix $\mathbf{W}$ records all infected cases with different resistance states in each column, and the individuals corresponding to these columns have different transmission rates $\beta=(\beta_{0},\beta_{\mathrm{INH}},\beta_{\mathrm{RIF}},\beta_{\mathrm{MDR}})^{\top}$. If we write $D_{\beta}$ for the diagonal matrix whose entries are from $\beta$, then the matrix $\mathbf{W}\,D_{\beta}$ is the infected population matrix $\mathbf{W}$ whose columns have been scaled by the entries of $\beta$ (the relevant transmission rates). The expression

[TABLE]

then gives the average transmission rate per susceptible individual. Since the population size $N$ is usually large and the time steps are short, the value for $p$ will nearly always be small. Accordingly, and to ensure that it does not exceed 1, we model the probability of transmission per susceptible individual as $\tilde{p}=\min\{1,p\}$.

At each time step the number $B$ of new infections is a random variable distributed as

[TABLE]

These $B$ new infections are then allocated across VNTR genotypes and resistance profiles according to the proportions represented by the matrix $\mathbf{W}\,D_{\beta}$. That is, a multinomial random sample distributes $B$ according to the existing infected population and their relative transmission rates, so that the resulting allocation is a $G\times 4$ matrix $\Delta_{\beta}$. Finally, as new infections are all assumed to be initially undetected, they are allocated to the untreated subpopulation, so that the matrix $\mathbf{U}$ is updated to $\mathbf{U}\rightarrow\mathbf{U}+\Delta_{\beta}$.

3.2 Cure, recovery and death

Infected individuals who are untreated (the population represented by the counts $\mathbf{U}$) recover or die at rate $\delta=\delta_{r}+\delta_{d}$ per case per time unit, where $\delta_{r}>0$ is the rate of recovery and $\delta_{d}>0$ is the rate of death due to any cause. The rate of cure due to successful treatment may vary according to resistance profile, so this rate is given by $\gamma_{k}$ for $k=0$, INH, RIF, MDR. The number of cures, recoveries and deaths in a time step is given by

[TABLE]

for the untreated population, where $U=\mathbf{1}_{G}^{\top}\,\mathbf{U}\,\mathbf{1}_{4}$ is the total number of all untreated cases, and

[TABLE]

for the treated population, where $T_{k}$ is the number of treated individuals with resistance profile $k$ (as defined at the start of this Section). The $R$ untreated recovered individuals are distributed across both VNTR genotypes and resistance profiles with a multinomial distribution according to the counts given in $\mathbf{U}$. These are recorded in the $G\times 4$ update matrix $\Delta_{\delta}$ (so that the sum of the entries in $\Delta_{\delta}$ is $R=\mathbf{1}_{G}^{\top}\,\Delta_{\delta}\,\mathbf{1}_{4}$). Similarly, the $C_{k}$ treated recovered individuals of resistance profile $k$ are distributed across the VNTR genotypes according to the distribution observed in $\mathbf{T}_{k}$. These recovered counts for all resistance profiles are recorded in the $G$ update matrix $\Delta_{\delta+\gamma},$ which is constructed from the column vectors of recovered treated counts for profile $k$ in the order $k=0$, INH, RIF and MDR. The matrices $\mathbf{U}$ and $\mathbf{T}$ are then updated to $\mathbf{U}\rightarrow\mathbf{U}-\Delta_{\delta}$ and $\mathbf{T}\rightarrow\mathbf{T}-\Delta_{\delta+\gamma}$ respectively. If the last instance of any genotype is removed by cure, recovery or death, the matrices $\bf U,T,V$ are adjusted by removing the rows corresponding to those genotypes, and the number of genotypes is updated with $G\rightarrow G-1$.

Similarly to the case of new infections (Section 3.1), we assume that the recovery rate due to treatment depends only on the number of drugs the infecting strain is resistant to. Specifically, this implies that $\gamma_{\mathrm{INH}}=\gamma_{\mathrm{RIF}}$.

3.3 Detection and treatment

In this model, the detection of cases and the commencement of treatment are combined as a single process. Detected cases are transferred from the untreated class to the treated class. We denote this combined detection and treatment rate, per case, per unit time, as $\tau>0$. With this rate we draw $D$ individuals to transfer between untreated and treated populations, where

[TABLE]

These $D$ individuals are then allocated across VNTR genotypes and resistance profiles according to the observed distribution of untreated cases, $\mathbf{U}$. As before, this results in a $G\times 4$ update matrix $\Delta_{\tau}$, which we use to update $\mathbf{U}\rightarrow\mathbf{U}-\Delta_{\tau}$ and $\mathbf{T}\rightarrow\mathbf{T}+\Delta_{\tau}$.

3.4 Acquisition of drug resistance

Individual treated cases are able to convert from one resistance profile to another through adaptive evolution. That is, under drug treatment, natural selection acts to favour increasing levels of resistance. As a result of this process, individuals may move from the $k=0$ resistance profile (sensitive to both drugs) to one of the other three resistance profiles: ${\mathrm{INH}}$, ${\mathrm{RIF}}$, or ${\mathrm{MDR}}$ (resistance to one or both drugs). Individuals may also move from resistance to exactly one of the drugs (${\mathrm{INH}}$ or ${\mathrm{RIF}}$) to the multiple drug resistance profile ${\mathrm{MDR}}$. We respectively denote the rate of acquisition of resistance to INH or RIF by $\rho_{\mathrm{INH}}$ and $\rho_{\mathrm{RIF}}$, and denote the rate of acquisition of resistance from individuals in the sensitive population to both drugs simultaneously by $\rho_{\mathrm{MDR}}$. These conversions and rates are illustrated schematically in Figure 2.

To model resistance acquisition, we select individuals to move between resistance profiles in the treated population i.e. between columns in the matrix $\mathbf{T}$. Acquiring resistance to the drug rifampicin will result in individuals moving from the column $\mathbf{T}_{0}$ to $\mathbf{T}_{\mathrm{RIF}}$, and from $\mathbf{T}_{\mathrm{INH}}$ to $\mathbf{T}_{\mathrm{MDR}}$ at a rate $\rho_{\mathrm{RIF}}$. Similarly, acquiring resistance to the drug isoniazid results in individuals moving from the column $\mathbf{T}_{0}$ to $\mathbf{T}_{\mathrm{INH}}$, and from $\mathbf{T}_{\mathrm{RIF}}$ to $\mathbf{T}_{\mathrm{MDR}}$ at a rate $\rho_{\mathrm{RIF}}$. Simultaneous acquisition of resistance to both drugs moves individuals from the column $\mathbf{T}_{0}$ to $\mathbf{T}_{\mathrm{MDR}}$ at the rate $\rho_{\mathrm{MDR}}$. These movements occur between columns but not across rows (infections do not change VNTR genotypes through this process).

Mechanistically, we can obtain the number of cases of genotype $i$ transitioning from resistance profile $k$ to resistance profile $k^{\prime}$, denoted $A_{i,k\rightarrow k^{\prime}}$, as

[TABLE]

where $A_{i,0\rightarrow*}=(A_{i,0\rightarrow\mathrm{INH}},A_{i,0\rightarrow\mathrm{RIF}},A_{i,0\rightarrow\mathrm{MDR}},A_{i,0\rightarrow 0})^{\top}$ is the vector of cases transitioning from sensitivity, $\mathbf{T}_{i,k}$ is the entry of the matrix $\mathbf{T}$ corresponding to the genotype $i$ and resistance profile $k$ (see Table 2), and $\rho_{0\rightarrow*}=(\rho_{\mathrm{INH}},\rho_{\mathrm{RIF}},\rho_{\mathrm{MDR}},1-\sum_{k}\rho_{k})^{\top}$ is the vector of probabilities of these events.

If we denote $\Delta_{k\rightarrow k^{\prime}}$ as column vectors of counts of movements from resistance profile $k$ to $k^{\prime}$ across all $G$ genotypes, we can then construct the overall $G\times 4$ update matrix

[TABLE]

from the column vectors $\Delta_{k}$, which denote the total population change for resistance profile $k$, where

[TABLE]

The population of treated cases is then updated to $\mathbf{T}\rightarrow\mathbf{T}+\Delta_{\rho}$.

3.5 Mutation of the marker

The set of $L=15$ VNTR loci constitute the genetic marker used to genotype bacterial isolates (see Section 2). Each genotype is a list of numbers of tandem repeat units at the $L$ loci. The states of all VNTRs in the infected population are given by the $G\times L$ matrix $\bf V$ with elements $V_{ij}$ describing the repeat number of locus $j$ in genotype $i$. Each locus mutates through a stepwise mutation process at rate $\mu$ per locus per case per unit time. When mutation occurs, the repeat number $V_{ij}$ at a locus $j$ of genotype $i$ changes by $+1$ or $-1$, each with probability 0.5. A repeat number of 1 is treated as an absorbing boundary (i.e. there is zero probability of the repeat number increasing from 1 to 2) because at state 1 there is no longer a genetic sequence that is tandemly repeated and no mechanism such as replication slippage acts to expand it from 1 to 2.

Mutation of the marker has the effect of moving cases between the rows of the matrix $\mathbf{W}$. We first identify the number of mutation events in the population, $M$, where $M\sim\mathrm{Binomial}(S,\mu)$ and $S=N-\mathbf{1}_{G}^{\mathrm{T}}\,\mathbf{W}\,\mathbf{1}_{4}$ is the size of the susceptible population (see Section 3). The $M$ cases are then distributed across the population of VNTR genotypes and resistance profiles, according to the entries of the matrices $\mathbf{T}$ and $\mathbf{U}$. Each individual case undergoing mutation corresponds to a specific entry in either $\mathbf{T}$ or $\mathbf{U}$. This entry is described by its VNTR genotype $\mathbf{V}_{i}=(V_{i,1},\ldots,V_{i,L})$ where $L=15$ for the Bolivian data, and its resistance profile, $k=0$, INH, RIF, MDR. The result of the mutation is a change to the VNTR genotype, which is represented by a change in the repeat number at a single locus, $V_{ij}$, by $\pm 1$. This may or may not result in a VNTR genotype that is already present in the population.

If the new VNTR genotype already appears as a row in the matrix $\mathbf{V}$ as an existing type in the data, then there is no change to $\mathbf{V}$. The matrix $\mathbf{T}$ or $\mathbf{U}$ on the other hand is changed by subtracting 1 from one entry and adding one to another entry in the same column (the resistance profile, $k$, does not change). In matrix terms, supposing the change is to a treated case, this can be described by updating $\mathbf{T}\rightarrow\mathbf{T}-e_{i,j}+e_{i,k}$, where $e_{i,j}$ is the matrix whose entries are zero except for a 1 in the $(i,j)$-th position, and where the VNTR genotype changes from row $j$ to row $k$.

If the new VNTR genotype does not already appear in the population, then the matrix $\mathbf{V}$ is expanded to include a new row describing the new genotype, so that $\mathbf{V}$ becomes a $(G+1)\times 15$ matrix. The update for $\mathbf{T}$ or $\mathbf{U}$ is the same as described above, except that now both matrices are $(G+1)\times 4$ dimensional. Subsequent to this update we increment $G\rightarrow G+1$. If mutation of a VNTR genotype removes the last instance of the original genotype from $\bf U$ and $\bf T$ the corresponding rows of matrices $\bf V$, $\bf U$ and $\bf T$ are deleted, requiring the update $G\rightarrow G-1$.

3.6 Initial conditions of the model

The model covers the period from when drugs are introduced at time $t=0$ to when sampling occurs. Since the main first-line anti-tuberculosis drugs were discovered/developed in the 1940s to early 1960s, we assumed treatment commenced around 1960 and ran the simulation for a period of 50 years. We assumed that both drugs, isoniazid and rifampicin, were introduced at the same time and are administered together in combination therapy. The standard course of treatment includes both drugs along with other first-line drugs [\citeauthoryearWHOWHO2015].

We assume that at the start of the process all cases are sensitive to both drugs and that the number of cases is at equilibrium in the absence of treatment and resistance. To compute this equilibrium state, we consider the differential equation describing the deterministic version of the model ignoring VNTR genotypes. Namely,

[TABLE]

where $S=N-U$ and $t$ indicates time. Setting $dU/dt$ to zero and solving for the dynamic variables we obtain equilibrium values of

[TABLE]

for $U>0$.

The basic reproduction number of a pathogen $R_{0}$ is defined to be the average number of new infectious cases caused by a single infection in a completely susceptible population. In our model, before there is any treatment, assuming all cases are doubly susceptible, a single case on average persists for $1/\delta$ years and generates $S\beta_{0}/N$ new cases per unit time but since $S=N$ in a wholly susceptible population then $R_{0}=\beta_{0}/\delta$.

All cases are initially untreated and sensitive. From time $t=0$ treatment in the population commences. To reintroduce into the model genetic variation at the marker loci, the initial distribution of genotype clusters is a random sample drawn from the infinite alleles model from population genetic theory \shortciteEwens.72.Sampling,Hubbell.01.unified,Luciani.08.Interpreting. The infinite alleles model depends on a single parameter, the diversity parameter, which we set to $2\hat{U}\mu L$ where $\hat{U}$ is the number of cases, taken from the equilibrium value described above, $\mu$ is the mutation rate per VNTR locus and $L$ is the number of VNTR loci used in genotyping isolates. To initialise the multi-locus VNTR genotypes, each genotype is a sequence of random integers, of length $L$, with each VNTR number $V_{ij}$ drawn from a discrete uniform distribution over $\{1,\ldots,10\}$.

Although the initial distribution of genotype clusters is set under the infinite alleles model, the mutation process for VNTRs brings the distribution in line with the stepwise model over time.

The initial conditions are a function of the parameters which are set according to the priors specified in Section 4.2.

4 Inference with approximate Bayesian computation

For the model in Section 3, when the data are only observed at a single point in time, the cost of evaluating the likelihood function is computationally prohibitive. This results from the “incomplete” nature of the observed data (see Section 2) in the sense that we only have access to a snapshot of the population, via the observed sample, at the time the study was conducted, with no direct measurements of the system as it progressed. Computing the likelihood then requires integrating over all potential trajectories the population could have gone through before reaching its final, observed state.

As such we adopt approximate Bayesian computation (ABC) methods as a means of performing Bayesian statistical inference for the unknown model parameters $\theta=(\beta_{0},\mu,\rho_{\mathrm{INH}},\rho_{\mathrm{RIF}},\rho_{\mathrm{MDR}})^{\top}$. As observed in other chapters in this Handbook, the ABC approximation to the true posterior distribution is given by

[TABLE]

where $\pi(\theta)$ is the prior distribution, $s=S(\mathbf{X})$ is a vector of summary statistics with $s_{obs}=S(\mathbf{X}_{obs})$, $p(s|\theta)$ is the computationally intractable likelihood function for the summary statistics $s$, and $K_{h}(u)=K(u/h)/h$ is a standard smoothing kernel with scale parameter $h>0$. In the following analyses we used the uniform kernel on $[-h,h]$ for $K_{h}(u)$. The quality of the ABC approximation depends on the information loss in the summary statistics $s$ over the full dataset $\mathbf{X}$, and the size of the kernel scale parameter $h$ with smaller $h$ producing greater accuracy and increased computational cost. Choice of both $s$ and $h$ are typically driven by the amount of expert knowledge and computation available for the analysis.

For the present analysis we implement a version of a simple ABC importance sampling algorithm, as outlined in the box. Given a suitable importance sampling distribution $q(\theta)$, the algorithm produces a set of weighted samples from the ABC approximation to the true posterior $(\theta^{(1)},w^{(1)}),\ldots,(\theta^{(\tilde{N})},w^{(\tilde{N})})\sim\pi_{ABC}(\theta|s_{obs})$. As with standard importance sampling, suitable choice of $q(\theta)$ is important to avoid high variance in the importance weights, and also to avoid needlessly generating datasets $s=S(\mathbf{X}^{(i)})$, $\mathbf{X}^{(i)}\sim p(\mathbf{X}|\theta)$ for which $s^{(i)}$ and $s_{obs}$ will never be close.

ABC Importance Sampling Algorithm

Inputs:

•

A target posterior density $\pi(\theta|\mathbf{X}_{obs})\propto p(\mathbf{X}_{obs}|\theta)\pi(\theta)$, consisting of a prior distribution $\pi(\theta)$ and a procedure for generating data under the model $p(\mathbf{X}_{obs}|\theta)$.

•

A proposal density $q(\theta)$, with $q(\theta)>0$ if $\pi(\theta|\mathbf{X}_{obs})>0$.

•

An integer $\tilde{N}>0$.

•

An observed vector of summary statistics $s_{obs}=S(\mathbf{X}_{obs})$.

•

A kernel function $K_{h}(u)$ and scale parameter $h>0$.

*Sampling:

*For $i=1,\ldots,\tilde{N}$:

Generate $\theta^{(i)}\sim q(\theta)$ from sampling density $q$.

Generate $\mathbf{X}^{(i)}\sim p(\mathbf{X}|\theta^{(i)})$ from the likelihood.

Compute the summary statistics $s^{(i)}=S(\mathbf{X}^{(i)})$.

Assign $\theta^{(i)}$ the weight $w^{(i)}\propto K_{h}(\|s^{(i)}-s_{obs}\|)\pi(\theta^{(i)})/q(\theta^{(i)})$.

*Output:

*A set of weighted parameter vectors $\{(\theta^{(i)},w^{(i)})\}_{i=1}^{\tilde{N}}$ $\sim$ $\pi_{ABC}(\theta|s_{obs})$.

To determine a suitable importance sampling distribution $q(\theta)$ we adopt a two stage procedure, following the approach of \citeNfearnhead+p12. In the first stage we perform a pilot ABC analysis using a sampling distribution that is diffuse enough to easily encompass the ABC posterior approximation obtained for a moderate value of the kernel scale parameter $h$. We specified $q(\theta)\propto\pi(\theta)I(\theta\in A)$ which is proportional to the prior, but restricted to the hyper-rectangle $A$. Here, $A$ is constructed as the smallest credible hyper-rectangle that we believe contains the ABC posterior approximation. As such, this $q(\theta)$ will identify the general region in which $\pi_{ABC}(\theta|s_{obs})$ is located. Specifically, for $\theta=(\beta_{0},\mu,\rho_{\mathrm{INH}},\rho_{\mathrm{RIF}},\rho_{\mathrm{MDR}})^{\top}$ we adopt $q(\theta)=\tilde{\pi}_{15}(\beta_{0})\times U(0,.005)\times U(0,.01)\times U(0,.005)\times U(0,.001)$, where $\tilde{\pi}_{15}(\beta_{0})$ is the prior $\pi(\beta_{0})$ for $\beta_{0}$ specified in Section 4.2, but truncated to exclude density above the point $\beta_{0}=15$.

For posterior distributions with strong dependence between parameters, defining $q(\theta)$ over such a hyper-rectangle may be inefficient as it will cover many regions of effectively zero posterior density. Accordingly we construct the sampling distribution for the second stage, with the lowest value of $h$, as a kernel density estimate of the previous ABC estimate of the posterior distribution: $q(\theta)=\sum_{i}w^{(i)}L(\theta|\theta^{(i)})$, where $L$ is a suitable kernel density (not to be confused with the kernel $K_{h}$). This approach follows the ideas behind the sequential Monte Carlo-based ABC samplers of \shortciteNScott2007 and others. At each stage the kernel scale parameter $h$ is decreased, and determined as the value which results in $\sim$ 2,000 posterior samples with non-zero weights, for the given computational budget.

To ensure greater efficiency at each stage we also performed a non-linear regression adjustment using a neural network with a single hidden layer (see \shortciteNPblum+f10,abc.package,Beaumont2002), as implemented in the R package abc. The adjustment used logistic transformations for the response.

For samples drawn from the final importance sampling distribution $q(\theta)$, the data generation procedure took on average $\sim 40$ seconds in R. This is computationally expensive from an ABC context, and could be reduced by recoding the simulator in a compiled language such as C, or by adapting the “lazy ABC” ideas of \shortciteNprangle16lazy to terminate early those simulations that are likely to be rejected. In this implementation we performed importance sampling from each distribution $q(\theta)$ in parallel on multiple nodes of a computational cluster.

4.1 Summary statistics

Considering the matrix structure of the observed data $\mathbf{X}_{obs}$ (see Section 2), we determine the information content in $\mathbf{X}$ as if it was the design matrix of a regression model and summarise it accordingly. Specifically, we define the summary statistics $s=S(\mathbf{X})$ to be the upper-triangular elements of the matrix

[TABLE]

where the vertical lines denote the addition of an extra column. The added columns of ones enriches the set of summary statistics by including the row and column totals of $\mathbf{X}$. Alternatively, these summary statistics can be described as:

1. i)

$g$: the number of distinct genotypes in the sample. 2. ii)

$n_{k}$: the number of isolates with resistance profile $k=$ [math], INH, RIF and MDR. 3. iii)

$c_{k,k^{\prime}}=(\mathbf{X}_{k})^{\top}\,\mathbf{X}_{k^{\prime}}$: the dot product between the resistance profiles of $k$ and $k^{\prime}$ within $\mathbf{X}$.

Note that these summary statistics are over specified in that $n_{0}+n_{\mathrm{INH}}+n_{\mathrm{RIF}}+n_{\mathrm{MDR}}$ equals the total number of isolates sampled from the population, which is known and equal to the number of isolates in the observed data sample (100 for the Bolivian data). Accordingly, and without loss of generality, we remove $n_{\mathrm{MDR}}$ as a summary statistic to avoid collinearity. In combination, this set of 14 summary statistics efficiently encapsulates the available information about the covariance structure of the original dataset $\mathbf{X}$, the distribution of the isolates among the different resistance profiles and the degree of diversity of isolates within the sample.

For the Bolivian dataset, there are $g=68$ distinct genotypes, $n_{0}=78$ sensitive isolates, $n_{\text{INH}}=8$ isolates resistant to isoniazid only, $n_{\text{RIF}}=0$ isolates resistant to rifampicin only and $n_{\text{MDR}}=16$ doubly resistant isolates (see Table 1). The remaining statistics, $c_{k,k^{\prime}}$, are computed as:

[TABLE]

Finally, in order to reduce the impact of summary statistics operating on different scales, we compare simulated and observed summary statistics within the kernel $K_{h}(\|s-s_{obs}\|)$ via the $L_{\frac{1}{2}}$ norm

[TABLE]

where $\dim(s)=14$ is the number of summary statistics. Alternative approaches could rescale the statistics via an appropriate covariance matrix (e.g. \shortciteNPLuciani.09.epidemiological,erhardt+s16) or use other norms, however the results in the following Section proved to be robust to more structured comparisons, so we did not pursue this further. In particular the following results were robust to these choices because of the use of a good (non-linear) regression adjustment, which greatly improves the ABC posterior approximation, and which has a larger impact on this approximation than the choice of metric $\|\cdot\|$.

4.2 Parameter specifications and prior distributions

Of the 13 model parameters (see Table 3), eight of these are known well enough for the purposes of our analysis to fix their values. Namely, the parameters $(\delta,\gamma_{0},\gamma_{\mathrm{INH}},\gamma_{\mathrm{RIF}},\gamma_{\mathrm{MDR}},N,\tau,c)^{\top}$ are set to these fixed values. We justify our choices for these values below. The remaining five parameters $\theta=(\beta_{0},\mu,\rho_{\mathrm{INH}},\rho_{\mathrm{RIF}},\rho_{\mathrm{MDR}})^{\top}$ are to be estimated, and require a prior distribution specification.

The rate of death or recovery, $\delta$, is fixed and set to be $\delta=0.52$ per case per year following \citeNDye.01.Will and \citeNCohen.04.Modeling. Similarly, following \citeNDye.01.Will, untreated individuals are detected and treated at rate $\tau=0.5$ per case per year. The rates of recovery due to treatment, $\gamma_{k}$, for resistance profiles $k=0$, INH, RIF and MDR, can be written in terms of the probability of treatment success

[TABLE]

We set the cure rates to be $\gamma_{0}=0.5,\gamma_{\rm INH}=\gamma_{\rm RIF}=0.25,\gamma_{\rm MDR}=0.05$, which, by using $\delta_{r}=0.2$ \shortciteDye.01.Will,Cohen.04.Modeling, corresponds to treatment success probabilities of approximately $p_{\rm 0}=0.69$, $p_{\rm INH}=p_{\rm RIF}=0.58$ and $p_{\rm MDR}=0.44$. These values are within the supported ranges in the literature, namely, $p_{\rm 0}=0.45-0.75$, $p_{\rm INH}=p_{\rm RIF}=0.3-0.6$ and $p_{\rm MDR}=0.05-45$ [\citeauthoryearBlower and ChouBlower and Chou2004]. We chose higher values within these ranges since \citeNBlower.04.Modeling explored a wide range of possibilities in models including epidemiologically pessimistic scenarios.

The fitness cost of drug resistance, $c$, was fixed and set to be $c=0.1$ based on estimates by \shortciteNLuciani.09.epidemiological. To set the total population size $N$ we first observe that because the sample of 100 isolates represents $\sim 1.1\%$ of the population, this implies that the infected population is 9091. We expect that the number of susceptible individuals who are exposed to disease is somewhat higher than this. Accordingly, we assumed that the total size of the population susceptible to tuberculosis is $N=10,000$. Larger total population sizes can be used, at the price of greater computational overheads for generating data under the model.

Previous work estimated rates of resistance acquisition by mutation to be around $0.0025-0.02$ per case per year \shortciteLuciani.09.epidemiological. The rate of mutation of the VNTR loci in M. tuberculosis was estimated to be around $10^{-3}$ per locus per case per year \shortciteReyes.10.Mutation,Aandahl.12.model,Ragheb.13.mutation but lower estimates have also been found \shortciteWirth.08.Origin,Supply.11.mutation. All of these mutation rates are much lower than 1. We treat these mutation rate parameters as probabilities and conservatively set the standard uniform distribution as a wide prior on each parameter. That is, for the acquisition of resistance to isoniazid or rifampicin (or both), we specify priors for the rates of resistance acquisition as $\rho_{\mathrm{INH}},\rho_{\mathrm{RIF}},\rho_{\mathrm{MDR}}\sim U(0,1)$. Similarly, for the mutation rate of the VNTR molecular marker, $\mu$, we use the prior $\mu\sim U(0,1)$.

The transmission parameter for doubly sensitive strains $\beta_{0}$ is given the shifted gamma prior

[TABLE]

where the parameters are chosen such that the resulting prior distribution of the basic reproduction number $R_{0}$ closely resembles the distribution obtained in a numerical analysis of tuberculosis dynamics by \shortciteNBlower.95.intrinsic. Note that the prior on $\beta_{0}$ is shifted in order ensure the realistic condition that $R_{0}>1$. A value of $R_{0}$ lower than unity would lead to extinction of M. tuberculosis.

We reiterate that we interpret the rate parameters as probabilities per time step and handle the parameters so that their values remain in (0,1). This approximation increases in accuracy as the time unit decreases. Here we divide the natural time unit of one year into new units of 1/12 year per time step.

5 Competing models of resistance acquisition

We estimate the rates of acquisition of drug resistance to rifampicin and isoniazid by fitting the model described in Section 3 to the Bolivian data \shortciteMonteserin.13.genotypes with the ABC method described in Section 4. Additionally, by constraining particular resistance-acquisition parameters $\rho_{k}$ to produce meaningful submodels of the full model, we are able to examine two specific biological questions. The relationships between the two submodels and the full model are illustrated in Figure 3. First, we ask whether it is possible for multidrug resistance to evolve directly from doubly sensitive bacteria or whether this direct conversion does not occur (i.e., $\rho_{\rm MDR}=0$: Submodel 1). Second, we ask whether differences between rates of mutation to rifampicin and isoniazid resistance are apparent at the epidemiological scale (i.e., $\rho_{\rm INH}=\rho_{\rm RIF}=\rho_{\mathrm{single}}$: Submodel 2).

Figure 4 illustrates the ABC marginal posterior density estimates of each parameter under the three different sets of model assumptions. Under the full model there is a clear visual difference between the rates of mutation of rifampicin and isoniazid resistance, with the latter occurring at a much higher rate. In contrast, the rate of simultaneous resistance acquisition appears to be higher than that for rifampicin alone. When eliminating the possibility of simultaneous acquisition of multiple drug resistance $\rho_{\mathrm{MDR}}=0$ (Submodel 1), $\rho_{\mathrm{INH}}$ and $\rho_{\mathrm{RIF}}$ both increase, relative to the full model, to compensate for the imposed restriction when fitting to the observed data (Figure 4b). Similarly, when we fix the identity $\rho_{\mathrm{INH}}=\rho_{\mathrm{RIF}}=\rho_{\mathrm{single}}$ (Submodel 2) to impose a single rate of resistance acquisition, the posterior density of this parameter moves to intermediate values compared to the two distinct rates of acquisition estimated under the full model (Figure 4c). The estimated posterior densities for the transmission ($\beta_{0}$) and mutation ($\mu$) parameters are visually similar across all models. ABC marginal posterior means and highest posterior density (HPD) credible intervals for all models are reported in Table 4.

5.1 Can resistance to both drugs be acquired simultaneously?

To determine whether resistance to both drugs can evolve directly from a double sensitive strain within an infection, we compare Submodel 1 ($\rho_{\mathrm{MDR}}=0$) against the full model. Formal standard Bayesian model comparison typically occurs through Bayes factors. In the ABC framework this task is complicated by the need to perform ABC with summary statistics that are informative for the model indicator parameter, in addition to those informative for the model specific parameters. Such summary statistics can not only be difficult to identify, but the resulting composite vector of summary statistics can be high dimensional, which may then produce more inaccurate inference than if each model was analysed independently. See e.g. \shortciteNrobert+cmp11, \shortciteNmarin+prr14 and \shortciteN[this volume]marin+pr17 for a discussion of these issues. A useful alternative is to consider posterior predictive checks or related goodness-of-fit tests (e.g. \shortciteNPthornton+a06,csillery+bgf10,Aandahl.12.model,prangle+bps14).

Figure 5 shows the posterior predictive distribution of the summary statistics $(n_{0},n_{\mathrm{INH}}+n_{\mathrm{RIF}},n_{\mathrm{MDR}})$ described in Section 4.1, for the full model (panel (a)) and Submodel 1 (panel (b)), where a darker intensity indicates higher density. This predictive distribution graphically illustrates each model’s ability to generate the observed summary statistics (78, 8, 16), indicated by the asterisks, which represent the number of individuals in the sample sensitive to both drugs ($n_{0}$), resistant to a single drug ($n_{\mathrm{INH}}+n_{\mathrm{RIF}}$) and resistant to both drugs ($n_{\mathrm{MDR}}$).

The predictive distributions for each model are diffuse, particularly for the full model. This variability is expected given that the sample size is small (100 isolates) and that the evolution of drug resistance from sensitivity is a relatively rare stochastic event. In the case of Submodel 1 (Figure 5 panel (b)) where we impose the condition $\rho_{\mathrm{MDR}}=0$, the density of samples is shifted away from the bottom-right corner which represents double resistance. This pattern is due to the lack of the direct route to multidrug resistance. The observed data (asterisk) is in a region of low posterior predictive density under Submodel 1, and so we conclude that this model is not particularly supported by the data. In contrast, the observed data lie more clearly within a moderately high density region of the posterior predictive under the full model (Figure 5 panel (a)). This analysis therefore suggests that of the two competing hypotheses, it is more likely that resistance to both drugs can be acquired simultaneously ($\rho_{\rm MDR}>0$) than otherwise. Note, however, that this direct route is not the only possible path to double resistance, which can still occur in stages through single resistance.

5.2 Is resistance to both drugs acquired at equal rates?

In order to determine whether the rates of acquisition of resistance to the two drugs are equal ($\rho_{\mathrm{INH}}=\rho_{\mathrm{RIF}}$), we compare Submodel 2 against the full model. Figure 6 depicts the posterior predictive distribution of $(n_{\mathrm{INH}},n_{\mathrm{RIF}})$ under each model – the number of cases resistant only to isoniazid ($n_{\mathrm{INH}}$) and the number of cases resistant only to rifampicin ($n_{\mathrm{RIF}}$) in the sample. The observed values of these summary statistics are $n_{\mathrm{INH}}=8$ for isoniazid and $n_{\mathrm{RIF}}=0$ for rifampicin, illustrated as the asterisk in Figure 6. As Submodel 2 does not favor any drug over the other, the predictive surface is symmetric with respect to the line $n_{\mathrm{INH}}=n_{\mathrm{RIF}}$. The extra flexibility provided by the full model shifts the predictive distribution towards the observed data. While the distribution under the full model comfortably accommodates the empirical point in a high density region, the predictive distribution under Submodel 2 is much more diffuse. This indicates that while the observed data is not unsupported under Submodel 2, it is far more likely to be observed under the full model. As a result, we conclude that the evidence favours the drugs being acquired at different rates; specifically, isoniazid resistance evolves faster than rifampicin resistance.

5.3 The relative contribution of transmission and treatment failure to MDR-TB

In addition to estimating the rates of acquisition of drug resistance and assessing whether rates differ, we may also consider where doubly resistant cases come from. That is, estimation of the relative contribution to multidrug resistant cases of transmission of existing MDR-TB strains compared to treatment failure leading to evolution of multidrug resistance. The posterior predicted samples generated under the full model provide a clear portrait of the relative contribution of the different paths to achieving double resistance (see e.g. \shortciteNPLuciani.09.epidemiological for an additional illustration of this procedure).

Table 5 shows the means, medians and the 95% HPD credible intervals for the predicted proportion of cases of double resistance from each potential source. These proportions are obtained conditionally on there being at least one case of double resistance in the predictive sample. Simulated samples of this nature account for 99.67% of all predictive samples. The predictive distributions of the proportions are highly asymmetric (not shown), making the median a more reliable point estimate than the mean.

In the overwhelming majority of posterior predictive samples, direct transmission was the main source of acquisition of double resistance, followed by conversion in a single step directly from a sensitive profile (from profile [math] to MDR) and conversion in two steps via a state of resistance to a single drug (from profile [math] to INH to MDR, or from [math] to RIF to MDR). This analysis corroborates the finding from Section 5.1 that $\rho_{\mathrm{MDR}}$ is most likely positive, and furthermore that this path is likely to be of even greater importance than conversion in two steps.

6 Conclusions

In this chapter we have estimated epidemiological parameters describing the acquisition of multi-drug resistance in M. tuberculosis from molecular epidemiological data \shortciteMonteserin.13.genotypes using approximate Bayesian computation. The underlying model is intended to capture essential processes that give rise to the data, namely, transmission of the disease, recovery or death, and within-host evolution giving rise to drug resistance and new genotypes at the molecular marker loci. From this analysis we may draw three major biological conclusions about the manner in which drug resistance arises.

First, there is an asymmetry in the acquisition of resistance to isoniazid and rifampicin. Specifically, isoniazid resistance occurs approximately an order of magnitude more frequently than resistance against rifampicin (see Table 4). This asymmetry in rates is consistent with in vitro (that is, through laboratory experiments) microbiological estimates of mutation rates per cell generation which find around 1 to 2 orders of magnitude difference between the two rates \shortciteDavid.70.probability,Ford.13.Mycobacterium.

Second, the analysis supports the occurrence of direct conversion from doubly drug sensitive to doubly resistant (MDR) infections. This may be initially unintuitive because under mutation alone, if mutation occurs at rate $\rho$ per gene per unit time, the rate of appearance of double mutants is $\rho^{2}$, which would be vanishingly small if $\rho$ is low. However, using a mathematical model, \shortciteNColijn.11.Spontaneous argued that direct conversion can occur surprisingly fast because resistant cells are sometimes present at low frequencies in a within-host population even before treatment commences. Our analysis of data at the epidemiological level is consistent with that theoretical result. This direct conversion to double resistance is epidemiologically important as it accelerates the accumulation of resistance, in that resistance evolution does not have to take place sequentially. Once double resistant mutants appear, transmission of these mutants further increases their prevalence in the population.

Third, the overwhelming majority of cases of multidrug resistant tuberculosis come from transmission of already multidrug resistant strains (see Table 5), a finding that is consistent with those of \shortciteNLuciani.09.epidemiological. This large contribution of transmission occurs despite the 10% transmission cost of each resistance which results in a $\sim 20\%$ cost for MDR-TB. This implies that in controlling drug resistance, although there is widespread concern about treatment failure leading to rising resistance, most resistant cases may be due to transmission. Therefore, although it is important to support treatment adherence, public health efforts may benefit from focusing more on preventing disease transmission. That is, control measures that reduce the incidence of new cases are likely to help reduce MDR-TB.

By developing epidemiological models with evolutionary processes we have been able to estimate parameters describing how drug resistance – particularly multidrug resistance – emerges in M. tuberculosis. Although there is existing knowledge of rates of mutation to resistant states in vitro, there is a need to assess the extent to which those rates translate to the epidemiological level. Large scale molecular epidemiological models, such as those presented here, are highly complex and multidimensional, and as such, likelihood-based analyses are not straightforward mathematically or computationally. In such cases, approximate Bayesian computation methods present a practical and viable approach to making statistical inferences, particularly as continually advancing molecular technologies require dynamical models to be extended and refined.

Acknowledgements

GSR is funded by the CAPES Foundation via the Science Without Borders program (BEX 0974/13-7). SAS is supported by the Australia Research Council through the Discovery Project Scheme (DP160102544). MMT is supported by grant DP170101917 from the Australian Research Council. This research includes computations using the Linux computational cluster Katana supported by the Faculty of Science, UNSW Australia.