# Seasonal and spatial heterogeneities in host and vector abundances impact the spatiotemporal spread of bluetongue

Maud Charron, Georgette Kluiters, Michel Langlais, Henry Seegers, Matthew Bailys, Pauline Ezanno

► To cite this version:

Maud Charron, Georgette Kluiters, Michel Langlais, Henry Seegers, Matthew Bailys, et al.. Seasonal and spatial heterogeneities in host and vector abundances impact the spatiotemporal spread of bluetongue. BMC Veterinary Research, 2013, 44 (1), pp.1-12. 10.1186/1297-9716-44-44. hal-00992664

HAL Id: hal-00992664

https://hal.science/hal-00992664

Submitted on 29 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH

Open Access

# Seasonal and spatial heterogeneities in host and vector abundances impact the spatiotemporal spread of bluetongue

Maud VP Charron<fnref n="1" />,<fnref n="2" />,³,⁴*, Georgette Kluiters⁵, Michel Langlais³,⁴, Henri Seegers<fnref n="1" />,<fnref n="2" />, Matthew Baylis⁵
and Pauline Ezanno<fnref n="1" />,<fnref n="2" />

## 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 wind-induced 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 vector-borne diseases of animals with a local to regional spread by vector diffusion.

## Introduction

There is significant concern regarding the resurgence or emergence of vector-borne diseases of animals with serious consequences for animal health and economics [1-3]. Climate change and socio-economic 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 non-contagious vector-borne 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<fnref n="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],

<footnote>* Correspondence: maud.charron@oniris-nantes.fr</footnote>

<footnote n="1">¹INRA, UMR1300 Biologie, Epidémiologie et Analyse de Risques en santé animale, CS 40706, F-44307 Nantes, France</footnote>

<footnote n="2">²LUNAM Université, Oniris, Ecole nationale vétérinaire, agroalimentaire et de l'alimentation Nantes-Atlantique, UMR BioEpAR, F-44307 Nantes, France
Full list of author information is available at the end of the article</footnote>
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 vector-borne diseases [8,13-16]. As hosts and vectors are mobile entities, a spatial component in vector-borne disease models should be taken into account to better consider the evolution of the biological system [17-19]. 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 vector-borne 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 [21-24]. 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) [21-24]. 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/\rho_{Ei}$, Table 1) was taken into account as it is of the same order as life expectancy. At the disease-free 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 density-dependant 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).

<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>
<table><caption>**Table 1 Parameters of the model of BTV8 spread in midge and cattle populations**</caption><thead><tr><th>Host parameters</th><th>Description</th><th>Value</th><th>References</th></tr></thead><tbody><tr><td>b<sub>H</sub></td><td>Birth rate (per day)</td><td>6,94.10<sup>-4</sup></td><td>a</td></tr><tr><td>m<sub>H</sub></td><td>Exit rate (mortality, selling, culling) (per day)</td><td>b<sub>H</sub></td><td></td></tr><tr><td>1/α<sub>I</sub></td><td>Duration of viremia (days)</td><td>60</td><td>[43]</td></tr><tr><td>c<sub>VH</sub></td><td>Probability of transmission from vector to host</td><td>0.92</td><td>[44,45]</td></tr><tr><th>Vector parameters</th><th>Description</th><th>Value</th><th>References</th></tr><tr><td>c<sub>HV</sub></td><td>Probability of transmission from host to vector</td><td>0.15</td><td>[46,47]</td></tr><tr><td>n</td><td>Biting rate (per day)</td><td>0.25</td><td>[45]</td></tr><tr><td>b<sub>V</sub></td><td>Fertility rate (per day)</td><td>6.1</td><td>[45]</td></tr><tr><td>m<sub>V</sub></td><td>Mortality rate (per day)</td><td>1/21</td><td>[34,45]</td></tr><tr><td>K(t)</td><td>Carrying capacity</td><td>*</td><td></td></tr><tr><td>k<sub>V</sub></td><td>Density-dependence mortality rate (per day)</td><td>(b<sub>V</sub>-m<sub>V</sub>)/K (t)</td><td></td></tr><tr><td>h</td><td>Maximum of K(t)</td><td>variable</td><td></td></tr><tr><td>d</td><td>Duration of favourable period (days)</td><td>243</td><td></td></tr><tr><td>Nb</td><td>Number of vectors during the unfavourable period</td><td>100</td><td></td></tr><tr><td>1/ρ<sub>E</sub></td><td>Duration of extrinsic incubation period (days)</td><td>10</td><td>[4,45]</td></tr><tr><td>D</td><td>Diffusion coefficient (km<sup>2</sup>/day)</td><td>1,25.10<sup>-2</sup></td><td></td></tr></tbody></table>

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.

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:

