Bluetongue Disease Risk Assessment Based on Observed and Projected Culicoides obsoletus spp. Vector Densities

Bluetongue is an arboviral disease of ruminants causing significant economic losses. Our risk assessment is based on the epidemiological key parameter, the basic reproduction number. It is defined as the number of secondary cases caused by one primary case in a fully susceptible host population, in which values greater than one indicate the possibility, i.e., the risk, for a major disease outbreak. In the course of the Bluetongue virus serotype 8 (BTV-8) outbreak in Europe in 2006 we developed such a risk assessment for the University of Veterinary Medicine Vienna, Austria. Basic reproduction numbers were calculated using a well-known formula for vector-borne diseases considering the population densities of hosts (cattle and small ruminants) and vectors (biting midges of the Culicoides obsoletus spp.) as well as temperature dependent rates. The latter comprise the biting and mortality rate of midges as well as the reciprocal of the extrinsic incubation period. Most important, but generally unknown, is the spatio-temporal distribution of the vector density. Therefore, we established a continuously operating daily monitoring to quantify the seasonal cycle of the vector population by a statistical model. We used cross-correlation maps and Poisson regression to describe vector densities by environmental temperature and precipitation. Our results comprise time series of observed and simulated Culicoides obsoletus spp. counts as well as basic reproduction numbers for the period 2009–2011. For a spatio-temporal risk assessment we projected our results from the location of Vienna to the entire region of Austria. We compiled both daily maps of vector densities and the basic reproduction numbers, respectively. Basic reproduction numbers above one were generally found between June and August except in the mountainous regions of the Alps. The highest values coincide with the locations of confirmed BTV cases.


Introduction
Veterinary authorities are generally interested in a spatiotemporal risk assessment to establish efficient protection and control measures when facing diseases with high economic impact such as caused by the Bluetongue virus (BTV). Such risk assessments are frequently based on one common epidemiological parameter, the basic reproduction number R 0 . It describes the number of secondary cases caused by a primary case in a completely susceptible population at the beginning of an epidemic. Thus, R 0 may be interpreted as a threshold parameter: a major disease outbreak may only occur for R 0 w1, while for R 0 v1 the disease will fade out [1].
The motivation for this study was the Bluetongue virus serotype 8 (BTV-8) epidemics in Europe, which emerged in 2006 at the border triangle of Belgium, Germany and the Netherlands [2,3]. Within a few years this arboviral disease rapidly spread across North-western and Central Europe [4,5] and caused substantial losses in agriculture amounting to millions Euro [6]. Beside movement restrictions and surveillance, the main veterinary measure was to vaccinate the susceptible populations beginning in 2008 [7]. In Austria, only a few cattle near the border to Germany were confirmed to be BTV-8 seropositive. Due to the German vaccination campaigns achieving a coverage of 78% [8] and preventive vaccination in Austria the eastwards spread of the Bluetongue disease could be stopped at the German-Austrian border ( Figure 1).
In the course of the European BTV-8 outbreak several methods for the calculation of R 0 have been proposed. Some of these methods are based on data from surveillance and monitoring programmes. For example, [9] estimated R 0 between herds for the Netherlands. Furthermore, a theoretical study on R 0 using a deterministic epidemic model was introduced by [10] to investigate the seasonal spread of BTV over several years and to evaluate the effectiveness of vaccination strategies. Most approaches presented so far considered temperature dependent vector parameters as well as host and vector densities. While [11] presented a general concept on the calculation of R 0 for which they discussed the influence of temperature on the transmission cycle, [12] and [13] generated spatial R 0 maps for the Netherlands and Switzerland, respectively. Furthermore, the uncertainty and sensitivity of such temperature-dependent models have been analysed and temperature turned out to be the most dominant factor for determining the magnitude of R 0 [14]. [15] estimated the possible impact of past and future climates on the basic reproduction number R 0 on an European scale. For their study they used input data from climate model runs applied to different emission scenarios of the Intergovernmental Panel on Climate Change (IPCC). The authors have used a similar approach for modelling another mosquito-borne disease outbreak, the Usutu virus epidemics in Vienna [16]. Here we present a predictive risk assessment method based on the basic reproduction number R 0 after [12] applied to a potential Bluetongue outbreak in Austria. Results comprise regions at risk on a spatio-temporal scale representative for a resolution of 10 km and 1 day.
In general, BTV circulates in a natural transmission cycle between vectors (Culicoides biting midges) and hosts (ruminants: mainly cattle, sheep, and goats). For risk assessment using the concept of R 0 the knowledge of both the vector and the host densities are of fundamental importance. While the population densities of farm animals are well documented in national veterinary databases or available worldwide on a regular grid [17], the vector density is generally not well known or rather unknown. Therefore, a central part of our study is concerned with the estimation of the vector population. We established a continuously operating daily Culicoides spp. monitoring at the University of Veterinary Medicine Vienna, Austria, to quantify the seasonal abundance for at least three years. This approach is in contrast to former -mainly weekly -monitoring programs, such as the large-scale entomological surveillance program throughout the European Union to investigate the occurrence and geographical distribution of Culicoides spp. ( [18]; results were published e.g. by [19,20]).

