Estimation of the basic reproductive number (R₀) for epidemic, highly pathogenic avian influenza subtype H5N1 spread

Summary

Three different methods were used for estimating the basic reproductive number (R₀) from data on 110 outbreaks of highly pathogenic avian influenza (HPAI) subtype H5N1 that occurred in village poultry in Romania, 12 May to 6 June 2006. We assumed a village-level infectious period of 7 days. The methods applied were GIS-based identification of nearest infectious neighbour (based on either Euclidean or road distance), the method of epidemic doubling time, and a susceptible–infectious (SI) modelling approach. In general, the estimated basic reproductive numbers were consistent: 2·14, 1·95, 2·68 and 2·21, respectively. Although the true basic reproductive number in this epidemic is unknown, results suggest that the use of a range of methods might be useful for characterizing epidemics of infectious diseases. Once the basic reproductive number has been estimated, better control strategies and targeted surveillance programmes can be designed.

INTRODUCTION

The basic reproductive number (R ₀) is the average number of secondary cases caused by one infectious individual during its entire infectious period in a fully susceptible population [1]. In epidemiological studies, the unit of interest might also be functional groups of individuals, such as villages, towns and cities, or herds and flocks. There is not an overall R ₀ for an infection – R ₀ is specific to populations. Estimating R ₀ for an infection in a population is critical to designing disease control and prevention strategies [1].

R ₀ consists of three components [1]: the rate of contact between the proportion (true mass action principle) or number (pseudo-mass action principle) of susceptible and infectious individuals (c), the probability of transmission on contact (p), and the duration of infectiousness (D). These components are related to the three principal factors determining infectious disease epidemiology: (1) natural history (‘course’) of infection in the individual; (2) transmission route; and (3) the environment and/or behavioural characteristics of the specific population. For the purposes of estimation from field data, these components can be combined into the number of contacts that occur during the infectious period (cD) and the number of ‘successful’ contacts (pc), also referred as the transmission coefficient (β).

The preceding expressions for R ₀ always contain some reference to the rate of contact between susceptible and infectious individuals (or epidemiological units). However, contact rate is difficult or impossible to estimate in the field for most infections, especially those transmitted directly via close contact or via multiple routes of infection [2]. Thus, other methods must be used to estimate R ₀. In general, these are based on empirical data, either at the invasion (epidemic) or endemic phases of the host–parasite relationship. During the epidemic phase, the full potential for disease spread (R ₀) is exhibited: there are no constraints on the growth rate of the epidemic. Reported data from epidemics are a valuable source for estimating R ₀. However, appropriate data is rarely available. Methods for estimating R ₀ from epidemic data rely on circumstances in which the population is naive (entirely susceptible) to the infection (e.g. foreign animal diseases) and the infectious disease spread is rapid and disease reporting occurs on a fine time-scale. Methods used to estimate R ₀ from epidemic data apparently have not been compared.

Highly pathogenic avian influenza (HPAI) virus subtype H5N1 is a threat to world health: it has caused numerous disease outbreaks in domestic poultry and wild bird populations, and there is a fear that it could become the next pandemic influenza strain [3]. Because exposure to sick or dead poultry is a strong risk factor for human disease caused by HPAI subtype H5N1 [4], the threat of pandemic influenza can be effectively reduced by controlling and preventing the spread of HPAI subtype H5N1 through national poultry flocks. Avian influenza virus infection is endemic in a range of free-living bird species worldwide [3], particularly species associated with water [5]. Waterfowl can be infected by all subtypes of type A influenza viruses, with few or no symptoms [6]. Although these species are capable of spreading influenza viruses between regions [7, 8], and HPAI outbreaks in poultry are sometimes assumed to occur from contact with wild avian species [9–11], the transport of infected poultry and contaminated poultry products have been blamed for spreading the disease within regions. To plan the most effective control programmes, quantitative estimates of the spread of the disease through populations of poultry are needed. Epidemic data provide a valuable insight into disease spread patterns.

The aim of the present study was to estimate R ₀ from a subset of data collected during the epidemic of HPAI subtype H5N1 that occurred in village poultry in Romania, 2005–2006. Tracing information from field investigation of the outbreaks was unavailable, therefore three different indirect estimation methods were used.

MATERIALS AND METHODS

Data source