$\vec{q}(x, y, t) = -D(x, y) \nabla SV(x, y, t)$, where $\vec{q}(x, y, t)$ is the diffusion flux of the vector population and $D(x, y)$, the positive diffusion matrix. Applying the conservation law (i.e. the variation of the amount of species in a volume is equal to the balance of entering and outgoing flux), we obtain

the second Fick's law: $\frac{\partial SV(x,y,t)}{\partial t} + div \vec{q}(x, y, t) = 0$. Combining these two laws, we obtain: $\frac{\partial SV(x,y,t)}{\partial t} = -div \vec{q}(x, y, t) = -div(-D(x, y) \nabla SV(x, y, t))$. We consider that the dispersion of vectors is homogeneous in space; therefore $D(x, y) = D$. Then, $\frac{\partial SV(x,y,t)}{\partial t} = D \div (\nabla SV(x, y, t)) = D \Delta SV(x, y, t)$.

By adding the reaction term, we obtain the following system of equations (Eq. 1) describing the spatiotemporal spread of the BTV8, for $(x, y) \in \Omega$ and $t > 0$:

$$ 
\left\{
\begin{aligned}
\frac{\partial SH}{\partial t} &= -c_{VH}n \frac{IV}{HP} SH + b_H HP - m_H SH \\
\frac{\partial IH}{\partial t} &= c_{VH}n \frac{IV}{HP} SH - (\alpha_I + m_H) IH \\
\frac{\partial RH}{\partial t} &= \alpha_I IH - m_H RH \\
\frac{\partial SV}{\partial t} &= D \Delta SV - c_{HV}n \frac{IH}{HP} SV - (m_V + k_V VP) SV + b_V VP \\
\frac{\partial EV}{\partial t} &= D \Delta EV + c_{HV}n \frac{IH}{HP} SV - (m_V + k_V VP) EV - \rho_E EV \\
\frac{\partial IV}{\partial t} &= D \Delta IV + \rho_E EV - (m_V + k_V VP) IV
\end{aligned}
\right.
\quad (1) $$

The initial time corresponds to the first day of the favourable season for the vector population. The spatial domain $(\Omega)$ is discretized into cells; each cell of surface area of 1 km². 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 non-negative 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).

$$ 
\left\{
\begin{array}{l}
SH(x, y, 0) = SH^0(x, y), SV(x, y, 0) = SV^0(x, y) \\
IH(x, y, 0) = 0 \quad \forall (x, y) \neq (21, 21); IH(21, 21, 0) = 1 \\
RH(x, y, 0) = EV(x, y, 0) = IV(x, y, 0) = 0 \\
\frac{\partial SH}{\partial n} = \frac{\partial IH}{\partial n} = \frac{\partial RH}{\partial n} = \frac{\partial SV}{\partial n} = \frac{\partial EV}{\partial n} = \frac{\partial IV}{\partial n} = 0, \forall (x, y) \in \partial\Omega, t > 0
\end{array}
\right.
\quad (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 semi-implicit 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).

<footnote>* K(t) = 1/[1 + α(t) * h * sin([π(365-t)/d])] + 1/[d+1,365](t) * Nb.</footnote>

<footnote>a: agricultural statistics.</footnote>
<table><caption>Table 2 Hypotheses of spatial heterogeneities in abundance and distribution of hosts and vectors</caption><thead><tr><th></th><th>H</th><th>Host</th><th>Vectors</th><th>Results figure</th></tr></thead><tbody><tr><td>Homogeneous in hosts and vectors (4 Scenarios)</td><td>H1</td><td>500 S/C 1 l in central cell Total number of hosts = 840 501</td><td>100 S/C h1=10<sup>6</sup>, or h2=10<sup>7</sup> or h3=10<sup>8</sup>, or h4= 10<sup>9</sup></td><td>5</td></tr><tr><td>Heterogeneous in hosts Homogeneous in vectors (16 scenarios)</td><td>H2</td><td>Total number of hosts S≈ 840 501 Four densities of occupied cells (Figure 2): 90% OC: 554 S/OC75% OC: 666 S/OC50% OC: 1000 S/OC25% OC: 2025 S/OC 1 l in central cell</td><td>H1</td><td>6 and 7</td></tr><tr><td>Heterogeneous in vectors Homogeneous in hosts (2 scenarios)</td><td>H3</td><td>H1</td><td>100 S/C Grid divided into four sub-areas of different maximum carrying capacities (h1, h2, h3 and h4) (Figure 3a)</td><td>8</td></tr><tr><td></td><td>H4</td><td></td><td>100 S/C 25% C: h1, 25% C: h2, 25% C: h3, 25% C: h4 (Figure 3b)</td><td>9 and 10</td></tr><tr><td>Heterogeneous in hosts and vectors (6 scenarios)</td><td>H5</td><td>Crossing hypotheses: H2-H4</td><td>Real area (Figure 4a)</td><td>10 and 11</td></tr><tr><td></td><td>H6</td><td>Real area (Figure 4b)</td><td>Multiplication of the number of trapped vectors by 100 or 1000</td><td>12</td></tr></tbody></table>

We set *H* hypothesis, *C* cell, *OC* occupied cell, *S* susceptible individual, *I* infected individual.

From this assumption, four levels of heterogeneity were analyzed and compared.

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).