Vector monitoring
For the Culicoides spp. monitoring a black light suction trap of Onderstepoort type, as described by [21], was placed next to an automatic weather station at the university campus (l~16:40 0 E geographical longitude and Q~48:25 0 N geographical latitude). The location in Vienna is characterised by a temperate, fully humid climate with warm summers, corresponding to Cfb climate following the Köppen-Geiger climate classification [22]. This standardized trap is commonly used in monitoring and surveillance programs, because it is very efficient compared to other suction light traps [23]. Following the guideline of [24], the trap was hung at a height of 1.5 m above the ground. The distance to the stables and paddocks with cattle, sheep, and horses was less than 20 m. A collection bottle rotator (model 1512, John W. Hock Company, FL, USA) with eight beakers was used to segregate collections at 24 h intervals. For species evaluation the catches were separated first in Culicoides spp. and other insects (bycatches) under a stereomicroscope. Afterwards the Culicoides species complexes were determined by the characteristic pattern and coloration of the wings according to [24][25][26].
During the monitoring period 2009-2011 a total of 18,952 Culicoides spp. were captured, in which the predominantly trapped complex was C. obsoletus (85.7%). The C. obsoletus complex comprises C. obsoletus, C. scoticus, C. dewulfi, and C. chiopterus. In addition we identified midges of the C. pulicaris complex (9.8%), amongst others (4.5%). In this study we used only the most frequent vectors, i.e. those of the C. obsoletus complex.

