Introduction

Highly pathogenic avian influenza virus (HPAI) H5N1 has led to the death of more than 200 million domestic birds during the period 1999–2004 alone (Capua and Alexander, 2004) and has been regularly detected in various species of wild waterbirds (Alexander, 2006). Although virus transmission from wild to domestic birds (and vice versa) is a major concern, very little is known about the prevalence, behavior, and spread of H5N1 and other avian influenza (AI) strains in wild waterbirds and in the poultry population. There is a need for research on the potentially high threat that the H5N1 virus might pose to domestic poultry or wild waterbirds. The study of HPAI H5N1 cases in wild waterbirds is an important step to improve our understanding of the spread in wild waterbird populations and the potential introduction of HPAI virus into the domestic poultry sector.

During the winter of 2005–2006, the HPAI H5N1 virus spread from Asia throughout Europe (Alexander, 2007) and reached Lake Constance, which is the most important inland wintering area for waterbirds in Europe (Dalessi et al., 2007). From January to April 2006, the H5N1 virus was found and confirmed in 82 dead wild waterbirds (Hofmann et al., 2008). The course of the disease in wild waterbirds is poorly understood but the presence of H5N1 in wild waterbirds in areas, such as Lake Constance, where no cases in domestic poultry have been reported suggests that the virus can be spread by wild waterbirds before they potentially succumb to the disease. Consequently, constant surveillance is essential to better assess the prevalence of HPAI in wild waterbirds and also in domestic poultry and thus lead to a better understanding of the risks posed by H5N1 and other HPAI to wider communities.

The various AI surveillance activities and research programs undertaken worldwide demonstrate that the prevalence of the virus in wild waterbird populations depends heavily on the host species, age group, geographic location, and seasonality (Olsen et al., 2006). A 7-year survey in eastern Alaska reported a prevalence of various AI strains (H3, H4, and H6) of 0.06% among waterbirds (Anatidae) and shorebirds (Chlaradriidae and Scolopacidae) (Winker et al., 2007). However, a similar study performed mostly among mallards (Anas platyrhynchos) in northern Europe found an overall mean prevalence of 12.5% with a high seasonal variability. In fact, this estimated prevalence ranges from a mean of 15% in autumn to 4% in spring (Wallensten et al., 2007). These different findings demonstrate a need for research concerning the distribution of AI virus in the diverse areas where wild waterbirds breed, molt, and rest during migration as well as during winter. We recognize that the above-mentioned low pathogenic avian influenza (LPAI) strains are a different disease entity compared with the highly pathogenic avian influenza (HPAI) presented in this paper and should not be compared directly. Better assessment of the wild waterbird population dynamic and the transmission pathways of the virus is required. The data of the H5N1 outbreak around Lake Constance during the winter of 2005–2006 and 12 years of waterbirds census data provide an ideal opportunity for such an assessment.

In this paper, we examine the transmission dynamics of HPAI H5N1 in wild waterbirds and estimate the basic reproductive ratio (R 0) for this outbreak by defining the most appropriate model for the dynamics of the susceptible wild waterbird population at Lake Constance and by investigating different transmission models based on the H5N1 cases observed in 2005–2006.

Materials and Methods

Available census data of birds during migration and wintering and data on the number of dead waterbirds found at Lake Constance in 2006 that tested positively for H5N1 are used to determine the most appropriate population dynamic and transmission model, respectively.

Waterbird Population at Lake Constance

Between September and April of every winter, a local ornithological working group (Ornithologische Arbeitsgemeinschaft Bodensee) performs monthly censuses of the number and species of waterbirds, including both resident and migrating waterbirds, on Lake Constance. These observations represent the best available data on local avian populations. In our model, we use the survey results for 12 years from 1995 to 2006 to determine the net population change at Lake Constance. The analysis is based on a subset of 12 species, which were selected based on their overall abundance, their susceptibility to H5N1, and their potential contact with domestic poultry. The species concerned include both migrant and waterbirds resident at the lake: Great Crested Grebe (Podiceps cristatus), Cormorant (Phalacrocorax carbo), Mute Swan (Cygnus olor), Whooper Swan (Cygnus cygnus), Gadwall (Anas strepera), Teal (Anas crecca), Mallard (Anas platyrhynchos), Red-crested Pochard (Aythya nyroca), Pochard (Aythya ferina), Tufted Duck (Aythya fuligula), Goosander (Mergus merganser), and Coot (Fulica atra). We use the total number of individuals in this group of species to model the total susceptible population. (We did not consider the dynamics of each species separately because we could not develop a multi-compartmental model for HPAI by host species because the risk of cross-species transmission is unknown).

