Elucidating Drivers for Variations in the Explosive Human Immunodeficiency Virus Epidemic Among People Who Inject Drugs in Pakistan.

BACKGROUND: Pakistan's explosive human immunodeficiency virus (HIV) epidemic among people who inject drugs (PWID) varies widely across cities. We evaluated possible drivers for these variations.
METHODS: Multivariable regression analyses were undertaken using data from 5 national surveys among PWID (n = 18 467; 2005-2017) to determine risk factors associated with variations in city-level HIV prevalence. A dynamic HIV model was used to estimate the population-attributable fraction (PAF; proportion of HIV infections prevented over 10 years when that risk factor is removed) of these risk factors to HIV transmission and impact on HIV incidence of reducing their prevalence.
RESULTS: Regression analyses suggested that city-level HIV prevalence is strongly associated with the prevalence of using professional injectors at last injection, heroin use in last month, and injecting ≥4 times per day. Through calibrating a model to these associations, we estimate that the 10-year PAFs of using professional injectors, heroin use, and frequent injecting are 45.3% (95% uncertainty interval [UI], 4.3%-79.7%), 45.9% (95% UI, 8.1%-78.4%), and 22.2% (95% UI, 2.0%-58.4%), respectively. Reducing to lowest city-level prevalences of using professional injectors (2.8%; median 89.9% reduction), heroin use (0.9%; median 91.2% reduction), and frequent injecting (0.1%; median 91.8% reduction) in 2020 reduces overall HIV incidence by 52.7% (95% UI, 6.1%-82.0%), 53.0% (95% UI, 11.3%-80.2%), and 28.1% (95% UI, 2.7%-66.6%), respectively, over 10 years.
CONCLUSIONS: Interventions should focus on these risk factors to control Pakistan's explosive HIV epidemic among PWID, including a concomitant expansion of high-coverage needle/syringe provision, opioid substitution therapy, and antiretroviral therapy.

Injecting drug use (IDU) is a global public health problem that is associated with a high burden of human immunodeficiency virus (HIV) [1, 2]. This is exacerbated in low- and middle-income country settings, where there is inadequate provision of harm reduction and treatment interventions among people who inject drugs (PWID) [3]. Although numerous countries have experienced steep declines in HIV among PWID when they have scaled up interventions [4], other countries are experiencing uncontrolled epidemics [5-7] or HIV outbreaks [8, 9]. These epidemics and outbreaks are generally associated with contextual factors that heighten the vulnerability of PWID, such as homelessness [9], conflict [7], economic recession [10], or high levels of criminalization preventing PWID from accessing services [11].

Pakistan has an expanding yet heterogeneous HIV epidemic among PWID, with the estimated national HIV prevalence increasing from 10.8% to 38.4% over 2005 to 2017 [5], and city-level HIV prevalences varying from 3.4% to 50.8% in 2016–2017 [12]. While reasons for this high burden of HIV among PWID are uncertain, previous epidemiological studies from Pakistan have suggested that risk behaviors such as high heroin and polydrug usage, high injecting frequency, sharing syringes and other injecting equipment, engaging in commercial sex, and widespread use of “professional injectors” who receive remuneration for injecting individuals may play a role [13-16]. This is compounded by a low coverage of needle and syringe provision (NSP) (<18 needles and syringes distributed to each PWID annually), lack of opioid substitution therapy (OST), and low coverage of antiretroviral therapy (ART) among HIV-infected PWID (<5% in 2017) [3, 17].

HIV in Pakistan is concentrated among PWID, but has been expanding in other key populations, including female sex workers and male/hijra sex workers, likely through infection between overlapping key populations [12, 18], as in other HIV epidemics in Asia [19]. Alongside this, HIV outbreaks in the general population have been occurring in Pakistan with increasing frequency in recent years [20]. This is thought to be due to high rates of therapeutic injections across Pakistan, which exposes the general population to increased risk of HIV bridging from key populations, most likely driven by the considerable HIV epidemic among PWID [21].

To better understand the transmission dynamics of HIV among PWID in Pakistan, and to optimize strategies to tackle this epidemic, it is important to understand what city-level risk factors or behaviors have contributed to the heterogeneous HIV epidemics among PWID across Pakistan. We consider this question using statistical and mathematical modeling techniques, utilizing data from 5 rounds of national integrated biological and behavioral surveillance (IBBS) surveys undertaken across 25 cities in Pakistan during 2005–2017. We then evaluate the potential impact of reducing the prevalence of these city-level risk factors to understand the importance of developing interventions to mitigate these risks.

METHODS

Data Sources and Analyses

