Abstract
Bluetongue (BT) can cause severe livestock losses and large direct and indirect costs for farmers. To propose targeted control strategies as alternative to massive vaccination, there is a need to better understand how BT virus spread in space and time according to local characteristics of host and vector populations. Our objective was to assess, using a modelling approach, how spatiotemporal heterogeneities in abundance and distribution of hosts and vectors impact the occurrence and amplitude of local and regional BT epidemics. We built a reaction–diffusion model accounting for the seasonality in vector abundance and the active dispersal of vectors. Because of the scale chosen, and movement restrictions imposed during epidemics, host movements and windinduced passive vector movements were neglected. Four levels of complexity were addressed using a theoretical approach, from a homogeneous to a heterogeneous environment in abundance and distribution of hosts and vectors. These scenarios were illustrated using data on abundance and distribution of hosts and vectors in a real geographical area. We have shown that local epidemics can occur earlier and be larger in scale far from the primary case rather than close to it. Moreover, spatial heterogeneities in hosts and vectors delay the epidemic peak and decrease the infection prevalence. The results obtained on a real area confirmed those obtained on a theoretical domain. Although developed to represent BTV spatiotemporal spread, our model can be used to study other vectorborne diseases of animals with a local to regional spread by vector diffusion.
Introduction
There is significant concern regarding the resurgence or emergence of vectorborne diseases of animals with serious consequences for animal health and economics [13]. Climate change and socioeconomic change are both believed to contribute to the emergence and spread of such diseases [4]. A recent example is the unexpected introduction of the serotype 8 of the bluetongue virus (BTV) in northern Europe in 2006. Bluetongue is a noncontagious vectorborne disease affecting domestic and wild ruminants with high direct and indirect economic consequences [5,6]. It spread for three years with an annual epidemic peak followed by the disappearance of clinical cases.
A better understanding of the temporal and spatial spread of BTV has direct consequences for the disease prevention and control. The recent incursion of BTV in Europe has been controlled using a massive vaccination. To propose alternative to such a massive strategy if BTV incursions were to occur, we need to better identify where and when targeted strategies should be implemented. Therefore, the occurrence and amplitude of both local (a few km^{2}) and regional epidemics should be more precisely predict.
The spatiotemporal heterogeneity in abundance and distribution of hosts and vectors generally has a strong impact on pathogen spread and persistence [7,8]. In a seasonal environment such as in Europe, bluetongue is characterized by strong seasonal variations in incidence related to the seasonality of the vector population [9,10], whose lifecycle largely depends on environmental factors, such as humidity and temperature. As a result, clinical cases almost disappear during the unfavourable season for the vector. In addition to the temporal heterogeneity in vector abundance, the heterogeneity in the spatial distribution of vectors and hosts may also impact bluetongue spread [11,12]. Such heterogeneities in vector and host abundance and distribution can differ between geographic areas. In livestock populations, they relate to the landscape structure as well as to farming practices as animal populations are managed by farmers.
Mathematical modelling is a relevant approach for investigating the spread of vectorborne diseases [8,1316]. As hosts and vectors are mobile entities, a spatial component in vectorborne disease models should be taken into account to better consider the evolution of the biological system [1719]. This spatial component is not only due to space structuring in terms of density and location of host and vector populations, but also to host and vector movements over space. Different methods of various levels of complexity exist to include this spatial component in epidemiological models. To study the spread of vectorborne diseases, spatially explicit models are generally preferred [13,20]. They permit to take into account both vector active and passive movements and host movements that occur at different scales. Moreover, such models have been used also to describe the velocity of travelling waves of epidemics [8,14,16].
Such a modelling approach has been used to represent the spatiotemporal BTV spread in specific areas [2124]. If the published models took into account the heterogeneity in distribution and abundance of livestock populations, they tend to assume uniform densities of vectors (i.e. the same number per farm) or uniform host: vector ratio (i.e. more midges on bigger farms) [2124]. Recently, it has been shown that these are probably unrealistic [25]. Therefore, models that take full account of vector heterogeneity both in space and time have to be developed.
Our objective was to assess, using a modelling approach, how spatiotemporal heterogeneities in abundance and distribution of hosts and vectors impact the occurrence and amplitude of local and regional BT epidemics. We first studied different hypothetical scenarios of spatial heterogeneity in host and vector populations, and then illustrated the theoretical results in a real geographic area.
Material and methods
Model description
Three actors are necessary for bluetongue spread: the virus (here BTV8), the vector (a midge, here Culicoides), and the host (a ruminant, here cattle). As the virus is not excreted, we have assumed all transmission is vectorial. The developed model is based on a more complex model of the seasonal temporal spread of bluetongue in cattle [26]. This model has been simplified and made spatial. Here, the vertical transmission (in utero) has not been taken into account, as this hypothesis has not been shown to influence the infection [26]. The temporal dynamics of the vector population has been modified using a more flexible function. We used a standard compartment model (Figure 1) to describe the transmission of a pathogen between a vertebrate host population (HP) and a vector population (VP). The parameters of this model are defined in Table 1. The host population (HP) is divided into three health states: susceptible (SH), infectious (IH), and immune (RH). It is assumed to remain constant: the entry rate (b_{H}) compensates the exit rate (m_{H}). The vector population (VP) is divided into three health states: susceptible (SV), exposed (EV), and infectious (IV). For vectors, a latency period (1/ρ_{E}; Table 1) was taken into account as it is of the same order as life expectancy. At the diseasefree state and with seasonality, the vector population was assumed to have a logistic growth with K(t) = (b_{V}m_{V})/k_{V}(t) (Table 1), the carrying capacity of the environment depending on the vector fertility (b_{V}), mortality (m_{V}) and densitydependant mortality (k_{V}) rates. In periods favourable for vectors, the K function is a sinusoidal function with a maximum h. In unfavourable periods, the K function is constant and equal to Nb. The vectorial transmission takes place when an infectious vector (IV) bites a susceptible host (SH) which becomes infectious (IH), or when a susceptible vector (SV) bites an infectious host (IH) and then becomes exposed (EV). Incidence functions were frequency dependant for hosts and vectors. The mean host viremia duration (i.e. the time spent in IH) was 1/α_{I}. With recovery, animals move from state IH to RH.
Figure 1. Conceptual model of BTV8 spread. Flow diagram describing the model used for BTV8 spread in midge and cattle populations. Squares represent the health states of hosts (H), circle those of vectors (V), with S for susceptible, E for latent, I for infectious, R for recovered. The descriptions, values and sources of all parameters in the epidemiological model are found in Table 1.
Table 1. Parameters of the model of BTV8 spread in midge and cattle populations
Let Ω be the square spatial domain. X = X (x, y, t) represents time dependant population densities in (x, y) ∈ Ω. During the epidemic, host movements are controlled, therefore the spatial spread of the epidemic is due to vector movements rather than host movements. Moreover, we focus on a local to regional scale, and therefore assume that the spatial spread of BTV8 is exclusively due to local movements of vectors. BTV8 having no detrimental impact on vectors, thereby the diffusion process is similar whatever the health state. Therefore, the diffusion process follows the first Fick’s law:
By adding the reaction term, we obtain the following system of equations (Eq. 1) describing the spatiotemporal spread of the BTV8, for (x, y)∈ Ω and t > 0:
The initial time corresponds to the first day of the favourable season for the vector population. The spatial domain (Ω) is discretized into cells; each cell of surface area of 1 km^{2}. Initially, all hosts and vectors are susceptible, except the primary case which corresponds to an infected host, placed in the centre of the grid, in a cell where there are hosts. We set nonnegative initial conditions (Eq. 2). The flow of individuals across the domain boundary is assumed to be zero, i.e. we do not consider immigration or emigration of individuals and set Neumann boundary conditions (Eq. 2).
To discretize the problem (Eq. 1 and 2) we used the finite difference method in space, and we converted the continuous time model into a discrete time one by using the semiimplicit Euler method, that we implemented in Scilab 5.1.
Hypotheses of spatial heterogeneities in hosts and vectors
Assumptions are described by increasing level of complexity in Table 2. A first reference hypothesis (H1) considers a homogeneous spatial distribution and abundance in hosts and vectors. Four scenarios were studied for four maxima of the carrying capacity in vectors (h1 to h4). From this assumption, four levels of heterogeneity were analyzed and compared.
Table 2. Hypotheses of spatial heterogeneities in abundance and distribution of hosts and vectors
First, a heterogeneous distribution of hosts was considered (H2), vectors being homogeneously distributed. We kept the total number of hosts on the grid constant. We tested five densities of occupied cells from 25% to 100%, where H1 is 100% of cells occupied (Figure 2). For each density, we randomly drew the occupied cells among the 1681 grid cells. For each scenario, three draws of occupied cells were compared. Results were identical whatever the distribution; therefore it did not influence the spatiotemporal spread of the virus on the grid. Thereafter, only one distribution was kept (Figure 2).
Figure 2. Spatial distribution of cells occupied by hosts for four density levels (hypotheses H2 and H5). Black cells correspond to occupied cells; white cells correspond to empty cells; the yellow and crossed cell corresponds to the cell wherein the primary case occurs.
Secondly, a heterogeneous distribution of vectors was considered (H3 and H4), hosts being homogeneously distributed. In H3, the domain (Ω) consisted of four subareas, each having a different maximum of the carrying capacity in vectors (h1 to h4, Figure 3a). In H4, we considered a random equiprobable distribution of the same four maxima of the carrying capacity in vectors on the grid (Figure 3b).
Figure 3. Spatial distribution of the four maximum carrying capacities in vectors. a) Hypothesis H3, b) Hypotheses H4 and H5. Yellow cells correspond to cells having a maximum carrying capacity in vectors given by h1 = 10^{6}_{,} orange cells correspond to cells having a maximum carrying capacity in vectors given by h2 = 10^{7}_{,} red cells correspond to cells having a maximum carrying capacity in vectors given by h3 = 10^{8}, brown cells correspond to cells having a maximum carrying capacity in vectors given by h4 = 10^{9}. The black cell corresponds to the cell wherein the primary case occurs.
Thirdly, we considered a heterogeneous distribution of hosts and vectors simultaneously (H5). This hypothesis crosses previous hypotheses H2 and H4.
Fourth, a last hypothesis (H6), based on real data, served to illustrate this theoretical work, in particular hypothesis H5.
Data
The Culicoides trap catches used for modelling were collected in the Welsh province of Bala, situated in Snowdonia National Park (for full methods, see [25]). The trapping farms were selected using a 6 × 6 km grid, whereby one farm was selected from each grid square (36 in total, Figure 4). Each farm was sampled for three nights using Onderstepoorttype down draught black light traps positioned as close to livestock as possible. Large collections were subsampled [27] and only females were considered in the analyses as males do not take blood meals or, consequently, transmit disease between vertebrates. The maximum trapcatch of Culicoides per farm, out of the three trapping nights, was used for modelling purposes (Figure 4a). Due to the nature of the terrain, two squares contained no properties. Culicoides counts for these grid squares were estimated using the models of [25]. The vector abundance is difficult to quantify. Therefore, we considered two scenarios, one multiplying the number of vectors per cell by 100, and another by 1000. The numbers of cattle per farm were recorded on questionnaires during data collection (Figure 4b). For farms with unknown numbers of cattle, values were interpolated using the known values of farms in adjacent grid squares.
Figure 4. Vector (a) and host (b) spatial distribution in the real area (in the Welsh province of Bala). From yellow to brown, increasing numbers of individuals (white: empty cells); the cell with the star corresponds to the cell wherein the primary case occurs.
Outputs
The date and the prevalence at the epidemic peak in each cell were analysed, as well as the total prevalence on the grid over time. Thereafter, these three outputs are respectively named peak date, local prevalence and total prevalence. The peak dates were compared among the cells located on the four lines between the central cell (the halfdiagonal), i.e. the cell of virus introduction, and the corners of the grid. This enabled us to numerically calculate an effective speed of the virus spread. We calibrated the diffusion coefficient (D) and the initial conditions (SH00^{0} and SV^{0}) to have an effective speed of the virus spread similar to the estimated velocity by Pioz et al. [28], for hypothesis H1 and a maximum of the carrying capacity in vectors equals to 10^{7}. The theoretical grid is a 41 × 41 km square, each halfdiagonal measuring about 29 km. For hypotheses H3 to H6, peak dates and local prevalences were studied for comparable cells, i.e. equidistant and having the same maximum of the carrying capacity in vectors.
Results
Homogeneous in abundance and distribution of hosts and vectors (H1)
For the lowest maximum of carrying capacity in vectors (h1 = 10^{6}), there is no epidemic. For the other values tested, there is an epidemic peak in all grid cells. A larger maximum of carrying capacity in vectors leads to an earlier peak and a faster speed of virus spread (Figure 5). For h2 = 10^{7}, five days are necessary for the virus to cover the halfdiagonal from the cell of the virus introduction, with the epidemic peak occurring in the central cell and in the cells at the end of the halfdiagonals at 96 and 101 days respectively after virus introduction. As expected in such a homogeneous environment, results are identical for the four halfdiagonals. For h3 = 10^{8} and h4 = 10^{9}, only two and one day, respectively, are necessary for the virus to cover the halfdiagonals, the epidemic peak occurring 36 days (19 days, respectively) after virus introduction. A larger maximum of carrying capacity in vectors leads to larger local prevalences (Figure 5). Local prevalences are almost constant over different cells of the grid and are worth 78%, 90% and 94% for h2, h3 and h4 respectively. The total prevalence on the grid over time confirms these results. A larger maximum of carrying capacity in vectors leads to an earlier and larger epidemic peak (Figure 6).
Figure 5. Peak date and prevalence as a function of the distance from the primary case. When hypothesis H1 holds. Black line corresponds to the peak date, green line corresponds to the prevalence. Solid lines correspond to a maximum carrying capacity in vectors given by h4 = 10^{9}, dashed lines correspond to a maximum carrying capacity in vectors given by h3 = 10^{8}, dasheddotted lines correspond to a maximum carrying capacity in vectors given by h2 = 10^{7}.
Figure 6. Total prevalence over time. In green when hypothesis H1 holds, in black when hypothesis H2 holds. Solid lines correspond to the case where 90% of cells are occupied, dashed lines correspond to the case where 75% of cells are occupied, in dasheddotted lines correspond to the case where 50% of cells are occupied and dotted lines correspond to the case where 25% of cells are occupied.
Heterogeneous in abundance and distribution of hosts and homogeneous in abundance and distribution of vectors (H2)
For all maxima of carrying capacity in vectors (except h1 whatever the density of occupied cells and h2 for a density of occupied cells equal to 25%, for which there is no epidemic) a smaller density of occupied cells leads to a later epidemic peak and a smaller local prevalence in all grid cells. For example, for a density of occupied cells equal to 75% and a maximum carrying capacity of h2, the epidemic peak was 37 days later than in H1, occurring at day 133 in the central cell and at day 139 in the cell at the end of the halfdiagonal. Indeed, we kept the total number of hosts on the grid constant, so the number of hosts per occupied cell increases when the density of occupied cells decreases. As incidence functions are frequency dependent for hosts and vectors (Eq. 1), this increase in the number of hosts per occupied cell has little effect on the frequency of infection. However, the number of unoccupied cells increases, and therefore the virus transmission slows. The local prevalence is almost constant in all grid cells but is decreased by 9% compared with H1. For a maximum carrying capacity of h3 and h4, the epidemic peak was respectively 9 and 2 days later than in H1 (Figure 7). Similarly, the local prevalence is decreased by about 2% and 1%, respectively, compared with H1 (Figure 7).
Figure 7. Peak date and prevalence as a function of the distance from the primary case. When hypothesis H2 holds. Green star correspond to the case where 100% of cells are occupied, lozenge correspond to the case where 90% of cells are occupied, circle correspond to the case where 75% of cells are occupied, triangle correspond to the case where 50% of cells are occupied and square correspond to the case where 25% of cells are occupied.
The total prevalence on the grid over time confirms these results. A delay of the epidemic is observed as the density of hostoccupied cells decreases. There is also a decrease of the total prevalence compared with H1 (Figure 6).
Heterogeneous in abundance and distribution of vectors and homogeneous in abundance and distribution of hosts (H3, H4)
Hypothesis 3: definition of four subareas
Compared with H1, an epidemic peak is observed in all grid cells, even for the lowest maximum of carrying capacity, h1 (Figure 8). The further the distance from the grid centre, the more the epidemic peak is delayed and the more the local prevalence decreases. For the carrying capacity h2, the observed epidemic peaks are earlier than for H1, including in the cell at the end of the halfdiagonal. The local prevalence slightly decreases with increasing distance from the grid centre, to reach an equilibrium equal to the observed infection prevalence at the epidemic peak for H1 (Figure 8). For the two largest maxima of carrying capacity in vectors, h3 and h4, the peak dates and the local prevalences are similar to H1 results (Figure 8). This hypothesis highlights the influence of diffusion in each subarea on the diffusion in others.
Figure 8. Peak date and prevalence as a function of the distance from the primary case and the vector maximum carrying capacity. When hypothesis H3 holds (to help for the comparison, crosses correspond to the homogeneous case, hypothesis H1). Black lines correspond to a maximum carrying capacity in vectors given by h1 = 10^{6}, very dark grey lines correspond to a maximum carrying capacity in vectors given by h2 = 10^{7}, dark grey lines correspond to a maximum carrying capacity in vectors given by h3 = 10^{8} and light grey lines to a maximum carrying capacity in vectors given by h4 = 10^{9}.
Hypothesis 4: variable maximum of carrying capacities in vectors
Compared with H1, an epidemic peak is observed in all grid cells, even for the lowest maximum of carrying capacity, h1 (Figure 9). However, the peak dates are delayed. Indeed, epidemic peaks in cells having the same maximum of carrying capacity and equidistant ranged between 96 and 134 days after virus introduction, with no effect of distance on the peak date, (Figure 3, Figure 9). Therefore, cells near the introduction cell can be infected later than cells more distant (Figure 9). While the distribution of the maximum of carrying capacity in vectors is random, clusters of peak dates and of local prevalence arise (Figure 9). The local prevalences vary between 63% and 82% whatever the maximum of the carrying capacity and the distance from the introduction cell (Figure 9). A balance is observed between the peak dates and the local prevalences (Figure 9).
Figure 9. Spatial distribution of peak dates and prevalence. When hypothesis H4 holds. The black cell corresponds to the cell wherein the primary case occurs; rounds correspond to cells having the same maximum carrying capacity. From yellow to brown: earlier peak date and increasing prevalence.
The total prevalence on the grid over time confirms these results. A delay of the epidemic and a decrease in the total prevalence are observed compared with H1 (Figure 10).
Figure 10. Total prevalence over time. In green results when hypothesis H1 holds (same curves as in Figure 6, given here to help for comparison among hypotheses); solid lines correspond to a maximum carrying capacity in vectors given by h4 = 10^{9}, dashed lines correspond to a maximum carrying capacity in vectors given by h3 = 10^{8}, dotted lines correspond to a maximum carrying capacity in vectors given by h2 = 10^{7}. In thick black results when hypothesis H4 holds, in black results when hypothesis H5 holds; solid lines correspond to the case where 90% of cells are occupied, dashed lines correspond to the case where 75% of cells are occupied, dasheddotted lines correspond to the case where 50% of cells are occupied and dotted lines correspond to the case where 25% of cells are occupied.
Heterogeneous in abundance and distribution of hosts and vectors (H5, H6)
Hypothesis 5: theoretical landscape
As for H2, although later, peak dates are delayed as the density of hostoccupied cells decreases (Figure 11). Likewise, the local prevalence decreases when the density of occupied cells in hosts decreases. However, for a density of occupied cells equal to 25%, there is no epidemic; and for a density of occupied cells equal to 50%, the epidemic is small. By comparison with H2 and H4, the same tendencies are observed, with a delay in peak dates and a decrease in local prevalences. The distance has no effect on the peak date, unlike the maximum of carrying capacity in vectors. By comparing the distribution of the maximum of carrying capacity in vectors (Figure 3) with the peak dates and the local prevalences (Figure 11), a balance is observed. Cells with the highest maximum of carrying capacity have an earlier peak date and a larger local prevalence, and vice versa for cells with the lowest maximum of carrying capacity. This is true whatever the distance to the primary case.
Figure 11. Spatial distribution of peak dates and prevalence. When hypothesis H5 holds. a) 90% of occupied cells, b) 75% of occupied cells, c) 50% of occupied cells. The black cell corresponds to the cell wherein the primary case occurs; rounds correspond to cells having the same maximum carrying capacity. From yellow to brown: earlier peak date and increasing prevalence.
The total prevalence on the grid over time confirms these results. A delay in the epidemic is observed as the density of hostoccupied cells decreases, as well as a decrease in the total prevalence compared with H1 (Figure 7, Figure 10). A balance is observed between H2H5 and H4H5. As for H2 and for a maximum of carrying capacity in vectors of h2, the lower the density of hostoccupied cells, the longer the peak date is delayed and the lower is the total prevalence (Figure 7, Figure 10). However, the total prevalence for the highest density of occupied cells is similar with hypothesis H4 (Figure 10). Therefore, as soon as the distribution and the abundance of vectors are heterogeneous, they strongly influence the global epidemic dynamics.
Hypothesis 6: application to a real landscape
As in the theoretical case, for realistic distribution and abundance of hosts and vectors on a grid, the peak dates and the local prevalences are affected by the maximum of the carrying capacity in vectors in each cell (Figure 12). For roughly equal numbers of hosts and for a high maximum of carrying capacity in vectors (compared to other cells) (cells: A3 and F1, Figure 4), the peak dates and the local prevalences are close (e.g. cells A3 and F1; Figure 12). On the contrary, for roughly equal number of hosts but different maximum of carrying capacity in vectors (e.g. cells C4 and C5; Figure 4), cells with the highest maximum of carrying capacity in vectors have an earlier peak date and a higher local prevalence (e.g. cells C4 and C5; Figure 12). These trends are observed even if the observed number of vectors is multiplied by 100 or 1000, although the peak dates are earlier and the local prevalences higher. Indeed, in our mathematical model, multiplying the observed number of vectors by 100 or 1000 is equivalent to multiplying the maximum carrying capacity of vectors in each cell; that is, to divide k_{V} by 100 or 1000 in Eq. 1. This has a scaling effect on the model: it does not modify the time periodicity in the population dynamics, but it does change the abundance of vectors. This causes an earlier increase in the number of available infectious vectors and, therefore, the peak dates occur earlier together with higher local prevalences.
Figure 12. Spatial distribution of peak dates and prevalence in a real area. When hypothesis H6 holds. Multiplication of the number of trapped vectors by: a) 100; b) 1000. The cell with the star corresponds to the cell wherein the primary case occurs. From yellow to brown: earlier peak date and increasing prevalence.
In addition, the number of hosts in each cell has an impact on the peak date and the local prevalence. For equidistant cells having the same maximum of the carrying capacity in vectors (e.g. cells A2 and C1; Figure 4a), the cell having the largest number of hosts (C1) shows a later peak date and a lower local prevalence than the cell having the lowest number of hosts (A2; Figure 12).
Discussion
A mathematical modelling approach allowed us to assess the impact of spatiotemporal heterogeneities in abundance and distribution of hosts and vectors on the spatiotemporal spread of BTV8. Individually and jointly, the heterogeneities in abundance and distribution of hosts and vectors delay the peak date and decrease the total infection prevalence. The different hypotheses of heterogeneity that we have tested allowed us to highlight the importance of the maximum of the carrying capacity in vectors and its influence on the spread of BTV8 within each cell. Indeed, cells next to the primary case can become infected later than more distant cells, if the maximum of carrying capacity is lower. Moreover, the density of cells occupied by hosts plays an important role in cases where the maximum carrying capacity in vectors is low (for homogeneous conditions for vectors) and when hosts and vectors are heterogeneous.
The spatial heterogeneity in host and/or vector abundances influences the infection frequency [7,17,29]. In our study, a decrease in the density of cells occupied by hosts results in a delay of the peak date and a decrease in the infection prevalence if vectors are homogeneously distributed and their population is large. In the case where hosts are homogeneously distributed, the same trend is observed for cells where the maximum of carrying capacity in vectors is large. However, for cells where this maximum is the lowest, an epidemic can occur, in contrast to the case with homogeneous vector and host populations. The coupling of heterogeneities in hosts and vectors increases the delay of the epidemic and decreases the prevalence.
Different models of the spatial spread of bluetongue have previously been published for specific geographic areas [2124]. Szmaragd’s model describes the BTV spread within and between farms in Great Britain via a generic kernel, which includes both animal and vector movements. Ducheyne’s model was calibrated with data from the BTV1 and BTV8 epidemics in Southern France. It describes the spatiotemporal BTV spread between farms assuming that the number of new cases per week is half attributable to local dispersion (active) of vectors, and half to longdistance dispersion (passive) of vectors by the wind. Graesboll’s model describes the BTV spread with a high spatial resolution, which includes both animal and vector movements and the seasonality of vectors. Turner’s model is a network model. It takes into account explicitly the spatial dispersal of both hosts and vectors and a seasonal vector to host ratio [24]. It studies the BTV spread between farms in England. Taking into account climatic and environmental data, all of these models consider the spatial heterogeneity of the landscape. Our model advances the field by representing spatial heterogeneity in both hosts and vectors. Here, we highlighted how such heterogeneities concretely impact BTV spread. Moreover, the seasonality of the vector population is managed by a simple function that can be easily related to observed data of vector abundance.
We chose to model the spatiotemporal BTV8 spread by a reaction–diffusion model. Several shapes of the transmission kernel are possible, but it is difficult to choose the most appropriate one to describe the dynamics of a pathogen spread. Indeed, Szmaragd et al. showed that a Gaussian kernel was the most appropriate to describe the BTV8 spread in northern Europe during 2006 [21]. If this kernel shape, comparable to reaction diffusion models, has been shown to underestimate the probability of the longdistance transmission, and thus is not appropriate to describe the pathogen spread on a larger scale [30], it can be used on a smaller scale. Graesboll et al. used a Gaussian kernel too, but coupled this approach with the wind dispersion [23]. In our study, the theoretical spatial domain is a 41 × 41 km square. The primary case is always located in the centre of the square, i.e. at 29 km from the domain edges. One limitation is that longdistance dispersal has been neglected; the wind dispersal responsible for the passive movements of vectors generating dispersal up to several hundred kilometres [31,32]. Coupling short and longdistance dispersal is necessary to study arbovirus spread in animal populations once the spatial scale is large enough that host movements and passive movements of vectors cannot be neglected anymore [3335].
Observational studies have been conducted to assess risk areas and to predict the spatiotemporal spread of BTV in BTVfree areas [11,36,37]. These studies have shown that the landscape heterogeneity, climatic conditions, the distribution and the density of host populations and the abundance of vectors were linked and influenced BTV spatiotemporal spread. Our model highlights similar findings but also allows us to distinguish between the impact of vector versus host heterogeneity. Maps of the basic reproduction number have highlighted the link between vector abundance and BTV spread [12]. However, the vector abundance is difficult to quantify. Hartemink et al. have used trapping data, multiplying the number of trapped Culicoides by 100 to obtain a local density of vectors [12]. Our real geographic area illustrates these differences in local abundance. On a small scale, large differences may exist between cells, whether they are occupied by hosts or not. Entomological studies identify and quantify the different vector species present in different geographic locations [9,10,38]. However, the real number of vectors remains difficult to approximate. As shown in our results, the abundance in vectors has a significant impact on the date and on the observed prevalence at the epidemic peak. However, by multiplying the number of vectors by 100 or 1000, we obtained the same qualitative findings on the real geographic area.
Modelling is a relevant approach to investigate the concept of spatiotemporal heterogeneities on the dynamics of virus spread. The distribution of hosts and vectors, and vector abundance strongly influence the dynamics of BTV spread. The application of our model on a real geographic area, although of limited size, allowed us to illustrate the conclusions drawn from a theoretical domain. The reaction–diffusion models are classically used in plant epidemiology [3941], with the modelled movements generally being for highly volatile entities, and the short versus longdistance movements being taken into account via different diffusion coefficients [40,42]. Although developed to represent BTV8 spatiotemporal spread at a local to regional scale, our model can be used to study other vectorborne diseases of animals and its extension to a larger area remains possible.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
MC and PE conceived of the study, carried out the model development and analysis, and drafted the manuscript. ML and HS conceived of the study, participated in its design and coordination, and drafted the manuscript. MB and GK provided data used in this research, and both commented on draft and final manuscripts. All authors read and approved the final manuscript.
Acknowledgements
Financial support for this research was provided by INRA, IRSTEA and BasseNormandie, Bretagne, Pays de la Loire and PoitouCharentes Regional Councils under SANCRE project, in the framework of “For and About Regional Development” programs. GK was supported by a BBSRC DTGfunded PhD studentship awarded to MB and Dr Jon Read.
References