Vector modelling
The aim of our vector modelling is to statistically describe the midge catches by operationally available meteorological quantities. Assuming that the relationship between meteorological quantities and C. obsoletus abundance is location-independent (but not independent on the host density), the model will be used to estimate the spatial distribution of C. obsoletus. Additionally, the C. obsoletus abundance may be projected into the future by using meteorological output parameters from numerical weather prediction models. The development of the statistical model was carried out in two steps. In the first step, we performed a cross correlation map (CCM) analysis to investigate which of the meteorological quantities correlates best with the observed C. obsoletus counts. In the second step, we used these meteorological quantities in a Poisson regression model. This procedure was done first for Vienna to demonstrate the performance of the statistical model by comparing it's output with observations and subsequently applied for the entire region of Austria.
The CCM analysis is an extension of the ordinary crosscorrelation by introducing a second time lag. This second time lag is used to average or accumulate meteorological quantities over a period beginning at time lag 1 and ending at time lag 2. Correlation coefficients are then calculated for these averaged quantities, in which the CCMs depict the averaging period resulting in the highest correlation coefficient ( Figure 2). CCMs were adopted for example by [27] to investigate the relationship between Ochlerotatus sollicitans abundance and environmental parameters. Further applications comprise studies on Aedes sollicitans [28] as well as Aedes sollicitans and Culex salinarius [29].
Using a more mathematical notation, the statistical analysis reads as follows: For n i representing a time series of the observed population of midges and X i{j,i{k a meteorological quantity with time index i, the CCM considers cor(log e (n i z1),X i{j,i{k ) for each time lag j and k, at which j §k. Generally, the population data were normalized using a log e (n i z1)-transformation and X i{j,i{k is averaged over the time interval ranging from i{j to i{k. The correlations for all possible time intervals are displayed in a CCM with the time lags j and k range along the two axes. Note, that for the special case of j~k (45 0 -line in a CCM), the correlation coefficients are equal to those of a cross-correlogram. We developed a source code for the calculation of CCMs using the R statistical computing environment [30].
Based on the correlations found by CCMs, a Poisson regression was applied to describe daily vector abundance by meteorological quantities. Again, the analysis was performed using the R statistical computing environment.

Risk assessment
For the different tasks we used vector densities, estimated from both observed and simulated counts of midges, to assess the risk of a potential BTV-8 outbreak in cattle and small ruminants. Therefore, we applied the basic reproduction number R 0 , which may be interpreted as a threshold value for the possibility of a major outbreak (R 0 w1). An analytical solution for R 0 is generally derived from epidemic differential equation models. Here the formula used in recent studies by [12], [15] and [31] is applied, which may be traced back to the historical developments by Ronald Ross and George Macdonald [32].
R 0 equations given by [12], [15] and [31] are formally equivalent but differ in the complexity of parameter formulation. The most comprehensive formulation by [31] provides a gamma distribution for the extrinsic incubation period, which reduces to equation (1) by setting the scale parameters for the gamma distributions to one. Model parameters are probabilities for an infection of cattle p C , small ruminants p S and midges p M , the removal rate (recovery or death rate) for cattle a C and small ruminants a S , as well as some functions describing the temperature dependent behaviour of midges. The latter comprise the biting rate of midges k(T), the virus reproduction rate in midges c M (T), i.e. the reciprocal of the extrinsic incubation period, and the mortality rate of midges m M (T). All constant parameters as well as temperature dependent parameter functions were taken according to [12] and are summarized in Table 1. Note that the temperature dependent parameter functions are not well defined for low temperatures. Therefore, we calculated R 0 only for threshold values of T §13 0 C. Note further that several other functions were proposed in the literature. A sensitivity study performed by [33] using alternative functions for midge mortality, extrinsic incubation period and biting rate results in similar R 0 values (see Text S1 for details). Additionally, the basic reproduction number depends on the host and vector densities. We distinguish between cattle N C , small ruminants N S and midges N M .

Simulated vector population
In a first step we investigated the correlations between time series of midges of the Culicoides obsoletus complex and environmental quantities measured by the automatic weather station beside the light trap. The latter comprise amongst others air  temperature, precipitation, wind speed and relative humidity. CCMs as discussed above were compiled. In accordance with the life cycle of midges we considered time frames of 120 days (4 months) preceding the catches. From all investigated quantities only temperature and precipitation were found to correlate quite well with the daily midge counts. Note that daily precipitation has a skewed frequency distribution and was normalized by logtransformation before CCMs were calculated (Figure 3).
We found a maximum correlation for midges vs. temperature averaged over the preceding 37 days of r(37,0)~0:734 (Figure 3 left). This correlation is slightly higher than the corresponding dayto-day correlation of r(0,0) = 0.697. Between midges and precip-itation we found no day-to-day correlation, but for the precipitation averaged over the period between the preceding days 100 to 16 a correlation of r(100,16) = 0.569 was calculated ( with the daily temperature T, the mean temperature T, and the mean logarithmic precipitation P (contribution of all coefficients significant with pv0:001). As a result, observed vs. predicted daily Culicoides counts for Vienna are depicted in Figure 4. According to that, temperature and precipitation explain about 48.1% of the daily variance of the C. obsoletus counts. We plotted also bi-weekly and monthly time series, although it is clear that the explained variances are higher for accumulated data, because they are usually appropriate for many epidemiological applications. For monthly data, for example, the explained variance increases up to 83.3%.

Dynamic Bluetongue risk assessment for Vienna
We calculated the daily reproduction number R 0 for Vienna by applying Eq. 1, in which vectors and hosts may be alternatively inserted as absolute numbers or densities. Assuming our university campus covers an area of 1 km 2 , values given in absolute numbers or densities are equal. We used an average host density of cattle and small ruminants of N C~NS~2 5 animals/km 2 . The true vector density N M , however, is unknown. In every case it is much higher than the number of midges caught at the location of the light trap. In the absence of any other estimate we applied the reasonable assumption proposed by [12]. Accordingly, the trapped midges were assumed to reflect 1% of the local vector population, i.e. our catches were multiplied by a factor of 100 before entering Eq. 1 (N M~n : 100). Additionally, the temperature dependent functions given in Table 1 are defined for temperatures w13 0 C, hence R 0 is calculated for only 46.3% of the days (mainly between April and October). Figure 5 depicts the time series of the daily reproduction number R 0 for Vienna. While in Figure 5a the observed numbers of midges were used, Figure 5b depicts R 0 values based on the simulated numbers of midges. Taking the first as a gold standard we evaluated the influence of simulated midges N M on the calculation of R 0 by using the two verification measures sensitivity and specificity. We calculated a sensitivity of 0.81 (probability that days with R 0 w1 were correctly realized using R 0 estimates based on simulated N M ) and a specificity of 0.53 (probability that days with R 0 v1 were correctly realized).

Spatio-temporal Bluetongue risk assessment for the entire region of Austria
Assuming that R 0 calculations in Vienna based on simulated N M values satisfy our requirements concerning the accuracy, we applied the R 0 estimation method to the entire region of Austria. Therefore we defined a grid with 10 km spacing, equidistant in geographical longitude (Dl~0:15 0 ) and latitude (DQ~0:09 0 ). Gridded data of temperature and precipitation were taken from the Austrian meso-scale numerical weather prediction model ALADIN [34]. Only data for the period 2010-2011 were used, because data for earlier years were not available. The advantage of ALADIN data over observations (in situ or satellite data) is that they are available on a higher temporal resolution, will be available in near future also on a much higher spatial resolution and gives us the opportunity to project R 0 values to the future. It should be mentioned that we found a few extraordinary high precipitation values within the ALADIN data leading to unrealistic high numbers of midges simulated by Eq. 2. To avoid such artefacts we applied a truncation for the accumulated logarithmic precipitation of Pw1:5 mm/day. This cut-off coincides with annual precipitation maxima of about 1200 mm/year as observed in high Alpine regions [35].
Host and vector densities used for the R 0 calculation over Austria are depicted in Figure 6. The constant densities of cattle N C (Figure 6a) and small ruminants N S (Figure 6b) were taken from the Austrian veterinary database. The spatial distribution of the vector density N M , a function of temperature, precipitation, N C and N S , is given in units of 10 3 midges/km 2 . Again we assumed that the simulated midges reflect 1% of the real vector population. Further we adjusted the vector density for host densities following the findings of [36]. Thus, the vector density at a grid point with geographical coordinates l and Q is calculated from the simulated vectors at this coordinate n(l,Q) multiplied by Note that the vector density varies temporally and spatially, because n(l,Q) is a function of temperature and precipitation. As an example, Figure 6c depicts the vector density for 11 July 2010.
Using these daily vector densities, the host densities as well as the parameters and temperature dependent parameter functions we compiled daily risk maps for the period 2010-2011. Movies of simulated daily vector densities and Bluetongue risk maps are provided as supporting information (see Videos S1 and S2) or on our website http://epidemic-modeling.vetmeduni.ac.at/btvmodel. htm. Figure 7 depicts mean daily R 0 maps for June, July, August, and September 2010. While red areas indicate a potential risk for a major outbreak (R 0 w1), green areas are associated with minor or no risk at all (R 0 v1). For example, the R 0 values for July 2010 are within the range 0:1{4:6, in which the maximum indicates that from one infected animal on average 4.6 animals may be newly infected with BTV. Interestingly areas with a high vector density coincided with high host densities, even without the adjustment after [36]. Verification of the R 0 maps is currently not possible, because so far no major BTV outbreaks occurred in Austria. However, the regions with BTV-8 seropositive cases at the Austrian-German border ( Figure 1

Discussion
We presented a risk assessment for a potential Bluetongue disease outbreak in Austria depicted by time series and maps of the basic reproduction number R 0 . Input data were constant host densities and fluctuating vector densities simulated by a Poisson regression model using temperature and precipitation fields from the Austrian numerical weather prediction model ALADIN. Although the methodology presented here is generally accepted and was recently applied to BTV by [13], [12] and [15], there are some uncertainties in estimating R 0 .
First of all, the applied R 0 formula is originated from an epidemic model, usually formulated by differential equations, by determining the largest eigenvalue of the so-called next-generation matrix. In the case of BTV it has never been verified by comparison with real outbreak data. Such a verification should demonstrate that the observed BTV-8 dynamics may be approximately reproduced by the underlying epidemic model with the parameter functions selected. Reasons for a missing verification may be due to the fact that it is extremely time-consuming, complex and expensive. However, other epidemic models for arboviral diseases were successfully verified against outbreak data, e.g. for the well-documented Usutu virus outbreak in Vienna [37].
Further, it is difficult to estimate vector densities, even if numbers of trapped midges from a monitoring programme are available. We noticed that the official monitoring data sampled at 54 Austrian locations once a week according to the regulation of the European Union [18] provide an excellent overview on the species composition and on the seasonal activity of our native midge populations. Both were hitherto unknown or at least remain undocumented. A quantitative interpretation of the monitoring data, however, is hardly applicable. This is documented by the following small statistical experiment, where we compiled weekly time series based on Monday catches (after [18]) and compared them with weekly time series based on all daily catches. An explained variance of r 2~0 :59 indicates that the numbers of midges sampled according to [18] were afflicted with significant undersampling errors. Further, the official sampling locations were neither representative for their vicinity nor comparable among each other. Also weather stations at short distances to the sampling locations were missing. Therefore, we established our own monitoring with continuously daily observations of midges and environmental parameters from 2009 to present. It allows us to calculate plausible R 0 values, although they should be interpreted with care. Due to the generally unknown magnitude of the vector density (we assumed numbers of midges per square kilometre to be two orders of magnitude higher than sampled at a point location), our R 0 values should be interpreted as relative rather than absolute measures.
Another uncertainty in the vector monitoring lies in the nature of the trapping method. [38] indicated that black light traps are adequate for measuring relative abundance, but not species composition. Further, we don't take into account that midges themselves can be transported by wind over considerable distances, although several models have been developed to investigate this wind-borne spread of BTV [39,40]. An alternative study on C. obsoletus modelling applied a negative binomial regression model (see also Text S1) to simulate the seasonal cycle of daily catches at several locations in England [41]. It distinguishes between the influence of seasonality and meteoro-  logical quantities (temperature, precipitation and wind), but does not consider accumulated quantities as defined by CCMs. Further differences comprise the fraction of trapped midges of the C. obsoletus complex, which is 86% in Vienna and 48% in England. Unfortunately, due to differences in visualization and verification an objective comparison between [41] and the results presented here is not possible. Both approaches underestimate the extremely high variance of the observations, which is a general problem in vector modelling.
Similar to domestic ruminants most species of wild ruminants, e.g. red deer, are susceptible to BTV infection. We neglected the influence of these potential host species due to their still unclear epidemiological role in BTV transmission [42].