We analyzed data from 5 cross-sectional IBBS survey rounds collected among PWID in Pakistan by the Canada-Pakistan HIV/AIDS Surveillance Project during 2005–2017 (n = 18 467), which covered 25 unique cities, with 13 cities surveyed in multiple rounds for a total of 57 city and round data points. The surveys included data on sociodemographic characteristics, sexual and injecting risk behaviors, uptake of services, recent HIV testing, and treatment. The HIV status of participants was determined using enzyme-linked immunosorbent assays performed on dried blood spot samples.

To examine which variables were associated with differences in city-level HIV prevalence (Supplementary Figure 1) among PWID over time, we restructured the dataset to give summary estimates for each variable for each city and round. Regression analyses determined associations between city-level HIV prevalence and population prevalences of different risk factors. Each survey variable was included separately in a univariable mixed-effects model with round nested within city and HIV prevalence as the outcome. All variables with a P value < .05 from univariable analyses were included in a multivariable model with the same structure. Variables in this model were then selected for the final multivariable model if they had a P value < .1. The same mixed-effects structure was also applied to the individual-level data using a multivariable logistic regression model to determine whether these variables were also associated with HIV at the individual level. All statistical analyses were conducted in Stata version 15.1 software.

PATIENT CONSENT STATEMENT

The study protocols for the IBBS surveys were reviewed and approved by the Pakistan Medical Research Council. Informed consent was obtained for all IBBS study participants. For this modeling study, patient consent was not required because it involved a retrospective analysis of de-identified data.

DYNAMIC MODEL STRUCTURE

The city-level multivariable regression analyses identified 3 factors associated with HIV prevalence at the city level: use of professional injector at last injection (“ProfInjUse”) vs no professional injector use, heroin use within the past month (“HeroinUse”) vs injecting other drugs (such as prescription opioids/tranquilizers; see Supplementary Figure 2), and frequent injecting (defined as injecting ≥4 times per day [“Inj4xpd”]) vs less frequent injecting. From these analyses, we constructed a dynamic, deterministic mathematical model of HIV transmission that incorporated states for PWID having 8 possible combinations of these 3 risk factors (Figure 1A). The HIV model dynamics (Figure 1B) assumed that individuals could either be susceptible to HIV infection or HIV infected, with HIV infection resulting in individuals progressing through acute, latent, and pre-AIDS stages of HIV disease, and transmissibility being heightened during acute infection and pre-AIDS HIV disease [22, 23]. We did not model ART as coverage is low (<5% in 2017) among PWID in Pakistan [3, 17].

Figure 1.

Schematic of dynamic model of human immunodeficiency virus (HIV) transmission among people who inject drugs (PWID) incorporating risk factors (A) and disease progression due to chronic HIV infection (B). Risk factors are assumed to be independent.

Individuals at model initiation are distributed into the 8 risk categories according to IBBS data on the proportion of PWID with each risk factor combination (Supplementary Materials). New individuals enter the model uninfected and without any risk factor. The model allows individuals to progressively transition to different risk categories at constant rates, and assumes that they then remain with that risk factor until they leave the modeled population. Susceptible individuals are infected with HIV at a rate dependent on their risk strata and prevalence of HIV for different disease stages. Individuals leave the model through cessation of IDU or due to natural or drug-related mortality, with the presence of HIV infection contributing heightened non-AIDS-related mortality [24] and additional AIDS-related mortality during the pre-AIDS stage. The initiation rate of injecting is set to balance these leaving rates, except for mortality due to HIV. Further model details are shown in the Supplementary Materials.

Baseline Model Parameterization and Calibration

City-level data for the model on injecting duration; prevalence of using professional injectors, heroin use, and frequent injecting; and HIV prevalence came from the IBBS surveys, with estimates being obtained for each city and round. Average injecting duration ranged from 2.1 to 10.8 years, with these durations remaining relatively stable for most cities across rounds (Supplementary Table 1; Supplementary Figure 3). The prevalence of using professional injectors and heroin use varied from <10% to 74% and <3% to >95%, respectively, while the prevalence of frequent injecting varied from <3% to 38%. Changes in the prevalence of risk factors across rounds were observed in some cities (Supplementary Table 1; Supplementary Figure 3). Across all cities and rounds, the proportions of PWID engaging in overlapping risk factors were consistent with each risk factor being independent of one another (eg, the proportion of PWID who use professional injectors and inject heroin is similar to the product of the proportions that have each risk factor across cities). HIV prevalence also exhibited heterogeneity across cities, with HIV generally increasing over time (Supplementary Table 1; Supplementary Figure 3).

