Biting midge dynamics and bluetongue transmission: a multiscale model linking catch data with climate and disease outbreaks
We have used a multi-scale modelling approach combining multiple modelling and inference types and techniques: generalised linear mixed-effect models (GLMMs), a mechanistic BTV transmission simulation model, marginal maximum likelihood inference, and direct calculation of resultant BTV risk predictions. We have parameterised our models with data from multiple existing sources (spatio-temporal climate time series, livestock density estimates, midge life process rates, BTV incidence time series), as well as data from our recent field study on the activity of midges at different latitudes in Europe34. The workflow and data sources of this study are summarised in Fig. 1.
Biting midge collection and identification
The method of biting midge sampling and identification was described in our earlier study34. In short, three habitat types were defined in which Culicoides specimens were collected: “farm”, “peri-urban”, and “wetland”. Traps were placed within a 50 m radius of cattle, a house, or waterbody, for farm, peri-urban, and wetland habitats respectively. Habitat types generally matched the classification of the CORINE European Land cover database46. Collections were performed in Sweden (surroundings of Linköping N58.410808, E15.621532), The Netherlands (surroundings of Wageningen N51.964795, E5.662898), and Italy (surroundings of San Benedetto del Tronto N42.949483, E13.878503). Onderstepoort Veterinary Institute black light traps were placed at three independent locations for each selected habitat type. Traps were at least 100 m apart to prevent overlap of the active trapping radius47. More details and exact trap locations can be found in Vogels et al.48 and Möhlmann et al.34. Collections were performed for six consecutive days in each month in all three countries, during the period from July 2014 to June 2015 except the winter months of December, January and February (and March for Sweden). Traps were emptied and repositioned at different locations every 24 h. Collections were sorted and stored in 70% ethanol at − 20 °C. Samples were identified to species level using the Interactive Identification Key for Culicoides (IIKC) developed by Mathieu et al.49.
Predictor variables for midge activity
Tinytag® meteorological data loggers (Gemini Data Loggers, Chichester, UK) were used to record local temperature and relative humidity every 30 min during the collection period from 17th July 2014 until 3rd July 2015. For each habitat in each country, one data logger was used (3 countries × 3 habitats). Additional meteorological data was collected from a number of sources: hourly wind velocities and local temperatures were available from the weather station closest to the trap location in Sweden (Swedish Meteorological and Hydrological Institute (SMHI) Linköping weather station), The Netherlands (Koninklijk Nederlands Meteorologisch Instituut (KNMI) weather stations “De Bilt”, “Deelen”, “Cabauw”, and “Volkel”), and Italy (San Benedetto del Tronto weather station). Finally, additional minimum, maximum, and mean daily temperatures along with precipitation and air pressure were sourced from the E-OBS European climate database50. This data was available at a daily temporal resolution and a spatial resolution of 0.25 degrees latitude and longitude. For climatic variables such as wind only one source per area could be used, whereas for temperature we had the opportunity to explore the most predictive of a number of data sources per country. Habitat effect was included by using the habitat types (farm, peri-urban, wetland) as a categorical predictor in the regression models.
BTV-seroprevalence survey data from Dutch sentinel herds
After the initial BTV outbreaks in 2006, the Dutch government decided to establish a sentinel network of dairy cattle herds in the winter of 2006/2007 to monitor the re-emergence of BTV-8 in 2007, using repeated milk ELISA testing35,36. For study purposes, The Netherlands was divided into 20 compartments based on geographic boundaries as proposed in Commission Decision 2005/393/EC. In each compartment, at least 10 randomly selected herds had to be sampled (with at least sixteen cows per herd) to obtain the required sample size. Herds were not necessarily completely BTV-8 seronegative at initial investigation, but cows designated for the sentinel program had to be BTV-8 seronegative at the moment of selection in May 2007. Therefore, dairy herds to act as sentinels for BTV incidence were selected, that had at least sixteen seronegative cows and at least 50 cattle in total. Monthly milk samples were collected from the sentinel cows in each herd unless prevalence had already reached 100%. The first round of monthly testing of sentinel cows was done in June 2007 and continued until January 2008. The monthly milk samples were tested at GD Animal Health, Deventer The Netherlands for antibodies to BTV-8 using a commercially available ELISA test. For further details on the sampling protocol and commercial ELISA see Santman-Berends et al.35,36.
Regression model for biting midge catches
A preliminary inference investigation using a hurdle negative binomial model to explain trap catches, inferred using Metropolis–Hastings MCMC, suggested that the excess of zero catches could be explained without zero inflation using the hurdle mechanism. We therefore used GLMM regression for its convenience and flexibility. Regression for the fixed effect coefficients and variance parameters of the random effects was performed via maximum likelihood using the Laplace approximation method implemented by the fitglme MATLAB® function. We followed the recommendation guidelines of Bolker et al.51 for using generalized linear models in the context of applied ecology, starting from a model with a full set of predictors and performed systematic model reduction using AICc to score model improvement (backwards model selection). AICc is a well-established information criterion for model selection since it is easy to calculate and interpret for GLMMs. The AICc difference between two models (ΔAICc) estimates the relative likelihood of the two models (with (sim {text{exp(}} – Delta {text{AICc}}/2)) as relative likelihood for the model with greater AICc). Moreover, the model selected by lowest AICc is also the model that would be selected by leave-one-out cross-validation52.
Candidate predictor variables for removal from the model were chosen by assessing which fixed-effect coefficients had the greatest P value (for the null hypothesis that the coefficient is zero) and which random effect coefficients had the smallest predicted standard deviation. This was followed by trialling the removal of either of the variables and removing the variable that reduced AICc the most, until no further reduction could be made. Several trials were made where we started from a number of different highly over-parameterised models, which all ended with the same best model.
For any given combination of predictor variables, the catches were assumed to be conditionally Poisson distributed, with the conditional mean for each collection defined by the log-link relationship,
$$lnleft( {mu_{lct} } right) = beta cdot X_{lt} + b_{l} cdot Z_{lt} + rho_{ct} + varepsilon_{lt}$$
(1)
where (beta) is the fixed effect regression coefficient that applies to all catches, with ({X}_{lt}) denoting the fixed effect predictors at trap location l on day t. ({b}_{l}sim Normal(0,Sigma )) are the location-grouped random effect regression coefficients with covariance matrix (Sigma), with ({Z}_{lt}) denoting the random effect predictors at location l on day t. The number of midges per catch was highly variable. Therefore, we included an independent random effect for each catch location and day ({varepsilon }_{lt}sim Normal(0,{sigma }_{varepsilon })). This corresponds to assuming that the number of midges per catch is Poisson log-normally distributed (a standard distribution for over-dispersed count data53). Spatial autocorrelation has been found in other midge catch regression studies38, so we also included an autocorrelation random effect grouped by country of trapping and day; ({rho }_{ct}sim Normal(0,{sigma }_{varrho })). This effect accounts for events influencing trap catch at a wider scale than just the location definition of each trap. See supplementary information for full model variables, regression coefficients and random effect variance estimates as well as AICc improvements from other models.
Connecting the regression model for midge captures to a midge biting prediction for herds
We connected the regression model for daily midge trap catch to a prediction of daily biting on cattle. In the BTV modelling literature it is common to assume that the vector-to-host ratio, and therefore the biting rate per livestock host, is proportional to the expected capture rate of a midge trap catch regression model28,54. In this study, we assumed that the biting rate per livestock host is proportional to the regression model presented in Results, with the major difference that we infer the scaling parameter that best fitted the serological data from Dutch sentinel herds. The location-group random effects in the regression model modelled unexplained spatial variation in average midge capture counts between trapping locations. We assumed that this unexplained spatial variation in midge abundance (as measured by trapping) mirrored the spatial variation in midge biting between different herds. Combining the proportionality assumption with the spatial variation assumption gave the biting rate among herds as,
$${text{Biting }},{text{rate }},{text{of }},{text{susceptible}},{text{ midges }},{text{per}},{text{ cattle }},{text{at }},{text{herd}},h,{text{on}},{text{ day}},t = xi mu_{ht} .$$
(2)
where ({mu }_{ht}) is the expected trap catch from the regression model given the climate condition local to herd h on day t and (xi) is a scaling parameter between the mean catch prediction and the biting rate prediction. In order to infer the proportionality factor between the catch model and the daily biting rate ((xi)), the outcomes of a dynamic transmission model for each herd were linked to the observed BTV seroprevalence data (see—“Mechanistic transmission model for BTV transmission within herds” section).
Mechanistic transmission model for BTV transmission within herds
We used a dynamic and mechanistic model of BTV transmission within herds, which was then matched to the data from the Dutch sentinel study. The dynamic BTV transmission model was formulated using disease compartments and rate-based transitions (see Keeling and Rohani55 for further details on this class of transmission model). In addition, it is in most respects similar to the model presented by Gubbins et al.27 in treating infectiousness amongst cattle, and latency amongst midges, as multi-stage processes that evolve deterministically. In particular, we follow Gubbins et al. in modelling the life processes of the infectious and latent midges (mortality, extrinsic incubation of BTV, and biting) as temperature dependent, and therefore varying daily with local background atmospheric temperature. (see Fig. 6 for a schematic diagram of the transmission model).
Figure 6
Schematic representation of the cattle herd level BTV transmission model. The population of cattle and infected biting midges are divided amongst discrete disease compartments. BTV latent midges (EM) enter the model at a rate proportional to the daily prediction of the catch model. The extrinsic incubation period for latent midges is modelled as a multi-stage process before midges become infectious (IM). Susceptible cattle (SC) become infectious cattle (IC) after a bite from infectious midges (IM), to become resistant cattle (RC) in time (also modelled as a multi-stage process; red box). Transitions are shown as solid lines, coloured according to their dependence on environmental variables: constant per-capita (black), daily mean temperature dependent (red), all predictor variables of capture model and the catch-to-bite scale parameter ({varvec{xi}}) (blue). Dotted lines indicate where the number of infected individuals in one species increases the incidence rate in the other species. Outcomes of the model are linked to observed cattle milk serology time series by the first two infectious stages for cattle with virus being undetectable by ELISA (blue box), whereas subsequent infectious stages and the recovered stage are detectable by ELISA (green box). The likelihood function for (xi) was inferred by marginalisation over the latent stochastic variables affecting model outcomes (e.g. herd-specific random effects, daily fluctuations in midge activity). Used images were available under open licence Creative Commons Deed CC0.
Full size image
The dynamic and mechanistic BTV transmission model used in this study describes the evolution of the numbers of susceptible, infected, and recovered cattle as well as latent and infectious midges for each herd (Fig. 6). Temperature-dependent midge bionomic rates were used for biting frequency of individual infectious midges56, the incubation rate of BTV within the midge57, and for the midge mortality58. The midge bionomic rates at each herd on each day were determined by the local mean temperature day according to the E-OBS climate dataset (see above—Predictor variables for midge activity). We modelled the incubation period of BTV within the midge vector as a ten-stage process, which is within the range of best fit models found in meta-analysis and laboratory studies of BTV incubation57 (see supplementary information for complete model details and literature estimates for rates).
The period during which BTV-infected cattle are detectable using an ELISA test (typically from 8–9 days post-infection onwards59) does not match the period during which the cattle are infectious (rapidly post-infection and then for an average of 20.6 days27). BTV-infected cattle can be in four states that are relevant to transmission modelling and their milk serology: 1) uninfected and susceptible to BTV, 2) infectious but undetectable by milk ELISA, 3) infectious and detectable by milk ELISA or, 4) non-infectious recovered from BTV but still milk ELISA positive. The BTV infectious period for cattle is usually modelled as a 5-stage process27, therefore it was convenient to model cattle in the first two stages of their infectious period (an average duration of 8.2 days) as infectious but undetectable. Cattle in the final three stages of the BTV infectious period are infectious and detectable (see Fig. 6 for a schematic representation of the BTV transmission and serology model).
The rate at which new susceptible midges arrived to bite cattle was assumed to be proportional to the expected number of midges from the regression model (({mu }_{ht})), with predictions of midge catches on each day and in each regional compartment of The Netherlands whilst the sentinel herd study was ongoing using the E-OBS historic climate records (see above—“Connecting the regression model for midge captures to a midge biting prediction for herds” section). The random effects in the catch model imply that ({mu }_{ht}) is a daily varying random variable, and that our transmission model is in the class of piecewise-deterministic Markov processes60. We assumed that the location-grouped random effects observed in the catch model became herd-grouped random effects for the biting model. In other words, we assumed that the high variance in midge catching between trapping location reflects high variance in midge biting between different cattle herd locations. Although an assumption, this would explain the highly variable intensity of BTV transmission observed between different herds in The Netherlands sentinel survey despite each herd experiencing a similar climate36. Because the herd locations were known only as geographic compartment occupancy, the daily local climate variables from the gridded E-OBS data used for predicting ({mu }_{ht}) were averaged over all spatial grids overlapping the herd’s geographic compartment. Accessing detailed and spatio-temporally resolved wind data across Europe was challenging, therefore we used the long-term average wind velocity of The Netherlands weather stations (see above—“Predictor variables for biting midge activity” section) as a constant predictor.
When simulating an outbreak of BTV within a herd, we first determined all relevant climatic predictors for the herd’s regional compartment and the daily temperature dependent midge bionomic rates. Second, we generated the herd-grouped and daily varying random effect coefficients which determined how biting at the herd from susceptible midges varied from a median prediction. Third, we solved the resultant deterministic BTV transmission model for each farm using the ode45 MATLAB® function (see Fig. 6 for an overview).
Inferring the catch-to-biting scale parameter from serological data
We inferred a maximum likelihood estimator for the trap-to-bite scale parameter by repeated simulation of the percentage of cattle detectable by milk ELISA tests herds in the sentinel herd network. For this we used the climatic conditions of The Netherlands in 2007, and at each simulation repeated redrawing the unobserved random effects for each herd and day. The average likelihood over many repeated simulations corresponds to a Monte Carlo estimate of the true marginal likelihood of the parameter (xi). Estimating the likelihood over a range of values of (xi) allowed the construction of a log-likelihood profile.
The stochastic elements of the piecewise-deterministic BTV transmission model were (i) the herd-grouped random coefficients (this modelled how biting varied from herd-to-herd) and (ii) the daily varying random effects (this modelled how biting varied from day-to-day). It is convenient to denote ({W}_{h}= ({{b}_{l}}^{(h)},{{rho }_{0}}^{(h)},{{rho }_{1}}^{(h)},{{rho }_{2}}^{(h)},…,{{epsilon }_{0}}^{(h)},{{epsilon }_{1}}^{(h)},{{epsilon }_{2}}^{(h)},…)) as the collection of all stochastic elements for herd h. For each simulation of the transmission model in each herd, we first drew ({W}_{h}) from their inferred distribution (see supplementary information for distribution parameters of best fitting model).
The likelihood of ({W}_{h}) and (xi) for each herd h was the chance of selecting the numbers of ELISA seroconverted cattle observed at the herd each month by The Netherlands sentinel study from the underlying distribution of ELISA detectable cattle implied by simulating the transmission model conditional on ((xi ,{W}_{h})),
$$L_{h} left( {xi ,W_{h} } right) = P({text{Serology }},{text{data }},{text{collected }},{text{at }},{text{herd }},h| xi ,W_{h} ).$$
(3)
Since we were not interested in inferring the specific values of ({W}_{h}) for each herd, they were treated as “nuisance” parameters. We inferred a maximum likelihood estimate, with confidence intervals, for (xi) by first estimating the marginal likelihood for (xi) (that is the likelihood after integrating over all possible values of the nuisance parameters) at each herd h,
$${L}_{h}(xi ) = int {L}_{h}left(xi ,wright)fleft(wright)dw.$$
(4)
where (f) is the density function for the distribution of random effects derived from the trapping model. The marginal log-likelihood function for the trap-to-bite scaling parameter, (l(xi )), for serological data over a number of herds, was then just the sum of the individual herd marginal log-likelihoods,
$$l(xi ) = sum_{h}log {L}_{h}left(xi right).$$
(5)
The herds we chose to contribute to the log-likelihood were those where BTV was found to be already present at the beginning of the study (see supplementary information for more details), to avoid making further assumptions about the introduction mechanism into the herd.
In practice, the log-likelihood was estimated for a profile of values of (xi) by simulating multiple realisations of ({W}_{h}) for each herd, that is we estimated (4) by Monte Carlo integration for (3) over a range of values of (xi), and interpolating between points with polynomial regression. The maximum likelihood estimator, ({xi }^{*}), was the maximizer of the marginal log-likelihood function presented in the main text along with confidence intervals derived by a standard comparison to the ({chi }^{2}) distribution (see methods section of King et al.61 for a brief but comprehensive introduction to maximum likelihood estimation using log-likelihood profiles in the context of inference for dynamical systems).
Calculating and mapping the herd reproductive ratio for bluetongue
The reproductive ratio for BTV will differ from day to day and across space. This reflects seasonality and variation in both climatic trends, and the population density of midges and livestock hosts. We approached estimating the reproductive ratio for BTV in the spirit of the case reproductive ratio62 using a technique already developed for midges spreading BTV37. That is, we calculated the expected number of secondary cases amongst hosts due to a host initially infected on each day t in each grid cell x whilst taking into account how the conditions for BTV transmission at location x changed after time t, and using the maximum likelihood model for midge biting. The size of each grid cell was determined by the resolution of the relevant climate data. We used the 0.25 degrees longitude and latitude grid resolution of the E-OBS climate datasets to map reproductive ratio estimates for Europe across space and time. Cattle and sheep densities across Europe were drawn from the livestock Geo-Wiki dataset63. The livestock Geo-Wiki datasets were more finely resolved (0.0083 degrees grid resolution) than the E-OBS climate datasets. We calculated the reproductive ratio for each grid cell x on each day t at the resolution of the E-OBS datasets. The cattle and sheep densities for this coarser grained grid were the average densities over the Geo-Wiki cells contained within the coarser grained grid.
The average number of secondary BTV cases amongst all hosts (cattle and sheep) given a host initially infected on day t and at grid cell x, is denoted ({R}^{(C)}(x,t)) for an initial infected cow and ({R}^{(S)}(x,t)) for an initial infected sheep. For both host species, the average number of secondary cases can be calculated by considering; how many days the host’s viraemia will last, the rate at which the host is bitten each day, the percentage of the biting midges that will become infected, how many of these biting midges are expected to survive their EIP to become actively infectious, and how many livestock will be successfully infected by those actively infectious midges. The methodology for combining these estimates using information about midge bionomic rates, EIP distribution, and the daily temperatures on each day after the initial host was infected has been developed by Brand and Keeling37.
In this study, we adapted the Brand-method for calculating the reproductive ratio to two species, and used the catch-to-biting scalar derived from comparison between the mechanistic transmission model and the herd sentinel serological survey. The cross-transmission between host species depends on how midge bites are distributed between cattle and sheep. We estimated the proportion of midge bites on cattle at grid cell x, ({phi }^{(C)}(x)), given the availability of sheep using a common relative preference model, e.g. Szmaragd et al.64,
$${phi }^{left(Cright)}left(xright)=frac{{N}^{left(Cright)}left(xright)}{{N}^{left(Cright)}left(xright)+ pi {N}^{left(Sright)}left(xright)}.$$
(6)
where ({N}^{(C)}(x)) and ({N}^{(S)}(x)) are, the local density of cattle and sheep at grid cell x. Parameter (pi) is a measure of the vector preference for sheep compared to cattle; (pi 1) preference for sheep. A relative biting study for sheep and cattle has revealed a preference for biting cattle30, from which we derived an estimate (pi =0.115) for use in this study. We combined ({R}^{(C)}(x,t)) and ({R}^{(S)}(x,t)) into a single reproductive ratio by calculating the leading eigenvalue of the next-generation matrix65,
$$Rleft(x,tright)= sqrt{{phi }^{left(Cright)}left(xright){R}^{left(Cright)}left(x,tright)+left(1-{phi }^{left(Cright)}left(xright)right){R}^{left(Sright)}left(x,tright)}.$$
(7)
An attractive feature of using the reproductive ratio as a measure of transmission intensity is its uncomplicated relationship with the persistence of transmission; if (Rle 1) then the infectious pathogen cannot persist. However, we expect that the biting rate, and therefore the reproductive ratio, will vary from herd-to-herd. From Eqs. (1) and (5) we see that the rate of biting from the susceptible midge population at each herd on each day depends on the random coefficients, ({{b}_{l}}^{(h)}), and daily varying random effects,({{rho }_{ct}}^{(h)}) and ({{varepsilon }_{t}}^{(h)}),
$$text{Biting rate of susceptible midges per cattle at herd }htext{ on day }tpropto expleft({{b}_{l}}^{left(hright)}cdot {Z}_{lt} +{{rho }_{ct}}^{left(hright)} + {{varepsilon }_{t}}^{left(hright)} right).$$
(8)
The daily varying random effects ((rho) and (epsilon)) are averaged over our estimates for ({R}^{(C)}) and ({R}^{(S)}) (this can be achieved analytically; see supplementary information for further details), and therefore our estimate of the reproductive ratio does not depend on daily fluctuations in midge activity. However, variation in the herd-grouped random coefficients indicated systematic differences in midge activity between herds that will not ‘average out’ over time. The distribution of ({{b}_{l}}^{(h)}) therefore implied a distribution of biting rates for herds within each grid cell on each day, and therefore a distribution of values of (R) for herds in each grid cell and on each day.
We present the distribution of (R) for herds by considering the reproductive ratio that would be calculated if the random variable in Eq. (8), ({{b}_{l}}^{(h)}cdot {Z}_{lt}), took its pth percentile value every day, denoting this reproductive ratio, ({R}_{p}(x,t)). ({R}_{p}(x,t)) estimates the reproductive ratio that p% of herd reproductive ratios are less than in grid cell x on day t. Also, we numerically invert the threshold relationship to find the percentage value, ({varvec{P}}), such that ({R}_{p}(x,t)) satisfies the threshold quantity,
$${varvec{P}}left(x,tright)={1-p in [mathrm{0,1}] | {R}_{p}(x,t) = 1 }.$$
(9)
({varvec{P}}(x,t)) is therefore an estimate of the percentage of herds that could have a multiplying BTV outbreak in grid cell x if BTV was introduced on day t.
To enable spatially extending risk predictions for BTV across Europe certain assumptions about how the midge biting model could be interpolated between the trap locations were necessary. The best regression model for midge catches found significantly different seasonality at the Italian catch locations compared to Sweden and The Netherlands. We assumed that day length was a determinant of midge seasonality and overwintering at different latitudes; for example, the first day of the year with a day length shorter than 9 h has been associated with the onset of overwintering in UK midges66. We found no midges on days with less than 8.5 h of daylight when trapping, although catching was not attempted outside of March-November so the number of samples was small. Therefore, at latitudes where no day is ever shorter than 8.5 h (lower latitudes than 46oN, which is more or less the border of Switzerland and Italy) we used Italian seasonality to predict midge biting. At latitudes that have at least one day shorter than 8 h (higher latitudes than 49oN) we used the Swedish and Dutch midge seasonality. Between these two latitudes, a linear interpolation between the predictions of the two seasonal models was applied.
All map images were generated using MATLAB® contourf function. Both the livestock density and E-OBS datasets included grid cells that contain only water, these cells were coloured white.
The MATLAB® code for generating the reproductive ratio estimates is available at a github repository (https://github.com/SamuelBrand1/BTV-European-Reproductive-Ratios). More
