Ascertaining the 2004-2006 HIV type 1 CRF07_BC outbreak among injecting drug users in Taiwan.

OBJECTIVE: To ascertain the explosive 2004-2006 outbreak of HIV-1 CRF07_BC among intravenous drug users (IDU) in Taiwan, which more than doubled the total number of reported HIV cases in less than 3 years, resulting in a 45-fold increase in cumulative IDU/HIV cases and a 40-fold increase in previously seldom-reported female IDU/HIV cases.
METHODS: A mathematical model was utilized to fit the monthly case data, in order to estimate the turning points (peak incidence) and the reproduction number R of the outbreak. Furthermore, correlation analysis was carried out to assess the correlation between infections among the male and female IDUs.
RESULTS: Model fit revealed a two-wave epidemic during April 2004-March 2007. The larger second wave started shortly after May 2005 and peaked in October 2005 before gradually subsiding. R was estimated to be 3.15 (3.14-3.16) and 27.21 (26.73-28.05) for the two respective waves. The time series of monthly differences in male and female case data were found to be most significantly correlated at lag 0 (i.e., r>0.7) with r=0.906 and 0.804, respectively in each direction. The Granger causality test indicated that the male time series caused the corresponding female time series with a lag of 2 months or less.
CONCLUSIONS: The modeling results revealed the presence of a small first wave in 2004, before an explosion of cases after May 2005. Furthermore, a harm reduction program implemented in August 2005 contributed to the downturn in the epidemic after October. Correlation results also suggest that the upsurge in male HIV cases led to the subsequent drastic surge in female cases.

Introduction

In Taiwan, where active HIV/AIDS surveillance has been in place since 1996, with 2.1–2.5 million annual screening tests and free antiretroviral therapy (ART) introduced since 1997,1, 2, 3, 4 a dramatic increase in HIV incidence occurred in 2004–2006, almost tripling the cumulative HIV case number from 5650 at the end of 2003 to 15 345 by October 2007, and over 20 000 by 2010. The major cause of this upsurge is attributable to an explosive increase in HIV infections among intravenous drug users (IDU), from 112 total cases (2.15% of all cases) at the end of 2003 to 5034 (38.42%) by the end of 2006, a 45-fold increase within 3 years.

During this time-period, an outbreak of HIV-1 CRF07_BC infections among IDUs, including many previously seldom-seen female IDU HIV-infected cases, resulted in more reported cases than all reported HIV cases among all risk groups combined in the previous 20 years since 1984, when the first AIDS case was reported. It has been speculated that the source of this outbreak was a drug-trafficking route to Taiwan from Yunnan Province via southeast China, Guangxi Province, and Hong Kong,6, 7, 8 from where a substantial amount of heroin was being smuggled into Taiwan. Moreover, five IDUs from southern Taiwan were diagnosed as the country's first HIV-1-seropositive cases infected with CRF07_BC in 2002.2, 3 The total number of female HIV/AIDS cases nearly tripled between 2003 and 2005, from 320 female HIV/AIDS cases by January of 2003 to 903 cases by January 2006. The male-to-female ratio of HIV/AIDS cases was 13:1 before 2003 but declined to 7:1 in 2005–2006, with most of the newly reported female cases being IDUs. The cumulative number of infected female IDUs increased by 40-fold, from 16 cases at the end of 2003 to 648 by the end of 2006.

It has been reported that the percentage of persons receiving a diagnosis of AIDS within 12 months of diagnosis of HIV infection dropped suddenly from 20% in 2003 to 8.3% in 2005, during the time when most of the newly diagnosed cases came from the IDU population, which implies that the detected IDU cases were in the early stage of HIV infection. It has also been reported that 96% of HIV cases infected through IDU diagnosed in 2003–2004 were infected with HIV subtype CRF07_BC, which is totally different from the previously common subtype B and subtype CRF01_AE in Taiwan.