Chomel BB, Belotto A, Meslin FX: Wildlife, exotic pets, and emerging zoonoses.
Emerg Infect Dis 2007, 13:611. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Jones KE, Patel NG, Levy MA, Storeygard A, Balk D, Gittleman JL, Daszak P: Global trends in emerging infectious diseases.
Nature 2008, 451:990993. PubMed Abstract  Publisher Full Text

Roche B, Guegan JF: Ecosystem dynamics, biological diversity and emerging infectious diseases.
C R Biol 2011, 334:385392. PubMed Abstract  Publisher Full Text

Gould EA, Higgs S: Impact of climate change and other factors on emerging arbovirus diseases.
Trans R Soc Trop Med Hyg 2009, 103:109121. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Velthuis AG, Saatkamp HW, Mourits MC, de Koeijer AA, Elbers AR: Financial consequences of the Dutch bluetongue serotype 8 epidemics of 2006 and 2007.
Prev Vet Med 2009, 93:294304. PubMed Abstract  Publisher Full Text

Velthuis AGJ, Mourits MCM, Saatkamp HW, de Koeijer AA, Elbers ARW: Financial evaluation of different vaccination strategies for controlling the bluetongue virus serotype 8 epidemic in The Netherlands in 2008.
PLoS One 2011, 6:e19612. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Smith DL, Dushoff J, McKenzie FE: The risk of a mosquitoborne infection in a heterogeneous environment.

