Impact of scaling up dolutegravir on antiretroviral resistance in South Africa: A modeling study

Background Rising resistance of HIV-1 to non-nucleoside reverse transcriptase inhibitors (NNRTIs) threatens the success of the global scale-up of antiretroviral therapy (ART). The switch to WHO-recommended dolutegravir (DTG)-based regimens could reduce this threat due to DTG’s high genetic barrier to resistance. We used mathematical modeling to predict the impact of the scale-up of DTG-based ART on NNRTI pretreatment drug resistance (PDR) in South Africa, 2020 to 2040. Methods and findings We adapted the Modeling Antiretroviral drug Resistance In South Africa (MARISA) model, an epidemiological model of the transmission of NNRTI resistance in South Africa. We modeled the introduction of DTG in 2020 under 2 scenarios: DTG as first-line regimen for ART initiators, or DTG for all patients, including patients on suppressive NNRTI-based ART. Given the safety concerns related to DTG during pregnancy, we assessed the impact of prescribing DTG to all men and in addition to (1) women beyond reproductive age; (2) women beyond reproductive age or using contraception; and (3) all women. The model projections show that, compared to the continuation of NNRTI-based ART, introducing DTG would lead to a reduction in NNRTI PDR in all scenarios if ART initiators are started on a DTG-based regimen, and those on NNRTI-based regimens are rapidly switched to DTG. NNRTI PDR would continue to increase if DTG-based ART was restricted to men. When given to all men and women, DTG-based ART could reduce the level of NNRTI PDR from 52.4% (without DTG) to 10.4% (with universal DTG) in 2040. If only men and women beyond reproductive age or on contraception are started on or switched to DTG-based ART, NNRTI PDR would reach 25.9% in 2040. Limitations include substantial uncertainty due to the long-term predictions and the current scarcity of knowledge about DTG efficacy in South Africa. Conclusions Our model shows the potential benefit of scaling up DTG-based regimens for halting the rise of NNRTI resistance. Starting or switching all men and women to DTG would lead to a sustained decline in resistance levels, whereas using DTG-based ART in all men, or in men and women beyond childbearing age, would only slow down the increase in levels of NNRTI PDR.

. It models the continuum of care -including NNRTI-based rst-line and PI-based second-line regimens -, the disease progression, acquisition and transmission of NNRTIresistance, and its impact on the ecacy of NNRTI-based regimen. The model was calibrated using dierent sources of data: 1) cohort data about more than 30,000 people living with HIV from IeDEA collaboration [1], 2) data from literature, and 3) general HIV estimates at the country scale produced by the Thembisa model. The Thembisa model is a compartmental model providing UNAIDS with estimates on the South African HIV epidemic [2]. We adapted this model to investigate the impact of the introduction of DTG-based regimens in South Africa from 2020. The changes include 1) incorporating DTG-based regimen into the continuum of care, 2) distinguishing between DTG-eligible and -ineligible individuals, and 3) adding a NRTI-resistance dimension.