As is typical of underreporting of HIV/AIDS among many hard-to-count high-risk populations for HIV/AIDS,9, 10 there is a significant discrepancy between reported and estimated HIV/AIDS cases in the IDU group because of the difficulty reaching IDUs. In Taiwan, where needle-sharing or apparatus-sharing behaviors have been found to be common among IDUs, the majority of newly diagnosed HIV/IDUs were detected through mandatory inmate screening upon entry to correctional facilities. However, mandatory HIV screening of persons under police custody due to violation of the Narcotics Control Act since late 2004 could also have partially contributed to the sharp increase in detection.

In response to the outbreak in IDUs, a harm reduction program, which involved a needle-syringe program (NSP) and substitution treatment, was implemented by the Taiwanese government in August 2005 as an intervention to the rapidly increasing HIV epidemic since 2004. There has been a steady drop in HIV incidence among IDUs in recent years, down from being the major mode of HIV transmission in Taiwan in 2004–2006 to below that of homosexual and heterosexual transmissions since 2007, but still significantly higher than its pre-2003 level. A recent molecular epidemiology study after the CRF07_BC outbreak in Taiwan concluded that while the percentage of CRF07_BC among all HIV infections decreased from 2007 to 2009, the percentage of subtype B actually increased. However, many questions remain regarding this sudden outbreak among the IDU population.

In this study the reported HIV/IDU case data and a simple mathematical model, the Richards model,12, 13, 14 were used to investigate the temporal progression of this epidemic among IDUs in Taiwan. In particular, the total case data are fit to the model, as well as the male and female case datasets separately, in order to ascertain the epidemic. Correlation analysis was performed in an attempt to determine the relationship between the male and female IDUs.

Materials and methods

Data

The data used here were extracted from the monthly reported HIV case data between April 2004 and March 2007, for a total of 36 months, made available by the Taiwan Centers for Disease Control and Prevention (TCDC) on the TCDC website. The data are provided for each risk group/factor and gender. However, the HIV/IDU case data by gender are only available after August 2004, the fifth month of the dataset. In what follows, for ease of illustration in the tables and figures, the months are numbered, namely, April 2004 is month 1 and March 2007 is month 36.

The Richards model

The Richards model, a logistic-type mathematical model was used in this study. The explicit solution of the Richards model is of the form:where C(t) is the cumulative number of deaths at week t and the prime ‘′’ denotes the time rate of change. K is the final outbreak size over a single wave of outbreak, r is the per capita growth rate of the cumulative case number, a is the exponent of deviation of the cumulative case curve, and t is the turning point of the epidemic (which signifies the moment of upturn or downturn for the increase in the cumulative case number).

The basic premise of the Richards model is that the incidence curve of a single wave of infections consists of a single peak of high incidence, resulting in an S-shaped cumulative case curve and a single turning point (or the inflection point of the cumulative case curve) of the outbreak. This turning point t i, which is defined as the point in time at which the rate of accumulation changes from increasing to decreasing, or vice versa, can easily be pinpointed via the Richards model.

When more than one wave of infection occurs, a variation of the S-shaped Richards model is proposed, which makes the distinction between two types of turning points. Other than the first turning point ending the initial exponential growth of the cumulative case number, a second type of turning point is present in a multi-wave epidemic where the growth rate of the cumulative case number begins to increase again, signifying the beginning of the next wave. For further illustrations, the readers are referred to Hsieh and Cheng and Hsieh and Chen, in which the incidence curves for the 2003 Great Toronto Area severe acute respiratory syndrome (SARS) and the 2007 Taiwan dengue outbreaks containing two peaks (or two turning points of the first type) and one valley (or a turning point of second type) are investigated.

