Understanding the effects of different HIV transmission models in individual-based microsimulation of HIV epidemic dynamics in people who inject drugs

Summary

We investigated how different models of HIV transmission, and assumptions regarding the distribution of unprotected sex and syringe-sharing events (‘risk acts’), affect quantitative understanding of HIV transmission process in people who inject drugs (PWID). The individual-based model simulated HIV transmission in a dynamic sexual and injecting network representing New York City. We constructed four HIV transmission models: model 1, constant probabilities; model 2, random number of sexual and parenteral acts; model 3, viral load individual assigned; and model 4, two groups of partnerships (low and high risk). Overall, models with less heterogeneity were more sensitive to changes in numbers risk acts, producing HIV incidence up to four times higher than that empirically observed. Although all models overestimated HIV incidence, micro-simulations with greater heterogeneity in the HIV transmission modelling process produced more robust results and better reproduced empirical epidemic dynamics.

INTRODUCTION

HIV transmission between serodiscordant sexual and injecting partners is determined in part by the per-act HIV transmission probability, and by the number of unprotected sex and/or syringe-sharing events within partnerships [1–5]. However, prior modelling work has found that the simple binomial model for per-partnership transmission probabilities (with a constant per-act transmission probability) may provide an unsatisfactory model of HIV transmission, and is unable to reproduce HIV incidence and prevalence estimates observed in most settings [6]. Two major concerns have been raised with the binomial model for HIV transmission between serodiscordant partners. First, substantial heterogeneity in HIV susceptibility/infectiousness has been observed for both sexual and parenteral transmission pathways [1, 5, 7–11]. Second, homogeneous transmission models (i.e. those that assume a constant per-act transmission probability) produce inadequate infectivity estimates as suggested by observed data [12], and are misrepresentative of true risks associated with heterosexual contact. In part, these discrepancies are a manifestation of the assumption of independence between transmission risk and number of exposures [12].

Despite the fact that binomial models of HIV transmission are known to be inadequate, few studies have been conducted to develop, examine, and validate more sophisticated models of HIV transmission associated with injecting behaviour. This suggests that further work is needed to better estimate transmission probabilities related to contaminated injections and their likely contribution to overall HIV transmission in people who inject drugs (PWID) [13]. Finally, HIV infectiousness is known to vary substantially during the course of HIV infection (being higher during primary infection and in advanced disease), and also among persons in clinical latency [14, 15]. Therefore, further investigation is required to improve mathematical modelling efforts to reproduce and understand the HIV transmission process, particularly in PWID.

Individual-based models are an increasingly common method to better understand HIV transmission dynamics and evaluate HIV treatment and prevention interventions [16–19]. In contrast to compartmental models, individual-based microsimulations permit the analysis of interacting heterogeneous individuals and their behaviours in diverse environments [16–21]. Despite their increasing popularity in HIV research, the core model structures and processes that produce meaningful and valid outputs remain poorly understood [22]. A recently published systematic review on individual-based modelling of HIV transmission concluded that there is need for more consistent evaluation and explicit comparison of HIV transmission models to increase confidence in existing and future modelling results [23]. Here, we compare the validity of increasingly sophisticated (i.e. more complicated) models for HIV transmission among PWID in a North American setting, and how different models affect quantitative understanding and simulation of the HIV transmission process.

MATERIAL AND METHODS

We used an individual-based model, previously calibrated to be representative of the adult population (i.e. aged 15–59 years) of the New York metropolitan statistical area (MSA) [24–34]. This model therefore reproduces an urban, mixed HIV epidemic in which multiple low- and high-risk groups [i.e. PWID, men who have sex with men (MSM)] interact, and in which several modes of HIV transmission (vaginal, anal, parenteral) are present. This urban population within an established mixed epidemic is well studied [18, 24, 35–37] and serves well for our investigation in accounting for heterogeneity in HIV transmission, by comparing the different groups at risk (i.e. PWID and non-drug users). In order to account for differences in transmission risk across HIV disease stages, the microsimulation proceeds through discrete monthly time steps, thus simulating transmission during acute, latent, and AIDS stages of the disease. In this analysis, we used the same parameter values as were coded in the existing model [28, 36, 38–46], but observed changes in estimated HIV incidence after incorporating increasing heterogeneity in the models for HIV transmission and engagement in agent risk behaviour. Key parameters and other relevant data sources have been described in detail previously [18, 20].

