The cost-effectiveness of oral HIV pre-exposure prophylaxis and early antiretroviral therapy in the presence of drug resistance among men who have sex with men in San Francisco

Background Poor adherence to either antiretroviral treatment (ART) or pre-exposure prophylaxis (PrEP) can promote drug resistance, though this risk is thought to be considerably higher for ART. In the population of men who have sex with men (MSM) in San Francisco, PrEP coverage reached 9.6% in 2014 and has continued to rise. Given the risk of drug resistance and high cost of second-line drugs, the costs and benefits of initiating ART earlier while expanding PrEP coverage remain unclear. Methods We develop an infection–age-structured mathematical model and fit this model to the annual incidence of AIDS cases and deaths directly, and to resistance and demographic data indirectly. We investigate the impact of six various intervention scenarios (low, medium, or high PrEP coverage, with or without earlier ART) over the next 20 years. Results Low (medium, high) PrEP coverage with earlier ART could prevent 22% (42%, 57%) of a projected 44,508 total new infections and 8% (26%, 41%) of a projected 18,426 new drug-resistant infections, and result in a gain of 43,649 (74,048, 103,270) QALYs over 20 years compared to the status quo, at a cost of $4745 ($78,811, $115,320) per QALY gained, respectively. Conclusions High PrEP coverage with earlier ART is expected to provide the greatest benefit but also entail the highest costs among the strategies considered. This strategy is cost-effective for the San Francisco MSM population, even considering the acquisition and transmission of ART-mediated drug resistance. However, without a substantial increase to San Francisco’s annual HIV budget, the most advisable strategy may be initiating ART earlier, while maintaining current strategies of PrEP enrollment. Electronic supplementary material The online version of this article (10.1186/s12916-018-1047-1) contains supplementary material, which is available to authorized users.