There were no estimates of non-AIDS-related mortality for PWID in Pakistan, and so estimates from other Asian settings (India [25] and Vietnam [26]) were used, suggesting a non-AIDS-related mortality rate for PWID of 3.22 per 100 person-years among HIV-negative PWID, with this rate being adjusted by a mortality rate ratio of 1.73 for HIV-positive PWID [24]. We also assumed an average time from HIV seroconversion to death of approximately 10 years (Table 1) based on studies conducted in India and other low- and middle-income countries [28]. Analyses of phylogenetic data from a 2014 survey were used to estimate the differing start dates for the HIV epidemics in each city, which were thought to be between 1980 and 1990 for Hyderabad and Karachi; 1990 and 1995 for Quetta, Peshawar, and Larkana; and 1980 and 1995 for all other cities (unpublished data; Supplementary Materials). Uncertainty was incorporated in all model parameters and calibration data (Table 1; Supplementary Table 1). Last, no data exist on the degree to which PWID transition from not using professional injectors, not injecting heroin, or not injecting frequently (lower risk) to engaging in these behaviors (higher risk), and so they were determined by fitting to the prevalence of each city-level risk factor in each city (see next section). Similarly, the HIV transmission risks among PWID who do or do not use professional injectors, inject heroin, and/or inject frequently were also estimated through model calibration as described in the next section.

Table 1.

Model Parameters With Associated Uncertainty Ranges/Distributions

Parameters or Initial Conditions	Symbol	Baseline or Fitted Value (Uncertainty Distribution/Range)	Source/Comment
Demographic parameters
Recruitment rate to initiating use of professional injectors	κ_P	…	Fitted to prevalence of using professional injectors, heroin use, and frequent injecting (ie, injecting ≥4 times per day) across cities in each round of the IBBS survey. Uncertainty is incorporated by sampling from normal distributions derived from conducting binomial trials on the data for each city.
Recruitment rate to initiating heroin use	κ_H	…
Recruitment rate to initiating frequent injecting	κ_I	…
Injecting duration	ν	Range: 2.1 y in Peshawar (2008) to 10.8 y in Mirpurkhas (2016–2017)	IBBS city-level injecting duration estimates used—reciprocal gives the injecting cessation rate. Uncertainty incorporated through sampling uniformly between 0.5 and 2.0 times the point estimate of injecting duration.
Non-AIDS-related PWID mortality rate	μ	Crude mortality rate for HIV-negative PWID: Log-normal: 0.0322 (95% CI, .0416–.0399) Mortality rate ratio HIV+ vs HIV–: Log-normal: 1.73 (95% CI, 1.15–2.61)	Crude mortality rate for non-AIDS-related death among PWID, stratified by HIV status, came from systematic reviews [24, 27].
Population recruitment rate	ϕ	Set ϕ=ν+μ	Assume inflow into model matches outflow due to non-AIDS-related deaths.
Epidemic parameters
HIV transmission rate per susceptible	β	Prior: (Uniform: 0–0.10)	Estimated in model calibration—see Methods
Relative risk of HIV infection due to use of professional injectorsa	χ_P	Prior: (Uniform: 1–8)	Estimated in model calibration—see Methods
Relative risk of HIV infection due to heroin usea	χ_H	Prior: (Uniform: 1–8)	Estimated in model calibration—see Methods
Relative risk of HIV infection due to frequent injectinga	χ_I	Prior: (Uniform: 1–8)	Estimated in model calibration—see Methods
Enhanced HIV transmission risk by disease progression stage
Acute phase	ψ_acu	Log-normal: 276 (95% CI, 131–509)	[22, 23]
Latent phase	ψ_lat	Log-normal: 10.6 (95% CI, 7.6–13.3)	[22, 23]
Pre-AIDS phase	ψ_pre	Log-normal: 76 (95% CI, 41.3–128)	[22, 23]
HIV progression parameters
Duration of acute to latent/chronic	1/σ	Triangular: 2.90 (95% CI, 1.23–6.00) mo	[22]
Duration of latent to pre-AIDS	1/γ	Calculate duration of chronic phase using duration from seroconversion to AIDS	Derived using the formula: ^1╱_γ=^1╱_π−^1╱_σ−^1╱_ξ
Duration of pre-AIDS to AIDS or AIDS-related death	1/ξ	Triangular: 9.00 (95% CI, 4.81–14.0) mo	[22]
Duration of seroconversion to AIDS	1/π	Triangular: 10.2 (95% CI, 9.7–10.5) y	[28]
Initial conditions
Initial prevalence of using professional injectors		Uniform: 0–73.5%	City-specific data used from IBBS sampled within a range that gives an odds ratio of 0.5 to 2.0 compared to the point estimate.
Initial prevalence of heroin use		Uniform: 0–99.0%
Initial prevalence of frequent injecting		Uniform: 0.8%–38.3%
Initial HIV prevalence		Uniform: 0–52.5%	Initial seeded HIV prevalence is sampled from a range between 0.1% and up to 1% of the HIV prevalence for that city at the first available data point. These HIV infections are distributed proportionately across the different risk categories and are assumed to be in the latent stage of HIV infection.
Starting year of HIV epidemic		Uniform: 1980–1995	Assumed to vary across cities and is sampled uniformly between 1980 and 1995 as informed by phylogenetic analyses (see Supplementary Materials for details).