To parameterize the individual-based model, we used previously collected empirical data [47], and an iterative indirect approach that has been described in detail previously [48]. First, we identified which real-world phenomena we were interested in reproducing (i.e. drug use prevalence, HIV prevalence/incidence, etc.), and successively developed a conceptual framework to guide the selection of processes and behaviours that would be modelled in the agent-based environment as described in Marshall et al. [20]. As a second step, we constructed a model that reflected known empirical and experimental evidence about these behaviours (see Tables 1–3). The parameter values in these tables represent point estimates (or in some cases, summary estimates) from relevant literature that were applied to the conceptual model structures for models 1–4. The relative ability of each conceptual model to reproduce estimates close to empirically observed HIV incidence and prevalence trajectories (providing information on which model structure has greatest validity) was of primary interest. Although extensive calibration procedures may theoretically enable each model type to reproduce empirical results for HIV incidence and prevalence, our goal was to determine which conceptual model structure was inherently more capable of reproducing these results, without extensive calibration procedures that may mask their limitations.

Table 1. Initial population distribution of the individual-based model (row percentages)

HM, Heterosexual male; HF, heterosexual female; MSM, men who have sex with men; PWID, people who inject drugs; PWUD, people who use drugs.

Proportions estimated from: [24–34].

Table 2. Initial parameter estimates and data sources for PWID (non-drug users) (PWUD) individuals

ART, Antiretroviral therapy; HAART, highly active antiretroviral therapy; HF, heterosexual female; HM, heterosexual male; MSM, men who have sex with men; NSP, needle-and-syringe exchange programme; PWID, people who inject drugs; PWUD, people who use drugs; SA, substance abuse.

* Defined as <100% correct condom use between individual dyads.

† Defined as <100% sterile syringe use with injecting partners.

‡ Defined as proportion of partners from preferred sexual orientation and same drug use strata.

§ Individuals who discontinue treatment at t = j can re-initiate treatment at some t > j with probability p = 0·18.¶ 60% of individuals achieve ⩾90% of adherence upon initiating ART (the remaining 40% are assigned to four other quartiles [0–29%, 30–49%, 50–69%, 70–89%] with probability 0·10).

¶ 60% of individuals achieve ⩾90% of adherence upon initiating ART [the remaining 40% are assigned to four other quartiles (0–29%, 30–49%, 50–69%, 70–89%) with probability 0·10].

Table 3. Per-act transmission probabilities according to type of transmission, sexual orientation, and stage of HIV disease

ART, Antiretroviral therapy; HF, heterosexual female; HM, heterosexual male; PWID, people who inject drugs; MSM, men who have sex with men, n.a., not applicable.

Agent network, characteristics, partnership formation and model structure

The individual ‘agents’ in the model are characterized by three time-varying drug-use categories [PWID, people who use non-injection drugs (PWUD), and non-drug users]. PWID agents represent individuals who have injected an illicit drug in the previous month, and PWUD represent recent users of an illicit drug (other than marijuana) through non-injection modes of consumption [109]. Agents are further stratified by sex (male, female) and sexual orientation – MSM, and heterosexual males (HM) and females (HF). Table 1 shows the population distribution specified when the individual-based model is initialized in the year 1992. During the first time step, sex, sexual orientation and drug-use status are attributed randomly to individuals to match the estimated proportion in the New York MSA in 1992, such that 6·5% are PWUD and 1·2% are PWID, of whom 7·3% and 6·5% are also MSM, respectively [25, 30, 31]. Among non-drug users, 2·2% are MSM [28, 31].

In the network, given an individual i, the number of (sexual and/or drug using) partnerships with other individuals at time step t (i.e. k _i,t ), follows a negative binomial distribution for all individuals per time step:

(1) $$\eqalign{& k_{i,t} \sim{\rm NB}\left( {\,p,r} \right) = \displaystyle{{\left( {k_{i,t} + r - 1} \right)!} \over {\left( {r - 1} \right)!k_{i,t} !}}p^r \left( {1 - p} \right)^{k_{i,t}}, \cr & k_{i,t} \in \aleph _0, \quad i = 1, \ldots, \; 20000, \quad t = 1, \ldots, \; 120,}$$