Model formulation
We extended our previously developed model [1,2] by considering the use of PrEP among male homosexual population in San Francisco which is decomposed into twelve categories (see Figure 1(a) in the main text and Appendix Figure 1): susceptible individuals without PrEP (S) and with PrEP (P ), untreated individuals infected with drug-sensitive strains at the primary stage (I U s1 ), chronic stage (I U s2 ) and AIDS stage (I U s3 ), untreated drug-resistant cases at the primary stage (I U r1 ), chronic stage (I U r2 ) and AIDS stage (I U r3 ), individuals treated with combination antiretroviral therapy (ART) but did not develop drug resistance (I T s1 ) and those who have entered the AIDS stage (I T s2 ), and ART-treated individuals with drug resistance before the AIDS stage (I T r1 ) and at the AIDS stage (I T r2 ). The variable's subscript identifies whether the infection is drug-sensitive (s) or drug-resistant (r); the superscript specifies whether the individuals are treated with ART (T ) or untreated (U ).
Denote the duration of the primary, chronic, AIDS stage for untreated drug-sensitive individuals as a p , d c , d A , and assume untreated drug-resistant individuals have a longer chronic stage d r (≥ d c ) due to weaker viral replication capacity (lower viral load in the absence of drug pressure) and thus a longer life expectancy [3,4].
We assume that the duration of the primary and AIDS stages did not differ with or without treatment and resistance, as in [1,5] infection and ART initiation is 2 years [6]). The parameter f r is the fraction of treated individuals who develop drug resistance and we assume that all of these drug-resistant cases use second-line drugs. Let t denote time and a denote the infection age. We assume that all of the infected individuals with the same infection age are homogeneous and have the same rates.
Let i U qj (a, t), i T q1 (a, t) and i T q2 (a, t) (where j = 1, 2, 3 and q ∈ {s, r}) be the respective density of infected individuals in I U qj , I T q1 and I T q2 classes at time t and infection age a. It follows that are the number of infected individuals in I U qj , I T q1 and I T q2 classes (j = 1, 2, 3 and q ∈ {s, r}), respectively, at time t ≥ 0. We denote the disease-induced mortality rates in the classes I U s2 , I U r2 , I T s1 , I T r1 and AIDS classes (including I U s3 , I U r3 , I T s2 , I T r2 ) as µ U s , µ U r , µ T s , µ T r and µ A , respectively. The probability that an infected individual in the I U qj , I T q1 and I T q2 class (where j = 1, 2, 3 and q ∈ {s, r}) still stays in the original class at infection age a [1] is given by where δ 1 is the progression rate to the chronic stage and δ U q , δ T q (q ∈ {s, r}) are the progression rates to the AIDS stage for untreated and treated individuals, respectively.
Let F q be the fraction of the untreated population that receives treatment [3]. It is given by Denote the transmission rate at the primary stage, chronic stage and AIDS stage for untreated drug-sensitive and drug-resistant individuals as β p s , β U s , β A s , and β p r , β U r , β A r , respectively. The transmission rate of a treated drug-sensitive (β T s ) or drug-resistant (β T r ) case is the infectivity of an untreated individual (β) multiplied by a constant, i.e., β T s = α s β U s and β T r = α r β U s where α s ≤ α r , i.e., the second-line drug effectiveness 1 − α r is assumed to be lower than the first-line drug effectiveness 1 − α s (see Appendix ) due to lower adherence [7].
Patients starting ART with higher baseline CD4 counts had longer life expectancies [5,[8][9][10][11]. The relationship between prior-treatment CD4+ count and infection age shown in Fig. 1(B) in [12] also suggested that a higher CD4+ count corresponded to an earlier infection stage. Thus, the earlier ART starts, the longer the patient is expected to live and vice versa. Similar to the assumption in [5], we assume that the duration from treatment initiation to death for treated individuals is a linear decreasing function of ART initiation timing a ART as follows (see Appendix Figure 3 in [2]): where L 0 q is the average maximum expectancies for those who are treated immediately after infection (i.e.,a ART = 0) and ξ T q is the slope. Therefore, we have We assumed the PrEP effectiveness against drug-sensitive strains (defined as the reduction in susceptibility to HIV infection upon exposure to an HIV-infected drug-sensitive partner) was ϵ s . We defined PrEP effectiveness to represent the product of PrEP's biomedical efficacy and PrEP adherence. We assumed PrEP effectiveness against resistant strains was ϵ r = δ r ϵ s (δ r is defined as relative effectiveness of PrEP, i.e., the ratio of PrEP effectiveness against drug-resistant strains ϵ r relative to drug-sensitive strains ϵ s ) [13][14][15].
We develop the complete dynamical model as follows where i U q10 (a) ∈ L 1 Here L 1 + is the space of functions that are nonnegative and Lebesgue integrable over the specified interval.
The number of persons living with HIV/AIDS at any time t is given by The total population size at any time t is given by The PrEP coverage among susceptible MSM population among any time t is given by Notice that we keep track of the numbers of newly diagnosed AIDS cases and AIDS deaths at time t using the following equations: , There are two ways to model the development of drug resistance. The first way is to model the progression The second is to model the consequence of treatment, assuming that if resistance is acquired then it happens fast enough such that amongst those who develop resistance, it is not important to consider that they move through a treated sensitive state. In this formulation, a fraction of treated patients will develop drug resistance after ART treatment and a fraction won't. In this paper, we adopted the second choice. We did this because it allows us to assume that a fraction (f r ) of treated individuals will acquire drug resistance, as shown in Figure 1(a) in the main text and Appendix Figure   1. With the first method, drug resistance is developed at a certain rate, and the proportion whoever develop drug resistance would also depend on other rates of the model, making it more difficult to match the model to available data.
A second advantage of using the second method is that we can track the timing of ART administration easily for both treated drug-sensitive and drug-resistant individuals using our infection-age-structured model. We used the same modeling approach in our earlier work [2] and adapt the parameter estimates from that paper here. If and the cumulative number of new drug-resistant infections over T years is calculated by where i U q1 (0, t), q ∈ {s, r} is the new infections at time t shown in Eq. (6).

Prevalence
The HIV prevalence at any time t is given by where I total (t) and N (t) are given by Eqs. (7) and (8) respectively.

The fraction of new infections that are drug resistant
The prevalence of transmitted drug resistance (TDR) among newly infected individuals (the fraction of new infections that are drug resistant) at any time t is given by where i U q1 (0, t), q ∈ {s, r} is the new infections at time t shown in Eq. (6).

QALYs
Denote the quality of life of classes S, P, I U qj , I T q1 , I T q2 , q ∈ {s, r}; j = 1, 2, 3 as Q S , Q P , Q U qj , Q T q1 , Q T q2 , and the life discount rate asr (see Appendix Table 1), then we have the total health benefits for the entire population measured in discounted quality-adjusted life years (QALYs) during the intervention duration T (see [16,17]): where , q ∈ {s, r}; j = 1, 2, 3 are given by Eq. (1). Here we assume that quality of life for drug-resistant individuals decreases by 5% relative to drug-sensitive individuals at the same stage [18,19], i.e.,

Costs
Denote annual costs of PrEP as c P , and annual HIV-related health care costs (exclude ART costs) for classes Here, we assume the HIV-related health care costs for drug-sensitive individuals are the same as these for drug-resistant individuals, i.e., c U sj = c U rj ; j = 1, 2, 3, c T s1 = c T r1 , and c T s2 = c T r2 . Let annual costs of ART (first-line drugs) for treated drug-sensitive individuals be c ART s and annual costs of second-line drugs for treated drug-resistant cases be c ART r (here we assume all drug-resistant individuals use second-line drugs). Assume ART was immediately following diagnosis so that the test rate ϕ can be obtained as ϕ = 1/a ART . For uninfected individuals, the total costs of testing (antibody test) and counseling are denoted as c neg . For infected individuals, the total costs of testing (antibody test) and counseling are denoted as c pos , and the costs of diagnosis is denoted as c diag . Denote the costs of genotype resistance test as c GT . See Appendix Table 1 for these costs. The cost discount rate is assumed as r. Then we calculate the total discounted costs for the entire population as the sum of annual health care costs for all individuals over the intervention's duration T as follows (see [16,17]): where ) dt, ) dt,  Table 2 for details.

ICERs
We calculated the incremental cost-effectiveness ratio (ICER) of each intervention strategy, relative to the status quo as follows (see [16,17,20,21]): In some cases, we also calculated the ICER of one intervention strategy A relative to another strategy B (for example, combination of high PrEP and earlier ART versus high PrEP alone).
According to the WHO standards [22][23][24], if ICER<1 per capita GDP, the strategy is regarded as very costeffective; if 1 per capita GDP<ICER<3 per capita GDP, the strategy is regarded as cost-effective; if ICER>3 per capita GDP, the strategy is regarded as not cost-effective.

Parameter estimation
We obtained data of the annual newly diagnosed AIDS cases and AIDS deaths from 1980 to 2014 in MSM population in San Francisco from the Department of Public Health HIV Epidemiology Section. Using the maximum likelihood estimation, we fit the model to the data between 1980 and 1995 to estimate the prior-treatment parameters (see more details in [2]): the recruitment rate of susceptible MSM is b = 4000 (95%CI:2295-5705) per year, the disease-induced death rate at the chronic stage for untreated drug-sensitive individuals is µ U s = 0.28 (95%CI:0.18-0.39) per year, and the transmission rate in this stage is β U s = 0.62 (95%CI:0.57-0.68) per year. The estimation of the transmission rate is the same as the value of 0.62 shown in [25] for MSM in San Francisco.
The values of other parameters are given as follows. The initial MSM population size is chosen as 69122 [26,27].
The expected duration of a sexual career in San Francisco is assumed to be about 47 years (age 18-65) as in [28].
Thus, we have the removal rate m = 1/47 = 0.021 per year. For untreated drug-sensitive individuals, we choose the duration of the primary stage as a p = 1.7 months [29], the duration of the chronic stage as d c = 7.5 years [30], and the duration of the AIDS stage as d A = 1.2 years (12 months to 20 months in [31]). Thus, we obtain the rate of progression to the asymptomatic stage is δ 1 = 1/a p = 1/(1.7/12) = 7.06 per year. Similarly, we have the rate of progression to the AIDS stage is δ U s = 1/d c = 1/7.5 = 0.13 per year and the disease-induced death rate in the AIDS stage is µ A = 1/d A = 1/1.2 = 0.83 per year. We assume the chronic stage d r for untreated drug-resistant cases is 25% longer than untreated drug-sensitive cases, i.e., d r = 1.25d c = 9.38 years, then the progression rate to the AIDS stage is δ U r = 1/d r = 1/9.38 = 0.11 per year. The transmission rates in the primary stage and AIDS stage are assumed to be 5.3 and 7 times more infectious than during chronic infection, respectively, i.e., β p s = 5.3β U s , β p r = 5.3β U r [29] and β A s = 7β U s , β A r = 7β U r [32] (see Appendix We used the data from 1996 to 2006 (because ART was widely available after 1995 [33]) to estimate treatmentrelated parameters. The treatment rate is estimated as η = 0.38 (95%CI:0.08-0.81) per year, i.e., the average time from infection to ART initiation is a ART = 2.8 years according to the relationship η = 1/(a ART − a p ).
The fraction of treated gay men in San Francisco is calculated as F s = 46.4% (95%CI:15.3%-64.9%) based on Eq. (3), which is close to the fraction in [3] where about 50% of HIV-infected MSM take ART. The diseaseinduced death rate in the post-treatment chronic stage is estimated as µ T s = 0.05 per year (95%CI:0.01-0.27) for drug-sensitive individuals and µ T r = 1.75µ T s = 0.088 per year for drug-resistant individuals [34]. In this fitting process, we chose a bigger recruitment rate b = 5600 per year to yield simulated prevalence, total infected individuals and population size simultaneously consistent with the prevalence data (Figure 2c in the main text), persons living with HIV/AIDS data and total MSM population size data respectively as closely as possible, which we did not fit directly (see [2] for more details). We also chose the relative transmissibility for treated drug-resistant individuals (β T r = 0.2β U s , i.e, the baseline second-line drug effectiveness was estimated as 80%) to match the prevalence data of transmitted drug resistance (Appendix Figure 3, see [2] for more details) under the assumption that the transmission rate for untreated drug-resistant individuals β U r was the average of that for treated drug-resistant β T r and untreated drug-sensitive individuals β U s based on their relationship β U s > β U r > β T r [3], i.e., β U r = (β T r + β U s )/2 = 0.6β U s . The fraction of treated cases that are drug resistant is chosen as f r = 25% (33% of MSM are virally unsuppressed [28] and of which 76% have drug resistance [35]).
The result of the HPTN 052 clinical trial [36,37] showed that treatment led to 96% reduction in infectivity.
Thus, the transmission rate in treated individuals without drug resistance is only 4% as transmissible as HIV positives not receiving ART (β T s = 0.04β U s ). After 2006, San Francisco had name-based HIV reporting and incorporated monitoring initial primary care visit into standard HIV public health investigation for newly diagnosed cases which improved the treatment rate and shortened the time to entry into HIV medical care [38]. So we used the data from 2006 to 2012 to estimate the growing treatment rate η = 0.7 per year and earlier ART initiating timing a ART = 1.6 years. All the other parameters are fixed in Appendix Table 1.
In 2012, the PrEP was approved by FDA and PrEP use has increased to 9.6% in 2014 in San Francisco [39]. We calibrated the model by choosing the PrEP receiving rate as ξ = 0.06 per year to match the 9.6% PrEP coverage in 2014 and the coverage will reach 25% in 2023 (low coverage) at this PrEP receiving rate. If we consider the cases that PrEP coverage will reach 50% (medium coverage) or 80% (high coverage) after 5 years (in 2023), then the rate that susceptible receive PrEP should increase to ξ = 0.20 or ξ = 0.68 per year respectively.
We assumed the PrEP effectiveness against drug-sensitive strains was ϵ s = 53% in the base case based on the meta-analysis results in [40] and against the resistant strains was ϵ r = 50% × ϵ s = 26.5% [13][14][15]. The rate discontinuing PrEP was chosen as w = 8% per year based on the data in [41] that approximately 8% of participants gave up PrEP due to side effect concerns and no perceived HIV risk in a one-year project in San Francisco. We considered how the above PrEP scenarios (low, medium, high PrEP coverage) interacted with the implementation of new even earlier ART guidelines (initiation at 1 year post infection averagely), in contrast with previous ART guidelines (initiation at 1.6 years post infection averagely), will affect the cost-effectiveness of PrEP over the next 20 years (from 2018 to 2038).
Next, we derive the relationship between extended life expectancies and infection age (more detail was shown in [2]). It can be seen from [10,11] that suppressed patients (HIV-1 RNA ≤ 400 copies/ml) who had CD4+ count < 200 or ≥ 350 cells/mm 3 at ART start can live mean 30 or 45 years after treatment, respectively. In addition, Fig. 1(B) in [12] shows that CD4+ count decreases with time since infection (infection age). Specifically, CD4+ count < 200 and ≥ 350 cells/mm 3 corresponds to the infection age of 7-9 years and 0-6 years, respectively. We Notice that 76% of treatmentfailed patients have resistance to one or more antiretroviral drugs [35]. Therefore, we assume that the above relationships between the extended life expectancy and ART initiating timing for unsuppressed and suppressed patients still hold for treated individuals who do or do not develop drug resistance. See [2] for more details. All the other parameters can be found in Appendix Table 1.

Supplementary Tables
Appendix Table 1 shows the parameters used in simulation. Appendix Table 2 gives more detail on the costs of expanding PrEP coverage and earlier ART, incremental to status quo (base case). Appendix