Rates are per year.

Abbreviations: CI, confidence interval; HIV, human immunodeficiency virus: IBBS, integrated biological and behavioral surveillance; PWID, people who inject drugs.

aIt is assumed that the risks of HIV infection due to use of professional injectors , heroin use , and frequent injecting are independent, so that the relative risk when multiple risk factors are present is multiplicative, eg, when 2 risk factors are present, the relative risks for their combinations are and when all 3 risks are present, the relative risk is .

Model Estimation of the Relative Risks due to Each City-Level Factor

In our modeling framework, we assumed that the observed associations between city-level HIV prevalence and the prevalence of using professional injectors, heroin use, and frequent injecting resulted from these city-level risk factors being associated with an increased risk of HIV transmission and acquisition. Our model calibration sought to estimate these increased risks through mimicking the observed city-level associations between HIV prevalence and the prevalence of the 3 risk factors. For simplicity, we assumed that the increased HIV transmission risks associated with using professional injectors, heroin use, or frequent injecting were constant across cities. We also assumed that these risks were independent. To mimic the associations, we did not calibrate closely to each city, but instead repeatedly simulated sets of 57 different HIV epidemics (using sampled values for the city and round prevalence of each risk factor and unknown baseline transmission risk, and city-level sampled values for duration of injecting and epidemic start date). For each set, we independently sampled values for the 3 relative risks (sampled from 1 to 8), seeking estimates that gave regression associations between the modeled data (using same mixed-effects multivariable regression model) that mimicked what we obtained from the empirical data. A good fit was defined as the modeled regression coefficients being within the 95% confidence intervals (CIs) of data-derived estimates, and the R2 value being within ±25% of the data estimate. The simulating of different epidemic sets was repeated until we obtained at least 1000 parameter sets that gave a good fit to the regression associations. In our analyses, a final 1080 parameter sets were obtained, which were defined as the “baseline model fits” and used to estimate the relative risks associated with using professional injectors, heroin use, and frequent injecting (median and 95% uncertainty interval [UI]). Further details are provided in the Supplementary Materials.

CITY-LEVEL ANALYSES

The parameter sets from the baseline model fits were used to undertake a more detailed investigation of the HIV epidemics in each city (n = 25). For the 13 cities with multiple rounds of data, these city-level models incorporated temporal changes in the prevalence of high-risk factors. These temporal changes were calibrated by varying the rates at which PWID initiate use of professional injectors, heroin use, and frequent injecting between survey rounds to fit to these changes (see details in Supplementary Materials). We then estimated a new baseline transmission risk (assumed constant over time) to recalibrate each baseline model fit to available data on overall HIV prevalence for that city across the IBBS rounds, accounting for uncertainty in HIV prevalence by sampling across the 95% CI for each city and round (Supplementary Table 1).

The resulting model fits for each city were used to estimate the contribution (population-attributable fraction, PAF) of using professional injectors (PAFP), heroin use (PAFH), and frequent injecting (PAFI) to new HIV infections over the next 10 years, assuming no future changes in the prevalence of these risk factors. The PAFs were estimated as the percentage reduction in new HIV infections over the next 10 years if the relative risk associated with each risk factor were set to one. We then estimated the overall national PAF for each risk factor, weighted by the estimated number of HIV-infected PWID in each city (details in Supplementary Materials and Supplementary Table 1).

Due to uncertainty in how individuals transition to and from using professional injectors, heroin use, and frequent injecting, we undertook sensitivity analyses to determine how the PAF estimates change if we assume that PWID either move back and forth between the risk categories or that PWID enter the model in specific risk categories and remain there (Supplementary Materials).

Last, the calibrated models were used to estimate the impact on the overall HIV incidence over 2020–2030 of reducing the prevalence of using professional injectors, heroin use, or frequent injecting to the lowest modeled levels in 2020 observed across all cities (2.8% [95% UI, 1.3%–4.8%], 0.9% [95% UI, 0.1%–2.4%], and 0.1% [95% UI, 0.01%–0.53%], respectively), corresponding to reductions in these city-level risk factors of 89.9% (95% UI, 85.1%–94.6%), 91.2% (95% UI, 85.5%–96.0%), and 91.8% (95% UI, 86.3%–96.3%), respectively. All modeling analyses were performed using MATLAB (version 2019a).