Adapted MARISA model
The adapted MARISA model is split in 5 dimensions: 1) care stages (15 levels), 2) disease progression, characterised by the CD4 counts (4 levels), 3) sex (2 levels), 4) NNRTI resistance (2 levels) and 5) NRTI resistance (2 levels). The rst dimension of the model accounts for the whole continuum of care (see Fig A). The rst three compartments model respectively HIV-infection of susceptible individuals and diagnosis (with a distinction between DTG-eligible and -ineligible women). We then considered the three dierent regimens -NNRTI-based, PI-based and DTG-based -, again with a distinction by DTG-eligibility for individuals on a NNRTI-based regimen. For each of the three regimens, three compartments are used to model treatment initiation ("Treat init.") with subsequent virological suppression ("Supp") or failure ("Fail"). "Treat init." compartments represent individuals who initiated treatment less than 3 months ago. Before 2020, all individuals receive a NNRTI-based rst-line regimen and switch to the second-line PI-based regimen in case of prolonged failure. From 2020, the DTG-based regimen is used as a rst-line regimen for all DTG-eligible individuals. From this time, DTG-eligible individuals who are currently on NNRTI-based regimen can transition to DTG-based regimen. PI-based regimen is still used as a second-line regimen, either for DTG-ineligible patients failing NNRTI-based ART, or for patients failing DTG-based ART. The second dimension splits individuals in 4 classes according to CD4 counts: 1) CD4 > 500 cells/µL, 2) 350 < CD4 < 500 cells/µL, 3) 200 < CD4 < 350 cells/µL and 4) CD4 < 200 cells/µL. The third dimension makes the distinction between male and female. The fourth and fth dimensions respectively model NNRTI-and NRTI-resistance. They each have two layers that distinguish between NNRTI-/NRTI-susceptible and -resistant individuals. We used the following indices to indicate a layer of a dimension: j for the second dimension (j = 1, 2, 3, 4), k for the third dimension (k = 1, 2), l for the fourth dimension (l = 1, 2) and m for the fth dimension (m = 1, 2).  Rates related to disease progression ν CD4 andν CD4 as well as rates related to continuum of care γ, which respectively model transition from one to another CD4 class and transition from one to another care stage, were estimated using observational cohort data from IeDEA-SA collaboration. Survival analyses were performed using information of more than 30'000 patients from South Africa. Mean estimates and 95% condence intervals (95%CI) are reported in Table A and B. Supplementary Material  Diagnosis rates depend on sex and CD4 class and treatment initiation rates depend on CD4 class. They have been described in details in [4], S1 File, Section 1.3. In the adapted MARISA model, we assumed that diagnosis rates are constant from 2016, while the treatment initiation rate has been adapted in order to model the impact of the Treat-All policy. We increased treatment initiation rates for the rst three CD4 count classes from 2017 to 2022 in order to have identical rates from 2022, irrespective of CD4 counts (see Fig B). To ensure a proportion p 1 of DTG-eligible women, two diagnosis rates are used γ k,elig I→D := p 1 · γ k I→D , γ k,inel I→D := (1 − p 1 ) · γ k I→D , in order to distribute women into the two DTG-eligibility classes. We rescaled the switching rates from unsuppressed NNRTI-based regimen to PI-based regimen γ k,elig/inel F1→T2 in order to reect PI-coverage in South Africa (∼ 4% in 2016 according to [5]). The new rates γ k,elig/inel F1→T2 can be found in Table B. 2.3 Resistance rates Two rates model the ow between the two NNRTI-resistance layers: the reversion rate σ rev and the rate of acquiring NNRTI-resistance σ N N RT I res . Reversion to wild-type occurs when no more drug pressure is exerted, i.e. in the "Infected" and "Diagnosed" compartments. An individual can acquire NNRTI-resistance when failing rst-line regimen. Both parameters σ rev and σ N N RT I res were collected from literature and can be found in Table  C. Supplementary Material Table C: Parameters collected from literature. As mortality estimates in the fourth CD4 class vary according to the proportion of people with CD4 < 50 cells/µL, lower and upper bounds are given (see [4] S1 File Section 1.2 for more details). The mortality risk µ j X in CD4 class j (i = 1, . . . , 4) and care stage X (X = I, D, T 1 , . . .) is given by: µ j X = µ 0 ·μ j X . For the sake of simplicity and in view of the scarcity of evidence on the impact of NRTI-resistance on DTGbased regimen, the dimension modelling NRTI-resistance has only two layers that distinguish between NRTIresistant and NRTI-susceptible individuals. NRTI resistance is dened as having both the K65R and the M184V mutations, which confers high level of resistance to tenofovir (TDF) and lamivudine/emtricitabine (3TC/FTC), respectively. In view of the low level of NRTI pre-treatment drug resistance (PDR) [12,23], we assume that NRTI resistance is not transmitted. The rate σ N RT I res models the process of acquiring NRTI-resistance, which occurs when individuals are failing rst-line NNRTI-based regimen. We calibrated σ N RT I res using results from a meta-analysis that estimates the prevalence of NRTI resistance mutation after 3 years on a failing NNRTI-based rst-line regimen [11]. This meta-analysis found that 75% of them had the K65R mutation and 73% the M184V mutation. Assuming no association between the two mutations as suggested by [24], we inferred σ N RT I res so that 54.8% (i.e. 75% · 73%) of individuals failing NNRTI-based regimen were resistant to NRTI after 3 years of ART. We found σ N RT I res = 1/40 months −1 (see Table C).

