## Abstract

Recent avian flu epidemics (A/H5N1) in Southeast Asia and case reports from around the world have led to fears of a human pandemic. Control of these outbreaks in birds would probably lead to reduced transmission of the avian virus to humans. This study presents a mathematical model based on stochastic farm-to-farm transmission that incorporates flock size and spatial contacts to evaluate the impact of control strategies. Fit to data from the recent epidemic in the Netherlands, we evaluate the efficacy of control strategies and forecast avian influenza dynamics. Our results identify high-risk areas of spread by mapping of the farm level reproductive number. Results suggest that an immediate depopulation of infected flocks following an accurate and quick diagnosis would have a greater impact than simply depopulating surrounding flocks. Understanding the relative importance of different control measures is essential for response planning.

## 1. Introduction

Highly pathogenic avian influenza (HPAI) outbreaks have been reported worldwide since surveillance started (Alexander 2000). Virus transmission to humans has been reported, leading to fears of a new human pandemic stemming from reports of outbreaks of avian influenza around the world (Saito *et al*. 2001; Gilbert *et al*. 2006; Wee *et al*. 2006). In general, avian influenza viruses do not replicate efficiently in humans, indicating that efficient direct transmission of avian viruses to humans would be an event as rare as human-to-human transmission (Beare & Webster 1991). Despite this, a dramatic increase in the number of reports of avian virus transmission to humans has been observed since 1997 (Webby & Webster 2003). Moreover, it appears that the number of human infections follow a similar trend, at a smaller scale, as the avian epidemic (Koopmans *et al*. 2004) suggesting that human infections would be reduced when an avian epidemic is quickly controlled. This concept highlights the importance of understanding the control strategies, which would have the greatest impact in diminishing the epidemic among poultry. Should policy makers focus on a quick destruction of infected flocks? Should poultry around infected cases be depopulated? How large should the depopulation radius be?

Mathematical models provide a framework to address these questions and determine the impact of different control strategies. Although mathematical models of pandemic influenza among humans have been developed to evaluate the impact of measures such as vaccination or antiviral drugs on controlling virus diffusion (Ferguson *et al*. 2005; Longini *et al*. 2005), this approach has been rarely applied to avian influenza. To date, little detailed data about avian influenza epidemics have been published with the exception of the Dutch outbreak in 2003 (European commission 2003). To evaluate the overall effectiveness of control measures, Stegeman *et al*. (2004) estimated that the end of the epidemic was due to extensive culling of the susceptible population rather than to other specific control measures from the retrospective epidemic data. Questions concerning control strategies which would have the greatest impact in reducing epidemic magnitude remain unanswered. To our knowledge, no other mathematical models have been developed to explore avian influenza dynamics in poultry and candidate control measures.

The aim of this paper is to identify the strategies which would have the greatest impact on controlling avian influenza epidemics. Exposure to a virus has a spatial component that influences the spatial and temporal spread (Lemenach *et al*. 2005). We developed a mathematical model, calibrated on the Dutch outbreak, taking into account both spatial component and farm size. We weighed the impact of different control strategies on the forecast cumulative epidemic and performed a multivariate sensitivity analysis.

## 2. Methods and data

### (a) Model: formulation and assumptions

A spatial farm-based model was developed to mimic the transmission dynamics of avian influenza and the implementation of control strategies. We assumed the farm to be an individual unit because of the rapid transmission of the virus and hence animals would be kept inside holdings if an epidemic occurred (Stegeman *et al*. 2004). The model includes two components: a stochastic one to assess the spread among farms and a deterministic one describing the changes in the farm disease status (figure 1).

Firstly, each farm was classified as susceptible (S). The infection process was modelled stochastically using a discrete-time model formulation with a one-day-time step using Monte Carlo sampling. On day *t*, the probability that a susceptible farm *i* becomes infected is given byThis hazard-based probability depends on the farm sizes (*N*_{i}, *N*_{j}) and on the distances between farms, which are taken into account through the contact rates. The contact rates are defined as follows: is the contact rate between farms within short-range distance (aerosol dispersion, movements of poultry over pasture within the protection area), within medium range distance (contacts with neighbours and movement of workers among farms with the surveillance area) and within long-range distance (truck transport of living poultry, eggs and manure). The use of the natural log allows for large differences in the number of birds per farm to be taken into account. This scale was considered most appropriate as differences in biosecurity risks between small farms and large ones are maintained but the differences between similar industrial farms with different sizes are smoothed. Assuming a log-scaled transmission did not impose specific mathematical constrains on likelihood optimization, since the optimization as explained below was performed numerically. Thus, this approach could easily be derived to other nonlinear scaling. We also assumed the risk of being infected to be equal within each of the three transmission buffers.