where partners are acquired with probability p until r suitable partners are found. Negative binomial regression was used because it is appropriate for overdispersed count data [110, 111]. This model for partnership formation and resulting partner distributions reflects empirical behaviour and data observed in previously conducted sexual and drug-using network studies (see Table 2 for sources). Partnership formation occurs when two individuals have the same sexual orientation or they are both PWID. The assortative mixing (i.e. sexual and/or needle-and-syringe (NS) sharing partnerships among people with similar risk for acquiring HIV) has been incorporated by weighting the probability of each contact to favour the formation of links between individuals with similar characteristics. For example, 90% of MSM agents who are not PWID interact exclusively with other MSM [85, 86], and the remaining 10% are linked randomly to other individuals with whom they engage in sexual intercourse. We assume that 80% of PWID have parenteral interaction only with individuals in their own subgroup [82, 83]. For PWUD, 60% are connected with other PWUD and 18% with PWID. More details about the assortative mixing are presented in Table 2.

In order to avoid overestimation of the sexual and/or NS sharing partnership turnover or underestimation of the partnership duration that may result from randomly re-assigning all links at each new time step, we developed the following algorithm that describes the process of partnership formation and dissolution. First, the number of target partners is randomly generated for each individual based on the negative binomial distributions described in equation (1), based on agent type (i.e. MSM, PWID). Second, agents whose partner numbers has decreased at the beginning of new month lose partnerships stochastically until the new partner number is reached. This process represents partnership dissolution (i.e. ‘break-ups’). Third, agents whose partner numbers increased at the beginning of new month are assigned new partners from the available pool of agents. This process proceeds iteratively through the agent population until all agents receive the targeted number of partners. Note that partnerships are also dissolved as a result of the death of one member of the pair. As the simulation proceeds, the number of partners for each agent varies; thus, both concurrent and sequential partnerships are possible.

In primary analyses, the model consisted of a population of N = 20 000 agents run over 120 monthly time steps (representing 10 years, 1992–2002). The model assumes that the characteristics (e.g. drug use, HIV disease stage) of each individual i, i = 1, … , 20 000 is updated on a discrete and monthly time step t, t = 1, … , 120, following pre-programmed rules and interactions with other individuals, 1, … , i − 1, i + 1, … , 20 000.

During the simulation HIV-uninfected agents can acquire HIV through unprotected sexual intercourse and/or NS sharing with HIV-infected agents. We defined unprotected intercourse between two agents as <100% correct and consistent condom use. Once infected, the natural history of HIV infection was modelled by considering three disease stages d [acute stage (AS), latent stage (LS), and AIDS], d = AS, LS, AIDS. Acute HIV infection is defined as the period immediately after infection during which the initial viraemia (and high infectiousness) occurs. On average, acute infection is considered to last 3 months following HIV acquisition [112, 113]; thus, we considered the first three time steps following infection as the acute phase.

At each time step, HIV-infected agents have a pre-specified probability of ‘accessing’ HIV testing; following which ART can be initiated (after 1996). Once an HIV-positive agent is assigned ART, the model randomly assigns one of the six adherence levels of ART (j, j = 0, … , 5), which varies between 0% and 100%, such that 60% of individuals achieve ⩾90% adherence [107]. These adherence levels correspond to a viral load (log₁₀ copies/ml) that decreases from 4·5 to 0·5.

HIV transmission probability models

HIV transmission within serodiscordant partnerships (i.e. dyads) is determined by the per-partner transmission probability, $\beta _p^{i,j,d,t} $ , which depends on per-act transmission probability, $\beta _a^{i,j,d,t} $ , and the number of unprotected sex, $n_{{\rm unp}}^{i,t} $ , and/or syringe-sharing acts, $n_{{\rm NS}}^{i,t} $ , within a partnership for each individual i (i = 1, … , 20 000), for every j (j = 0, … , 5) adherence levels, in disease stages d (d = AS, LS, AIDS), at time step t (t = 1, … , 120). We assume that the per-partnership transmission probability follows a Binomial distribution [1], where transmission is sexual or parenteral, respectively in equations (2) and (3):

(2) $$\eqalign{\beta _{\,p,{\rm unp}}^{i,j,d,t} & \sim {\rm Bin}\left( {1,n_{\rm unp}^{i,t}, \beta _{a,{\rm unp}}^{i,j,d,t}} \right) \cr & = \beta _{a,{\rm unp}}^{i,j,d,t} \left( {1 - \beta _{a,{\rm unp}}^{i,j,d,t}} \right)^{1 - n_{\rm unp}^{i,t}},}$$

