Correction to: Drivers of HIV-1 drug resistance to non-nucleoside reverse-transcriptase inhibitors (NNRTIs) in nine southern African countries: a modelling study

1 Additional data description 2 1.1 Systematic review of NNRTI PDR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Initial systematic review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Update to the systematic review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.3 Included papers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Data on four key indicators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1 Additional data description 1.1 Systematic review of NNRTI PDR

Initial systematic review
Survey data on PDR in adults was obtained from a systematic review conducted by [1]. The authors searched for studies in PubMed, Embase and for conference abstracts published between Jan 1, 2001 and Dec 31, 2016, and complemented the systematic review with additional unpublished datasets from WHO-supported surveys of drug resistance. From the items included in the original systematic review, we selected 56 surveys of NNRTI PDR conducted in one of our 9 countries of interest: Botswana, Lesotho, Malawi, Mozambique, Namibia, South Africa, Eswatini, Zambia and Zimbabwe.

Update to the systematic review
Search strategy. The original search strategy not available in the online appendix of the journal article [1]. Therefore, we reproduced the search according to the Methods section (Methods > Search strategy and selection criteria) of the published article, without any language restrictions. We reproduced the search in "All fields" allowing a mapping to thesaurus terms (Mesh and Emtree). We limited time range to the interval January 1, 2017 to July 31, 2019 for PubMed and to August 9, 2019 for Embase. We concentrated on 9 countries of interest. Identified research items. The search resulted in the identification of 525 items in Pubmed and 913 items in Embase (tables 1 and 2). After deduplication this led to 1,039 items to screen. Screening and inclusion. These 1,039 items were independently screened with title and abstract by two reviewers, CD and JR. The objective was to identify relevant studies according to our inclusion criteria: studies reporting original results on NNRTI PDR in adults in one of the included countries. This step led to the identification of 32 items that were then read by the two reviewers to confirm the inclusion criteria. This second step led to the inclusion of 8 additional studies.
Data extraction. Extracted variables included, for each study i, the number of ART initiators with NNRTI resistance mutations E i , the sample size N i and the median year of sampling T i .

Included papers
In total, we included 64 surveys of NNRTI PDR in adults conducted in 9 countries between 2000 and 2018, for a total of 14,567 individuals. All data used in the analysis is available to download from https: //github.com/jriou/nnrtipdr_sa.

Data on four key indicators
For each of the nine countries, we obtained yearly estimates from UNAIDS regarding the number of adults living with HIV (A t ), the number of adults living with HIV on ART (B t ), and the number of AIDS-related deaths among adults (C t ) for each year t ∈ {2000, . . . , 2018} [2]. This data constitute model-based estimates (using Spectrum or Thembisa [3,4]) and not direct measurements. Data on the size of the adult population (>15) of each country (D t ) was obtained from the World Bank website [5]. Data for the year 1999 (A 0 and D 0 ) was also extracted to inform on initial conditions.

Pubmed n = 525
Embase n = 913 Deduplication Unique items n = 1, 039 Title and abstract screening Screened n = 32 Full-text screening 2 Additional model description