For the computation of the basic reproduction number R 0, the formula R 0  =  exp (rT) was used, where T is the generation interval of the disease or the average interval from onset of one individual to the onset of his/her contacts. It has been shown mathematically that, given the growth rate r, the expression R 0  =  exp (rT) provides an upper bound for the basic reproduction number, regardless of the assumed distribution of the generation interval. We noted that in this instance, the estimate obtained is not the basic reproduction number, but the effective reproduction number R, since in Taiwan HIV is endemic among the IDU population and multiple intervention measures have already been in place for some years.

The model parameters of epidemiological importance are K, r, and the turning point t i of the epidemic. The cumulative death data can be fitted to the Richards model to obtain estimates of these model parameters, using any standard software with a least-squares approximation tool, e.g., SAS, MATLAB, etc. More applications of the Richards model on other infectious disease outbreaks such as dengue can also be found in Hsieh and Chen and Hsieh and Stefan.

Statistical methods

Test for stationarity (unit root test, ADF)

We first examined whether the times series were ‘stationary’, in order to avoid spurious regression, which could possibly result in a biased and inconsistent estimator. A stationary time series means that its statistical characteristics do not change in time. In the event of a non-stationary time series, it can be suitably transformed to achieve stationarity. The augmented Dickey–Fuller (ADF) test was employed to verify if the random variables were indeed stationary series. Three equations were used to test if the series process has a non-stationary character:

Without drift and trend terms:

Including drift term without trend term:

Both drift and trend terms:

The null hypothesis is non-stationary or unit root, i.e. H 0  :  γ 1  = ⋯ =  γ  = 0. The series {Yt} has unit root if we cannot reject the hypothesis.

Correlation coefficient, r

To determine the correlation between the epidemic among the male and female IDUs through their time series of monthly case numbers, a distributed lag model (DLM) was employed to describe the relationship between the male and female time series. A DLM is a regression model that includes current and lagged values of one or more explanatory variables. This model allows the determination of what the effects are for a change in a time series. The resulting correlation coefficient, r, is a useful measure of linear strength between two random variables. The mathematical formula for computing r is:where n is the number of pairs of data. The value of r is such that −1 ≤ r ≤ + 1. The ‘+’ and ‘−’ signs are used for positive linear correlations and negative linear correlations, respectively. If there is no linear correlation or a weak linear correlation, r is close to 0. A value near zero means that there is a random, nonlinear relationship between the two variables. ‘|r| ≤ 0.4’ means ‘low correlation’, ‘0.4 < |r| ≤ 0.7’ means ‘moderate correlation’, and ‘|r| > 0.7’ means ‘high correlation’. In other words, a correlation greater than 0.7 is generally described as strong, whereas a correlation less than 0.4 is generally described as weak.

Granger causality test