(3) $$\eqalign{\beta _{p,{\rm NS}}^{i,j,d,t} & \sim{\rm Bin}\left( {1,n_{{\rm NS}}^{i,t}, \beta _{a,{\rm NS}}^{i,j,d,t}} \right) \cr & = \beta _{a,{\rm NS}}^{i,j,d,t} \left( {1 - \beta _{a,{\rm NS}}^{i,j,d,t}} \right)^{1 - n_{{\rm NS}}^{i,t}}.}$$

We constructed and analysed four HIV transmission models, each based on the previous iteration but incorporating more heterogeneity as described below.

Model 1: constant model

The simplest model assumed that the per-partnership transmission probability varies by ART adherence (for agents assigned ART) and HIV disease stage, but does not depend on the number of unprotected sex and the NS sharing acts, which are considered to be constant across all partnerships per time step. Although agent-based models of HIV transmission typically incorporate behavioural differences in engagement in risk behaviour at the individual level, we assess this model structure, as it is similar in concept to some deterministic models that have evaluated the preventive benefits of ART [114]. The per-act transmission probability is equal to c ^j,d , varying by ART adherence and HIV disease stage:

$$\eqalign{&{n_{{\rm unp}}^{i,t}} = n_{{\rm unp}} \,\; {\rm and}\,\beta _a^{i,j,d,t} = \beta _a^{\,j,d} = c_{{\rm unp}}^{\,j,d}, \quad j = 0, \ldots, \; 5;\cr & \qquad d = {\rm AS},{\rm \; LS}, \;{\rm AIDS\;} \,\left[ {{\rm sexual \;transmission}} \right].}$$

$$\eqalign{&{n_{\rm NS}}^{i,t} = {n_{\rm NS}}\,{\rm and}\,\beta_a^{i,j,d,t} = \beta_a^{j,d} = c_{{\rm NS}}^{j,d} ,\quad j = 0, \ldots ,\; 5; \cr & \qquad d = {\rm AS},\; {\rm LS},\; {\rm AIDS}\,\left[ {\hbox{parenteral transmission}} \right].}$$

Thus, depending on the type of HIV transmission, the per-partnership transmission probability is given by:

$$\eqalign{&{\beta _p^{i,j,d,t}} = c_{{\rm unp}}^{\,j,d} \left( {1 - c_{{\rm unp}}^{\,j,d}} \right)^{n_{{\rm unp}} - 1}, \quad j = 0, \ldots, \; 5;\cr & \qquad d = {\rm AS},\; {\rm LS},\; {\rm AIDS}\,\left[ {{\rm sexual \; transmission}} \right].}$$

$$\eqalign{&{\beta _p^{i,j,d,t}} = c_{{\rm NS}}^{j,d} \left( {1 - c_{{\rm NS}}^{j,d}} \right)^{n_{{\rm NS}} - 1}, \quad j = 0, \ldots, \; 5;\cr & \qquad d = {\rm AS},\; {\rm LS},\; {\rm AIDS}\,\left[ {{\rm parenteral \; transmission}} \right].}$$

Table 3 shows the values for the per-act transmission probabilities for syringe-sharing events [13, 115, 116] and unprotected sex acts [116–121]. Also shown are the monthly transition probabilities for progressing to AIDS [66, 122, 123].

Model 2: random number of acts

Previous studies have shown differences in condom usage by serostatus and by partners of PWID, and should thus be incorporated into HIV epidemic models [76]. The second model is more complicated by introducing stochasticity to the number of unprotected sex [equation (4)] and the NS sharing [equation (5)] acts. This model structure has been used to evaluate HIV prevention services (including NS programmes) for people who inject drugs [124, 125]. Specifically, these values are assigned within a partnership for each individual and from Poisson distributions [126]:

(4) $$\eqalign{&{n_{{\rm unp}}^{i,t}} \sim{\rm Poi}\left( {\lambda _{{\rm unp}}} \right) = \displaystyle{{\lambda _{{\rm unp}}^{n_{{\rm unp}}^{i,t}} {\rm e}^{ - \lambda _{{\rm unp}}}} \over {n_{{\rm unp}}^{i,t} !}},\; \quad n_{{\rm unp}}^{i,t} \ges 0;\cr & i = 1, \ldots, N;\quad t = 1, \ldots, \; 120,}$$