RESULTS

Association of HIV Prevalence With City-Level Risk Factors

Our multivariable regression model showed that the proportion of PWID in a city who used a professional injector at their last injection (regression coefficient = .25 [95% CI, .10–.40]), used heroin in the last month (= .19 [95% CI, .11–.26]), or injected frequently (≥4 times per day) in the last month (= .47 [95% CI, .23–.71]) were all strongly associated with HIV prevalence in a city (Table 2). Scatterplots of each association can be seen in Figure 2. No other variables were associated with HIV prevalence at the city level; for example, variables related to syringe or injection equipment sharing, housing status, engaging in commercial sex, and injecting drugs other than heroin were not found to be associated with HIV infection. The final multivariable model produced the following equation:

Table 2.

Univariable and Multivariable Coefficients for City-Level Human Immunodeficiency Virus Status, From Mixed-Effects Regression Models With City and Round as Random Effects

Variable	Univariable		Multivariable (Selection)		Multivariable (Final)
	Coefficient (95% CI)	P Value	Coefficient (95% CI)	P Value	Coefficient (95% CI)	P Value
Injecting duration (years)	0.01 (–.02 to .03)	.649	…		…
Injected 4 times per day last month	0.60 (.28–.92)	<.001	0.47 (.22–.72)	<.001	0.47 (.23–.71)	<.001
Used a used syringe, last time	0.10 (–.06 to .27)	.229	…		…
Same injection equipment used by others	0.01 (–.36 to .38)	.966	…		…
Currently lives on the street/lane	0.09 (–.12 to .30)	.412	…		…
Paid female sex worker for sex	–0.12 (–.33 to .09)	.258	…		…
Paid man/hijra for sex	–0.19 (–.42 to .04)	.100	…		…
Exchanged sex for money	–0.00 (–.20 to .19)	.965	…		…
Uses Restoril (temazepam) (capsule)?	0.06 (–.14 to .26)	.559	…		…
Uses diazepam?	–0.22 (–.34 to –.09)	.001	–0.06 (–.19 to .07)	.353	…
Uses heroin?	0.26 (.17–.35)	<.001	0.11 (–.01 to .23)	.063	0.19 (.11–.26)	<.001
Uses Pentazegon (pentazocine)?	–0.07 (–.31 to .17)	.570	…		…
Uses Phenergan (promethazine)?	0.04 (–.19 to .27)	.728	…		…
Uses Sosegon (pentazocine)?	–0.14 (–.31 to .03)	.105	…		…
Uses Marzine (cyclizine) (tablet)?	–0.03 (–.25 to .20)	.816	…		…
Uses Tamgesic (buprenorphine)?	–0.20 (–.33 to –.07)	.002	–0.05 (–.19 to .08)	.433	…
Uses other drugs?	0.19 (–.03 to .40)	.084	…		…
Did not always use clean syringe, last month	0.06 (–.07 to .19)	.382	…		…
Tried to clean used syringe/needle, last time	0.24 (.04–.44)	.016	0.02 (–.12 to .16)	.795	…
Injected outdoors, last time	0.19 (.02–.36)	.026	0.05 (–.08 to .18)	.433	…
Injected with others, last time	0.31 (.10–.51)	.003	0.12 (–.03 to .28)	.110	…
Injected by professional, last month	0.19 (.08–.31)	.001	0.02 (–.11 to .15)	.745	…
Injected by professional, last time	0.43 (.26–.60)	<.001	0.19 (–.01 to .39)	.066	0.25 (.10–.40)	.001

Variables from the univariable analyses with P < .05 were entered into the final multivariable model.

Abbreviation: CI, confidence interval.

Figure 2.

Scatterplot and corresponding data regression lines (solid black line showing the median, with dashed lines showing 95% confidence intervals [CIs]) and modeled regression lines (gray lines/shading) for the 57 city-level estimates of human immunodeficiency virus (HIV) prevalence against prevalence of using professional injectors (ProfInjUse), heroin use (HeroinUse), or frequent injecting (Inj4xpd). Note that the regression lines do not track the data points as closely as one would expect because the regression coefficients come from a model accounting for all 3 independent factors and the graphs shown are the regression lines projected onto the 2-dimensional plane with HIV prevalence and each of the risk factors separately.