HPAI H5N1 Cases

The winter of 2005–2006 saw a total of 82 confirmed cases of H5N1 in dead wild waterbirds found at Lake Constance. These positive cases included Little Grebe (Tachybaptus ruficollis), Common Pochard (Aythya ferina), Tufted Duck (Aythya fuligula), and some birds only identified as “duck” (Rutz et al., 2007) (Happold et al., 2008). In subsequent winters, no cases were detected despite more intensive and active surveillance and monitoring of both dead and live birds (BVET, 2008).

SIR Models

Various systems of coupled differential equations are examined. Each equation is based on the compartmental SIR models considering wild birds divided into three different groups, defined as S = susceptible, I = infectious, and R = removed (deceased, recovered). Each equation also accounts for both the population dynamics of the wintering waterbirds and the transmission of HPAI.

There are various simplifying assumptions to this modeling approach. These simplifications are pooling of all wild waterbird species as potentially susceptible to HPAI, modeling total number of waterbirds, as opposed to densities, assuming homogeneous distribution and contact of these populations over the study area and assuming density dependent transmission so that potential spread of H5N1 increases with density as opposed to frequency dependent transmission. Furthermore, environmental effects on yearly population dynamics could not be considered.

Population Dynamics

To determine the dynamics of population growth and migration of waterbirds at Lake Constance, four simple models (A–D) are considered. In these models, the population dynamics are represented by differential equations and take the general form,

$$ \frac{dS}{dt} = \Upomega (S,t) + \chi_{1} \sin \left( {\chi_{2} \Upgamma t + \chi_{3} } \right) $$
(1)

S(t) (birds) represents the total number of waterbirds at Lake Constance at time t (weeks) and represents the susceptible compartment for the transmission models described in the following section. The function \( \Upomega (S,t) \) takes the form of one of the four population models (A–D) under investigation. The remainder of equation (Akaike, 1974) is a sinusoidal function to account for the cyclic-seasonal dynamics, growth, and migration of S(t). The constant \( \Upgamma \) is given by

$$ \Upgamma = {\frac{2\pi }{365/7}}, $$
(2)

indicating a yearly cycle of (365/7) weeks per year. The parameters \( \chi_{1} \) (birds/weeks), \( \chi_{2} \) (per year), and \( \chi_{3} \) (dimensionless) are fit to the survey population data. The model given by (Akaike, 1974) represents both the yearly net growth (births and deaths) as well as the migration dynamics of the 12 considered species, which represent the susceptible compartment of the SIR model.

Four separate models (A–D) are considered, each with slightly different contributions to the net change in total population. Model A describes the rate of change of S(t) to be dependent on the sinusoidal function only,

$$ \Upomega (S,t) = \, 0. $$
(3)

Model B includes a linear growth term, model C an exponential growth term, and model D a logistic growth term. The corresponding analytic solutions for model A is given by:

$$ S(t) = S_{0} - {\frac{{\chi_{1} }}{{\chi_{2} \Upgamma }}}\left( {\cos \left( {\chi_{2} \Upgamma t + \chi_{3} } \right) - \cos \left( {\chi_{3} t} \right)} \right). $$
(4)

The population models A–D are fit to survey count data with t = 0 being the first survey data point in mid September 1995. In each case the number of fitting parameters depends on the model choice. The fit for each model to observed counts of observed waterbirds is performed using regression analysis, minimizing sums of squares by assuming both normally distributed and log-normally distributed errors. For validation, both the analytic solutions and the differential equations of models A–D are fitted and S(t) set to always be greater than zero, ensuring physically realistic results. The model producing the best fit was determined using the Akaike’s Information Criterion (AIC), a measure of the goodness of fit of an estimated statistical model (Akaike, 1974).