(5) $$\eqalign{&{n_{{\rm NS}}^{i,t} \sim{\rm Poi}\left( {\lambda _{{\rm NS}}} \right) = \displaystyle{{\lambda _{{\rm NS}}^{n_{{\rm NS}}^{i,t}} {\rm e}^{ - \lambda}_{{\rm NS}}}} \over {n_{{\rm NS}}^{i,t} !}},\; \quad n_{{\rm NS}}^{i,t} \ges 0;\cr & i = 1, \ldots, N;\quad t = 1, \ldots, \; 120,}$$

As in model 1, the per-partnership transmission probabilities also vary by partnership and time, and depend on highly active ART adherence levels (for agents on ART) and HIV disease stage.

Model 3: individual viral load

In the third model, we varied the per-act transmission probability by individual plasma HIV RNA viral load. We used each HIV-infected agent's viral load (stochastically assigned as described below) to calculate unique per-act risks of HIV transmission for each serodiscordant dyad.

Fraser et al. [127] found significant heterogeneity in asymptomatic (or set-point) viral load, which varies between and within the stages of the disease. Previous cross-sectional surveys have also shown that, during untreated latent stage, HIV viral load varies between 4 and 5 log₁₀ copies/ml [128–130]. Therefore, in model 3, we assumed a viral load of 4·5 log₁₀ copies/ml as the baseline, and assigned a viral load equal to 0 or 4·5 log₁₀ copies/ml to all the individuals before year 1996. It is also known that immediately after exposure and transmission, the viral load is undetectable in plasma and this generally lasts 7–21 days [131, 132]. Thus, when an individual becomes HIV infected, the adherence levels correspond to a log₁₀ viral load copies/ml that decreases from 4·5 to 0·0 during the first month. In the following 2 months the individual will be allocated a log₁₀ viral load between 6·9 and 0·5 copies/ml. During the acute and AIDS stage all agents were assumed to have 7·0 log₁₀ viral load copies/ml [14, 113, 133].

Each individual's per-act risk of HIV transmission was calculated based on results of a study by Baggaley et al. [119], in which the relationship between probability of HIV transmission and plasma HIV RNA viral load (copies/ml) was determined. The vaginal per-act probability of HIV transmission, $\beta _{{\rm a},{\rm vag}}^i $ , was best estimated as a function of an individual's viral load as follows:

$$\eqalign{& \beta _{{\rm a},{\rm vag}}^{i,j,d,t} = 1 - \left( {1 - \displaystyle{{0 \!\cdot\! 317\left( {{\rm VL}^{i,j,d,t}} \right)^{1 \cdot 02}} \over {\left( {{\rm VL}^{i,j,d,t}} \right)^{1 \cdot 02} + \;13\,938^{1 \cdot 02}}}} \right)^{\textstyle{1 \over {83 \!\cdot\! 17544}}}, \cr & i = 1, \ldots, \; 20\,000;\quad\!\! j = 0, \ldots, \; 5; \!\!\!\quad d = {\rm AS},\; {\rm LS},\; {\rm AIDS};\cr & t = 1, \ldots, \; 120,} $$

and the corresponding anal and parenteral HIV per-act probability, $\beta _{{\rm a},{\rm anal}}^{i,j,d,t} $ and $\beta _{{\rm a},{\rm inj}}^{i,j,d,t} $ , are respectively given by:

$$\eqalign{&\beta _{{\rm a},{\rm inj}}^{i,j,d,t} = 4.67 \times \beta _{{\rm a},{\rm vag}}^{i,j,d,t}, \,\; \; \; \; i = 1, \ldots, \; 20\,000;\cr & j = 0, \ldots, \; 5;\quad d = {\rm AS},\; {\rm LS},\; {\rm AIDS};\quad t = 1, \ldots, \; 120,}$$

$$\eqalign{&\beta _{{\rm a},{\rm anal}}^{i,j,d,t} = 3.50 \times \beta _{{\rm a},{\rm vag}}^{i,j,d,t}, \,\quad i = 1, \ldots, \; 20\,000;\cr & j = 0, \ldots, \; 5;\quad d = {\rm AS},\; {\rm LS},\; {\rm AIDS};\quad t = 1, \ldots, \; 120,}$$

where 4·67 and 3·50 represent the estimated increase in per-act risk of HIV transmission for parenteral and anal intercourse compared to vaginal intercourse, respectively [119, 134].