where = –.05 (95% CI, –.11 to .01), with , , and having P values < .001 and an R2 value of 0.60. This regression model implies that for every 10-percentage point increase in the prevalence of using professional injectors, heroin use, or frequent injecting in a city, we would expect a 2.5 (95% CI, 1.0–4.0), 1.9 (95% CI, 1.1–2.6), and 4.7 (95% CI, 2.3–7.1) percentage point increase in HIV prevalence, respectively. At the individual level, HIV infection was associated with heroin use (odds ratio [OR], 1.21 [95% CI, 1.09–1.34]) and frequent injecting (OR, 1.60 [95% CI, 1.45–1.76]), with both P values < .001, but not use of professional injectors.

Estimation of Elevated HIV Transmission Risk Associated With City-Level Risk Factors

Using the final 1080 baseline fits from our dynamic model yields estimates that using professional injectors, heroin use, and frequent injecting are associated with a 2.3 (95% UI, 1.1–5.4), 1.9 (95% UI, 1.1–3.8), and 2.9 (95% UI, 1.2–6.7) increased relative risk of HIV infection, respectively (Figure 2; Supplementary Figure 4).

CITY-LEVEL SCENARIO ANALYSES

When the model was fit to data from multiple survey rounds for each city, modeled projections were within the 95% CI of the data for most points across cities and risk groups (Supplementary Figure 5; Supplementary Table 2). The modeled HIV epidemic trajectories suggest that overall HIV prevalence has increased from 24.7% (95% UI, 17.4%–34.4%) to 41.1% (95% UI, 33.2%–51.3%) between 2000 and 2020, coinciding with increases in the prevalence of using professional injectors and heroin use, stabilizing thereafter (Figure 3), with considerable variability between cities (Supplementary Figure 5). The model projected that HIV prevalence increased by >10 percentage points over 2000–2020 in 10 cities, increased by <10 percentage points in 14 cities, and decreased in 1 city (Hyderabad, 12.5 percentage point decrease). Overall HIV incidence peaked in 2011 at 22.4 (95% UI, 16.7–29.9) per 100 person-years, remaining relatively stable thereafter.

Figure 3.

Modeled overall human immunodeficiency virus (HIV) prevalence and HIV incidence, as well as overall prevalence of the 3 risk factors (use of professional injectors, heroin use, and frequent injecting), across all cities over time compared to published data estimates. The model for each city was calibrated to the 5 rounds of integrated biological and behavioral surveillance (IBBS) data on HIV prevalence at the city level. Both the modeled HIV prevalence projections and the IBBS data for each city were weighted by the estimated number of HIV infections among people who inject drugs in each city to give overall HIV prevalence projections and data estimates. The weighted data estimates only include cities with available data for the corresponding IBBS round, whereas the weighted model estimates include projections at each timepoint across all 25 cities. Overall HIV prevalence estimates were not fit to. HIV incidence is reported as new HIV infections per 100 person-years (PY).

For each city, Figure 4 shows that using professional injectors or heroin use will each contribute nearly half of new HIV infections over the next 10 years (PAFP = 45.3% [95% UI, 4.3%–79.7%] and PAFH = 45.9% [95% UI, 8.1%–78.4%]), frequent injecting will contribute nearly one-quarter (PAFI = 22.2% [95% UI, 2.0%–58.4%]), and all 3 combined will contribute 88.7% (95% UI, 74.9%–95.4%). There is considerable variability across cities, with the median PAF of using professional injectors, heroin use, and frequent injecting varying from 3.2% to 84.0%. However, the combined PAF showed less variability, varying from 78.3% to 96.0%.

Figure 4.

Population-attributable fraction (PAF) of using professional injectors (ProfInjUse), heroin use (HeroinUse), frequent injecting (Inj4xpd), or all 3 risk factors combined to new human immunodeficiency virus (HIV) infections over 10 years from 2020 for each city and overall. The lines show the median, limits of boxes give 25th and 75th percentiles, and whiskers give the 95% uncertainty interval across the final model runs. The weighted estimate is obtained by first calculating the weighted average across cities for each model run (weighted average by estimated number of HIV infections in each city) and then calculating the median, interquartile range, and 95% uncertainty interval across these average estimates.

Sensitivity analyses suggest that the 10-year PAF estimates were generally robust to different assumptions on how PWID transition between different risk groups, with median estimates from the sensitivity analyses varying <30% from the base case for PAFP, PAFH, and PAFI, except for the scenario where PWID enter the different risk categories when they start injecting and stay there, in which PAFH was reduced by around 40%. The combined PAF varied <3% relative to the base case across all sensitivity analyses (Supplementary Figure 6).

