- Research
- Open
- Published:

# Estimating front-wave velocity of infectious diseases: a simple, efficient method applied to bluetongue

*Veterinary Research***volume 42**, Article number: 60 (2011)

## Abstract

Understanding the spatial dynamics of an infectious disease is critical when attempting to predict where and how fast the disease will spread. We illustrate an approach using a trend-surface analysis (TSA) model combined with a spatial error simultaneous autoregressive model (SAR_{err} model) to estimate the speed of diffusion of bluetongue (BT), an infectious disease of ruminants caused by bluetongue virus (BTV) and transmitted by *Culicoides*. In a first step to gain further insight into the spatial transmission characteristics of BTV serotype 8, we used 2007-2008 clinical case reports in France and TSA modelling to identify the major directions and speed of disease diffusion. We accounted for spatial autocorrelation by combining TSA with a SAR_{err} model, which led to a trend SAR_{err} model. Overall, BT spread from north-eastern to south-western France. The average trend SAR_{err}-estimated velocity across the country was 5.6 km/day. However, velocities differed between areas and time periods, varying between 2.1 and 9.3 km/day. For more than 83% of the contaminated municipalities, the trend SAR_{err}-estimated velocity was less than 7 km/day. Our study was a first step in describing the diffusion process for BT in France. To our knowledge, it is the first to show that BT spread in France was primarily local and consistent with the active flight of *Culicoides* and local movements of farm animals. Models such as the trend SAR_{err} models are powerful tools to provide information on direction and speed of disease diffusion when the only data available are date and location of cases.

## Introduction

Understanding the spatial dynamics of an infectious disease is critical when attempting to predict where and how fast the disease will spread. The use of simulation modelling for estimating the spread of infectious animal diseases has now become common [1–4]. However, most of these modelling approaches are complex, based on spatially explicit, stochastic, state-transition or network models [2, 5–7] or on diffusion equations [8, 9]. While such modelling approaches require precise knowledge on the model parameters, available data on emerging animal diseases are often restricted to case reports providing only date and location. These limited data nevertheless can provide useful information on the direction and speed of disease diffusion and be used to estimate front-wave velocity of an infectious disease. We illustrate an approach using a trend-surface analysis (TSA) model combined with a simultaneous autoregressive spatial error (SAR_{err}) model to estimate the speed of diffusion of bluetongue (BT).

BT is a non-contagious, infectious, viral disease of ruminants caused by bluetongue virus (BTV) and transmitted by biting midges (*Culicoides*) [10]. 24 serotypes of BTV have now been described. Until recently, BT was thought to be restricted to tropical regions and southern Europe where competent *Culicoides* species vectors are present [11]. The large-scale, BTV serotype 8 (BTV-8) epidemic in north-western Europe in 2006-2008 consequently surprised the veterinary community and caused major economic losses [12]. Contrary to what was previously thought, the abundant local *Culicoides* species in north-western Europe are vector competent [13, 14]. Disease progression was thus rapid: 2 000 infected farms across Belgium, Germany, The Netherlands, France, and Luxembourg in 2006 [12], more than 30 000 farms in 2007 with an expansion in range to Denmark, the United Kingdom, Switzerland, and the Czech Republic [12, 15].

While active flight of infected *Culicoides* is responsible for local propagation of BT, the movement of viremic animals and passive flight of infected *Culicoides* carried by the wind are responsible for long distance (>100 km) dissemination of the infection [16]. Quantification of the respective importance of local spread and long distance dissemination of BT currently is lacking. Gerbier et al. [17] suggested that the effect of wind was probably negligible for BT diffusion in 2006 and that local spread of BTV-8 can be explained mainly by active flight of *Culicoides*. Moreover, Mintiens et al. [18] showed that control measures implemented on animal movements failed to stop further spread of BTV because it is impossible to limit vector movements. The authors, however, noted that an absence of control measures probably would have resulted in an even wider and faster spread.

In a first step to gain further insight into the spatial transmission characteristics of BTV-8, we used 2007-2008 clinical case reports in France and TSA combined to spatial error modelling to identify the major directions and speed of diffusion of the BT epidemic.

## Materials and methods

### Compilation of BT clinical cases

We used BT cases recorded by the Direction Générale de l'Alimentation of the French Ministry of Agriculture, Food and Fishing in 2007-2008 to assess front-wave velocity. A case was defined as a bovine herd or an ovine or goat flock in which BT was clinically suspected and later confirmed by serological or virological analyses. Our analysis was performed on a municipality basis (the smallest administrative unit in France). Cases with missing date of record or location data were discarded, leaving 33 042 cases in 12 620 municipalities belonging to 82 departments (French administrative unit of a median surface area of 5 985 km^{2}). In 6.9% (*n* = 2 279) of these cases, the date of report was missing, however, the date of serological or virological confirmation was available. We were able to include these cases by extrapolating the date of report from the date of serological or virological confirmation. We did so by considering the 30 263 clinical cases with both date of report and date of confirmation, and verifying that the date of report and date of serological or virological confirmation were strongly correlated (*r*^{2} = 0.99) with a mean delay of 6 days (0-152), and 95% of the delay being lower than 19 days.