Model 4: two groups

Previous studies have suggested that observed heterogeneity in HIV transmission rates may be partially explained by differential risk behaviour across different types of partnerships. Earlier work has generally divided partnership types into two groups (i.e. regular partner, primary, and non-regular partners, casual) [5, 135]. In general, condom use is more frequent with casual partners than with a steady partner [136–138]. Similar discrepancies have also been observed between steady and non-regular partners for NS sharing (i.e. higher sharing rates with primary partners) [55, 139, 140].

Thus, the fourth and most complex model builds on the third, but considers two groups of partnerships, with either a higher (primary partner), ${}_{}^{{\rm primary}} n_{{\rm unp}}^{i,t} $ and ${}_{}^{{\rm primary}} n_{{\rm NS}}^{i,t} $ , or lower [casual partner(s)], ${}_{}^{{\rm casual}} n_{{\rm unp}}^{i,t} $ and ${}_{}^{{\rm casual}} n_{{\rm NS}}^{i,t} $ , number of unprotected sex and NS sharing acts, respectively. Stratifying the population into ‘low’ and ‘high’ risk groups is a common feature of previously published HIV transmission models [5, 55, 135–140].

Model outcomes and sensitivity analyses

The individual-based model was coded in Python^™ version 2.7.2, an open-source programming language [141], and the simulations conducted on a supercomputer at the Brown University Center for Computation and Visualization. The simulation was conducted using Monte Carlo methods to account for uncertainty in model outputs arising from the many processes and behaviours that are stochastic, by repeating each scenario 100 times. Mean estimates and 95% confidence intervals (CIs) for annual HIV incidence in the populations of non-drug users and PWID (and other outputs of interest) were obtained for each HIV transmission model. We then compared the trends of HIV prevalence and annualized incidence among PWID (per 100 persons), comparing models 1–4 to the empirical estimates for HIV prevalence [27], and incidence [142], observed in New York among PWID, from 1992 to 2002. Percent relative bias between model estimates and the HIV incidence rates observed empirically were calculated, and Pearson's χ ² test statistic [143] was used as a measurement of goodness of fit to determine which model produced more satisfactory estimates. A lower value of Pearson's χ ² test statistic indicates a better fit of the model to the observed data. As a secondary analysis, we also investigated the number of NS sharing acts per month outputted from the four different HIV transmission models.

Sensitivity analyses were performed to describe the extent to which changing the HIV transmission models and their respective parameter values affected the primary results. For each model (models 1–4), we considered five scenarios for the number of unprotected sex acts and syringe-sharing, respectively, and the per-act transmission probabilities, where the reference values were increased by 25% and 50%, and decreased by 25% and 50%, separately and then simultaneously.

RESULTS

Trends in HIV infections in PWID

Table 2 presents the values for monthly risk behaviour probabilities [28, 36, 38–46], NS sharing acts, and other model parameters (sources provided in table). In Figure 1, we show the trends in HIV prevalence and incidence obtained from each model, in addition to the empirical estimates. The most basic model (model 1, in which per-act transmission probabilities and numbers of unprotected sex and NS shared acts were constant), did not approximate the observed incidence and prevalence among PWID between 1992 and 2002, with Pearson's χ ² test statistic equal to 0·1248 and 0·1412, respectively. We were also unable to reproduce observed HIV epidemic trajectories when stochasticity was incorporated in the number of sexual and parenteral acts (model 2), resulting in the worst fit to the data compared to other models (i.e. the Pearson's χ ² test statistic was 0·2401 and 0·2559, respectively for HIV incidence and prevalence). Both models substantially overestimated HIV prevalence and incidence, particularly in later years of the simulation. In contrast, models 3 and 4, in which agents were assigned individual viral load values and also engaged in two types of partnerships, respectively, better re-produced HIV prevalence and incidence trends that were observed historically in PWID. For model 3 we obtained the following Pearson's χ ² test statistic: 0·0281 and 0·0203, respectively, for HIV incidence and prevalence. Model 4 presented better fit to the HIV incidence and prevalence in PWID, with Pearson's χ ² test statistic equal 0·0209 and 0·0107, respectively.