Using the city-level models, we found that reducing the prevalence of either using professional injectors or heroin use to the lowest modeled values from 2020 could halve the overall HIV incidence (52.7% [95% UI, 6.1%–82.0%] for using professional injectors, 53.0% [95% UI, 11.3%–80.2%] for heroin use; Supplementary Figure 7) across the 25 cities by 2030, with a smaller reduction if the prevalence of frequent injecting were reduced to the lowest value (28.1%; 95% UI, 2.7%–66.6%). Reducing the prevalence of all 3 risk factors concurrently would reduce overall HIV incidence by 87.4% (95% UI, 77.0%–94.9%). Smaller reductions are seen in HIV prevalence (Supplementary Figure 8).

DISCUSSION

Our analyses suggest that use of professional injectors, heroin use, and frequent injecting are 3 key risk factors that contribute substantially to HIV transmission in Pakistan, with these risk factors explaining nearly 90% of ongoing HIV transmission occurring across Pakistan. Use of professional injectors and heroin use may double the risk of HIV infection, while frequent injecting may triple HIV risk, with these factors contributing 45.3%, 45.9%, and 22.2% of all new HIV infections in Pakistan, respectively, and 88.7% combined. The smaller contribution of frequent injecting is due to its lower prevalence compared to the prevalence of using professional injectors and heroin use, which both showed increasing trends over time (Figure 3), with these changes likely explaining the increase in weighted HIV prevalence over time in Pakistan from 16.3% to 44.2% over 2005–2017. At the city level, projections suggest that changes in the prevalence of these risk factors coincide with observed changes in HIV prevalence. The prevalences of the 3 key risk factors, and HIV prevalence, showed a high degree of variability across cities and rounds, and even within the same city over time. It will be important to investigate the reasons for such variability in more detail. For example, in the IBBS dataset, we found that non-heroin-injecting PWID primarily injected other opioids (eg, buprenorphine [Temgesic], pentazocine [Sosegon]) and psychoactive drugs, namely, diazepam (Supplementary Figure 2), which were not associated with HIV infection. Cities reporting lower prevalence of heroin use in the past decade (eg, Pakpattan, Haripur) tended to be relatively smaller cities, which may have limited the availability of heroin.

Interestingly, the association between using professional injectors and HIV infection was observed at the city level, but not at the individual level. This may suggest that using professional injectors acts as a contextual factor increasing the overall risk of HIV at the community level, not just among PWID who use professional injectors. Alternatively, it is also possible that using professional injectors is a transient behavior, with IBBS data suggesting that PWID tend to use them earlier in their injecting career (Supplementary Materials), but less so later on, diluting any association between current use and HIV infection. The risk associated with using professional injectors is likely linked to their provision of needles and syringes, which may be reused between different clients, thus increasing HIV transmission [15]. This is possibly compounded by professional injectors being used more frequently by less experienced PWID who have less awareness of the risks associated with IDU. Indeed, large changes in the prevalence of using professional injectors over time (eg, >50% increase in Karachi and Larkana between 2006 and 2016) coincided with a reduction in the average age of initiating drug use (Karachi/Larkana: 27.7/28.6 years in 2006 vs 24.4/23.3 years in 2016). Meanwhile, our findings corroborate the well-established association between high injecting frequency and higher HIV prevalence [8]. The reason why heroin use is associated with greater HIV infection risk is less clear but could be due to heroin being more addictive than other drugs [29], or the prevalence of heroin use in cities being related to closer proximity to the opium trade route from Afghanistan to Turkey and then Europe [30]. This trade route may act as a channel along which HIV is disseminated.

Strengths and Limitations

To our knowledge, this is the first study to explore how city-level differences in the prevalence of risk factors or behaviors among PWID relate to the heterogeneous HIV epidemic in Pakistan. Strengths include utilizing detailed national datasets from 25 cities over 12 years (n = 18 467) and linking dynamic HIV transmission modeling for these cities to determine key risk factors that have been driving the expansion of HIV across Pakistan, which has allowed us to estimate the contribution of these key risk factors to incident HIV infections over the next 10 years, including the resulting chain of transmission.