Impact of NNRTI-resistance on NNRTI
Unlike the previous MARISA model, the adapted MARISA used two parameters α 1 and α 2 to model the impact of NNRTI resistance on NNRTI treatment response. Both parameters increase the previously estimated rates of failure γ T1→F1 , γ S1→F1 and decrease the suppression rates γ T1→S1 and γ F1→S1 for NNRTI-resistant individuals, but at dierent treatment stages. While α 1 represents the impact of NNRTI resistance among individuals having just started treatment (less than 3 months), α 2 models this impact at the later stage of treatment. In order that the MARISA model achieves the same suppression level with these modied rates as estimated from IeDEA-SA cohort data, we used a third scaling parameter α 4 which increases the overall suppression rates and decreases the failing rates. The dierent failing and suppression rates according to CD4 class j and NNRTI-resistance status l are given in Eq 1-8. The rates γ T1→F1 , γ T1→S1 , γ S1→F1 and γ F1→S1 represent the overall suppression and failure rate for NNRTI-based ART, as estimated with IeDEA cohort data (see Tables A and B).
The three parameters α 1 , α 2 and α 3 were simultaneously calibrated using two dierent kinds of data. To identify α 1 and α 2 , we used estimates from two studies that compared level of NNRTI failure between NNRTIsusceptible and NNRTI-resistant individuals. Both studies reported a hazard ratio (HR) of ART failure between NNRTI-resistant and NNRTI-susceptible people of 3.13. To identify α 4 , we used the overall suppression level of 88% for NNRTI-based regimen, as estimated from IeDEA cohort data. The values of the three parameter estimates are given in Table C. The higher estimated value of α 2 compared with α 1 (α 1 = 1.97, α 2 = 3.24) reects the long-term impact of NNRTI-resistance on the NNRTI-treatment response.

DTG-ecacy and impact of NRTI-resistance on DTG
In this updated MARISA model, we also model the potential impact of NRTI-resistance on DTG-based regimen. We used the same suppression and failure rates for DTG-based regimen as for NNRTI-based one, but replace the scaling factor α 4 by α 5 , to take into account the dierence in treatment ecacy between NNRTI and DTG. The scaling factor α 5 was calibrated to reect results of the NAMSAL study [18], which observed a crude odds ratio (OR) of failure of 1.46 between NNRTI-and DTG-based regimens, after 48 weeks of treatment. To do so, we tted the OR calculated by the MARISA model to the OR observed in NAMSAL studies, taking into account the dierent baseline characteristics (distribution of CD4 counts and level of baseline NNRTI-resistance) of the NNRTI-and DTG-groups. After these adjustments, we found an OR of 1.02 between the two groups, assuming that they are both susceptible to their respective ART regimen (i.e. no NNRTI-resistance). This decrease in OR after adjustment is due to the fact that the NNRTI-group in the NAMSAL had lower baseline CD4 counts and that part of them had baseline NNRTI-resistance. Other ecacies of DTG-based regimens, corresponding to ORs of 2 and 5, were investigated in the sensitivity analysis (see Section 5). As a simplifying assumption, all individuals that transitions to DTG-based regimen are considered to have received a NNRTI-drug combined with TDF and 3TC/FTC and to keep this NRTI-backbones combination after transitioning to DTG-based regimen. This assumption is motivated by the expected reluctance of clinicians to prescribe zidovudine (AZT) for TDF-experienced individuals transitioning to DTG due to its side eects. In the case where NRTI backbones would be adapted when transitioning to DTG, the model might overestimate the impact of NRTI resistance on DTG-based regimen. We applied the same approach to model the impact of NRTI resistance on DTG-based regimen as we did to model the impact of NNRTI resistance on NNRTI-based regimen. The suppression and failure rates for NRTI-resistant individual starting a DTG-based regimen are respectively divided and multiplied by a factor α 3 . The dierent failing and suppression rates according to CD4 class j and NRTI-resistance status m are given in Eq 9-16.