The data from the 33 042 clinical cases were then reduced to the first report of a case for each municipality. In the northern and central part of France, the French entomological surveillance reported only 7, 2 and 31 *Culicoides* captured in the 39, 34 and 46 traps set in January, February and March 2008, respectively (T. Balenghien, unpublished). Based on these entomological data and literature [13, 19], we considered that *Culicoides* activity was negligible from January to March, and thus excluded these three months from the subsequent analyses. As clinical signs of BT can be missed [16, 20], the disease may pass unnoticed in an infected herd over several weeks. To limit the bias due to unnoticed cases, we discarded municipalities with a very late first clinical case report in relation to the first clinical case report in the department. To do so, we estimated the mean length of a French department to be 105 km. Based on a mean velocity of 10 km per week [17], BT would pass through a French department within 75 days under optimal conditions. We added 75 more days to allow infected *Culicoides* to reach less accessible municipalities within the department. One hundred fifty days would thus be the restrictive delay to consider. However, as BT transmission depends on *Culicoides* density and activity, we took into account the month of the first case report of the department to balance the period of time during which we included municipalities depending on the month of the first case report in the department. The period of time during which we included municipalities with their first cases in relation to the month of the first case report in their department are presented Figure 1. The curve represents the distribution of the clinical cases reported per month over 2007-2008. The horizontal lines show the period of time during which we included the municipalities in each department depending on the month of the first clinical case report in the department. When the first clinical case report of the department occurred during the less favourable periods for vector activity (i.e., in April and October), we included municipalities for a longer period of time than when it occurred during the months of the peak of *Culicoides* activity (i.e., in July and August). In all, 1 627 municipalities in 51 departments were considered as having a late first case report and were discarded.

Our study consequently was based on 10 994 municipalities in 82 French departments which reported at least one case of BT (Figure 2).

### Trend-surface analysis

TSA is a least squares regression method used to study diffusion processes in space and time [21]. A surface pattern can be constructed by mapping the specific timing of events at each (X, Y) coordinate in two dimensions. The general procedure is described in Unwin [22]. Briefly, the method uses a model with power series polynomials, fitting linear, quadratic, cubic, and higher order trend-surfaces to the data [21]. The shape and flexibility of the trend surface are determined by the order of the polynomial chosen as the model [23]. A first-order polynomial restricts the trends to a plane through the data. Second-order polynomial models allow for curvature over the entire data set, while higher-order models allow for much more local curvature in the fitted surface [24]. In health research, TSA is a simple method which aims to capture the generalized direction(s) and speed(s) of the propagating wave of an infectious disease [24]. It previously has been used to assess the front-wave velocity of rabies and plague and to identify the pattern of disease diffusion [8, 23–25].

In our study, a polynomial surface was fitted to the set of spatially distributed times of first BT clinical case detection across the 10 994 municipalities. The (X, Y) coordinates of the municipality centroids were uploaded into ArcGIS software v.9.1 (ESRI Inc.). Original geographical coordinates were expressed in meters, *Lambert II étendu*. They were then translated into (X, Y) coordinates with the origin adjusted to the area of BT introduction. This area, referred to as the Area of First Infection (AFI), was identified in Gerbier et al. [17] as the border area between The Netherlands, Germany, and Belgium (latitude = 2 669 013 m; longitude = 3 726 278 m). Retrospective preliminary reports on the first BTV-8 outbreaks in Belgium and Germany indicated that the first clinical signs appeared between 17 July and 5 August 2006 [16, 26]. We therefore considered 17 July 2006 as the date of BT introduction in the AFI.

