Assessing the danger of self-sustained HIV epidemics in heterosexuals by population based phylogenetic cluster analysis

Abstract

Assessing the danger of transition of HIV transmission from a concentrated to a generalized epidemic is of major importance for public health. In this study, we develop a phylogeny-based statistical approach to address this question. As a case study, we use this to investigate the trends and determinants of HIV transmission among Swiss heterosexuals. We extract the corresponding transmission clusters from a phylogenetic tree. To capture the incomplete sampling, the delayed introduction of imported infections to Switzerland, and potential factors associated with basic reproductive number R0, we extend the branching process model to infer transmission parameters. Overall, the R0 is estimated to be 0.44 (95%-confidence interval 0.42-0.46) and it is decreasing by 11% per 10 years (4%-17%). Our findings indicate rather diminishing HIV transmission among Swiss heterosexuals far below the epidemic threshold. Generally, our approach allows to assess the danger of self-sustained epidemics from any viral sequence data.

Keywords: HIV; basic reproductive number; concentrated vs. generalised epidemic; epidemiology; global health; heterosexual; infectious disease; microbiology; molecular epidemiology; transmission; virus.

Figure 1.. Overall basic reproductive number R0 and R0 per subtype from stratified analysis.

The dark gray point indicates the overall basic reproductive number $R_0$ estimate (by neglecting the transmission chain subtypes) and the corresponding $95\%$-confidence interval is shown with the dark gray line and the gray-shaded band. The analogous results from the per-subtype stratified analysis are represented by colored points and lines, each color corresponding to one of the subtypes (B, C, CRF01_AE, CRF02_AG or A) or the group of subtypes (other).

Figure 2.. Time trends for R0.

The upper smaller panels show the time trends for $R_0$ from the subtype-stratified analyses, in which the $log\left(R_0\right)$’s were modeled as linear functions of establishment date (i.e., for each subtype the time trend rate was assumed to be constant). The colored shaded-bands correspond to the $95\%$-prediction bands. The (best-fitting) nonlinear time trend for $R_0$ from the overall analysis is displayed in the lower panel (dark gray curve) together with the $95\%$-prediction band (gray-shaded area). The black points represent the $R_0$ estimates from the per establishment year stratified analyses and the gray vertical lines the corresponding $95\%$-confidence intervals.

Figure 3.. Effect of different factors on the basic reproductive number R0 from the multivariate model with only linear factor terms.

The black square and the black line show the reference basic reproductive number $R_0$ and its $95\%$-confidence interval (for a transmission chain of subtype B which started on 1.1.1996, and in which the index case was diagnosed 3 years after the infection, was 32 years old upon infection, never reported on having sex with occasional partner and had the earliest CD4 cell count of 350 cells per μL). The vertical gray line separates the factors associated with lower $R_0$ (left; effect factor $<1$) and from the factors contributing to higher $R_0$ (right; effect factor $>1$). The black points on this line refer to the reference transmission chain. The colored and dark gray lines represent the effect sizes from multivariate model (black circles depicting the estimates) for different factors and their $95\%$-confidence intervals. The corresponding $p$-values are shown in the rightmost column. FUP, follow-up visit.

Figure 4.. Final multivariate model’s profile plots of factors associated with the basic reproductive number R0.

The vertical dotted lines depict the reference transmission chain (of subtype B, started on 1.1.1996, in which the observed index case did not report having sex with occasional partner and was diagnosed after 3 years after the infection). The left $y$-axis represents the basic reproductive number whereas the right $y$-axis corresponds to the relative values of $R_0$ as compared to the baseline $R_0$. The $R_0$ as the function of specific factor (with the other factors held fixed at the reference value) is displayed by the colored (for HIV-1 subtype) and the dark gray (establishment date, sexual risk behavior and time to diagnosis) lines. The vertical bars and the shaded bands, respectively, correspond to the $95\%$-confidence intervals.

Figure 5.. Graphical representation of our phylogeny-based statistical approach.

(i): HIV transmission among heterosexuals in Switzerland (white arrow) has never led to a self-sustained epidemic. However, the unknown potential of imported infections (black arrows) either from abroad or from other transmission groups in Switzerland remains a large concern. (ii): The HIV transmission chains corresponding to Swiss heterosexuals (depicted in red) were identified from the phylogenetic tree containing the SHCS and background viral sequences. (iii): Our mathematical model is based on the discrete-time branching process with nodes of three different types: sampled Swiss infection (red), unsampled Swiss infection (light red) and foreign infection infected by a Swiss index case before moving to Switzerland (green). (iv): Our method for inferring $R_0$ accounts for both imperfect sampling and modified transmission potential of the index case. (v): Moreover, it includes the baseline transmission chain characteristics to assess the determinants of $R_0$.

Appendix 1—figure 1.. Sensitivity analysis regarding the index case relative transmission potential.

Panel (i) shows the sensitivity of the $R_0$ estimates from baseline model and panel (ii) the sensitivity of the time trend factor. The colored lines represent the subtype-stratified analyses, while the results from the overall models are shown in gray. In the first sensitivity analysis, the ${\rho }_{\text{index}}$ of Swiss-originating transmission chains was held at $1$ and the ${\rho }_{\text{index}}$ of non-Swiss origin varied (solid lines). In the second analysis, the ${\rho }_{\text{index}}$ of Swiss and non-Swiss origin was the same (dashed lines). The dotted lines show the results from the sensitivity subanalysis including only the transmission chains of non-Swiss origin. The vertical and horizontal lines depict the parameters and estimates from the main analysis, respectively.

Appendix 1—figure 2.. Sensitivity analysis regarding the sampling density.