Backer JA, Nodelijk G: Transmission and control of African horse sickness in The Netherlands: a model analysis.
PLoS One 2011, 6:e23066. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Foxi C, Delrio G: Larval habitats and seasonal abundance of Culicoides biting midges found in association with sheep in northern Sardinia, Italy.
Med Vet Entomol 2010, 24:199209. PubMed Abstract  Publisher Full Text

Ander M, Meiswinkel R, Chirico J: Seasonal dynamics of biting midges (Diptera: Ceratopogonidae: Culicoides), the potential vectors of bluetongue virus, in Sweden.
Vet Parasitol 2012, 184:5967. PubMed Abstract  Publisher Full Text

Guis H, Tran A, de la Rocque S, Baldet T, Gerbier G, Barragué B, BiteauCoroller F, Roger F, Viel JF, Mauny F: Use of high spatial resolution satellite imagery to characterize landscapes at risk for bluetongue.
Vet Res 2007, 38:669683. PubMed Abstract  Publisher Full Text

Hartemink NA, Purse BV, Meiswinkel R, Brown HE, de Koeijer A, Elbers ARW, Boender GJ, Rogers DJ, Heesterbeek JAP: Mapping the basic reproduction number (R_{0}) for vectorborne diseases: a case study on bluetongue virus.
Epidemics 2009, 1:153161. PubMed Abstract  Publisher Full Text