A model of the form t = β_{0} + β_{1}X + β_{2}Y + β_{3}X^{2} + β_{4}XY + β_{5}Y^{2} + ε was used to estimate linear and quadratic surfaces by least squares. In this model, t is the number of days to the BT introduction in 2006, β_{i} are the fitted parameters, X and Y are the geographic coordinates of the municipality centroids adjusted to the area of BT introduction in Europe, and ε represents the error term. The main caveat for significance testing in TSA is that spatial autocorrelation in the residuals is almost always present to a certain extent by the nature of the spatial data [22]. In our study, residuals of the TSA models showed autocorrelation (Moran's I statistic = 0.3894, *p* < 0.001). Among the numerous methods available to deal with spatial autocorrelation, autoregressive models incorporate spatial autocorrelation using neighbourhood matrices. These matrices specify the relationship between the response values (in the case of Conditional Autoregressive CAR models) or residuals (in the case of Simultaneous Autoregressive SAR models) at each location *i* and those at neighbouring locations *j* [27]. SAR models can take three different forms depending on where the spatial autoregressive process is believed to occur [27]. The spatial error model (SAR_{err}) assumes that the autoregressive process occurs only in the error term, and neither in response nor predictor variables. The lagged response model (SAR_{lag}) assumes that the autoregressive process occurs in the response variable, and the lagged-mixed model (SAR_{mix}) assumes that the spatial autocorrelation affects both response and predictor variables. Kissling and Carl [28] tested the performance of the three different SAR model types (SAR_{err}, SAR_{lag}, and SAR_{mix}) and ordinary least squares (OLS) regression. They showed that SAR_{err} models were the most reliable SAR models and that they performed well in all cases (independent of the kind of spatial correlation induced and whether models were selected by minRSA, R^{2} or AIC), whereas OLS, SAR_{lag} and SAR_{mix} models showed weak type I error control and/or unpredictable biases in parameter estimates. Based on this conclusion we thus chose SAR_{err} model to account for the spatial autocorrelation observed in the model residuals.

The neighbourhood relationship is formally expressed in a *n* × *n* matrix of spatial weights (W), with elements (w_{ij}) representing a measure of the connection between locations *i* and *j*. Details on SAR models are provided elsewhere [see [27–29]]. Briefly, the usual ordinary least squares regression model Y = β X_{i} + ε is complemented by a term (λWμ), which represents the spatial structure (λW) in the spatially dependent error term (μ) [27]. The SAR_{err} model thus takes the form Y = β X_{i} + λWμ + ε where λ is the spatial autoregression coefficient. We combined TSA with a spatial error simultaneous autoregressive (SAR_{err}) model to account for spatial autocorrelation in the residuals [27, 28]. This combination, hereafter referred to as trend SAR_{err} model, leads to a model that takes the form

We estimated the best trend SAR_{err} model using R software v.2.10.1 [30] and the spdep package [31]. We used a row standardized spatial weights matrix with a neighbourhood distance of 80 km. The neighbourhood distance of 80 km was chosen based on biological hypothesis as recommended by Dormann et al. [27]. The neighbourhood distance reflects the maximal distance at which the date of the first BT clinical case in a municipality influences the date of the first BT clinical case in surrounding municipalities. Based on the Commission Regulation No 1266/2007 the BT restricted zone, which represents the maximal distance at which contamination can occur from a BT infected farm, was defined as a 70-km radius around the infected farms. 70 km could thus be chosen as the neighbourhood distance to use in the spatial weights matrix. However, the minimal distance at which all municipalities of our dataset had at least one neighbour was 80 km. To avoid the problem of some areal entities having no neighbours [32] we chose a neighbourhood distance of 80 km, which is close to the 70-km radius of the restricted zone (see Additional file 1 for a description of the spatial weights matrix). Moran's I value, a measure of autocorrelation, and correlograms were calculated with the functions moran.test from the spdep package and correlog from the ncf package [33], respectively. Correlograms, which plot Moran's I values on the y-axis against geographic distance in the x-axis, allow the assessment of the spatial autocorrelation pattern with increasing distance. As recommended by Kissling and Carl [28], we computed two model selection statistics, the Akaike's Information Criterion (AIC) [34] and the minimum residual autocorrelation (minRSA). The latter was obtained by summing up the absolute Moran's I values in the first 80 distance classes of the correlogram [28]. The final model was chosen to minimize both AIC and minRSA based on backward selection. We considered that two nested models differing by less than 2 AIC points received identical support from the data. In such a situation, the model with less parameters was preferred [34].

To evaluate the selected model, a deviance-based pseudo-*R*^{2} (in the following simply referred to as *R*^{2}) was calculated as the squared Pearson correlation between predicted and observed values [28]. *R*^{2} provides a measure of goodness of fit of the model. Observed Moran's I values, computed using the residuals from the trend SAR_{err} model, provide a measure of the spatial autocorrelation in the model residuals [32, 35].

The model predicted a date of the first clinical case occurrence in each contaminated municipality. These dates were then contoured using a contouring and splining function in the Spatial Analyst extension of ArcGIS software v.9.1 (ESRI Inc.). We used 30-day intervals for the contour lines so that each band would represent a month, which would facilitate comparison of disease spread between periods of time.

The velocities of the travelling waves of BT were calculated using the best-fit surface model. This model was used to derive partial differential equations ∂t/∂X and ∂t/∂Y [24]. Substituting the specific (X, Y) coordinates observed for a case into these equations allowed the X and Y geometric vectors contributing to the slope of the travelling wavefront (i.e., its magnitude) to be described [36]. For each municipality, we summed the X and Y geometric vectors to obtain a resultant slope vector with magnitude and direction. The inverse of the slope is the velocity or speed of diffusion of the epidemic at each municipality location in kilometers per day. The larger the velocity vector, the faster the speed of diffusion [24].

## Results

### Data description

The date of the first clinical case of BTV-8 was considered as relevant for 10 994 municipalities in 82 departments. Descriptive statistics are provided in Table 1. The average number of municipalities per department was 415 (range: 102-894) with an average number of contaminated municipalities per department of 134 (1-533). The percentage of contaminated municipalities per department varied from 0.2 to 94% with a median value of 25.5% (Figure 3).

### The trend SAR_{err} model

Based on backward selection from the full quadratic trend SAR_{err} model we selected the model with the lowest AIC and minRSA (Table 2). All the trend SAR_{err} models had lower AIC values and less spatial autocorrelation in the residuals (lower minRSA) than the simple TSA model, which did not account for spatial autocorrelation. Adding the spatial autocorrelation term to the TSA model also improved the model fit. Models m0 and m5 were similar in terms of AIC and minRSA (Table 2), but m5, which did not include a first order parameter for X, had the highest *R*^{2} and the lowest Moran's I. Moreover, the values of the coefficients of m0 and m5 were similar (see Additional file 2), differing by less than 0.008, and the value of the estimate of the first order parameter for X in m0 was -0.0690, standard error = 0.0388, with a statistically non-significant *p*-value = 0.076. Based on these results we selected the model m5. This model provided an increase in *R*^{2} from 0.85 to 0.96. Correlograms for the residuals of the TSA model and the selected trend SAR_{err} model are showed in Figure 4. Although positive autocorrelation at small distances (less than 40 km) was still present in the residuals of the trend SAR_{err} model, autocorrelation was greatly reduced in comparison to the TSA model. Moran's I values and minRSA decreased by 93% and 78%, respectively (Table 2).

The residuals of the final trend SAR_{err} model had a mean non significantly different from zero (-2.1 e-12, 95% Confidence Interval (CI): -0.605-0.605) and a bell-shaped distribution (Figure 5a). To identify areas of earlier or later than predicted diffusion, the residuals of the final model were categorized into classes of residual values and plotted on the map of France (Figure 5b). An area of negative residual values (faster than predicted diffusion) in dark blue followed by an area of positive residual (slower than predicted diffusion) values in red are visually apparent in the centre of France. Overall, the difference between the predicted and observed date of the first clinical case occurrence was less than 30 days for 79% of the municipalities.

We used the final trend SAR_{err} model to estimate the day of the first clinical case occurrence in each contaminated municipality, based on the geographic coordinates of the municipality centroids (Table 3). The contoured trend-surface map (Figure 6) of the predicted time of occurrence illustrates the general direction and movement of the diffusion of BT in France in 2007-2008, month by month. Contour lines that are far apart indicate that the disease spread rapidly through an area while lines close together indicate slow progression. The direction of diffusion is given by the front of the contour lines [24]. Overall, BT spread from north-eastern (introduction from Belgium and Germany) to south-western France. The disease spread slower during the less favourable months for *Culicoides* activity (Figure 6): November-December and April to July. The same pattern of disease spread was observed in 2007 and 2008 with the faster spread of BT occurring in August-October.

### Estimated velocity from the trend SAR_{err} model

The average estimated velocity across the country was 5.6 km/day. However, velocities varied according to areas and time period (Figure 6). Overall, estimated velocities of BT spread varied between 2.1 and 9.3 km/day (Figure 7). For 9192 municipalities (84%), trend SAR_{err}-estimated velocities were lower than 7 km/day. 2% of the municipalities had an estimated velocity strictly lower than 3 km/day, which is the value commonly admitted of maximal active dispersal of *Culicoides* by flight [37].

## Discussion

We used a trend-SAR_{err} model to assess the velocity of BTV-8 spread across France between 2007 and 2008 and compare the respective importance of local spread and long distance dissemination. The average estimated velocity of BT was 5.6 km/day, or 39.2 km/week. In line with our results, the velocity of BT spread in Sardinia during the Italian BTV-2 outbreak was about 30 km/week [38]. Our estimation is higher than the findings of Gerbier et al. [17], who found a 10-15 km/week rate of spread, during the early stage of the BTV-8 epidemics in The Netherlands, Germany, and Belgium. The difference could be due to the methods used in the two studies: whereas we corrected for spatial autocorrelated residuals, Gerbier et al. did not. Interestingly, we found a mean velocity of 11.9 km/week, a similar value than the findings of Gerbier et al. [17], when we did not corrected for spatial autocorrelation in the model residuals. The difference in the average estimated velocity of BT could also be related to the heterogeneity of landscape, host density, and altitude of our study area, as well as the broad range of period of time considered in our analysis. Indeed, we studied the speed of BT diffusion across France during two years whereas Gerbier et al. [17] focused only on the early stage of the epidemics, from 17 August to 15 September 2006, in a limited homogeneous rectangular area of 90 × 75 km.

The range of the estimated velocities suggests that BT transmission in France was primarily local: 84% of the contaminated municipalities had an estimated velocity lower than 7 km/day, and none above 10 km/day. Regarding long-range dissemination, *Culicoides* can be passively dispersed overseas over long distances (> 100 km) by prevailing wind [39–41]. Hendrickx et al. [39] developed wind models which suggested that long-distance spread over land of *Culicoides* by prevailing wind could have allowed the BT epidemics to cover large distances at a rate of 15 km/day. Here, with no velocity above 10 km/day, such long-range dispersal of infected *Culicoides* by wind was probably uncommon and had a negligible effect on BT spread. Considering that the flight range of *Culicoides* is short and most species disperse only at about 3 km from their breeding sites [37, 42], active flight of *Culicoides* from farm to farm can not be the only factor explaining the local spread of BT. Indeed, the estimated velocity was less than 3 km/day for only 222 (2%) contaminated municipalities. Consequently, other factors must have facilitated the local diffusion of BT in France.

Several hypotheses can be proposed to explain these findings. The first hypothesis is that wind, besides facilitating the long-range dispersal of *Culicoides*, could have also increased the local dispersal of infected vectors [37, 40, 43]. Flight and orientation depends on the speed and direction of wind [43], but studies have shown that *Culicoides* prefer to shelter and cease almost all activity at wind speeds above 3 m/s or 10.8 km/h [37]. Interestingly, Bishop et al. [43] showed in Australia that high frequencies of wind speed above 8 km/h slowed the dispersion of *Culicoides brevitarsis*. Consequently, by stopping *Culicoides* activity, wind speeds above 3 m/s decrease local dispersal of infected vectors, whereas lower wind speeds could slightly increase the distance of local dispersal. Further analysis should investigate the wind patterns observed over the study area.

A second hypothesis is that dispersal distances may differ between *Culicoides* species. The maximal distance of active dispersal have only be estimated for North-American species: 3.2 km in 24 h for *Culicoides mississippiensis* [44], 4.8 km in 24 h for *Culicoides mohave* [45] and 5.6 km in 24 h for *Culicoides variipennis* [42]. Higher average distance of dispersal of the European species compared to the North-American ones could explain the velocities of BT spread estimated in our study. It will be important to determine whether dispersal capacities of European species are greater than the North-American species to support or invalidate this hypothesis.

Thirdly, local movements of infected animals may have facilitated BT diffusion around the contaminated farm units and have slightly increased the velocity of BT spread. Livestock are regularly moved within farms between pastures. However, little is known in literature on these movements. In France cattle farms have a mean surface area of 0.77 km^{2}, which results in a mean radius of 0.5 km if we consider the farm as a circle. Brunschwig et al. [46] reported that in France dairy cows are moved in a 0.5 km radius around the farm, while beef cattle are regularly moved in a 5 km radius, and occasionally as far as 20 km from the farm. Regarding farm-to-farm cattle movements, in Great Britain Mitchell et al. reported a mean straight-line distance of 58 km with 43% of the movements occurring within less than 20 km [47], in Portugal 80% of the cattle trade movements were local (< 40 km) [48], and in Sweden 87% of the cattle movements were to farms within 100 km [49]. Overall, these studies showed that a large part of the cattle movements between farms occurs at distance less than 100 km. Restrictions on farm animals movements were implemented in France following the Directive 2000/75/EC, which defined the restricted zone for BT as a 150 km radius around the BT contaminated farm. In October 2007 the Commission Regulation No 1266/2007 reduced the restricted zone to a 70-km radius around the contaminated farms. While regulations on animal transport prevented any movements from the restricted zone to the non-restricted zones, movements between farms could occur within the restricted zone and may thus partly explain the range of observed velocities.

Finally, wild ungulates may also play a role in BT diffusion because i) high seroprevalence have been reported in cervid species, ii) similar patterns were reported in domestic herd and wild ungulates, and iii) experimental infections suggested that wild ungulates could be infectious to *Culicoides*. BTV antibodies have been reported in various wild cervid species in Europe [50–54]. In France, in 2008-2009 a mean BTV seroprevalence rate of 41% and 1% was found in red deer (*Cervus elaphus*) and roe deer (*Capreolus capreolus*), respectively [55]. In red deer seroprevalence rates varied from 8 to 70% according to the study areas. In Spain Ruiz-Fons et al. [50] found a seroprevalence rate of 22% in red deer and 5% in roe deer and Garcìa et al. [54] reported a seroprevalence rate of 66% in red deer. Furthermore, in Spain similar spatial and temporal BTV patterns were observed in red deer and livestock [50, 54]. Moreover, the detection of BTV in skin samples of experimentally infected red deer suggests that they could be infectious to *Culicoides* [53]. Overall, red deer, in which BTV-infection does not induce a significant mortality [51], may serve as reservoir hosts. Movements of wild ungulate species could thus play a role in BT spread [54].