Secondly, a heterogeneous distribution of vectors was considered (H3 and H4), hosts being homogeneously distributed. In H3, the domain (Ω) consisted of four

sub-areas, 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).

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

<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

**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.

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 Onderstepoort-type down draught black light traps positioned as close to livestock as possible. Large collections were sub-sampled [27] and only females were considered in the analyses as males do not take blood meals or, consequently, transmit disease between vertebrates. The maximum trap-catch 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.

## 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 half-diagonal), 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 half-diagonal 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 half-diagonal 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 half-diagonals at 96 and 101 days respectively after virus introduction. As expected in such a homogeneous environment, results are identical for the four half-diagonals. For $h3 = 10^8$ and $h4 = 10^9$, only two

<figure><img src="image_6.png" /></figure>

**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.
<figure><img src="image_7.png" /></figure>

**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$, dashed-dotted lines correspond to a maximum carrying capacity in vectors given by $h2 = 10^7$.

and one day, respectively, are necessary for the virus to cover the half-diagonals, 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).

## 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

<figure><img src="image_8.png" /></figure>

**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 dashed-dotted lines correspond to the case where 50% of cells are occupied and dotted lines correspond to the case where 25% of cells are occupied.
<figure><img src="image_9.png" /></figure>

equal to 75% and a maximum carrying capacity of *h*2, 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 half-diagonal. 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 *h*3 and *h*4, 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).

The total prevalence on the grid over time confirms
these results. A delay of the epidemic is observed as the

density of host-occupied 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
distances from the grid centre, the more the epidemic peak is delayed and
the more the local prevalence decreases. For the carrying
capacity *h*2, the observed epidemic peaks are earlier than
for H1, including in the cell at the end of the half-
diagonal. The local prevalence slightly decreases with
increasing distance from the grid centre, to reach an equi-
librium equal to the observed infection prevalence at the
epidemic peak for H1 (Figure 8). For the two largest max-
ima of carrying capacity in vectors, *h*3 and *h*4, the peak
<figure><img src="image_10.png" /><figcaption>**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 *h*1 = 10⁶, very dark grey lines correspond to a maximum carrying capacity in vectors given by *h*2 = 10⁷, dark grey lines correspond to a maximum carrying capacity in vectors given by *h*3 = 10⁸ and light grey lines to a maximum carrying capacity in vectors given by *h*4 = 10⁹.</figcaption></figure>

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.

**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, *h*1 (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).

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).

**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 host-occupied 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.

The total prevalence on the grid over time confirms these results. A delay in the epidemic is observed as the density of host-occupied cells decreases, as well as a decrease in the total prevalence compared with H1 (Figure 7, Figure 10). A balance is observed between H2-H5 and H4-H5. As for H2 and for a maximum of carrying capacity in vectors of *h*2, the lower the density of host-occupied 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.
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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.