Raffy M, Tran A: On the dynamics of flying insects populations controlled by large scale information.
Theor Popul Biol 2005, 68:91104. PubMed Abstract  Publisher Full Text

Tran A, Raffy M: On the dynamics of dengue epidemics from largescale information.
Theor Popul Biol 2006, 69:312. PubMed Abstract  Publisher Full Text

Lewis M, Rencławowicz J, van den Driessche P: Traveling waves and spread rates for a West Nile virus model.
Bull Math Biol 2006, 68:323. PubMed Abstract  Publisher Full Text

Maidana NA, Yang HM: Spatial spreading of West Nile Virus described by traveling waves.
J Theor Biol 2009, 258:403417. PubMed Abstract  Publisher Full Text

Ostfeld RS, Glass GE, Keesing F: Spatial epidemiology: an emerging (or reemerging) discipline.
Trends Ecol Evol 2005, 20:328336. PubMed Abstract  Publisher Full Text

Adams B, Kapan DD: Man bites mosquito: understanding the contribution of human movement to vectorborne disease dynamics.
PLoS One 2009, 4:e6763. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chen Q, Han R, Ye F, Li W: Spatiotemporal ecological models.
Ecol Inform 2011, 6:3743. Publisher Full Text

Ferreira LC, Udia , Pulino P, Yang , Hyun , Takahashi L: Controlling dispersal dynamics of Aedes aegypti.
Math Popul Stud 2006, 13:215236. Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ducheyne E, Lange M, Van der Stede Y, Meroc E, Durand B, Hendrickx G: A stochastic predictive model for the natural spread of bluetongue.
Prev Vet Med 2011, 99:4859. PubMed Abstract  Publisher Full Text