Fig. 1. (a) Estimated HIV prevalence and (b) annualized incidence in people who inject drugs (PWID) in New York, from 1992 to 2002, obtained from a Monte Carlo simulation of an individual-based model, considering four models for HIV transmission probability. In each panel, four HIV transmission probability models for sexual and parenteral transmission are presented: model 1, constant probabilities; model 2, random number of sexual and parenteral acts; model 3, viral load is individual assigned; model 4, two groups of partnerships (low and high risk). Red, dotted line indicates the empirical estimates of HIV (a) prevalence [27], and (b) incidence [142] observed in New York in PWID, from 1992 to 2002. In panel (b) the grey area represents the ‘burn-in’ period, since HIV incidence was annualized, and there are no estimated data before 1993. ART, Antiretroviral therapy.

The projected HIV prevalence in 2002 in the New York MSA was estimated to be 34% (95% CI 25–42), 40% (95% CI 32–49), 27% (95% CI 19–34), and 24% (95% CI 18–31) in PWID in models 1–4, respectively. HIV incidence in 2002 was estimated to be 2·4 (95% CI 0·0–5·0), 5·0 (95% CI 1·1–8·8), 1·3 (95% CI 0·0–3·0), and 1·0 (95% CI 0·0–2·7) per 100 person-years in models 1–4, respectively. In Figure 2, we show the relative bias in projected HIV incidence in 2002 compared to that which was actually observed in New York PWID [142]. Although all models overestimated HIV prevalence by the end of the simulation, the relative bias was smallest in model 4. Similar patterns were observed for relative bias of estimated HIV incidence in 2002 (Fig. 2b ). Once again, all models overestimated HIV incidence, but the relative bias was smallest in models 3 and 4.

Fig. 2. Relative bias in (a) HIV prevalence and (b) annualized HIV incidence in people who inject drugs (PWID) in New York, in 2002, obtained from a Monte Carlo simulation of the individual-based model, considering four models for HIV transmission probability. In each panel are represented estimates considering four HIV transmission probability models for sexual and parenteral infections: model 1, constant probabilities; model 2, random number of sexual and parenteral acts; model 3, viral load is individual assigned; model 4, two groups of partnerships (low and high risk). The percentage relative biases were calculated relative to the observed HIV prevalence and incidence observed in 2002, respectively to panel (a) and (b) HIV prevalence [27, 142].

Sensitivity analyses

In a series of sensitivity analyses (see Fig. 3), our simulation showed that the constant probabilities and random number of acts models were highly sensitive to changes in the number of sexual and/or parenteral acts per time step, in non-drug users and PWID, respectively. In scenarios where the number of sexual and/or parenteral acts were scaled up 25% and 50% (compared to those used in the primary analyses), unrealistically high HIV incidence rates in 2002 were projected in PWID: 11·8% (95% CI 5·2–18·5) and 22·2% (95% CI 12·6–31·8%) per 100 person-years in the constant probabilities model, and 8·9% (95% CI 3·3–14·6) and 15·2% (95% CI 6·3–24·2) per 100 person-years in the random number of acts model, respectively. HIV transmission models with individual viral load and the two-group partnership model were significantly less sensitive to changes in number of sexual and/or parenteral acts, with the individual-based model producing acceptable output (compared to the highest incidence observed in New York in PWID, from 1992 to 2002) [142] HIV incidence estimates in 2002 among PWID across all sensitivity analysis scenarios (see Fig. 3a ). Similar results were obtained when per-act transmission probabilities were increased (Fig. 3b ), and then number of monthly risk acts and per-act transmission probabilities were scaled up simultaneously (Fig. 3c ).

Fig. 3. Projected HIV incidence (per 100 person-years) in 2002 in people who inject drugs (PWID) in New York, considering four models for HIV transmission probabilities, for three different sensitivity analyses: (a) changes in number of risk acts, (b) per-act transmission probability, and (c) number of acts and per-act probabilities. In (a–c) for each HIV transmission probability model (1–4), we considered four scenarios for the number of unprotected sex acts and needle and syringe sharing, respectively, and the per-act transmission probabilities, where the baseline parameter values are increased by 25% and 50%, and decreased by 25% and 50%. In each panel, four HIV transmission probability models are presented: model 1, constant probabilities; model 2, random number of sexual and parenteral acts; model 3, viral load is individual assigned; model 4, two groups of partnerships (low and high risk). Red-dotted line indicates the highest HIV incidence (in log scale) observed in New York in PWID, from 1992 to 2002, as estimated in Des Jarlais et al. [142].