Other parameters: HIV transmission and mortality
The MARISA model accounts for both heterosexual and homosexual HIV-transmission, with dierent risks of transmission per intercourse. We also assumed that undiagnosed individuals have a more risky behaviour. Parameters related to HIV-transmission were either collected from literature (see Table C) or estimated using results from Thembisa model (see Table D). We assumed that mortality depends on both the CD4 counts and the treatment stage. Relative mortality estimates were collected from literature (see Table C) and a scaling parameter, representing the mortality risk among suppressed invidual with CD4>500 copies/ml, was tted to HIV-mortality estimate provided by the Thembisa model. More information about HIV-transmission and mortality can be found in S1 Text of [4]. 1. DTG only used in rst-line regimen of ART-initiators and, as second-line, in patients failing NNRTI-based ART (γ D→T 3 and γ F 1→T 3 in Fig A), 2. DTG used as initial rst-line regimen (for ART-initiators), with all patients on NNRTI-based regimens being switched to a DTG-based regimen (all the red arrows).
Within these 2 scenarios, 4 sub-scenarios investigated the impact of dierent percentages p 1 of DTG-prescription for women: a) no women (0%), b) women outside reproductive age (17.5%), c) women outside reproductive age or using contraception (63%), and d) all women (100%). Percentage in b) is calculated with the help of IeDEA-SA cohorts [1] which estimated that 17.5% of adult women under ART were older than 49. Percentage in c) is calculated with the help of both IeDEA-SA estimate and World Bank [25], which estimated that 54.6% of women aged 15-49 in South Africa were using any contraception method in 2015: As no information about contraceptive prevalence in South African adult women on ART have been found, we approximated it by the contraceptive prevalence in the general South African adult women population (see 17). By denition, the percentage p 0 of DTG-eligible men is 100%.

Sensitivity analysis
In the sensitivity analysis, we perturbed eight parameters 200 times using a Latin Hypercube Sampling method (see Table E). Table E displays the main values of the eight parameters, which were informed from literature and lower and upper bounds, chosen to reect plausible values of the parameters. As varying the transmissionrelated parameters may modify the overall transmission rate, an adjustment is made to have a transmission rate similar to the baseline model. We ran the sensitivity analysis for each prospective scenario (13 dierent scenarios in total).

Model ODEs
The rates γ represent the transition between care stages, ν CD4 the transition between CD4 stages and µ ij the mortality. The rate σ rev represents reversion of NNRTI-resistance when no more drug pressure is exerted, while σ N N RT I res and σ N RT I res represents the rates of acquiring NNRTI resistance and NRTI resistance, respectively, when an individual is failing NNRTI-based treatment. To model new infections, we used β u and β d the respective monthly number of sexual contacts among undiagnosed and diagnosed individuals, ρ k,k the assumed proportion of heterosexual individuals within men and women and ν k,k the probability of HIV transmission per sexual act. Finally, we also used a function δ(x), which is given by : S jklm 1,inel (t) = −ν S1,j−1 CD4 · S jklm 1,inel (t)1 j≥2 +ν S1,j CD4 · S (j+1)klm 1,inel (t)1 j≤3 − γ jl S1→F1 · S jklm 1,inel (t) + γ jl T1→S1 · T jklm 1,inel (t) + γ jl F1→S1 · F jklm 1,inel (t) − µ j S1 · S jklm 1,inel (t), Supplementary Material