Transmission Model

Given the best model for the population dynamics, two additional model compartments were added to account for the transmission of H5N1. I(t) and R(t) were introduced as the number of infectious waterbirds at time t and the cumulative number of confirmed H5N1 waterbirds found at the Lake, respectively. The following general SIR model is proposed:

$$ \frac{dS}{dt} = \chi_{1} \sin \left( {\chi_{2} \Upgamma t + \chi_{3} } \right) - \beta SI, $$
(5)
$$ \frac{dI}{dt} = \beta SI - \delta I - \mu I + f_{1} \left( {t,I} \right), $$
(6)
$$ \frac{dR}{dt} = f_{2} \left( t \right)\mu I, $$
(7)

S (birds), I (birds), and R (birds) represent the susceptible, infectious, and dead/removed waterbirds, respectively. β (birds−1 weeks−1) is a simple contact rate between S and I, 1/δ (weeks) is the infectious period before recovery, 1/μ (weeks) infectious period before death from H5N1, and the function Γ is given by equation (2). The model equations given by f1(t,I) and f2(t) for the three transmission models are summarized in Table 1 and reflect assumptions for the three transmission models (model I–III) considered. Model I assumes that all dead (H5N1-positive) waterbirds were found, model II accounts for only a proportion of the dead waterbirds being found, and model III includes a migration term for the infectious compartment. Our initial time, t*, is 1 week before the first H5N1 case was found on January 30, 2005.

Table 1 Model Equations for f 1(t,I) and f 2(t) for the SIR Model
Full size table

The three models (I–III) describing the transmission of HPAI are fit by comparing changes in R(t) for each week to the known number of dead, HPAI-infected waterbirds found in each week from January to April 2006. We use Matlab for solving the differential equations (ode package in Matlab) and for fitting the model to observed data. We use regression analysis assuming normally distributed errors (minimizing sums of squares) and Poisson distributed errors (maximizing Poisson log-likelihood). We check results obtained with Poisson distributed errors because the observed data are count values. The parameters fit in each model I–III are β, μ, and δ. Additionally, we fit for parameter γ in model II and parameter χ4 in model III. An additional fitting parameter I 0 is required when we investigate, not assuming the initial infective population to one. The basic reproductive ratio (R 0) (parameter that indicates the average number of secondary infections from one primary case throughout its infectious period) (Heffernan et al., 2005), the simple contact rate β between susceptible and infectious, and the peak incidence are calculated for the different models. R 0 is therefore a value that indicates whether an infection will spread or die out: if R 0 > 1, the disease can spread; if R 0 < 1, the transmission will fade out (Dietz, 1993). For reference, the derived R 0 for models I and II is given by

$$ R_{0} = {\frac{{S_{0}^{*} \beta }}{{\left( {\delta + \mu } \right)}}}, $$

and for model III the derived R 0 is

$$ R_{0} = {\frac{{S_{0}^{*} \beta }}{{\left( {\delta + \mu + {\frac{{\chi_{4} }}{{I_{0} }}}\sin \left( {\chi_{2} \Upgamma t^{*} + \chi_{3} } \right)} \right)}}}, $$

where S 0* denotes initial susceptible population level.

Results

Population Model

Model D has the smallest residual sum of squares (RSS), i.e., the smaller dispersion, but Akaike’s Information Criterion (AIC) indicated that model A produces a better fit (Table 2). However, there is minimal difference between models A–D. Additionally, the F-statistic comparing model A to models B, C, and D indicates that additional parameters to model A do not significantly improve the fit (not shown).

Table 2 Residual Sums of Squares (RSS) and the Akaike’s Information Criterion (AIC) for Models A–D
Full size table

Considering these results, we accept model A as appropriate to represent the population dynamics of the susceptible compartment in subsequent transmission models. The results of model A, with the survey data points for comparison, are shown in Fig. 1. Our results indicate that seasonal population changes due to the migration of waterbirds are apparent, and a general trend analysis of the data indicates a slight linear increase of population size.

Figure 1
figure 1