The Granger approach is used to ascertain how much of the current values of time series y can be explained by past values or some lagged values of time series y. The commonly used software EViews was developed originally by economists for use in economics applications, but can also be useful in other statistical applications. EViews version 5.0 (http://www.eviews.com/) was used to analyze the data. In general, it is better to use more lags rather than fewer lags, since the theory is couched in terms of the relevance of all past information. It is advisable to pick a lag length, l, which corresponds to the reasonable beliefs about the longest time over which one of the variables could help to predict the other.

EViews performs bivariate regressions of the form:for all possible pairs of series in the group. The reported F-statistics are the Wald statistics for the joint hypothesis:for each equation. The null hypothesis is that series x does not Granger-cause series y in the first regression and that y does not Granger-cause x in the second regression.

Results

The monthly time series data of reported HIV cases for male IDUs, female IDUs, and all IDUs in Taiwan were fit to the Richards model as in Figure 1 and Table 1 . For all IDUs, the model fit resulted in a two-wave epidemic from April 2004 to March 2007, separated by the month of May 2005. The turning points for the two waves were April 2005 and October 2005, respectively. For male IDUs, the data fit a two-wave epidemic between August 2004 and March 2007, separated also by the month of May 2005, with the turning points at March 2005 and October 2005. For female IDUs, the fitted model also yielded a two-wave epidemic between August 2004 and June 2006, separated again by the month of May 2005, with the turning points at April 2005 and October 2005. Note that the female dataset does not converge for any period after June 2006.

Figure 1

(a) Total IDU population in Taiwan modeled by a two-wave Richards model, with the two waves of months 1–14 (2004.04–2005.05) and 14–36 (2005.05–2007.03); (b) male IDU population in Taiwan modeled by a two-wave Richards model, with the two waves of months 5–14 (2004.08–2005.05) and 14–36 (2005.05–2007.03); (1c) female IDU population in Taiwan modeled by a two-wave Richards model, with the two waves of months 5–14 (2004.08–2005.05) and 14–27 (2005.05–2006.06).

Table 1

Estimation results for the HIV-infected IDU population in Taiwan during April 2004 to March 2007

Time period	Turning point ti (95% CI)	Growth rate r (95% CI)	R (95% CI)
			T = 5	T = 6
All cases
Months 1–14 (2004.04–2005.05)	11.25 (10.91, 11.58)	0.191 (0.179, 0.203)	2.60 (2.59, 2.61)	3.15 (3.14, 3.16)
Months 14–36 (2005.05–2007.03)	4.29 (3.03, 5.56)	0.551 (−0.346, 1.447)	15.69 (14.85, 16.53)	27.21 (26.37, 28.05)
Male
Months 5–14 (2004.08–2005.05)	6.94 (6.17, 7.72)	0.194 (0.178, 0.211)	2.64 (2.63, 2.66)	3.21 (3.19, 3.22)
Months 14–36 (2005.05–2007.03)	4.30 (3.02, 5.58)	0.497 (−0.243, 1.238)	12.03 (11.33, 12.72)	19.78 (19.08, 20.47)
Female
Months 5–14 (2004.08–2005.05)	7.05 (6.06, 8.05)	0.260 (0.226, 0.293)	3.66 (3.64, 3.69)	4.75 (4.72, 4.78)
Months 14–27 (2005.05–2006.06)	4.19 (2.42, 5.96)	0.368 (−0.316, 1.053)	6.31 (5.70, 6.91)	9.11 (8.51, 9.72)

IDU, injection drug user; CI, confidence interval; R, reproduction number.

The effective reproduction number R was computed for each wave. However, it is unclear what the generational interval is for HIV. It has been proposed that the rates of partner change for homosexuals and heterosexuals tend to be of the order of 1 year, but it is unclear what the rate of needle-syringe sharing is among IDUs,22, 23, 24 although it is most likely of shorter length than that of sexual transmission of HIV. Due to the lack of a reliable estimate for the HIV generation time among IDUs in the literature, a generational interval of 5 or 6 months was assumed, based on an estimate of doubling time for AIDS cases among IDUs in the northeastern USA early in the epidemic. It should be noted that the main purpose for estimating R in this study was to compare the transmissibility estimated using different datasets and to ascertain the relative temporal change in transmissibility that occurred in each wave of infection during the course of the epidemic.

To further ascertain the role of female and male IDUs in the spread of HIV, correlation analysis was carried out for the time series of reported HIV cases of male and female IDUs. The time series used in the correlation analysis covers the time periods of September 2004–March 2007 for the monthly reported cases and October 2004–March 2007 for the monthly differences. To test for stationarity of the time series under investigation, results of the ADF unit root test are given in Table 2 , where the time series of monthly reported cases and the difference of two successive monthly reported cases were found to be stationary for both males and females. Note that we use ‘Male’ and ‘Female’ to denote the respective time series of monthly reported male and female HIV cases, and ‘△Male’ and ‘△Female’ to denote the respective time series of differences of monthly male and female HIV cases.

Table 2

Results of the augmented Dickey–Fuller (ADF) unit root test, where the difference of monthly reported male cases at month t is △Male = Male − Male − 1. Similarly, the difference of monthly reported females cases at month t is △Female = Female − Female – 1. The ADF unit root test is t-statistic, where p-values are in parentheses

Test type	Male	Female	△Male	△Female
Type 1None	−0.418	−0.561	−10.896	−6.859
	(0.523)	(0.465)	(<0.001)a	(<0.001)a
Type 2Constant	−1.875	−3.153	−10.700	−6.727
	(0.339)	(0.033)a	(<0.001)a	(<0.001)a
Type 3Constant, linear trend	−3.662	−3.139	−4.060	−5.667
	(0.041)a	(0.116)	(0.020)a	(0.001)a

p-Value of <0.05 indicates that the null hypothesis of non-stationarity can be rejected.

Next, the correlation between the monthly reported case data (Male and Female) and monthly differences in case data (△Male and △Female) were analyzed. The correlation coefficient r is a useful measure of the linear strength between two random variables. Applying a univariate model of y  =  β 0  +  βx t − lag  +  ɛ for ‘x causes y’ (or x → y) with lag, Male → Female and △Male → △Female were found to be most significantly correlated (i.e., r  > 0.7) with r  = 0.906 and 0.804 with lag 0, respectively. The correlation plots for the correlation are given in Figure 2 , which indicates the most significant correlations between the male reported cases and female reported cases, as well as their differences, are consistently at lag 0 (in red).

Figure 2

Cross-correlation between Case and △Case for the time series of reported HIV cases among male and female IDUs between September 2004 and March 2007.

The Granger causality test was subsequently carried out between the time series of Male and Female and between the time series of △Male and △Female for time lags up to 5 months. Test results, also given in Table 3 with the causal direction for each pair indicated with an arrow, indicate that both male time series caused the corresponding female time series after a lag of 2 months.

Table 3

p-Values of the Granger causality test between the time series of Male and Female and between the time series of △Male and △Female, with <0.05 indicating significant causality

Causal direction	Lags
	1	2	3	4	5
Male → Female	0.848	0.007a	0.016a	0.049a	0.110
Female → Male	0.199	0.331	0.452	0.560	0.699
△Male → △Female	0.071	0.032a	0.068	0.162	0.329
△Female → △Male	0.127	0.274	0.437	0.628	0.734

p-Values of <0.05 indicate that the null hypothesis of non-stationarity can be rejected.

Discussion

The timelines of the epidemic for all cases, as well as for the male and female case numbers, are illustrated in Figure 3 , indicating good agreement among estimates of temporal progression of the epidemic using the three datasets, and pinpointing May 2005 as the month that separated the two waves obtained from all three datasets. It is further concluded that the turning points for the first of these two waves differ slightly at 16.25 (95% confidence interval (CI) 15.91, 16.58), 15.94 (95% CI 15.17, 16.72), and 16.05 (95% CI 15.06, 17.05), respectively, for the total, male, and female case data (see Table 1), but agree at month 16 or March 2005 when rounded off to the nearest integer. More precisely, the months 16.25 and 16.05 signify the eighth and second days, respectively, of April 2005, while 15.94 denotes the next to last day of March. Hence the use of the three datasets pinpoints this turning point (indicating the peak of the earlier wave) at the end of March or beginning of April in 2005.

Figure 3

Timeline of the HIV epidemic among the IDU population in Taiwan, April 2004–March 2007.

More importantly, the turning point for the second and much larger wave of these two waves, which signifies the peak of the 2004–2007 epidemic, is pinpointed at between the sixth and ninth day of October 2005 by using all three datasets – at month 22.29 (18 + 4.29), 22.30 (18 + 4.30), and 22.19 (18 + 4.19), for all cases, males, and females, respectively.

It is interesting to note that in a recent study, phylogenetic tree analysis of 451 HIV-infected inmates with an IDU history in Taiwan was employed to demonstrate that there were two waves of HIV-1 CRF07_BC infection from mainland China to Taiwan. Although no timeline of the two waves was available from the molecular study, it does corroborate our modeling results.

The estimates for R indicate that the transmissibility of HIV among the IDU population increased from the first wave to the second wave according to all three datasets, culminating in a peak around October 2005, when a downward trend ensued. It is interesting to note that the harm reduction program that was implemented by the Taiwanese government in August 2005 could have impacted the downturn in case numbers after October, as revealed by our modeling results. However, only a trial harm reduction program, which included a needle-syringe program (NSP) and substitution treatment in four counties, was established at that time. After 1 year of the pilot study, data indicated that the HIV incidence in cities with an NSP decreased from 13.9 to 13.3 per 100 000 persons compared to an incidence increase from 11.5 to 15.3 per 100 000 persons in cities without an NSP. Subsequently, the harm reduction program was expanded to the whole of Taiwan in July 2006. Therefore, the country-wide downturn in case numbers after October 2005 may only be partially attributable to the harm reduction program.

R for females was higher than that of males in the first wave, but lower in the second wave and ending earlier in June 2006, perhaps reflecting the more rapid initial increase in female incidence in the early stages of the epidemic. However, this drastic upsurge in female HIV-infected IDUs was relatively more difficult to sustain in the second wave, as indicated by the shorter length of this wave when compared to that of the males. Therefore the epidemic impacted female IDUs more drastically initially, but the overall magnitude of infection was still less than that of the males, which was the group most affected by the epidemic.

Correlation analysis of the time series of male and female cases indicates that the outbreak among female IDUs was most likely driven by the infections among male IDUs, who were far more numerous at the beginning of the epidemic. In the Granger causality test, the lag length corresponds to the longest time over which one of the variables could help to predict the other. Therefore, the minimum p-value at a lag of 2 months for both the time series of case numbers and the first differences of the time series suggests that the strongest causality relationship, of male infections causing female infections, was within 2 months, a rather immediate relationship that led to the sudden upsurge in reported female HIV/IDU cases shortly after the upsurge of the male IDU cases in early 2004. It also indicates that there is likely some needle-sharing or apparatus-sharing among and between male and female IDUs, although one cannot rule out the possibility of sexual contact between the male and female IDUs, since one study in Taiwan found that 99.1% of the HIV-infected IDUs participating in the study were heterosexual. Unfortunately no data on any relationships that might exist between the reported cases were available. The study further demonstrates how quickly a disease traditionally endemic in a male population can develop into an epidemic in the female population given appropriate circumstances; hence transmission across gender is an important aspect of disease surveillance.

In summary, the abrupt outbreak among IDUs in Taiwan in 2004–2007, which was caused by the recently introduced CRF07_BC recombinant, led to two waves of infection with increasing transmissibility as measured by the effective reproduction number R during each wave, suggesting the presence of a small first wave in 2004 before the explosion of cases after May 2005. This further demonstrates the future potential of real-time modeling and analysis of disease data as part of a disease surveillance system, which could conceivably detect and alert the authorities of a possible herald wave before the arrival of a major outbreak. The waves ended by March 2007, which could be attributable to a timely and effective harm reduction program implemented in August 2005; this has been essential in preventing further occurrences of wide-spread infections among IDUs in Taiwan since 2007. Furthermore, the infections among male IDUs led to the epidemic among the female IDUs, a population that had previously been mostly devoid of HIV infections.

The limitation of this modeling study arises mainly from the nature of the HIV surveillance data, which typically consists of a longer period of time due to the long HIV incubation period, and hence is highly dependent on temporal changes in testing, reporting, and interventions over the years. In this study, the data that were used spanned a period of roughly 3 years, during which time the explosive outbreak among IDUs emerged; this subsequently led to interventions that included wider testing of IDUs in Taiwan and the harm reduction program. Moreover, while the simple mathematical model that was employed was able to reveal the temporal progression, culmination, and conclusion of the outbreak, it was unable to further pinpoint the exact impact of these intervention measures on the outbreak, which would require much more detailed data on the reported cases as well as a far more complicated mathematical model.