<figure><img src="image_12.png" /><figcaption>**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 $h_4 = 10^5$, dashed lines correspond to a maximum carrying capacity in vectors given by $h_3 = 10^5$, dotted lines correspond to a maximum carrying capacity in vectors given by $h_2 = 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, dashed-dotted lines correspond to the case where 50% of cells are occupied and dotted lines correspond to the case where 25% of cells are occupied.</figcaption></figure>

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.

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 [21-24]. 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 long-distance 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 long-distance transmission, and thus is not appropriate to describe the pathogen spread on a larger scale [30], it can be

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>
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 long-distance 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 long-distance 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 [33-35].

Observational studies have been conducted to assess risk areas and to predict the spatiotemporal spread of BTV in BTV-free 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 [39-41], with the modelled movements generally being for highly volatile entities, and the short *versus* long-distance 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 vector-borne

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 Basse-Normandie, Bretagne, Pays de la Loire and Poitou-Charentes Regional Councils under SANCRE project, in the framework of "For and About Regional Development" programs. GK was supported by a BBSRC DTG-funded PhD studentship awarded to MB and Dr Jon Read.

**Author details**

¹INRA, UMR1300 Biologie, Epidémiologie et Analyse de Risques en santé animale, CS 40706, F-44307 Nantes, France. ²LUNAM Université, Oniris, Ecole nationale vétérinaire, agroalimentaire et de l'alimentation Nantes-Atlantique, UMR BioEpAR, F-44307 Nantes, France. ³Institut de Mathématiques de Bordeaux, UMR 5251, Université de Bordeaux, F33076 Bordeaux, France.

⁴CNRS, IMB, UMR 5251, F-33400 Talence, France. ⁵Institute of Infection and Global Health, Liverpool University Climate and Infectious Diseases of Animals (Lucinda) Group, Leahurst Campus, University of Liverpool, Neston CH64 7TE, UK.

Received: 6 February 2013 Accepted: 7 June 2013

Published: 19 June 2013

**References**

1. Chomel BB, Belotto A, Meslin FX: Wildlife, exotic pets, and emerging zoonoses. *Emerg Infect Dis* 2007, 13:6–11.

2. Jones KE, Patel NG, Levy MA, Storeygard A, Balk D, Gittleman JL, Daszak P: Global trends in emerging infectious diseases. *Nature* 2008, 451:990–993.

3. Roche B, Guegan JF: Ecosystem dynamics, biological diversity and emerging infectious diseases. *C R Biol* 2011, 334:385–392.

4. Gould EA, Higgs S: Impact of climate change and other factors on emerging arbovirus diseases. *Trans R Soc Trop Med Hyg* 2009, 103:109–121.

5. 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:294–304.

6. 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.

7. Smith DL, Dushoff J, McKenzie FE: The risk of a mosquito-borne infection in a heterogeneous environment. *PLoS Biol* 2004, 2:1957–1964.

8. Backer JA, Nodelijk G: Transmission and control of African horse sickness in The Netherlands: a model analysis. *PLoS One* 2011, 6:e23066.

9. 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:199–209.

10. 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:59–67.

11. Guis H, Tran A, de la Rocque S, Baldet T, Gerbier G, Barragué B, Biteau-Coroller F, Roger F, Viel J-F, Mauny F: Use of high spatial resolution satellite imagery to characterize landscapes at risk for bluetongue. *Vet Res* 2007, 38:669–683.

12. 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₀) for vector-borne diseases: a case study on bluetongue virus. *Epidemics* 2009, 1:153–161.

13. Raffy M, Tran A: On the dynamics of flying insects populations controlled by large scale information. *Theor Popul Biol* 2005, 68:91–104.
14. Tran A, Raffy M: On the dynamics of dengue epidemics from large-scale information. *Theor Popul Biol* 2006, 69:3-12.

15. Lewis M, Renclawowicz J, van den Driessche P: Traveling waves and spread rates for a West Nile virus model. *Bull Math Biol* 2006, 68:3-23.

16. Maidana NA, Yang HM: Spatial spreading of West Nile Virus described by traveling waves. *J Theor Biol* 2009, 258:403-417.

17. Ostfeld RS, Glass GE, Keesing F: Spatial epidemiology: an emerging (or re-emerging) discipline. *Trends Ecol Evol* 2005, 20:328-336.

18. Adams B, Kapan DD: Man bites mosquito: understanding the contribution of human movement to vector-borne disease dynamics. *PLoS One* 2009, 4:e6763.

19. Chen Q, Han R, Ye F, Li W: Spatio-temporal ecological models. *Ecol Inform* 2011, 6:37-43.