DISCUSSION

In this simulation study, we found that HIV transmission models with less heterogeneity were unable to reproduce HIV epidemic trajectories observed historically in PWID, and were highly sensitive to changes in key parameters (e.g. number of risk acts within serodiscordant partnerships). In contrast, models which incorporated more heterogeneity – for example, transmission dependent on individually assigned viral loads – produced satisfactory estimates of HIV prevalence and incidence, and showed less sensitivity to changes in key parameters. Our results support the work of previous mathematical modelling studies that have demonstrated the need for representing heterogeneity in behaviour and determinants of transmission probabilities to appropriately capture HIV transmission dynamics in non-drug using populations [144–146]. Furthermore, our study could serve as an important guideline for critical model structures and processes that should be included in future individual-based models of HIV transmission in dynamic networks.

The HIV transmission model that considered two types of partnerships (primary and casual) was able to better reproduce observed HIV prevalence and incidence trajectories than a model without partner-level differences in risk behaviour rates. This finding supports the work of Kretzschmar & Wiessing [147] who demonstrated that reducing rates of NS sharing with strangers is more effective than reducing the overall sharing frequency to decrease HIV transmission at the population level.

Some limitations to this simulation study should be noted. First, the modelled population size (20 000) was relatively small, and larger simulations are recommended to confirm our findings. However, we do not expect simulations with a larger population size to result in substantially different findings, given that we were able to qualitatively reproduce the HIV epidemic trajectories observed historically. Moreover, we parameterized the model with detailed data regarding risk behaviours among non-drug users, PWID, and PWUD (further stratified by sex and sexual orientation), which represents the MSA of New York. Finally, comparisons of the output from the various model structures were not limited by lack of precision in the primary estimates of interest (i.e. HIV incidence and prevalence in 2002). Second, although the effect of interventions such as provision of ART, NS programmes, and opioid substitution therapy were coded in the model, we did not incorporate other behavioural and structural interventions (e.g. serosorting). Furthermore, as our model sought to reproduce historical HIV epidemic trajectories, we did not incorporate interventions such as pre-exposure prophylaxis (PrEP) or post-exposure prophylaxis (PEP). Given the potential for PrEP and PEP to reduce the acquisition of HIV infection in MSM and PWID [148, 149], future research should be conducted to take into account the impact of these interventions in this setting. Third, although we did increase HIV transmission risk during the 3 months following infection, more sophisticated models are needed to understand the effect of varying durations of the acute infection phase in HIV transmission models. Fourth, further investigation is needed to lead to a better understanding why incidence and prevalence are overestimated by all four models. Friedman et al. [150] suggested that, in stable high-prevalence epidemic contexts, the relatively small size of sub-networks of linked seronegatives may limit infectious outbreaks. This process, known as the ‘firewall effect’ [150, 151], captures an important barrier to HIV propagation in PWID risk networks. Future work will investigate whether this effect is observed in the simulated network of PWID in our model. Fifth, a recent study indicated that current estimates of HIV-1 acute-phase infectivity relative to chronic infection may be substantially overestimated [152]. Future work will be conducted to determine whether reducing the relative acute-phase infectivity results in improved model fit to HIV incidence and prevalence trends. We note, however, that the effect of changing the relative infectivity in acute phase would not influence which model performed best, and consequently the interpretation of our current results. Sixth, our model did not incorporate more complex behavioural processes including mixing within and across age groups, age of sexual/injection debut, changing sexual/injection behaviour with age or time, and the effect of migration, given that there is not enough data to parameterize the model processes for some of these behavioural complexities. Future work will improve the individual-based model, with the goal of better understanding the HIV transmission dynamic process. Finally, as in all mathematical models, parameter values are subject to error. We conducted Monte Carlo simulations and sensitivity analyses to assess the degree to which uncertainty in parameter estimates affected our results, but nonetheless our findings should be interpreted within the context of this important limitation.

CONCLUSIONS

Overall, HIV transmission models with less heterogeneity were more sensitive to changes in numbers of sexual and parenteral acts, producing HIV incidence rates up to four times higher than those empirically observed. Although all models overestimated HIV prevalence and incidence, microsimulations with greater heterogeneity in the HIV transmission process, specifically transmission determined by individual viral loads, produced more robust results and better reproduced empirical epidemic dynamics among PWID in New York from 1992 to 2002.