Contact rate values were used at two time lags: before the virus was known to circulate (*k*=1, pre-control period: *T*_{1}) and afterwards when control strategies, such as a national ban (i.e. all exports of live poultry and eggs were blocked and movement of live poultry was banned in the Netherlands), were put into place (*k*=2, control period: *T*_{2}). For the ease of computation, we included a constant *c* to ensure estimated parameters were within reasonable bounds.

Once a farm was infected, the farm was classified as latent (L, infected but not contagious); after *δ* days it was classified as not reported (NR, infectious but not reported to authorities because of poor clinical expression). After *γ* days, the farm was classified as reported (R, infectious with explicit clinical signs reported to authorities); after *τ*_{t} days the farm was depopulated (D, removed from transmission chain), where *τ*_{t} decreases with time (table 1).

In addition to their disease status, farms were classified according to their control status. Farms were classified as: ‘infected premises’ (IPs) for reported farms targeted for depopulation; ‘ring-culled premises’ (RCPs) for farms in the vicinity (within a certain radius size, *rpolicy*) of IPs and targeted for the neighbourhood control policies, such as depopulation or vaccination; or ‘dangerous contacts’ (DCs) for farms not concerned by neighbourhood control policies but epidemiologically in contact with IPs. We let *ρ*_{t} denote the ratio of pre-emptive depopulated farms (RCPs and DCs) to infected farms (IPs) increasing with time. We simulated a set of DC farms, as this information was not available for the Dutch epidemic. We assumed the probability of a farm being classified as a DC to be slightly different from the risk of being infected and to decrease with distance. Since the DC status is established while investigating previous contacts of an infected farm, it does not depend as strongly on key radius zones as infected holdings. For each farm, a set of DCs was sampled using a function 1/distance.

### (b) Parameter estimation and initial values

To obtain values for the spatial contact rates, (*β*_{short, k}, *β*_{medium, k} and *β*_{long, k}), we used maximum likelihood estimation (MLE) at two times: the pre-control and control periods. The likelihood of an individual farm *i* over the observed epidemic length *T* was calculated as the probability of remaining susceptible on day *t* multiplied by the probability of being infected the day infection was noted according to the observed data (*tr*_{i}). In other words, . If the farm remained susceptible over the entire epidemic duration, then *tr*_{i}=*T*+1. If the farm is pre-emptively depopulated, *tr*_{i}=date of depopulation+1. The individual probabilities are multiplied together to achieve a single likelihood value.

The Simplex method was used to optimize the likelihood values against observed available data. To prevent the occurrence of local minima, 1000 sets of six initial parameters were generated and the resulting optimal set was used for simulations.

All other parameters were taken from the published literature, experimental studies or previous observed epidemics (table 1). The radius distances corresponding to the three transmission buffers (less than 3 km, between 3 and 10 km and more than 10 km) and to the zone concerned by the control policy were chosen according to what is usually applied in Europe when an infected farm (Council of the European Communities 1992). The same radius distances were employed during the Dutch outbreak.

### (c) Data

The model was fit to the 2003 A/H7N7 Dutch outbreak. We focused on the two most affected areas, Gelderland and Utrecht, where 185 commercials farms were infected starting from February 28 2003.