The epidemiological characteristics of the 2005–2006 Romanian epidemic of HPAI subtype H5N1 have previously been described [12], and some spatial risk factors of outbreak occurrence have been investigated [13]. Briefly, the epidemic occurred between 7 October 2005 and 6 June 2006, and affected a total of 161 villages. Three phases of the epidemic were characterized, based on detailed geostatistical and spatial statistical analysis of the data [12]: October–December (23 outbreaks, 14·3%), January–March (28 outbreaks, 17·4%) and May–June (110 outbreaks, 68·3%). Risk mapping suggested that outbreaks first appeared in eastern and southern Romania. During the autumn and winter of 2005 (epidemic phases 1 and 2), the environment and landscape (specifically the Danube River Delta) played a critical role in the introduction and initial spread of HPAI subtype H5N1 [13]. The transport of poultry and poultry products might have introduced the infection into central Romania during the spring and summer of 2006, where the disease spread rapidly during a 26-day period, presumably via direct or indirect contact between villages. The current study focused exclusively on this final phase of the epidemic, which was temporally and spatially distinct from preceding epidemic phases that had occurred in Romania [12]. For each outbreak, its location (x,y-coordinates) and reported date of first detection – but not its probable infectious source – was available.

Data analysis

As is the case in many disease epidemics, detailed and reliable tracing information from field investigations of the outbreaks was unavailable. Thus, the basic reproductive number for between-village spread of HPAI subtype H5N1 was indirectly estimated from the dataset using three different methods – nearest infectious neighbour, a susceptible–infectious (SI) epidemic model, and epidemic doubling time.

Nearest infectious neighbour method