HIV transmission dynamics
We developed a deterministic, population-based HIV-1 transmission model using ordinary differential equations (ODEs) (figure 2). Our objective was to jointly model four key indicators of the local HIV epidemic (HIV prevalence, ART coverage, AIDS-related mortality, and population size) and the level of PDR for the adult population of a given country during the period 2000-2018.
The adult population of each country is divided into six compartments: susceptible to HIV-1 (S), infected with a NNRTI-sensitive strain of HIV-1 and untreated (I), infected with a NNRTI-sensitive strain and treated with first-line ART (T ), infected with a NNRTI-resistant strain and untreated (J), infected with a NNRTIresistant strain and treated with first-line ART (virological failure, U ), and infected with a NNRTI-resistant strain and treated with second-line ART (V ). The model considers the following dynamics: Transitions between compartments reflect the following mechanisms: • Transmission. Susceptible individuals may acquire a NNRTI-sensitive strain of HIV-1 from I, or a NNRTI-resistant strain from J or U , according to a transmission rate β. This assumes that (1) individuals infected with a NNRTI-sensitive strain and treated with first-line ART (compartment T ) or individuals infected with a NNRTI-resistant strain and treated with second-line ART (compartment V ) do not transmit the virus at a meaningful level [6] and (2) there is no fitness cost associated with the transmission of NNRTI-resistant strains [7].
• First-line, NNRTI-based treatment. The progressive roll-out of NNRTI-based ART from the early 2000s is modeled using a logistic growth function. This feature is implemented as follows: a maximum treatment rate τ is constrained by an indicator variable f (t, ν, ξ): which takes as input ν the time until treatment reaches 50% of its maximal value, and ξ the steepness of the logistic curve. This rate applies to both I and J compartments, assuming that individuals infected with a NNRTI-sensitive or NNRTI-resistant strain of HIV-1 are equally likely to initiate ART as drug resistance tests are not performed in routine.
• Emergence of de novo NNRTI resistance. Individuals infected by a NNRTI-sensitive strain and under first-line ART (compartment T ) may acquire resistance to NNRTI with a rate ω in response to the selective pressure created by the treatment. As the model does not differentiate between virological suppression, treatment failure or non-adherence (because no data on these aspects could be found for each country over the period), ω is interpreted as a compound measure of all phenomenons leading to the occurrence of drug resistance mutations under treatment including poor ART adherence, low patient retention, quality of ART service delivery and unreliable supply chains.
• Second-line, NNRTI-free ART. Adults who fail first-line ART due to NNRTI resistance (compartment U ) may be switched to second-line ART regimens that do not include with a rate κ.
• Demography. New susceptible adults enter the population at rate η. All adults can die from background mortality at rate µ, which we fix to the inverse of the life expectancy in adults in each country. HIV-1 infected adults without treatment (compartments I and J) or who fail firt-line ART without being switched to second line (compartment U ) may die from AIDS-related mortality at rate δ. We select t 0 = 1999 as the starting date. Initial conditions are set-up using data for the year 1999 (population size D 0 and HIV prevalence A 0 ). An additional parameter ι is introduced to represent the initial proportion of individuals infected with a resistant strain: The ODE system thus relies upon seven parameters related to the natural history of HIV in a country θ 1 = {β, τ, ν, ξ, κ, η, δ} (µ being fixed to a chosen value) and two parameters related to resistance to NNRTI θ 2 = {ω, ι}. Upon solving the ODE system with numerical methods [8], we obtain as output the number of individuals in each compartment S(t), I(t), J(t), T (t), U (t) and V (t) for each year t ∈ {2000, · · · , 2018}.
In addition to these, we introduce three dummy compartments that do not have any influence on the system dynamics but only record the cumulative incidence of AIDS-related deaths C D : and the cumulative incidence of new treatment initiations of individuals infected a with NNRTI-sensitive (C T ) or a resistant strain (C U ): Outputs from the ODE system are further transformed into a set of values {a t , . . . , e t } matching the structure of the data, as shown in Table 3.
E t /N t Proportion of ART naïve adults infected with an NNRTI-resistant strain in a survey of size conducted at year t

Hierarchical structure
The previous section described the model in one country. In order to extend the model to 9 countries, we consider 9 parallel and hierarchically-structured ODE systems. In country j, parameters related to the natural history of HIV θ 1,j are considered as independent across countries. We impose a hierarchical structure by country on the parameters related to NNRTI resistance θ 2,j , so that: where θ 3 = {µ ω , σ ω , µ ι , σ ι } are hyperparameters governing the between-country distribution of the "randomeffects" ω j and ι j . We can consider the full ODE model as a function with as inputs the parameters and as outputs the model-predicted HIV-1 prevalence, ART coverage, AIDS-related mortality, population size and NNRTI PDR in each country for the period 2000-2018:

Statistical inference
The outputs of the ODE system {a j,t , . . . , e j,t } at time t in country j are fitted to data in order to estimate the parameters . For the four key indicators of the HIV-1 epidemic, we use a multivariate normal likelihood after a variance-stabilizing square-root transformation [9], which can be expressed as: Here, Σ j is a diagonal covariance matrix: and {σ 2 j,a , · · · , σ 2 j,d } are variance parameters, with covariance parameters being fixed to 0. Survey data on NNRTI PDR has a different structure. Considering survey i, data consists of the number of ART-naïve adults with NNRTI resistance E i , the sample size N i , the year of sampling T i and the country of sampling J i . We use a binomial likelihood of the following form: The full likelihood can be expressed as: Rate of treatment failure due to NNRTI resistance Gamma(µ ω , σ ω ) y −1 ι j Initial proportion of NNRTI resistance in 1999 N (logit(µ ι ), σ ι ) y −1 Hyperparameters µ ω Location of the distribution of ω j across countries Expon(5) y −1 σ ω Scale of the distribution of ω j across countries Expon(20) y −1 µ ι Location of the distribution of ι j across countries Beta (1,9) y −1 σ ι Scale of the distribution of ι j across countries Expon(5) y −1 Variance parameters Σ j Covariance matrix, reparametrized into {σ a , · · · , σ d } --{σ j,a , · · · , σ j,d } Scales for the multivariate normal distribution Expon(1) -

Prior predictive checks
We selected weakly-informative priors distributions for all parameters (table 4). We conducted prior predictive checks [10] to ensure that the chosen priors limited the range of explored parameter space to sensible values. For this, we generated fake data from the prior distributions and checked that the resulting values were compatible with actual data, without constraining parameter inference to specific values. The results are shown in figures 3 and 4. The large range of values generated from the prior distributions corresponds to the possibilities implied by our choice of priors. This is an indication that our assumptions regarding the priors did not influence parameter inference, and that they can be qualified as "weakly-informative".

Implementation and code
The model was implemented in Stan 2.18.2 [11]. The joint posterior distribution of the parameters was explored with the No-U-Turn Sampler (NUTS), a variant of Hamiltonian Monte Carlo [12]. The model code, as well as posterior samples and R scripts are available on https://github.com/jriou/hiv_nnrtipdr_sa.