The Emergence of HIV Transmitted Resistance in Botswana: “When Will the WHO Detection Threshold Be Exceeded?”

Background The Botswana antiretroviral program began in 2002 and currently treats 42,000 patients, with a goal of treating 85,000 by 2009. The World Health Organization (WHO) has begun to implement a surveillance system for detecting transmitted resistance that exceeds a threshold of 5%. However, the WHO has not determined when this threshold will be reached. Here we model the Botswana government's treatment plan and predict, to 2009, the likely stochastic evolution of transmitted resistance. Methods We developed a model of the stochastic evolution of drug-resistant strains and formulated a birth-death Master equation. We analyzed this equation to obtain an analytical solution of the probabilistic evolutionary trajectory for transmitted resistance, and used treatment and demographic data from Botswana. We determined the temporal dynamics of transmitted resistance as a function of: (i) the transmissibility (i.e., fitness) of the drug-resistant strains that may evolve and (ii) the rate of acquired resistance. Results Transmitted resistance in Botswana will be unlikely to exceed the WHO's threshold by 2009 even if the rate of acquired resistance is high and the strains that evolve are half as fit as the wild-type strains. However, we also found that transmission of drug-resistant strains in Botswana could increase to ∼15% by 2009 if the drug-resistant strains that evolve are as fit as the wild-type strains. Conclusions Transmitted resistance will only be detected by the WHO (by 2009) if the strains that evolve are extremely fit and acquired resistance is high. Initially after a treatment program is begun a threshold lower than 5% should be used; and we advise that predictions should be made before setting a threshold. Our results indicate that it may be several years before the WHO's surveillance system is likely to detect transmitted resistance in other resource-poor countries that have significantly less ambitious treatment programs than Botswana.


The mean-field dynamics
In Fig.1 we show the flow diagram of the compartmental Markov chain model that we used to describe the dynamics of the emergence of acquired and transmitted drug-resistant HIV in a population. The model identifies four population groups: sexually-active susceptible adults X, treatment-naive adults infected with wild-type HIV strains Y U S , treated adults infected with wild-type HIV strains Y T S , and adults infected with drug-resistant HIV strains Y R . The drug-resistant group could be subdivided into treated and treatmentnaive classes. In the limit of large population numbers present in all four compartments, the dynamics (1.7)

The parameters: definition and estimation
The parameters that appear in Eqs. (1.1)-(1.7) are defined as follows: π is the inflow rate of newly sexually active adults; µ −1 is the average length of time an adult acquires new sex partners; c is the average number of new sex partners a susceptible individual acquires per year; σ s is the per capita treatment rate per year; g s is the yearly proportion of patients suspending treatment and r −1 is the average time required for a patient to develop drug resistance when on treatment. The disease progression rate to AIDS for the wildtype treatment-naive, wild-type treated and drug-resistant populations are respectively described by the parameters v U S , v T S and v R ; the transmissibility (i.e., fitness) coefficients for the wild-type treatment-naive, wild-type treated and drug-resistant populations are respectively described by the parameters β U S , β T S and β R .
We assume that the expected time that an adult acquires new sex partners in Botswana is on average 34 years. This is exactly the difference in years between the upper (49 years) and the lower (15 years) adult age group that we consider. The inflow of at-risk susceptible adults, π, is chosen and calculated such that in the short term the total population is stationary in the absence of HIV. On the assumption that the population growth can be ignored over the short term, we obtained that π is 49,300 per year. The contact rate was found to be 1.76 per year on the assumption that HIV prevalence is constant over the short term.
In a similar way the contact parameter, c, is chosen and calculated such that in the absence of treatment (g s = 0 and σ s = 0) Y U S is stationary in the short term. We calculated the per capita treatment rate per year (σ s ) from the projected number of patients that the Botswana ART program plans to treat over the next four years. The program plans to treat 85,000 patients by 2009 [1] (see Fig. 2). The treatment rate was found by assuming that it remains constant such that by the year 2009, Botswana will have 85,000 patients that have received ART. This procedure gave a constant per capita treatment rate of 0.050 per year. For Botswana, we assumed that µ −1 is 34 years (from 15 to 49 years old). We used the values of 10, 18 and 12 years respectively for (v U S ) −1 , (v T S ) −1 and (v R ) −1 representing the average progression time to AIDS [2,3]. Thus, individuals that are on treatment remain sexually active for longer than individuals that are not on therapy, on average, for an additional 8 years. Since our predictions are valid for the short term to 2009, our results are robust to the values chosen for the average progression time to AIDS. We assumed that the yearly proportion of treated cases suspending ART, g S , in Botswana is 0.1 [4] and that, on average, drug resistance would develop in 3 to 5 years [5,6,7,8,9]  β T S ) were calculated to be 0.12 and 0.04 respectively [10,11,12]. Our values for β U S and β T S are consistent with those used by Wilson et al. [13] for South Africa for a similar drug regiment. We assume that β U S lies in the range 0.1 to 0.15 [11]. We then calculate β T S from β U S using the following expression where α T is the reduction in HIV transmissibility (i.e., fitness) due to treatment-induced viral load reduction and p S is the fraction of treated individuals infected that achieve viral suppression on ART. Following Wilson et al., the calculation of the drug-resistant transmissibility is obtained using the following expression where α R is a proportionality constant. Relatively little is known about the fitness of resistant strains of HIV in vivo. In our paper we assumed α R to be 0.25, 0.5 and 1.0 in different scenarios.
An alternative way of obtaining the values for transmissibility is to use the relationship that each logarithmic increase in viral load (w) is associated with an increase in the risk of transmission by a factor of 2.45 [12]. Explicitly, this relationship is formalized as where β(w) and β(ν) are transmissibilities as a function of viral loads w and ν.