Reported outbreaks were mapped (ArcGIS™ 9.0; ESRI Inc., Redland, CA, USA) using a shape file (Dealul Piscului 1970 Datum, Stereographic 70 Projection) of Romanian administrative units (counties). For each outbreak, outbreaks reported up to 7 days previously were selected and Euclidean distances (km) were calculated (Hawth's Analysis Tools for ArcGIS; http://www.spatialecology.com/htools/). The closest outbreak was classified as the infection source of the outbreak-of-interest. This process was repeated until all outbreaks (except the first outbreak reported) had a designated source.

The frequency distribution and average number of outbreaks attributable to each source was calculated; the latter was an estimate of R ₀. The spatial distribution of R ₀ estimates was characterized by Moran's autocorrelation statistic (ArcGIS™ 9.2, Spatial Statistics), using an inverse distance weight. Local clusters of the R ₀ estimates were investigated using Anselin's local indicator of spatial autocorrelation (LISA) statistic (ArcGIS™ 9.0, Spatial Statistics). Spatial weights were defined by the inverse of Euclidean distance, and were globally standardized.

Outbreak locations were projected on to a map of primary and secondary Romanian roads, and village locations were assigned to the nearest road. The above procedure was repeated, except that the source for each outbreak-of-interest was identified as the outbreak that had occurred up to 7 days previously and that was closest (km) by road distance. In addition, the sensitivity of the estimated R ₀ to assumed period of infectiousness (7 days) was assessed by repeating the above estimation procedures, but assuming a period of infectiousness of 14 days.

SI modelling method

R ₀ was estimated using a method described by Stegeman et al. [14]. The transmission of HPAI subtype H5N1 between villages was estimated from the relationship between the number of newly infected villages each day (case villages, C _t) and the number of infectious villages each day (infectious villages, I _t), assuming an infectious period of 7 days. Assuming that all newly infected villages were infected by an infectious village during this phase of the epidemic, a simple deterministic SI model was used to describe the transmission of HPAI subtype H5N1 between villages. The number of susceptible villages was the total number of villages in Romanian counties affected during this phase of the epidemic (17 counties containing 7872 villages). All villages were assumed to have poultry and be at risk of an outbreak of HPAI subtype H5N1. Thus, we assumed that at the beginning of the epidemic, 7872 villages were at risk. We assumed that all outbreaks were detected and reported. Furthermore, we assumed that during the epidemic, the number of susceptible villages decreased only via the depopulation of infected village poultry populations.

We assumed [14] that the rate at which susceptible villages became infected (C) depended only on the proportions of susceptible (S/N) and infectious (I/N) villages in the population and the transmission coefficient, β. Thus, β was estimated for each epidemic day as (NC)/(SI), and R ₀ was estimated as βT. The period of infectiousness (T) was assumed to be 7 days.

Beginning 5 May 2007 (7 days prior to the first reported outbreak in this epidemic phase), the number of newly infected villages (C _t), the number of susceptible villages (S _t) and the number of infectious villages (I _t) were calculated for each epidemic day, and β was estimated as described above (Table 1).

Table 1. Estimates of the number of villages newly infected (C), the number of infectious villages (I), the number of susceptible villages (S) and the daily infection rate parameter (β) (N=7872)

Epidemic doubling time method

During the exponential phase of an epidemic (introduction of an infectious individual into a susceptible population), the number of secondary disease cases increase at an exponential rate. Each infection gives rise to R ₀ new infections per generation of infection, assuming a constant doubling time. The doubling time (t _d, the time period in which the number of outbreaks doubles) can be approximated as (ln2D)/(R ₀−1), so that an estimate of R ₀ is 1+[(ln2D)/t _d], where T is the duration (days) of infectiousness of an outbreak [2]. The average doubling time was calculated for all possible combinations occurring during the ascending phase of the epidemic (days 1–13). Assuming that T was 7 days, R ₀ was estimated as 1+(7/t _d) ln2.

Results

The epidemic of outbreaks (n=110) between 12 May and 6 June 2006, is shown in Figure 1. The estimated R ₀ using the nearest infectious neighbour (Euclidean and road distance), the epidemic doubling time, and the SI modelling methods were 2·14, 1·95, 2·21 and 2·68, respectively. Assuming a period of infectiousness of 14 days and using the nearest neighbour method based on Euclidean and road distances, the estimated R ₀ (2·10 and 1·91, respectively) was not substantially different.

Fig. 1. Epidemic curve of highly pathogenic avian influenza subtype H5N1 outbreaks (n=110) in Romania, 12 May–6 June 2006.

Using Euclidean distance to define infectious sources, the mean distance and mean time elapsing between source villages and outbreak villages was 23·36 km (95% CI 18·19–28·53) and 3·08 days (95% CI 2·72–3·45), respectively. Using road distance to define infectious sources, these parameters were estimated to be 33·37 km (95% CI 25·17–41·56) and 3·09 days (95% CI 2·74–3·45), respectively. The distance and time period between source villages and outbreak villages was not significantly (P>0·5) correlated using either of the methods to define proximity. However, both distances (0·8925) and time periods (0·8168) estimated with the two methods were significantly (P<0·001) correlated. The frequency distributions of R ₀ are shown in Figure 2; the distribution defined by road distance was more skewed than that defined by Euclidean distance (skewness statistics 3·07 and 1·30, respectively).

Fig. 2. Frequency distribution of estimated secondary outbreaks of highly pathogenic avian influenza subtype H5N1 outbreaks in village poultry populations, Romania, 12 May–6 June 2006, using either Euclidean distance (■) or road distance (□) to define the spatial relationship between outbreaks and an assumed period of infectiousness of 7 days.

Significant spatial autocorrelation was not detected for the source locations, whether these source locations were identified using Euclidean distance (I=0·0125, P=0·5730) or road distance (I=–0·0773, P=0·7610). No significant (P>0·05) local clusters of spatial autocorrelation were detected.

The product of the proportions of susceptible and infectious villages (SI/N) increased quickly between epidemic days 1 and 11 (Fig. 3), peaking on day 15 and sharply decreasing to day 24. The estimated daily transmission coefficients decreased (0·8723 per day) significantly (P=0·0011) throughout the epidemic (Fig. 4). The estimated mean epidemic doubling time (t _d) was 2·89 days (95% CI 2·42–3·35).

Fig. 3. The proportion of susceptible and infectious villages (SI/N) during an epidemic of highly pathogenic avian influenza subtype H5N1 in village poultry populations, Romania, 12 May–6 June 2006 (epidemic days 1–26), assuming a period of infectiousness of 7 days.

Fig. 4. Estimated daily transmission coefficient (β) during an epidemic of highly pathogenic avian influenza subtype H5N1 in village poultry populations, Romania, 12 May–6 June 2006 (epidemic days 1–26), assuming a period of infectiousness of 7 days.

DISCUSSION

Estimates of R ₀ using the three methods of nearest (Euclidean or road distance) infectious neighbour, SI modelling, and epidemic doubling time ranged from 1·91 to 2·68. Although the true R ₀ in this epidemic is unknown (because the actual source of each outbreak is unknown), a range of methods might be useful for characterizing epidemics of infectious diseases. These methods of indirectly estimating R ₀ have application in the common field situation in which tracing information is unavailable or incomplete. Once R ₀ has been estimated, better disease control strategies can be designed.

R ₀ depends on contact patterns between susceptible and infectious individuals (or groups of individuals, such as village poultry populations), and this parameter is specific for a given population during a given time period. Thus, determining which method of estimating R ₀ in field situations is the most accurate is virtually impossible because the source of each infection can rarely be verified. For this reason, a comparison of R ₀ values estimated from field data using different methods is useful to determine if such methods provide consistent results. Consistency can provide a greater level of confidence in the estimate made. In addition, if results are consistent then a method that is technically easy to perform, demands few restrictive assumptions, and can be applied to the imperfect data that is often the only source of information available in disease outbreaks, would be preferred. In the present study, the epidemic doubling method might the method of choice, if the preceding characteristics are considered important.

Several assumptions are implicit in the approaches used to estimate R ₀. Most importantly, in all of the methods used it was assumed that a population of village poultry was infectious for 7 days and that this period was constant during the epidemic. This latter assumption appears reasonable, given the very short during of the epidemic (26 days) compared to the assumed period of infectiousness (7 days) and the fact that the epidemic did not spread across jurisdictions.

The incubation period of avian influenza viruses in individual poultry is generally 3–5 days (but may be longer) (http://www.oie.int/eng/avian_influenza/A_Fiches_IA.pdf). Other investigators have assumed incubation periods of 2–6 days [15] and 1–4 days [16]. However, accurate estimates of the period of infectiousness at the aggregated (village) level are unavailable. Assuming that poultry are infectious during the incubation period, the period of infectiousness might range from 1 to 6 days. The village-level period of infectiousness is likely to be longer, because we assumed that several birds with clinical signs would need to be observed before the outbreak was reported, investigated, and control measures that prevented further village-to-village transmission from occurring were applied. Although this parameter is uncertain, the focus of this study was the comparison of estimates from the different methods used, rather than precise estimation of R ₀ using only one method. In addition, estimation of R ₀ (using the nearest neighbour method based on Euclidean and road distances), assuming a period of infectiousness of 14 days, resulted in no substantial differences in the estimated R ₀ (2·10 and 1·91, respectively). Other investigators have assumed the flock-level period of infectiousness, based on assumed incubation period and the delay between clinical signs, reporting and flock depopulation, to be 11·3–12·3 [15] and 8–13·8 days [17]. Overestimating the period of infectiousness will result in underestimating R ₀, and vice versa. Better estimates of the period of infectiousness, and particular information on factors that influence this parameter, would be useful in future studies of the epidemiology of HPAI subtype H5N1.

An assumption made using the nearest infectious neighbour method was that the nearest (based on either Euclidean or road network distance) neighbouring village within the infectious window period of 7 days was a source for each outbreak of interest. Thus, it was assumed that proximity in both time and space are the strongest factors driving the spread of a contagious disease such as HPAI. Clustering of disease events in time and space provides compelling evidence of disease spread, however the same patterns can arise from the presence of other risk factors that are spatially correlated [18]. It is thought that HPAI virus can be introduced into populations of domestic poultry from wild waterfowl [19]. Although the mechanism of spread between populations of domestic poultry in districts is unclear, it is likely to be via the movement of live birds: illegal movements of poultry and other avian species have been documented as contributing to H5N1 spread in several instances [20, 21]. Live bird markets have also been the presumed source of avian influenza outbreaks in domestic poultry [21–24]. We assumed that the closest infectious village to each outbreak served as the source of that outbreak. Information on the contact structure of villages, with respect to the movement of poultry, poultry products and potentially contaminated fomites, would improve the validity of this method of estimating R ₀. However, in most situations it is unlikely that such information (with sufficient spatial and temporal detail) can ever be collected. The assumption that each village was only infected by one other village might not be realistic, but the chances of a village being infected by more than one other village simultaneously would seem to be remote. If such infections occurred commonly in this epidemic, then we might have underestimated R ₀.

In the SI modelling method, we assumed that at all villages located in counties in which one or more outbreaks occurred were at risk of HPAI subtype H5N1 virus infection. Three of the counties included in this current study had previously reported a total of 11 outbreaks (between 5 December and 10 March), all in areas peripheral to the third epidemic phase. Removing villages located in these three counties would have reduced the initial susceptible population from 7872 to 7165, but the estimated R ₀ would have increased only from 2·212 to 2·213. Estimated R ₀ is insensitive to the number of susceptible villages in this model because the number of infectious villages at any point in time is only small (<1%) compared to the number of susceptible villages [14]. In addition, there were no recovered (and therefore resistant) villages included in the model (all village poultry populations in which HPAI subtype H5N1 was detected were depopulated). Finally, it was assumed that the number of new outbreaks on each day of the epidemic depended only on the density of infectious and susceptible villages [14], the ‘true mass action’ formation of an SIR model [25]. In general, heterogenous mixing increases the potential for a pathogen to invade and persist in the population, so that we might have underestimated R ₀ using the SI modelling approach.

The highest estimate of R ₀ in this study (2·68) was made using the epidemic doubling time method. We assumed that each outbreak produced R ₀ outbreaks per generation of infection (7 days), and the distribution had a constant doubling time. This approach also assumes that the proportion of resistant individuals in a population is negligible compared to the population at risk. It approximates well the initial phases of an epidemic. This case study involved the invasion of HPAI subtype H5N1 into the central part of Romania, where the disease had not previously been reported. It is possible that other subtypes of avian influenza (e.g. low pathogenic strains) could exist in this part of the country and provide some level of protection against HPAI subtype H5N1. Also, the epidemic curve (Fig. 1) shows a reduction in reported outbreaks between days 8 and 11, but then an increase on days 12 and 13. These fluctuations may be real, or could represent reporting delays. In the present study, doubling time was estimated for the first 13 days of the epidemic. If this period is restricted to the first 8 days (data not shown), the mean doubling time is 1·73 days and the estimated R ₀ would be 4·05. The doubling method relies on accurate reporting of field data. Whilst it is a robust method, relying on few assumptions, it might not be appropriate in all situations.

There are few estimates of the R ₀ of HPAI in poultry populations. Using a generalized linear model, Mannelli et al. [15] estimated between-farm transmission parameters for the 1999–2000 HPAI subtype H7N1 epidemic in northern Italy. Between-farm R ₀ for the regions of Lombardy and Veneto, during the initial phase of the outbreak, were estimated to be 1·8 and 1·5, respectively. These estimates assumed average periods of infectiousness of 12·3 and 11·3 days, respectively [15]. If a period of infectiousness of 7 days is assumed (as for the present study), these estimates would be 1·05 and 0·91, respectively. The much lower R ₀ estimates, compared to the present study, might reflect underreporting of case flocks or could be a real effect as a result of lower contact rates between commercial poultry farms vs. villages in our study. In a study of an epidemic of HPAI subtype H7N7, occurring in The Netherlands in 2003 [17], the between-flock R ₀ before detection of the first outbreak was estimated for two affected regions (Gelderse Vallei and Limburg) to be 6·5 and 3·1, respectively (assuming periods of infectiousness of 13·8 and 8 days, respectively). If a period of infectiousness of 7 days is assumed, then these estimates would be 3·29 and 2·73, respectively. Thus, the estimates of R ₀ (1·95–2·68) in the present study lie in between those made from the 1999–2000 Italian epidemic of H7N1 (0·91–1·05) and those made from the 2003 Dutch epidemic of H7N7 (2·73–3·29). An important consideration when interpreting R ₀ estimates for between-flock transmission of HPAI is the relative density of the populations studied [17].

Surprisingly, few estimates of the R ₀ for pandemic influenza in human populations have been made. Massad et al. [26] examined the 1918 pandemic H1N1 outbreak in the city of Sao Paulo. Using a mathematical model, R ₀ was estimated to be 2·68. Mills et al. [27] estimated an R ₀ value of 2–3 for the 1918 influenza epidemic by fitting a deterministic SEIR model to pneumonia and influenza death curves from 45 US cities. R ₀ values ranging from 1·6 to 2·4 were considered in stochastic simulation models used to investigate pandemic influenza spread in rural Southeast Asia [28] and the US population [29]. These estimated or assumed values for R ₀ are consistent with the R ₀ values (1·95–2·68) estimated in the present study. Thus, the analysis of epidemiological data from outbreaks of HPAI subtype H5N1 in poultry might further assist with the development of control programmes for pandemic influenza.

Efforts to prevent between-flock spread, such a quarantine, depopulation and disinfection, if effective (49–63%, based on R ₀ estimates in this study and critical threshold theory), would provide a high pay-off in disease control programmes. Our understanding of the epidemiology of HPAI subtype H5N1 will be increased by further studies utilizing epidemic data, and thus will greatly assist in the design of better disease control programmes.