However, there are limitations to our approach. First, we only considered factors arising from our multivariable regression analyses to develop the dynamic transmission model. It is possible that unmeasured factors not included in the model may be important in determining levels of HIV transmission. Moreover, we did not incorporate wider contextual-level factors that may be associated with HIV but were not asked about in the IBBS surveys, including homelessness [31] and recent incarceration [32]. Similar analyses could also include data from other sources on these and other measures, such as the levels of conflict, economic recession, or disparities (Gini coefficient), to further our understanding of how wider contextual-level factors influence the HIV epidemic in Pakistan. Second, we focused on city-level effects because we wanted to explore factors that are associated with the different HIV epidemics occurring across Pakistan. Individual-level differences were accounted for in the model simply by allowing variation in the baseline HIV transmission parameter. Incorporating individual-level determinants of HIV infection in greater detail may reveal more complex patterns of transmission within cities. In addition, we assumed random mixing among the various risk groups, whereas different mixing patterns by risk (eg, like-with-like) may affect the epidemic dynamics and should be considered in future studies. Third, we did not incorporate other complexities into the model, such as gender or sexual HIV transmission from other key populations. However, surveys suggest that few PWID in Pakistan are female (<2%) [16], even if they are undersampled, and HIV prevalence is much lower among other key populations [12], so the exclusion of these added factors is unlikely to affect our findings. Fourth, we did not model the effects of ART and harm reduction interventions (eg NSP, OST), which may play an important role in determining the HIV epidemic trajectory. The current coverage of these interventions is either nonexistent (OST) or extremely low among PWID (ART: <5%; NSP: <18 needles/syringes per year) [3, 17] and so their scale-up will be crucial in forthcoming years. This is the focus of ongoing work using the insights from this analysis. Fifth, to obtain national-level estimates for the PAF and changes in HIV incidence and HIV prevalence, we weighted the effects in each of the 25 cities by the estimated number of HIV infections in each city. This depended on empirical PWID population size estimates, which can be inaccurate. Fortunately, though, the geographic and network mapping exercise undertaken for the Pakistan IBBS [12] is considered to be a reliable and representative method of estimating key populations, and so our national aggregated estimates should be sufficiently accurate.

Comparison With Other Studies

Previous studies have undertaken statistical analyses of the national IBBS data in Pakistan. Archibald et al used data from the IBBS in 2006 and 2011 to evaluate multivariable associations with HIV prevalence by province, finding that HIV was associated with a longer injecting duration in Punjab, but not in Sindh [13]. Our analyses build on this previous work by pooling data from all 5 IBBS rounds to determine what risk factors are associated with HIV prevalence across cities. Another mathematical modeling study in Pakistan parameterized a simple HIV epidemic model using the same IBBS data, but focused on projecting the HIV prevalence and incidence over time among different key populations without considering what factors drive HIV transmission [33]. Our model considers a different question to investigate how city-level risk factors or behaviors determine the varying levels of HIV transmission across different cities in Pakistan. Last, another model considered reasons for the limited HIV epidemic among PWID in Rawalpindi before 2009 [23]. The HIV epidemic in Pakistan has expanded considerably since then, meaning this older model has little relevance now.

CONCLUSIONS

Our results suggest that future HIV programming in Pakistan needs to address 3 important risk factors shown in our analysis to be contributing nearly 90% of ongoing HIV transmission among PWID. In settings with low usage of professional injector services, targeted interventions need to prevent uptake of these services, possibly through awareness campaigns to communicate the risks associated with this practice. Conversely, in settings where the use of professional injectors is more prevalent, strategies such as supplying professional injectors with clean needles and syringes could reduce HIV transmission in the short term while longer-term prevention initiatives are developed and implemented. This should include expanding the availability of NSP, OST, and ART, which in other settings has controlled similarly high-prevalence HIV epidemics [8, 34, 35] and is highly cost-effective [36, 37]. For instance, OST encourages the cessation of opioid injecting [38] and reduces injecting frequency [39], while NSP reduces the HIV infection risk [40] associated with frequent injecting and could reduce the need for professional injectors if also linked with education on how to inject safely. The current coverage of these interventions in Pakistan is poor (<5% of HIV-positive PWID on ART in 2016–2017 [3, 17], and 18 clean needles and syringes provided per PWID per year [3]) or nonexistent (OST [3]); fulfilling these basic needs is an urgent priority that could dramatically decrease the transmission of HIV among PWID in Pakistan. This needs to be paired with increased engagement with law enforcement and civil society to reduce the stigma and harassment that PWID encounter when accessing harm reduction and treatment services. In addition, improved access to OST and ART will reduce drug and HIV-related morbidity and mortality among PWID, with OST having numerous other benefits [41] that are crucial for stabilizing the injecting and opioid epidemic in Pakistan. Our findings also suggest that shifting to injecting drugs other than heroin (eg, prescription opioids), which are not associated with HIV infection, or preventing transition to heroin injecting may lead to lower HIV transmission. Last, although HIV has been concentrated in PWID, it is expanding to other key populations, with multiple recent HIV outbreaks also occurring among the general population [20]. By intervening among PWID, added benefits will also be seen in these other groups due to their overlapping network interactions and bridging populations [18, 33].

Supplementary Data

Supplementary materials are available at Open Forum Infectious Diseases online. Consisting of data provided by the authors to benefit the reader, the posted materials are not copyedited and are the sole responsibility of the authors, so questions or comments should be addressed to the corresponding author.

Click here for additional data file.