Although our data does not allow us to distinguish the effects of flight of *Culicoides* from the effects of short range movements of farm animals or wild ungulates (or any combination of the three) on BT spread, our results clearly show that BT diffusion in France was primarily local. Few measures were available in France to control the spread of the infection until the summer of 2008. While BTV-8 vaccines became available in the spring of 2008, and progressively were used in the second part of the year to protect livestock, regulations on animal transport probably were the most valuable measure in limiting the spread of BT through long-range dissemination [11].

TSA is one of numerous methods that can be used in the analysis of change over space. It mainly has been used in geological studies but we encourage researchers who are interested in a simple method to estimate the direction and speed of an infectious disease spread based solely on the geographical location and date of cases. However, investigators must be aware of the limitations of the method. A disadvantage of TSA is that predictions or extrapolations outside the area and time of study are not accurate and should not be made [21]. In our study, residual values of the TSA model showed high level of autocorrelation (Figure 4), but by combining the TSA model with a spatial error model, the observed Moran's I value decreased by 93%. The residuals of the trend SAR_{err} model remained slightly positively autocorrelated at small distances (< 40 km). This may be due to i) the large dataset or ii) the omission of spatially patterned explanatory variables [27, 32, 56, 57]. Indeed, Koenig stated [58]: "With large datasets, statistically significant values can be obtained even though the absolute degree of correlation is small; whether such low spatial autocorrelation is biologically significant or not must be considered on a case-by-case basis". Moreover, positive autocorrelation in the residuals at small distances can be observed when environmental factors associated to small-scale variations in the explained variable are omitted [56]. This paper shows that with simple epidemiological data (geographic coordinates and date of clinical cases), the TSA method enables to estimate and map the spread of a disease [24, 36]. This first step can then be completed by modelling the effect of a wide range of environmental variables on the estimated velocities. Overall, both our large dataset (10 994 municipalities) and the use of geographic coordinates may explain the slightly autocorrelated residuals of the trend SAR_{err} model. However, because i) TSA was specifically developed to fit spatial data, and ii) the remaining autocorrelation is small, we believe that it did not bias the estimated velocities. Unwin [22] noted that if the sample size is sufficiently large, and provided one avoids using a higher-order surface model, the TSA method is robust enough to allow some violation of the least-squares regression assumptions. Furthermore, the model was not use as a predictive tool to extrapolate the date of the first clinical case occurrence in other areas or at different periods of time. Overall, with a *R*^{2} value of 0.96, and 79% of the municipalities for which the difference between the predicted and observed date of the first clinical case occurrence was less than 30 days, our model fitted the data well.