Graesboll K, Bodker R, Enoe C, Christiansen LE: Simulating spread of Bluetongue Virus by flying vectors between hosts on pasture.
Sci Rep 2012, 2:863. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Turner J, Bowers RG, Baylis M: Modelling bluetongue virus transmission between farms using animal and vector movements.
Sci Rep 2012, 2:319. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kluiters G, Sugden D, Guis H, McIntyre KM, Labuschagne K, Vilar MJ, Baylis M: Modelling the spatial distribution of Culicoides biting midges at the local scale.

Charron MVP, Seegers H, Langlais M, Ezanno P: Seasonal spread and control of Bluetongue in cattle.
J Theor Biol 2011, 291:19. PubMed Abstract  Publisher Full Text

Van Ark H, Meiswinkel R: Subsampling of large light trap catches of Culicoides (Diptera: Ceratopogonidae).
Onderstepoort J Vet Res 1992, 59:183189. PubMed Abstract

Pioz M, Guis H, Calavas D, Durand B, Abrial D, Ducrot C: Estimating frontwave velocity of infectious diseases: a simple, efficient method applied to bluetongue.
Vet Res 2011, 42:60. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Kilpatrick AM, Daszak P, Jones MJ, Marra PP, Kramer LD: Host heterogeneity dominates West Nile virus transmission.
Proc Biol Sci 2006, 273:23272333. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