Sensitivity analysis and additional results
The main sensitivity analysis reects the uncertainty about eight parameters related to resistance and HIVtransmission. The 95% sensitivity ranges in 2040 for each scenario are shown in Fig 3 of the main manuscript. The evolution of uncertainty over time is represented in Fig C, which displays the 95% sensitivity ranges from 2005 to 2040 for each scenario. The dierence in NNRTI TDR levels over time between the dierent scenarios of DTG-introduction and the scenario where DTG is not introduced is displayed in Fig D with the 95% sensitivity ranges. Fig E diplays the predicted percentage of women failing NNRTI-based ART after 1 and 2 years of ART in 2025 and 2035, depending on the scenario of the rollout of DTG-based ART. Finally, we simulated three sensitivity analyses in order to investigate 1) the impact of the Treat-All policy, 2) the impact of treatment interruption, and 3) the impact of NRTI-resistance and higher ecacy of DTG.  q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Eect of no Treat-All policy
We previously assumed that the Treat-All policy increased the treatment initiation rates for people with CD4 > 200 cells/µL from 2017 to 2022. Here, we investigated the scenario where the Treat-All policy does not have any impact on the treatment initiation rates (which is equivalent to assuming no Treat-All policy, see FigF). Globally, assuming a Treat-All policy increases the levels of NNRTI PDR for each scenario, but does not change our conclusion (see FigG).  Calendar year NNRTI TDR (in %)

C No Treat−All policy
All on NNRTI Only men on DTG All men and women not of childbearing age on DTG All men and women not at risk of pregnancy on DTG All men and women on DTG Fig G: Levels of NNRTI resistance when assuming increase in treatment initiation rates due to the Treat-All policy ("Baseline Model") and when assuming identical treatment initiation rates from 2017 ("No Treat-All policy"). Dolutegravir is introduced in 2020 under three scenarios: DTG as rst-line regimen for ART-initiators (panel A), DTG for all patients (panel B), DTG for all patients, assuming an OR of failure of 2 when having NRTI-resistance (panel C), and with dierent eligibility criteria for women (colors).

Eect of treatment interruption
We introduced treatment interruption rates for the three ART regimens. Table G shows these rates estimated from IeDEA-SA data [1]. The introduction of treatment interruption did not substantially change the results (see Fig H). Calendar year NNRTI TDR (in %)

C Treatment interruption
All on NNRTI Only men on DTG All men and women not of childbearing age on DTG All men and women not at risk of pregnancy on DTG All men and women on DTG Fig H: Levels of NNRTI resistance using the baseline model ("Baseline Model") and when including treatment interruption ("Treatment interruption"). Dolutegravir is introduced in 2020 under three scenarios: DTG as rstline regimen for ART-initiators (panel A) or DTG for all patients (panel B), DTG for all patients, assuming an OR of failure of 2 when having NRTI-resistance (panel C), and with dierent eligibility criteria for women (colors).

Eect of NRTI-resistance and higher ecacy of DTG
We assessed the impact of NRTI resistance and higher ecacy of DTG-based regimen on the level of NNRTI PDR. We investigated three dierent scenarios regarding the impact of NRTI-resistance: 1) no impact (i.e OR of failure between NRTI-resistant and NRTI-susceptible individuals equals 1), 2) OR=2, and 3) OR=5. A meta-analysis comparing DTG-monotherapy with DTG-dual therapy found an odds ratio of failure of 13.9 after 48 weeks (8.9% vs 0.7% of failure, respectively) [17] However, we expect lower dierence between NRTI-resistant and -susceptible individuals, as some activity of the NRTI-backbones are observed even in resistant individuals [15]. Another study comparing DTG-ecacy according to the presence of specic NRTI-mutations found a HR of 3.23 (95%CI: 0.27-38.40) when having the K65R mutation and a HR of 0.99 (95%CI: 0.19-5.21) when having the M184V, suggesting low impact of NRTI-resistance on DTG-failure [16]. We also investigated three dierent scenarios regarding the ecacy of DTG compared with NNRTI: 1) OR of failure between NNRTI-and DTG-based regimen of 1.02, 2) OR=2, and 3) OR=5. The rst scenario refers to the results of the NAMSAL study after the adjusting for CD4 counts (see Section 2.5). The two other scenarios were investigated in view of the higher ecacy of DTG compared with NNRTI found in some studies [26].