For each of the 156 municipalities in Gelderland and Utrecht in 2002, we had the number of farms (*n*=1136 in both regions) and the number of animals for any avian species (broiler 15.4%, layers hens 75.4%, ducks 2.2%, turkey 4.8% and other 2.2%; Statistics in Netherlands 2004). As data were not available on exact location of farms, the geographical coordinates of each farm and the number of animals within a farm were simulated at the municipal level. We first built municipal borders around its gravity centres using Voronoi tessellation, (Duyckaerts & Godefroy 2000) and then farm locations were randomly distributed within each municipality (average surface area of 40.8 km^{2}, figure 2*b*). Suggested by data (figure 2*a*), the number of animals in each farm was assumed to follow a right-skewed distribution, i.e. a Poisson distribution parameterized byData on infected farms were collected through the European Commission reports published throughout the epidemic (European commission 2003). The daily incidence of infected farms and the delay between report of infection and depopulation for each of the 185 farms were provided. The dates for infection, suspicion, laboratory confirmation and depopulation were indicated as well. Data also indicate to which town, municipality and province the infected farm belongs, type and number of animals within the farm and the source of infection. The control measures were fully described as well: a national ban on poultry movement the day the virus is known to be circulating; date of implementation of the surveillance and control buffer areas; date at which pre-emptive depopulation starts (within a 1 km radius from the beginning of March 2003 and then 3 km radius from the beginning of April 2003); and date at which depopulation of dangerous contacts starts.

### (d) Model calibration

We simulated 1000 epidemics over a period of 100 days beginning with a unique index case located in Renswoude municipality where the epidemic began. The results of model simulations were then compared, by calculating the correlation coefficient, to the 5-day moving average incidence of reported cases (figure 3*a*). Simulations also provide spatial-kriged risk maps for the pre-control and control periods evaluating *R*_{t} values, that is the average number of secondary infected flocks by one primary infected flock (Anderson & May 1991). For each farm, *R*_{t} values were estimated by simulating infection of an initial flock and then counting the average number of secondary cases generated for the length of the infectious period. We assumed the infectious period was 10 days during the pre-control period and 6 days during the control period.

### (e) Multivariate sensitivity and uncertainty analysis

We performed a multivariate sensitivity analysis to assess the impact of control strategies and potential delays in diagnosis, reporting and implementation of control measures: (i) delay between onset of infectiousness and report to authorities *γ*, (ii) delay between report and depopulation *τ*, (iii) ratio of pre-emptive culled farms to infected farms *ρ*, (iv) delay before the virus is confirmed in the country and first control measures are put into place *T*_{1}, and (v) radius of neighbourhood culling policy (*rpolicy*). These five parameters, constant over the entire duration of the epidemic, were sampled using a Latin hypercube sampling scheme (Blower & Dowlatabadi 1994) based on the assumption of statistical independence between input parameters. The distribution of infections in flocks was computed using 500 sets of parameters drawn from uniform continuous or discrete distributions, with bounds described in table 1. The partial rank correlation coefficients (PRCC) were also calculated using the 500 values for each parameter and the 500 mean cumulative number of predicted cases over the 1000 simulations (Blower & Dowlatabadi 1994).

As the Simplex method provides only point estimation values, an uncertainty analysis was performed for the contact rates using the same Latin hypercube sampling scheme, explained above. The effect of variation in contact rates values on the precision of the simulated cumulative number of infected farms was explored.

Finally, the sensitivity of the model to distances used for transmission buffers was tested by calculating the likelihood for several couple of values around the historical values (3 and 10 km) given in the literature.

## 3. Results

### (a) Parameter estimation

Before control strategies were put to place, estimated contact rates exhibited, as expected, greater values at short distances (table 1). During the pre-control period, similar contact levels were noted for virus transmission within the protection (less than 3 km) and surveillance area (less than 10 km), but long-range contact were five times lower.