De Koeijer A, Boender G, Nodelijk G, Staubach C, Meroc E, Elbers A: Quantitative analysis of transmission parameters for bluetongue virus serotype 8 in Western Europe in 2006.
Vet Res 2011, 42:53. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Sellers RF: Weather, Culicoides, and the distribution and spread of bluetongue and African horse sickness viruses.
Bluetongue, African horse sickness, and related orbiviruses: Proceedings of the Second International Symposium 1992, 284290.

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 northwestern Europe bluetongue epidemic, 2006.
Prev Vet Med 2008, 87:162181. PubMed Abstract  Publisher Full Text

Murray MD: Akabane epizootics in New South Wales: evidence for longdistance dispersal of the biting midge Culicoides brevitarsis.
Aust Vet J 1987, 64:305308. PubMed Abstract  Publisher Full Text

Mellor PS, Boorman J, Baylis M: Culicoides biting midges: their role as arbovirus vectors.
Annu Rev Entomol 2000, 45:307340. PubMed Abstract  Publisher Full Text

Finlaison DS, Read AJ, Kirkland PD: An epizootic of bovine ephemeral fever in New South Wales in 2008 associated with longdistance dispersal of vectors.
Aust Vet J 2010, 88:301306. PubMed Abstract  Publisher Full Text