The index case relative transmission potential parameter ${\rho }_{\text{index}}$ was the same as used in the main analyses, while the sampling densities varied ($x$-axis). In the pooled analysis (larger plots) the sampling density was the same for all transmission chains. Panel (i) shows the corresponding estimates of the basic reproductive number $R_0$ and the time trend factor estimates are displayed in panel (ii). The dotted vertical lines depict the sampling densities used for each subtype in our study (subtype-stratified plots) and the mean sampling density over all transmission chains (overall plots). The horizontal dotted lines represent the estimates from the main analysis.

Appendix 1—figure 3.. Conservative (with respect to ongoing transmission) maximum number of completed transmission degrees by a given date.

The red lines show the date ($y$-axis) by which at least a certain number (red numbers) of transmission degrees have been completed for a transmission chain with a specific establishment date ($x$-axis). The diagonal dotted gray lines depict the number of years since the establishment date, and the horizontal blue line represents the last sampling date.

Appendix 1—figure 4.. Relative bias due to ongoing transmission.

The upper panel shows the relative bias of the basic reproductive number $R_0$ from the baseline model and the lower panel the relative bias of the linear time trend factor from the corresponding generalized linear model. The proportion of active transmission chains over time is represented by the black line. The relative bias associated with overestimation and underestimation is displayed with green and red bars-points, respectively. Absence of bias is depicted by the horizontal gray lines.

Appendix 1—figure 5.. Sensitivity analysis regarding the stuttering transmission chains assumption.

The Q-Q plots compare the hypothetical transmission chain size distributions ($y$-axis showing their empirical permilles) with the transmission chain size distribution (empirical permilles on the $x$-axis) inferred from the phylogeny. The upper left plot compares the distribution of the simulated transmission chain sizes based on the estimated $R_0$ with the (from the phylogeny) observed transmission chain sizes and thus verifies the $R_0$ estimate. The remaining plots compare the simulated transmission chain size distributions against the extracted transmission chain sizes for $R_0$ closer to $1$ to justify the subcritical transmission assumption. Each point represents a permille, hence the darker points indicate more overlapping permilles.

Appendix 1—figure 6.. Comparison of effect sizes in the multivariate model with linear terms only for different sexual risk behavior definitions of a transmission chain.

The thick lines with black circles show the original effect sizes (where the index case determined the sexual risk behavior of the transmission chain) and their $95\%$-confidence intervals. The empirical distribution of the effect sizes where a random individual in a transmission chain determines its sexual risk behavior is displayed by the shaded areas. The thinner horizontal double sided arrows with the filled circles correspond to the effect sizes and their $95\%$-confidence intervals for the transmission chain level fraction of follow-up visits (FUPs) with reported sex with occasional partner by any of the infected individuals from the transmission chain. The vertical dotted gray line depicts the reference $R_0$ from the original model, i.e., using the index case to define the sexual risk behavior.

Appendix 1—figure 7.. Comparison between the Poisson and the negative binomial offspring distribution baseline model R0 estimates.

The dark gray and colored lines show the estimates from the model with Poisson offspring distribution, while the black lines correspond to the negative binomial distribution. The index case relative transmission potential parameter ${\rho }_{\text{index}}$ was fixed to $1$ and the sampling density ($x$-axis) varied. In the overall analysis the sampling density was the same for all transmission chains regardless of their subtype. The vertical gray lines depict the sampling densities used for each subtype in our study (above panels) and the mean sampling density in the overall analysis (bottom panel).

Appendix 1—figure 8.. Sensitivity analysis regarding the transmission cluster definition.

The upper panel (i) compares the estimated $R_0$ with the original cluster definition (brighter lines) with the $R_0$ estimated based on the relaxed cluster definition (darker lines) from the overall analysis (in gray) and subtype-stratified analyses (in colors). Similarly, the bottom panel (ii) shows the comparison between the estimated time trend factors obtained from the transmission chain sizes based on different cluster definition thresholds.

Appendix 1—figure 9.. Subanalysis for the transmission chains with available follow-up information about sex with occasional partner of the index case compared to the main analysis with imputed data.

The effect sizes from the subanalysis are shown in brighter colors and those from the main analysis in dark. In the main analysis, the missing data were replaced by never reporting sex with an occasional partner.

Appendix 1—figure 10.. Empirical distribution of maximum likelihood (ML) estimator and the Wald-type confidence intervals (CI) coverage rates.

Each plot represents a single parameter from a single model (see Appendix 1—table 1 for the parameters overview including their values), where the number in the lower left corner denotes the parameter’s consecutive parameter number. The light gray-shaded area represents the proportion of the Wald-type $95\%$-CIs from the parametric bootstrap simulations which contained the true value (depicted by the vertical orange line), while the green-shaded area corresponds to those CIs from the simulations that missed the true value. The numbers in the upper left corners are the coverage rates from the parametric bootstrap. The original Wald $95\%$-CIs used in our study are displayed with the light orange-area. The dark blue and gray lines show the empirical distribution of ML estimators from the parametric bootstrap samples and the normal approximation based probability density function, respectively. The horizontal red lines depict the target coverage rate of $95\%$.

Appendix 1—figure 11.. Comparison of different types of 95%-confidence intervals (CI) with the normal approximation based Wald-type 95%-CIs.

Each column corresponds to a different type of CIs, namely the profile likelihood based CIs, the basic nonparametric bootstrap CIs and the basic parametric bootstrap CIs. Each row represents a single parameter (the overview of the parameters is provided in Appendix 1—table 1). The colorful lines show the specific CIs compared to the corresponding Wald-type CIs, namely their relative widths and positions. The gray-shaded areas represent the Wald-type $95\%$-CIs.