Another limitation of the TSA method is that interpretation at the edges of the map should be made with caution [24]. Residuals at the edges of the study area presented extreme positive and negative values (Figure 5b). The same pattern of extreme residual values also occurred around a geographical area in the centre of France where no or very few municipalities were included in our dataset, which induced an edge effect in the centre of the study area. The absence of BT infected municipalities in this area was due to the absence of, or to incomplete data on, clinical cases linked with low livestock densities. The trend SAR_{err} model had problems handling such areas with no or few data and this led to areas with extreme residual values.

The pattern of autocorrelation in residual values can however be used to indicate heterogeneities in disease spread. Residuals should always be examined and lead to questions of why the disease occurred at a particular location earlier or later than predicted by the model [24]. In our case, most areas with extreme residual values in the centre of France (in red and dark blue in Figure 5b) correspond to the geographic areas where BT spread was very low immediately prior to and after the vector-free period (January to March 2008).

The 10 994 municipalities included in the study did not represent an exhaustive sample of municipalities with BT-infected animals. Indeed, some municipalities with clinical cases but incomplete data were excluded from the analysis, and we did not consider municipalities where they were no clinical cases but where seropositive farm animals were reported. However, while the case detection system based on clinical suspicion underestimated the real impact of the epidemic, it indicated a correct spatial trend [59].