Calvete C, Estrada R, Miranda MA, Borrás D, Calvo JH, Lucientes J: Ecological correlates of bluetongue virus in Spain: predicted spatial occurrence and its relationship with the observed abundance of the potential Culicoides spp. vector.
Vet J 2009, 182:235243. PubMed Abstract  Publisher Full Text

Sanders CJ, Shortall CR, Gubbins S, Burgin L, Gloster J, Harrington R, Reynolds DR, Mellor PS, Carpenter S: Influence of season and meteorological parameters on flight activity of Culicoides biting midges.
J Appl Ecol 2011, 48:13551364. Publisher Full Text

Balenghien T, Fouque F, Sabatier P, Bicout DJ: Horse, bird, and humanseeking behavior and seasonal abundance of mosquitoes in a West Nile virus focus of southern France.
J Med Entomol 2006, 43:936946. PubMed Abstract  Publisher Full Text

Gilligan CA: Modelling soilborne plant pathogens: reaction–diffusion models.
Can J Plant Pathol 1995, 17:96108. Publisher Full Text

Burie JB, Calonnec A, Langlais M: Modeling of the invasion of a fungal disease over a Vineyard. In Mathematical Modeling of Biological Systems. Volume 2. Boston: Birkhäuser; 2008::1121. [Modeling and Simulation in Science, Engineering and Technology]