20. Ferreira LC, Udia, Pulino P, Yang, Hyun, Takahashi L: Controlling dispersal dynamics of *Aedes aegypti*. *Math Popul Stud* 2006, 13:215-236.

21. 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.

22. 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:48-59.

23. 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.

24. Turner J, Bowers RG, Baylis M: Modelling bluetongue virus transmission between farms using animal and vector movements. *Sci Rep* 2012, 2:319.

25. 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. *J Appl Ecol* 2012, 50:232-242.

26. Charron MVP, Seegers H, Langlais M, Ezanno P: Seasonal spread and control of Bluetongue in cattle. *J Theor Biol* 2011, 291:1-9.

27. Van Ark H, Meiswinkel R: Subsampling of large light trap catches of *Culicoides* (Diptera: *Ceratopogonidae*). *Onderstepoort J Vet Res* 1992, 59:183-189.

28. Pioz M, Guis H, Calavas D, Durand B, Abrial D, Ducrot C: Estimating front-wave velocity of infectious diseases: a simple, efficient method applied to bluetongue. *Vet Res* 2011, 42:60.

29. Kilpatrick AM, Daszak P, Jones MJ, Marra PP, Kramer LD: Host heterogeneity dominates West Nile virus transmission. *Proc Biol Sci* 2006, 273:2327-2333.

30. 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.

31. Sellers RF: Weather, Culicoides, and the distribution and spread of bluetongue and African horse sickness viruses. In *Bluetongue, African horse sickness, and related orbiviruses: Proceedings of the Second International Symposium*; 1992:284-290.

32. 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.

33. Murray MD: Akabane epizootics in New South Wales: evidence for long-distance dispersal of the biting midge *Culicoides brevitaris*. *Aust Vet J* 1987, 64:305-308.

34. Mellor PS, Boorman J, Baylis M: Culicoides biting midges: their role as arbovirus vectors. *Annu Rev Entomol* 2000, 45:307-340.

35. Finlaison DS, Read AJ, Kirkland PD: An epizootic of bovine ephemeral fever in New South Wales in 2008 associated with long-distance dispersal of vectors. *Aust Vet J* 2010, 88:301-306.

36. 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:235-243.

37. 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:1355-1364.

38. Balenghien T, Fouque F, Sabatier P, Bicout DJ: Horse-, bird-, and human-seeking behavior and seasonal abundance of mosquitoes in a West Nile virus focus of southern France. *J Med Entomol* 2006, 43:936-946.

39. Gilligan CA: Modelling soil-borne plant pathogens: reaction-diffusion models. *Can J Plant Pathol* 1995, 17:96-108.

40. Burie J-B, Calonnec A, Langlais M: Modeling of the invasion of a fungal disease over a Vineyard. In *Mathematical Modeling of Biological Systems, Modeling and Simulation in Science, Engineering and Technology, Volume 2*. Boston: Birkhäuser; 2008:11-21.

41. Sapoukhina N, Tyutyunov Y, Sache I, Arditi R: Spatially mixed crops to control the stratified dispersal of airborne fungal diseases. *Ecol Model* 2010, 221:2793-2800.

42. Zawolek MW, Zadoks JC: Studies in focus development - an optimum for the dual dispersal of plant pathogens. *Phytopathology* 1992, 82:1288-1297.

43. Singer RS, MacLachlan NJ, Carpenter TE: Maximal predicted duration of viremia in bluetongue virus-infected cattle. *J Vet Diagn Invest* 2001, 13:43-49.

44. Baylis M, O'Connell L, Mellor PS: Rates of bluetongue virus transmission between *Culicoides sonorensis* and sheep. *Med Vet Entomol* 2008, 22:228-237.

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

46. Gerry AC, Mullens BA, Maclachlan NJ, Mechem 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:197-209.

47. 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:73-78.

doi:10.1186/1297-9716-44-44

Cite this article as: Charron et al.: Seasonal and spatial heterogeneities in host and vector abundances impact the spatiotemporal spread of bluetongue. Veterinary Research 2013 44:44.

Submit your next manuscript to BioMed Central and take full advantage of:

* Convenient online submission

* Thorough peer review

* No space constraints or color figure charges

* Immediate publication on acceptance

* Inclusion in PubMed, CAS, Scopus and Google Scholar

* Research which is freely available for redistribution

Submit your manuscript at
www.biomedcentral.com/submit

BioMed Central