In conclusion, TSA modelling procedure, combined with a spatial error model to limit the effect of autocorrelated residuals, is a powerful tool to provide information on the direction and speed of disease diffusion when the only data available are case date and location. Our study was a first step in describing the diffusion process for BT in France and showed that the spread of BT was primarily local. The next step in the spatial analytic process is to understand why BT spread the way it did. This will be achieved by using logistic regression to investigate the correlation between ecological factors, i.e., landscape pattern, farm animal density, meteorological conditions, and velocity of spread.

## References

- 1.
Keeling MJ: Models of foot-and-mouth disease. Proc Biol Sci. 2005, 272: 1195-1202. 10.1098/rspb.2004.3046.

- 2.
Russell CA, Smith DL, Waller LA, Childs JE, Real LA: A priori prediction of disease invasion dynamics in a novel environment. Proc Biol Sci. 2004, 271: 21-25. 10.1098/rspb.2003.2559.

- 3.
Guitian J, Pfeiffer D: Should we use models to inform policy development?. Vet J. 2006, 172: 393-395. 10.1016/j.tvjl.2006.03.001.

- 4.
Hosseini PR, Dhondt AA, Dobson AP: Spatial spread of an emerging infectious disease: conjunctivitis in house finches. Ecology. 2006, 87: 3037-3046. 10.1890/0012-9658(2006)87[3037:SSOAEI]2.0.CO;2.

- 5.
Harvey N, Reeves A, Schoenbaum MA, Zagmutt-Vergara FJ, Dubé C, Hill AE, Corso BA, McNab WB, Cartwright CI, Salman MD: The North American Animal Disease Spread Model: a simulation model to assist decision making in evaluating animal disease incursions. Prev Vet Med. 2007, 82: 176-197. 10.1016/j.prevetmed.2007.05.019.

- 6.
Smith DL, Lucey B, Waller LA, Childs JE, Real LA: Predicting the spatial dynamics of rabies epidemics on heterogeneous landscapes. Proc Natl Acad Sci USA. 2002, 99: 3668-3672.

- 7.
Smith GC, Wilkinson D: Modelling disease spread in a novel host: rabies in the European badger

*Meles meles*. J Appl Ecol. 2002, 39: 865-874. 10.1046/j.1365-2664.2002.00773.x. - 8.
Maidana NA, Yang HM: Spatial spreading of West Nile Virus described by traveling waves. J Theor Biol. 2009, 258: 403-417. 10.1016/j.jtbi.2008.12.032.

- 9.
Lewis M, Rencławowicz J, den Driessche P: Traveling waves and spread rates for a West Nile virus model. Bull Math Biol. 2006, 68: 3-23. 10.1007/s11538-005-9018-z.

- 10.
Mellor PS, Baylis M, Mertens P: Bluetongue. 2009, Elsevier, Academic Press, London, UK

- 11.
De Koeijer A, Elbers ARW: Modelling of vector-borne diseases and transmission of bluetongue virus in North-West Europe. Emerging pests and vector-borne diseases in Europe. Edited by: Takken W, Knols BGJ. 2007, Wageningen Academic Publishers, Wageningen, The Netherlands, 99-112.

- 12.
Carpenter S, Wilson A, Mellor PS:

*Culicoides*and the emergence of bluetongue virus in northern Europe. Trends Microbiol. 2009, 17: 172-178. 10.1016/j.tim.2009.01.001. - 13.
Meiswinkel R, Baldet T, de Deken R, Takken W, Delécolle JC, Mellor PS: The 2006 outbreak of bluetongue in northern Europe - the entomological perspective. Prev Vet Med. 2008, 87: 55-63. 10.1016/j.prevetmed.2008.06.005.

- 14.
Takken W, Verhulst N, Scholte E-J, Jacobs F, Jongema Y, van Lammeren R: The phenology and population dynamics of

*Culicoides*spp. in different ecosystems in The Netherlands. Prev Vet Med. 2008, 87: 41-54. 10.1016/j.prevetmed.2008.06.015. - 15.
Wilson A, Mellor P: Bluetongue in Europe: vectors, epidemiology and climate change. Parasitol Res. 2008, 103: S69-S77. 10.1007/s00436-008-1053-x.

- 16.
Saegerman C, Mellor P, Uyttenhoef A, Hanon J-B, Kirschvink N, Haubruge E, Delcroix P, Houtain J-Y, Pourquier P, Vandenbussche F, Verheyden B, De Clercq K, Czaplicki G: The most likely time and place of introduction of BTV8 into belgian ruminants. PLoS One. 2010, 5: e9405-10.1371/journal.pone.0009405.

- 17.
Gerbier G, Baldet T, Hendrickx G, Guis H, Mintiens K, Elbers ARW, Staubach C: Modelling local dispersal of bluetongue virus serotype 8 using random walk. Prev Vet Med. 2008, 87: 119-130. 10.1016/j.prevetmed.2008.06.012.

- 18.
Mintiens K, Méroc E, Faes C, Abrahantes JC, Hendrickx G, Staubach C, Gerbier G, Elbers ARW, Aerts M, De Clercq K: Impact of human interventions on the spread of bluetongue virus serotype 8 during the 2006 epidemic in north-western Europe. Prev Vet Med. 2008, 87: 145-161. 10.1016/j.prevetmed.2008.06.010.