Sapoukhina N, Tyutyunov Y, Sache I, Arditi R: Spatially mixed crops to control the stratified dispersal of airborne fungal diseases.
Ecol Model 2010, 221:27932800. Publisher Full Text

Zawolek MW, Zadoks JC: Studies in focus development  an optimum for the dual dispersal of plant pathogens.
Phytopathology 1992, 82:12881297. Publisher Full Text

Singer RS, MacLachlan NJ, Carpenter TE: Maximal predicted duration of viremia in bluetongue virusinfected cattle.
J Vet Diagn Invest 2001, 13:4349. PubMed Abstract  Publisher Full Text

Baylis M, O’Connell L, Mellor PS: Rates of bluetongue virus transmission between Culicoides sonorensis and sheep.
Med Vet Entomol 2008, 22:228237. PubMed Abstract  Publisher Full Text

Epidemiological analysis of 2006 bluetongue virus serotype 8 epidemic in northwestern Europe 2007. 2007.
apart from the link http://www.efsa.europa.eu/fr/efsajournal/doc/34r.pdf webcite

Gerry AC, Mullens BA, Maclachlan NJ, Mecham JO: Seasonal transmission of bluetongue virus by Culicoides sonorensis (Diptera: Ceratopogonidae) at a southern california dairy and evaluation of vectorial capacity as a predictor of bluetongue virus transmission.
J Med Entomol 2001, 38:197209. PubMed Abstract  Publisher Full Text

Carpenter S, Lunt HL, Arav D, Venter GJ, Mellor PS: Oral susceptibility to bluetongue virus of Culicoides (Diptera: Ceratopogonidae) from the United Kingdom.
J Med Entomol 2006, 43:7378. PubMed Abstract  Publisher Full Text