The initial conditions
Our model is integrated starting from the initial conditions that no one is on treatment (Y T S (0) = 0). Therefore HIV drug resistance has not developed (Y R (0) = 0). We assume that the HIV-infected wild-type population is at equilibrium with the susceptible population for our initial conditions. These equilibrium population values are given by: . (1.12) 2 Stochastic formulation for the drug-resistant population dynamics

Stochastic vs. Deterministic
The ODEs formulated in section 1.1 can be numerically integrated to obtain the mean-field population evolution. As long as the populations present in each compartment are large, relative stochastic fluctuations are negligibly small. From the discussed initial conditions (Section 1.3), we are assuming that both the treated wild-type and the drug-resistant population begin their dynamics from zero. This means that relative stochastic fluctuations are important in their initial short-term dynamics. However, for the chosen parameters the initial inflow of individuals into the treated wild-type population is much larger than the inflow into the drug-resistant population. Thus, the time for which relative stochastic fluctuations are important is much greater for the drug-resistant population than for the treated wild-type population. Hence, the evolution of the wild-type population can be approximated by its mean-field dynamics at a much earlier time than the evolution of the drug-resistant population. For this reason, the forecast of the short-term dynamics of the drug-resistant population benefits from a stochastic analysis. This analysis provides not only the mean field evolution of the population numbers, but also the evolution of its variability. Our main interest is in the variability present in the prevalence and incidence patterns of the drug-resistant population. We therefore proceed by treating the drug-resistant population stochastically.
A description of the stochastic dynamics is often obtained via computational integration by means of Monte Carlo methods [18]. This would involve identifying the stochastic processes affecting the drug-resistant population dynamics and simulating the evolution discretely starting from the initial conditions. Then, in order to obtain the variability in the dynamics one would have to average over many independent ensemble realizations of this evolution. An alternative method is to try to solve the stochastic dynamics analytically by formulating and solving a Master equation [15,16]. This is an attractive method for proceeding as it produces exact equations to predict the variability in the dynamics. In problems in mathematical epidemiology, this can not always be done as transition rates present in the Master equation formulation are nonlinear functions of the population size. By assuming that the stochastic fluctuations have a predefined distribution (e.g., Gaussian or Log-Normal) this problem can be overcome to give an approximate solution via the method of moment closure [19]. However, in some problems the nonlinear transition rates can be approximated to linear transition rates that remain valid for the time interval of interest. Our stochastic formulation makes use of this latter approximation and we proceed in obtaining expressions for the variability of the dynamics given by time dependent functions of the variance and skewness of the process.