Plot of the results of the population model A (line) and the total counts (1995–2006) (star) of the wild bird species present at Lake Constance.

Full size image

Transmission Model

The three transmission models (I–II–III) are fitted assuming the best appropriate population dynamics model: model A. The fitted contact rates \( \beta \) for models I–III are given in Table 3, along with the 95% confidence intervals, the resulting reproductive ratio R 0, an estimate of the length of infection (1/(δ + μ)) (in days), and the peak incidence of infectious per 100,000 susceptible waterbirds.

Table 3 Fitted and Calculated Parameters (with 95% Confidence Intervals) for the Transmission Models I–III, Contact rate β, Length of Infection, Basic Reproductive Ratio (R 0), and Peak Incidence
Full size table

By assuming that the initial infective population (I 0) has a value of 1 or treating it as a fitted parameter, the reproductive ratio for model I is estimated at approximately 1.6 and the length of infection at approximately 3 days. When I 0 = 1, the peak incidence is approximately 90 infectious per 100,000 susceptible waterbirds, but in the case of fitting I 0, the incidence differs between assuming normal and Poisson distributed errors from 55 to 16 per 100,000 susceptible waterbirds, respectively. This is due to changes in the resulting fitted values of I 0 and other parameters (results not shown). For model II, by setting I 0 to be unity and varying γ from 5 to 100%, the reproductive ratio ranged between 1.62 and 1.68 (Table 3). The basic reproductive ratio for model III is estimated to be approximately 1.5 or 1.6, slightly lower than predicted by the previous models. The peak incidence of infectious per 100,000 susceptible waterbirds is 47 and 175, depending on the assumption of normal or Poisson errors, respectively, and the length of the infectious period is estimated at 2–3 days.

As an example only, Figs. 2 and 3 (fitted dead found and model results of dead found with infectious population, respectively, with normal and Poisson errors) illustrate the plotted results for model II when only 5% of H5N1 cases are assumed to be actually found, which is more realistic than model I where all infected individuals are collected. Further plotted results for model II assuming other proportions of H5N1 cases are actually found are not included here, but summarized values for peak incidence and R 0 are given in Table 3. As can be assessed from these figures (and other plots not shown), the epidemic is self-limiting because I(t) dies out as the number of susceptible waterbirds becomes lower. Essentially, the epidemic dies out coinciding with outward migration of the wintering waterbirds (as predicted by the population model) causing a decrease in transmission chance among the remaining waterbirds. Analytically, we can show that the number of infectious birds eventually decreases and the epidemic dies out, when the rate of change for I (Eq. 6) is <0. Namely, when the number of susceptible is less than

$$ S(t) = {\frac{{\left( {\delta + \mu - {{f_{1} (t,I)} \mathord{\left/ {\vphantom {{f_{1} (t,I)} I}} \right. \kern-\nulldelimiterspace} I}} \right)}}{\beta }}, $$

which is true as the population dynamics cause S(t) to decrease. Figure 3 shows that the number of infectious do decrease when S(t) falls below this threshold and consequently the epidemic dies out (time for S(t) is below threshold is marked by X on the infectious curves for model II).

Figure 2
figure 2

Plot of predicted R*(t) with observed data for model II (normal and Poisson errors) when I 0 is unity and γ = 5%.

Full size image
Figure 3
figure 3

Plot of predicted I(t) and R*(t) for model II (normal and poisson errors) when I 0 is unity and γ = 5%.

Full size image

We are unable to determine the best fitting transmission model when fit using regression assuming Poisson distributed errors. However, we can compare the RSS and AIC for models fit assuming normally distributed errors. In summary, model III (which includes out migration of infective) appears to have the lowest RSS but the AIC along with RSS indicates that model I (no out migration and assuming all dead birds are found) may be the best fitting model when we assume the initial number of infected is one. Model II, which includes no out migration of infective, appears to be the next best fitting model when we assume only 20% of the dead waterbirds are recovered.

Discussion