Once control strategies were implemented, the values of all contact rates decreased. The largest decreases were for long-range contacts (75% decrease beyond the surveillance area) and for short-range contacts (more than 50% in the protection area). Reduction of medium range contacts was smaller (less than 50% decrease in the surveillance area (table 1).

### (b) Model fitting and spatial risk maps

Our simulation results are in agreement with the observed dynamics of 2003 (correlation coefficient=0.94, figure 3*a*). In the 2003 Dutch epidemic, 185 commercial farms were reported as infected and the 5-day moving average of reported cases exhibited a peak of 6.8 at day 27. The number of simulated cumulative cases was of 184, 95% empirical confidence interval (34; 294) with a peak of 7.1 reported cases at day 29.

Simulated cases were located mainly in the Gelderse valley as was the case in 2003 (figure 3*b*). The most affected municipalities in the simulations were Ede and Barneveld, as was also observed during the Dutch epidemic.

Estimation of the reproductive number identifies and quantifies high risk regions for virus spread. During the pre-control period, in the Gelderse valley, *R*_{t} values were above 4 on an average and in the north or west of Gelderland *R*_{t} values were between 2 and 5 (figure 4*a*). Over all regions, *R*_{t} was estimated to be 5.2 (4; 6.9) during the pre-control period. During the control period, estimated *R*_{t} values strongly decreased but remained above one in the Gelderse valley indicating the control measures alone were inadequate to halt transmission, and so the virus may have continued spreading (figure 4*b*). On an average, *R*_{t} was estimated to be 1.5 during the control period (1; 2.5); although the measures were inadequate to halt the spread, they may have limited the size of the final epidemic.

### (c) Sensitivity and uncertainty analysis

The results of the multivariate sensitivity analysis on control parameters demonstrated that the most influential parameters were first, the delay between detection and depopulation (*τ*, PRCC=0.85) followed by the time needed to detect an infected farm when poor clinical signs are noticed (*γ*, PRCC=0.73). The greater the delay between diagnosis, reporting and implementation of depopulation, larger is the epidemic (figure 5*a*). The duration of the pre-control period would play a minor role (*T*_{1}, PRCC=0.55) but still, increasing the delay between virus importation and disease confirmation may lead to larger epidemics. Increasing the ratio of pre-emptive depopulated farms to infected ones (*ρ*, PRCC=−0.50) or the neighbourhood control policy radius (*rpolicy*, PRCC=−0.62) would both lead to fewer cases, but with less impact. The monotonic relationships between the outcome and the input variables, assumed by the sensitivity analysis, have been verified (figure 5*a*).

Considering the range of spatial contact rates explored in the uncertainty analysis, the average size of epidemics was 233 cases, the 0.25 and 0.75 quantiles were equal to 188 and 277, respectively (figure 5*b*) and the standard-deviation was equal to 60.

The model results depend on the radius lengths that define short, medium and long transmission. The lengths used here correspond to the radius lengths implemented during the Dutch epidemic, so they have a practical interpretation. By allowing the lengths to vary, we found that a best optimization was obtained with radius equal to 1.6 km for the protection area and to 10.9 km for the surveillance area, respectively. Nevertheless as the likelihood obtained for this couple (likelihood=1015) was close to the one obtained for the couple 3 and 10 km (likelihood=1020), we decided to keep the historical values.

## 4. Discussion

We evaluated key control strategies in an HPAI epidemic used to control the course of the epidemic among birds and thus reduce potential transmission of the virus to humans. Unlike previous studies, our model takes into account farm size and spatial contacts.

We first estimated spatial contact rates from data describing the 2003 A/H7N7 Dutch outbreak. A sensitivity analysis was performed to evaluate the impact of control strategies and potential delays in diagnosis, reporting and implementation of control measures on the cumulative number of forecast cases.

The sensitivity analysis points out that immediate depopulation has the greatest impact in reducing avian influenza cases, as soon as a flock is reported infected (PRCC=0.85). To be effective, this strategy must be enacted quickly with a rapid detection of infection within flocks (PRCC=0.73). Further, the importance of this is evidenced by the high reproductive number during the pre-control period. These strategies have also been confirmed by field surveys (Marangon *et al*. 2003) and should be applied especially in layer type farms (Thomas *et al*. 2005). These two measures are likely to have more impact than increasing neighbourhood depopulation in HPAI epidemics (PRCC=−0.62). However, it may be necessary to apply proximity depopulation, as shown in the Netherlands where a 3 km depopulation radius was efficient in controlling the epidemic when combined with pre-existing control strategies. As illustrated in figure 4, *R*_{t} drastically decreased after the implementation of control measures. Despite the persistence of regions with *R*_{t}>1 even during the control period, the epidemic died out. This is probably a consequence of a saturation effect: no more susceptible farms were available to be infected.

Estimated spatial parameters clearly showed effective decreases in contact rates once when the disease was confirmed illustrating the impact of the control measures. Since the optimization method used here furnishes only point estimates, we explored the variation in the cumulative number of cases for contact rate values around those estimated through the uncertainty analysis. As shown in figure 5*b*, the range in epidemic size may be large (from 96 to 397 cases), but the likelihood, re-calculated for each set of sampled contact rates values, was set to its optimum for the estimated points. Thus, even though contact rates have an impact on the total number of cases, here we used the estimated values from data that provides the best optimization.

There are several limitations that require discussion. First, as exact data on farm location was not available; the model was based on radius buffers. Individual farms were assigned to a random location within the municipality for the purposes of the simulations. The individual-based model with random locations was preferred to a municipality metapopulation model (i.e. where the farms in a municipality were well-mixed) since transmission between farms and actual control were explicitly spatial. We tested the random placement assumption by repeating the simulations in which farms were assigned to different random locations within a municipality. Similar results were obtained for parameter estimation and epidemic curve shape, but the use of geographical information systems (GIS) to obtain more precise data concerning farms' location would help in moving to a continuous distance model. Second, we did not consider species other than domestic poultry. Waterfowl probably introduced the virus into domestic flocks in the Netherlands, but they were not presumed to play a role in the course of the epidemic afterwards. From 85 carcasses of free-living birds collected in the Gelderse valley, all were tested negative to A/H7N7 (European commission 2003). But this feature may not be relevant for all avian influenza epidemics, as during H5N1 epidemics many wild animals, birds (Ellis *et al*. 2004) or mammals (Kuiken *et al*. 2004), have been infected and may have taken part in virus spread. Even though pigs may play a major role in avian influenza dynamics through virus reassortment and transmission (Choi *et al*. 2004), the swine species did not seem to participate in the Dutch epidemic (European commission 2003). Therefore, no other species but domestic poultry were implemented in our model.

Detailed data on other epidemics and agricultural structures would be helpful to evaluate further the impact of other control strategies: new transmission buffer may be put into place and thus parameters could be re-estimated using the same estimation method as described in the present paper. Our results were obtained focusing only on the Dutch epidemic: what would be the most important control strategies if an epidemic occurred in another country? Vaccination using the differentiating infected from vaccinated animals (DIVA) strategy (Capua *et al*. 2003) coupled with a wider complex territorial strategy is recommended when there is a risk of major spread. These strategies have been successfully used in Italy between 2000 and 2003 (Marangon *et al*. 2003). In addition to the field studies (Lee *et al*. 2004), the potential impact of vaccination could be evaluated through mathematical models and used to provide important information for an A/H5N1 epidemic (Capua & Marangon 2004).

Avian influenza outbreaks cause dramatic damage in the poultry industry leading to bans on exportation and to intensive depopulations (Henzler *et al*. 2003). However, these epidemics also have an important impact in human societies. The primary threat is the emergence of a new pandemic after reassortment of an avian virus in humans or in other species, such as pigs. To prevent this threat, control strategies for avian outbreaks should focus on both high and low pathogenic avian influenza viruses as the virus may acquire additional virulence following antigenic drift or shift (Suarez *et al*. 2004). Detection of avian to human virus transmission, beyond just rare events, should be examined more thoroughly using epidemiological tools such as cluster size distribution (Ferguson *et al*. 2004). An efficient monitoring response in animals will help to reduce the number of infected humans, diminishing the risk of a pandemic (Webster & Hulse 2005).

The approach used here provides a framework for both mapping the risk of avian influenza and evaluating the impact of various control strategies. The model can be applied to retrospective data, as was the case here, or in theory used in real time with adequate information available on farm structures and transmission risk and with estimation of new contact parameters. The results of these simulations suggest that avian influenza risk is higher with diagnostic or depopulation delays, thereby providing important observations with clear public health implications to prevent avian virus transmission to humans.

## Acknowledgments

We gratefully acknowledge the Institut National de la Recherche Agronomique de Clermont-Ferrand and the Agence Française de la Sécurité Sanitaire des Aliments de Ploufragan for their support in reviewing the manuscript.

## Footnotes

- Received April 2, 2006.
- Accepted April 30, 2006.

- © 2006 The Royal Society