The Master equation
Eq. (1.4) can be used to formulate a birth-death-immigration Master equation [15,16] for the drug-resistant population: where P k (t) denotes the probability that Y R (t) = k. The set of equations described by (2.13) can be used to derive other equations describing the evolution of the moments of the probability distribution of Y R (t).
By the standard approach, we multiply Eq.

∂K(θ, t)/∂t
14) The solution to the master equation allows us to extract cumulants of Y R (t): Y R (t) n = ∂ n K(θ, t)/∂θ n | θ=0 ; (2.15) the mean ( Y R (t) ), variance ( Y R (t) 2 )and the third central moment ( Y R (t) 3 ) are obtained for n = 1, n = 2 and n = 3 respectively. The skewness can be obtained from the third central moment The transition rate β R cX(t)/N (t) is not constant and depends on the evolution of Y R (t). Therefore, as in most epidemic models this makes this transition rate in the master equation a nonlinear function of the population size. To proceed we note that our interest is in the early stages of the dynamics when the susceptible population and the HIV-infected wild-type population have only marginally changed from their initial equilibrium values. Therefore, for the chosen parameters, the approximation that X(t)/N (t) stays constant at its equilibrium value is acceptable during the early stage of the evolution.
This approximation is expressed as X(t)/N (t) ∼ (X) * /[(X) * + (Y U S ) * ] = (X/N ) * 1 . Making use of this approximation allows us to solve for Y T S (t) exactly: 1). (2.16) where α = −(g s + r + v T s + µ) and I = rσ S (Y U S ) * /α. Making use of this solution and using X(t)/N (t) ∼ (X/N ) * = (v U S + µ)/(cβ U S ) we can solve Eq. (2.14) analytically. By using the method of characteristics [17] we obtain three equations: where β = β R c(X/N ) * and γ = (v R +µ). Solving for θ and letting our initial constant characteristic f = θ(0) we obtain which can be re-written in the form Substituting this expression in the equation forK (in (2.17)) gives us This integral can be evaluated by appropriate substitutions. Finally, we use our initial condition of P k (0) = δ k,YR(0) which gives K(θ, 0) = f Y R (0) and we obtain our solution as Here, where F [a, b, c, d] is the hypergeometric function.

Solution for the Mean, Variance and Skewness
By using Eq. (2.15), Eq. (2.21) and our initial condition of no drug-resistant cases (Y R (0) = 0) we obtain all the cumulants of Y R (t). The mean is given by our second central moment (variance) is given by and our third central moment is given by

Solution for the Cumulative Incidence
From the solution of the Master equation we can also extract the cumulant generating function for the cumulative incidence of transmitted drug resistance (T (θ, t)). To obtain this we need to solve the following partial differential equation that relates the rate of change in transmitted incidence to the prevalence of drug resistance in the population: Therefore, the generating function for the cumulative transmitted drug resistance incidence between times t 1 and t 2 is given by: (2.28) By taking the first, second and third derivatives of (2.28) with respect to θ and evaluating these at θ = 0 we obtain expressions for the dynamics of the mean, variance and skewness of the cumulative transmitted drug-resistant population in the interval t 2 − t 1 . We use these expressions in our main article to produce the probabilistic forecasts in the growth of the newly infected transmitted drug-resistant population for different parameter sets. In our manuscript we used time intervals t n+1 − t n each equal to three months.

Solution of the Master Equation versus the numerical integration of the ODEs
We mentioned in Sec. 2.1 that in our stochastic formulation we make use of an approximation whereby the nonlinear transition rates describing the infectious process is replaced by a linear transition rate. This approximation will stay valid initially and only for a short period. where S stochastic represents the mean-field solution of the dynamic variable S obtained from the master equation and S numerical represents the corresponding numerical solution. Fig. 3 shows a plot for Er[Y R (t)] (shown in blue) and Er[Y R (t)/N (t)] (shown in red). Therefore, we compute the percentage difference of both the mean-field solution for Y R (t) and Y R (t)/N (t) to that obtained by numerical integration of the ODEs. The numerical integration was done using