- 19.
Clausen P-H, Stephan A, Bartsch S, Jandowsky A, Hoffmann-Köhler P, Schein E, Mehlitz D, Bauer B: Seasonal dynamics of biting midges (Diptera: Ceratopogonidae,

*Culicoides*spp.) on dairy farms of Central Germany during the 2007/2008 epidemic of bluetongue. Parasitol Res. 2009, 105: 381-386. 10.1007/s00436-009-1417-x. - 20.
Szmaragd C, Wilson AJ, Carpenter S, Wood JLN, Mellor PS, Gubbins S: A modeling framework to describe the transmission of bluetongue virus within and between farms in Great Britain. PLoS One. 2009, 4: e7741-10.1371/journal.pone.0007741.

- 21.
Moore DA, Carpenter TE: Spatial analytical methods and geographic information systems: use in health research and epidemiology. Epidemiol Rev. 1999, 21: 143-161.

- 22.
Unwin DJ: An introduction to trend surface analysis. Concepts and Techniques in Modern Geography. 1975, Geo Abstracts Ltd., University of East Anglia, Norwich, 1-40.

- 23.
Lucey BT, Russell CA, Smith D, Wilson ML, Long A, Waller LA, Childs JE, LA R: Spatiotemporal analysis of epizootic raccoon rabies propagation in Connecticut, 1991-1995. Vector Borne Zoonotic Dis. 2002, 2: 77-86. 10.1089/153036602321131878.

- 24.
Moore DA: Spatial diffusion of raccoon rabies in Pennsylvania, USA. Prev Vet Med. 1999, 40: 19-32. 10.1016/S0167-5877(99)00005-7.

- 25.
Ball FG: Front-wave velocity and fox habitat heterogeneity. Population dynamics of rabies in wildlife. Edited by: Bacon PJ. 1985, Academic Press, New-York, 255-289.

- 26.
Mintiens K, Méroc E, Mellor PS, Staubach C, Gerbier G, Elbers ARW, Hendrickx G, De Clercq K: Possible routes of introduction of bluetongue virus serotype 8 into the epicentre of the 2006 epidemic in north-western Europe. Prev Vet Med. 2008, 87: 131-144. 10.1016/j.prevetmed.2008.06.011.

- 27.
Dormann CF, McPherson JM, Araújo MB, Bivand R, Bolliger J, Carl G, Davies RG, Hirzel A, Jetz W, Kissling WD, Kühn I, Ohlemüller R, Peres-Neto PR, Reineking B, Schröder B, Schurr FM, Wilson R: Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography. 2007, 30: 609-628. 10.1111/j.2007.0906-7590.05171.x.

- 28.
Kissling WD, Carl G: Spatial autocorrelation and the selection of simultaneous autoregressive models. Global Ecol Biogeogr. 2008, 17: 59-71.

- 29.
Haining R: Spatial data analysis: theory and practice. 2003, Cambridge University Press, Cambridge

- 30.
The R Foundation for Statistical Computing: R statistical software. 2009

- 31.
Bivand R: spdep: spatial dependence: weighting schemes, statistics, and models. R package version 0.5-24. 2005

- 32.
Bivand RS, Pebesma EJ, Gómez-Rubio V: Applied spatial data analysis with R. 2008, Springer, New-York

- 33.
Bjørnstad ON: ncf: spatial nonparametric covariance functions. R package version 1.1-3. 2005

- 34.
Burnham KP, Anderson DR: Model selection and Multi-Model Inference, a practical information-theoretic approach. 2002, Springer-Verlag, New-York

- 35.
Pfeiffer DU, Robinson TP, Stevenson M, Stevens KB, Rogers DJ, Clements ACA: Spatial analysis in epidemiology. 2008, Oxford University Press, Oxford

- 36.
Adjemian JZ, Foley P, Gage KL, Foley JE: Initiation and spread of traveling waves of plague,

*Yersinia pestis*, in the western United States. Am J Trop Med Hyg. 2007, 76: 365-375. - 37.
Mellor PS, Boorman J, Baylis M:

*Culicoides*biting midges: their role as arbovirus vector. Annu Rev Entomol. 2000, 45: 307-340. 10.1146/annurev.ento.45.1.307. - 38.
Calistri P, Giovannini A, Conte A, Nannini D, Santucci U, Patta C, Rolesu S, Caporale V: Bluetongue in Italy: Part I. Vet Ital. 2004, 40: 243-251.

- 39.
Hendrickx G, Gilbert M, Staubach C, Elbers A, Mintiens K, Gerbier G, Ducheyne E: A wind density model to quantify the airborne spread of Culicoides species during north-western Europe bluetongue epidemic, 2006. Prev Vet Med. 2008, 87: 162-181. 10.1016/j.prevetmed.2008.06.009.

- 40.
Ducheyne E, De Deken G, Bécu S, Codina B, Nomikou K, Mangana-Vougiaki O, Georgiev G, Purse BV, Hendrickx G: Quantifying the wind dispersal of

*Culicoides*species in Greece and Bulgaria. Geospat Health. 2007, 2: 177-189. - 41.
Sellers RF: Possible windborne spread of bluetongue to Portugal, June-July 1956. J Hyg. 1978, 81: 189-196. 10.1017/S0022172400025018.

- 42.
Lillie TH, Marquardt WC, Jones RH: The flight range of

*Culicoides variipennis*(Diptera: Ceratopogonidae). Can Entomol. 1981, 113: 419-426. 10.4039/Ent113419-5. - 43.
Bishop AL, Barchia IM, Spohr LJ: Models for the dispersal in Australia of the arbovirus vector,