This study is the first to attempt quantifying the transmission of HPAI H5N1 virus at Lake Constance during the outbreak of winter 2005–2006 and is a step toward the improvement of the knowledge of virus transmission among wild waterbirds. The population dynamics model of wintering waterbirds (resident and migrating), consisting of a sinusoidal function to describe net growth and migration, was chosen to represent the susceptible compartment. Three transmission models, developed to consider transmission dynamics during this outbreak, predicted a basic reproduction ratio (R 0) with values of approximately 1.6. Given that the basic reproduction ratio found in this work is greater than one, the outbreak can be described as a small epidemic (MacDonald, 1957). This work has shown that HPAI in 2005–2006 was of epidemics levels in wild waterbirds but that the spreading of the virus was self-limited due to the migration of susceptible wild waterbirds at the end of the winter, thus reducing the chance of transmission by any of the remaining infectious waterbirds. Wild birds at Lake Constance migrate between Eastern Europe, Russia, and the Balkans as well as the Mediterranean and West Africa depending on the species. Lake Constance is hence connected to the Euro-Asian interface of wild bird migration, which could be a future source of reintroduction.

Because of the unknown waterbird population size and the difficulty of collecting all positive cases, transmission dynamics models are rare in wild waterbird populations. In contrary, transmission analyses have regularly been performed in the context of HPAI outbreaks in the domestic poultry sector (Garske et al., 2007). Depending on the phase of the outbreak, the basic reproductive ratio of the H7N1 outbreak in Italy in 1999–2000 ranges between 0.6 and 1.8 (Mannelli et al., 2007). The calculation of this parameter allows one to assess the quality of the spreading and of the remedial measures. In the 2003 outbreak in the Netherlands, the basic reproduction rate depended on local farm density (Boender et al., 2007) and the effectiveness of the control measures could be estimated by observing the resulting decrease of R 0 (Stegeman et al., 2004).

Simplifications and assumptions in our models were necessary given the limited understanding of transmission, but primarily given the lack of data to validate more complex models. These simplifications included the pooling of all wild waterbird species, thus not accounting for individual species dynamics, and assuming homogenous distribution of both susceptible and infectious waterbirds even though the behavior of different species with regards to feeding, preferred habitat, and migration strategies varies strongly. Additionally, an explicit model for the arrival of infectious waterbirds was not included, but instead the initial concentration was either set as fixed with a value of unity or was allowed to be a fitting parameter.

To improve the model, it would be important to include environmental factors, such as temperature at the lake, temperature of the water, ice cover, or water level, which could strongly influence the contact rates between infectious and susceptible, the population migration, or the arrival of infectious waterbirds. The inclusion of environmental factors would allow assessment and prediction of outbreaks, especially when comparing the outbreak in the winter of 2005–2006 (Rutz et al., 2007) to the absence of an outbreak in the following winter (BVET, 2008). The behavior, seasonality, ecology, and contact rate of wild waterbirds varies strongly between species. Therefore, the frequently encountered registration of HPAI cases as “ducks” or “wild duck” only is not sufficient to elaborate complex transmission models (Yasue et al., 2006). Waterbirds often live as mixed-species flocks. Some of these species could have different infectiousness or spreading potential and therefore might not play the same role in the transmission of the virus (Stallknecht and Shane, 1988). An ideal epidemic model could therefore incorporate heterogeneous contacts between individuals (Volz and Meyers, 2007) but would necessitate good information on the waterbird populations at the outbreak site. It is important to highlight that any improvements to the models of this study would add more realism and more complexity but also will increase the number of parameters to be estimated. Without more data this work would unfortunately be more uninformative.

This work is not only driven by the aspiration to understand transmission dynamics in wild waterbirds but also to enable one to quantify the potential risk of the virus spreading to domestic poultry and its subsequent high economic impact in the private and public sector (Capua and Alexander, 2004). The results of the studies of Munster et al. (2005) and Campitelli et al. (Campitelli et al., 2004) highlight the potential transmission of AI virus between wild waterbirds and domestic poultry. By combining the transmission dynamics information estimated in this study with the still-unknown contact rate between wild waterbird and domestic poultry, it would be possible to identify adequate surveillance needs, depending on the different contact rates near and away from waterbodies. Risk-based surveillance could therefore be improved and resources be allocated more effectively and efficiently to protect not only the wildlife population but also the health of livestock and consumers (Stark et al., 2006).