*Culicoides brevitarsis*Kieffer (Diptera: Ceratopogonidae). Prev Vet Med. 2000, 47: 243-254. 10.1016/S0167-5877(00)00175-6. - 44.
Lillie TH, Kline DL, Hall DW: The dispersal of

*Culicoides mississippiensis*(Diptera: Ceratopogonidae) in a salt marsh near Yankeetown, Florida. J Am Mosq Control Assoc. 1985, 1: 465-468. - 45.
Brenner RJ, Wargo MJ, Stains GS, Mulla MS: The dispersal of

*Culicoides mohave*(Diptera: Ceratopogonidae) in the desert of southern California. Mosq news. 1984, 44: 343-350. - 46.
Brunschwig G, Josien E, Bernhard C: Contraintes géographiques et modes d'utilisation en élevage bovin laitier et allaitant. Fourrages. 2006, 185: 83-95. (in French)

- 47.
Mitchell A, Bourn D, Mawdsley J, Wint W, Clifton-Hadley R, Gilbert M: Characteristics of cattle movements in Britain - an analysis of records from the cattle tracing system. Anim Sci. 2005, 80: 265-273.

- 48.
Baptista FM, Nunes T: Spatial analysis of cattle movement patterns in Portugal. Vet Ital. 2007, 43: 611-619.

- 49.
Nöremark M, Håkansson N, Lindström T, Wennergren U, Sternberg Lewerin S: Spatial and temporal investigations of reported movements, births and deaths of cattle and pigs in Sweden. Acta Vet Scand. 2009, 51: 37-10.1186/1751-0147-51-37.

- 50.
Ruiz-Fons F, Reyes-Garcia AR, Alcaide V, Gortazar C: Spatial and temporal evolution of bluetongue virus in wild ruminants, Spain. Emerg Infect Dis. 2008, 14: 951-953. 10.3201/eid1406.071586.

- 51.
Linden A, Grégoire F, Hanrez D, Mousset B, Massart AL, De Leeuw I, Vandemeulebroucke E, Vandenbussche F, De Clercq K: Bluetongue virus in wild deer, Belgium, 2005-2008. Emerg Infect Dis. 2010, 16: 833-836.

- 52.
Rodríguez-Sánchez B, Sánchez-Cordón PJ, Molina V, Risalde MA, Pérez de Diego AC, Gómez-Villamandos JC, Sánchez-Vizcaíno JM: Detection of bluetongue serotype 4 in mouflons (

*Ovis aries musimon*) from Spain. Vet Microbiol. 2010, 14: 164-167. - 53.
López-Olvera JR, Falconi C, Férnandez-Pacheco P, Fernández-Pinero J, Sánchez MÁ, Palma A, Herruzo I, Vicente J, Jiménez-Clavero MÁ, Arias M, Sánchez-Vizcaíno JM, Gortazar C: Experimental infection of European red deer (Cervus elaphus) with bluetongue virus serotypes 1 and 8. Vet Microbiol. 2010, 145: 148-152. 10.1016/j.vetmic.2010.03.012.

- 54.
García I, Napp S, Casal J, Perea A, Allepuz A, Alba A, Carbonero A, Arenas A: Bluetongue epidemiology in wild ruminants from Southern Spain. Eur J Wildl Res. 2009, 55: 173-178. 10.1007/s10344-008-0231-6.

- 55.
Rossi S, Gibert P, Bréard E, Moinet M, Hars J, Maillard D, Wanner M, Klein F, Mastain O, Mathevet P, Bost F: Circulation et impact des virus de la fièvre catarrhale ovine (FCO) chez les ruminants sauvages en France. Bull Epidémiol AFSSA/DGAL. 2010, 35: 28-32. (in French)

- 56.
Diniz-Filho JAF, Bini LM, Hawkins BA: Spatial autocorrelation and red herrings in geographical ecology. Global Ecol Biogeogr. 2003, 12: 53-64. 10.1046/j.1466-822X.2003.00322.x.

- 57.
Lichstein JW, Simons TR, Shriner SA, Franzreb KE: Spatial autocorrelation and autoregressive models in ecology. Ecol Monogr. 2002, 72: 445-463. 10.1890/0012-9615(2002)072[0445:SAAAMI]2.0.CO;2.

- 58.
Koenig WD: Spatial autocorrelation of ecological phenomena. TREE. 1999, 14: 22-26.

- 59.
Méroc E, Faes C, Herr C, Staubach C, Verheyden B, Vanbinst T, Vandenbussche F, Hooyberghs J, Aerts M, De Clercq K, Mintiens K: Establishing the spread of bluetongue virus at the end of the 2006 epidemic in Belgium. Vet Microbiol. 2008, 131: 133-144. 10.1016/j.vetmic.2008.03.012.

## Acknowledgements

We thank the French Ministry of Agriculture for having funded this research and provided access to the national BT dataset, as well as to the data on entomological surveillance. We are grateful to T Balenghien for providing data on *Culicoides* abundance and we thank the two anonymous referees who provided valuable comments on an earlier draft of this paper.

## Author information

## Additional information

### Competing interests

The authors declare that they have no competing interests.

### Authors' contributions

MP participated in the design of the study, performed the statistical analysis and drafted the manuscript. HG, BD and DC participated in the design of the study and DA participated in the statistical analysis. CD conceived of the study, and participated in its design and coordination. All authors read, amended and approved the final manuscript.

## Rights and permissions

## About this article

#### Received

#### Accepted

#### Published

#### DOI

### Keywords

- Spatial Autocorrelation
- Spatial Weight Matrix
- Spatial Error Model
- Wild Ungulate
- Clinical Case Report