Assessing the risk of bluetongue to UK livestock: uncertainty and sensitivity analyses of a temperature-dependent model for the basic reproduction number

Since 1998 bluetongue virus (BTV), which causes bluetongue, a non-contagious, insect-borne infectious disease of ruminants, has expanded northwards in Europe in an unprecedented series of incursions, suggesting that there is a risk to the large and valuable British livestock industry. The basic reproduction number, R0, provides a powerful tool with which to assess the level of risk posed by a disease. In this paper, we compute R0 for BTV in a population comprising two host species, cattle and sheep. Estimates for each parameter which influences R0 were obtained from the published literature, using those applicable to the UK situation wherever possible. Moreover, explicit temperature dependence was included for those parameters for which it had been quantified. Uncertainty and sensitivity analyses based on Latin hypercube sampling and partial rank correlation coefficients identified temperature, the probability of transmission from host to vector and the vector to host ratio as being most important in determining the magnitude of R0. The importance of temperature reflects the fact that it influences many processes involved in the transmission of BTV and, in particular, the biting rate, the extrinsic incubation period and the vector mortality rate.


INTRODUCTION
Bluetongue virus (BTV) is a double-stranded RNA virus (Reoviridae: Orbivirus) that causes bluetongue (BT), a non-contagious, insect-borne infectious disease of ruminants. All ruminant species are able to be infected, but improved breeds of sheep and some species of deer tend to develop particularly severe clinical signs of disease. Owing to this and its ability to spread rapidly and without apparent warning, BT has been cited as one of the most important of the Diseases Notifiable to the OIE. Between 1998 and 2005, the virus expanded northwards in Europe in an unprecedented series of incursions, causing the most severe outbreak of this disease ever recorded and resulting in the deaths of over 1.5 million sheep Purse et al. 2005). In August and September 2006, the first outbreaks of BT were reported in northern Europe in Belgium, France, Germany and The Netherlands (Albers et al. 2007).
This dramatic expansion has been attributed to the northward movement of the major Old World vector of the disease, Culicoides imicola, probably under the influence of climate change, and to the involvement of more northerly Culicoides species whose distribution overlaps that of C. imicola and extends across the whole of central and much of northern Europe, including UK (Mellor & Leake 2000;Purse et al. 2005;Carpenter et al. 2006). The extension of BTV further north in Europe than ever before, and the involvement of vectors which also occur throughout UK, suggests that there is a risk to the large and valuable UK livestock industry ( Wittmann & Baylis 2000;Purse et al. 2005). Indeed, it has been shown already that Culicoides present in some areas of UK are capable of replicating BTV to transmissible levels (Carpenter et al. 2006). To date, however, the level of risk according to other epidemiological variables and how this risk varies temporally and geographically across UK has yet to be assessed.
The basic reproduction number, R 0 , provides a powerful tool when assessing the risk of disease invasion (Anderson & May 1986Heffernan et al. 2005). This quantity is defined as the average number of secondary cases arising from the introduction of a single infected individual to an otherwise susceptible population. A disease is able to invade a host population (i.e. increase in frequency) only if R 0 O1 and, consequently, R 0 provides a means to assess the level of risk posed by a disease. Vector-borne diseases present a particular problem when computing R 0 , owing to the need to consider the biology of the host, the disease agent and the vector. Previous studies have derived R 0 for a number of vector-borne diseases including malaria (MacDonald 1957), trypanosomiasis (Rogers 1988), African horse sickness (Lord et al. 1996), dengue (Newton & Reiter 1992;Luz et al. 2003) and West Nile virus ( Wonham et al. 2004;Cruz-Pacheco et al. 2005), but, to date, the basic reproduction number for BTV has not been considered.
In this paper, we developed a model for the transmission of BTV in a mixed population of cattle and sheep and used this to compute R 0 . Estimates for each parameter which influences R 0 were obtained from the published literature. Wherever possible estimates were derived to reflect the UK situation and, in particular, UK vector species. Moreover, explicit temperature dependence was included for those parameters for which it had been quantified. Uncertainty and sensitivity analyses were carried out to identify the most influential parameters in R 0 and, hence, those for which the most robust estimates are required.

TRANSMISSION MODEL FOR BTV
To capture the transmission dynamics of BTV, it is essential to consider the dynamics of infection in the ruminant host species (mainly cattle and sheep) and the Culicoides vector. It is important to include both cattle and sheep as hosts: cattle harbour the infection usually without displaying clinical disease and, hence, provide a covert reservoir of BTV for the vector, while sheep are often afflicted by severe clinical signs of disease.
The cattle and sheep populations are assumed to be constant (H i ), except for disease-associated mortality, and are subdivided into the proportions susceptible (i.e. uninfected), infected and recovered, denoted by x (i ) , y (i ) and z (i ) , respectively, where the superscript i indicates cattle (C ) or sheep (S ). To allow for a more general distribution for the duration of viraemia, the infected host population, y (i ) , is subdivided into a number of stages, with newly infected hosts entering the first stage and then passing through each successive stage. If the mean time spent in each stage is 1/nr, the total length of time spent in the n classes follows a gamma distribution, with mean 1/r and variance 1/nr 2 (Lloyd 2001). For an appropriate choice of n, this approach includes the common epidemiological assumption of an exponentially distributed infectious period (nZ1) and a fixed infectious period (n/N). The vector population (N ) is subdivided into the number of adult female midges that are susceptible (i.e. uninfected), latent (i.e. infected, but not infectious) and infectious, denoted by S, L and I, respectively. To allow for a more general distribution for the extrinsic incubation (i.e. latent) period (EIP), the latent class is subdivided into a number of stages in a similar approach to that described above for the duration of host viraemia. Once infectious, midges remain so for life. Adult males and immature (larval and pupal) stages are not considered as they do not blood feed and, hence, do not transmit BTV.
The dynamics of BTV infection in the cattle and sheep populations are described by the following set of linked differential equations: where i denotes host species and j the infection stage.
The duration of viraemia in each host species follows a gamma distribution with mean 1/r i and variance 1/n i r i 2 (i.e. there are n i stages each of mean duration 1/n i r i ) and d i is the rate at which infected animals die of severe clinical disease. Host reproduction and natural host mortality are ignored because the natural lifespan is much longer than the duration of an outbreak or the lifespan of the vector. The force of infection for each host species, l i , is given by which is the product of the probability of transmission from an infected midge to a host (b), the biting rate on the species (a i ), the ratio of vectors to hosts (N/H i ) and the proportion of bites which are from infectious midges (I/N ). The biting rate on species i can be decomposed such that a i Zaf i , where a is the reciprocal of the time interval between blood meals and f i is the proportion of bites on the species. The proportion of bites on each host species (f i ) is given by where a i is a measure of vector preference for host species i (Sota & Mogi 1989;Kelly & Thompson 2000).
In the case of two host species (i.e. cattle and sheep), equation (2.3) can be simplified to 2)) and the proportion of bites on cattle is fZm S /(m S Csm C ) (cf. equation (2.4)).
The dynamics of BTV infection in the vector population are described by ð2:5Þ where the subscript j denotes the stage of the EIP; r is the rate of recruitment to the adult female midge population; the EIP follows a gamma distribution with mean 1/n and variance 1/kn 2 (i.e. there are k stages each of mean duration 1/kn); and m is the vector mortality rate. The force of infection for vectors, l V , is given by which is the product of the probability of transmission from an infected host to a midge (b) and the biting rate on infectious hosts of each species. Point estimates or plausible ranges for each parameter were derived from the literature (table 1). Point estimates for parameters relating to the duration of host viraemia were obtained by fitting gamma distributions to published data. Ranges for the reciprocal of the time interval between blood meals (a), the vector mortality rate (m) and the EIP (1/n) reflect the dependence of these parameters on temperature and, furthermore, temperature-dependent functions have been derived for these parameters Mullens et al. 2004;table 1). Ranges for the remaining parameters reflect uncertainty in their values.

BASIC REPRODUCTION NUMBER FOR BLUETONGUE
In the mathematical theory of epidemics, the basic reproduction number, R 0 , is defined as the dominant eigenvalue, r (K ), of the next-generation operator, K (Diekmann et al. 1990;Heffernan et al. 2005). For the transmission model described by equations (2.1)-(2.6), the basic reproduction number is given by where F is a matrix reflecting the rate at which new infections arise (i.e. transmission from vector to host or host to vector) and V is a matrix reflecting the rate at which individuals enter or leave infection classes (either due to completion of an infection stage or death). The methods used to compute R 0 are presented in the electronic supplementary material, appendix A, though it is not possible to derive an expression in terms of the parameters for the general model. For the special case in which the duration of viraemia in cattle and sheep and the EIP follow exponential distributions (i.e. n C Zn S ZkZ1), the basic reproduction number can be derived explicitly (see electronic supplementary material, appendix A). In this case, it is given by where the parameters are summarized in table 1 (cf. Lord et al. 1996). Heuristically, expression (3.1) can be understood as follows. For each host species, an infected individual remains infectious for 1/(r i Cd i ) days, during which time it is bitten by susceptible midges on average am i f i times per day, a proportion, b, of which results in an infected midge. Thus, on average a single infected host will result in bam i f i /(r i Cd i ) infected midges. The probability of a vector surviving the EIP to become infectious is n/(mCn). It survives for 1/m days, during which time it will bite susceptible hosts at a rate af i , a proportion of which, b, will result in an infected host. Thus, on average a single infected midge will result in baf i (n/mCn)/m infected hosts. The basic reproduction number (3.1) is the geometric mean of the number of infected midges per infected host and the number of infected hosts per infected midge.

UNCERTAINTY AND SENSITIVITY ANALYSES
For the uncertainty and sensitivity analyses, replicated Latin hypercube sampling (LHS) was used to explore the parameters influencing the basic reproduction number, R 0 (Blower & Dowlatabadi 1994;Sanchez & Blower 1997;Luz et al. 2003). The LHS results were used to compute the median and maximum value for R 0 , and the proportion of samples yielding R 0 O1. The distributions for each parameter which resulted in values of R 0 !1 or R 0 O1 were plotted to identify parameters associated with R 0 O1. The sensitivity of R 0 to changes in each parameter was assessed by calculating the partial rank correlation coefficients (PRCCs). Full details of these methods are presented in the electronic supplementary material, appendix A. The analyses were implemented for three types of farm depending on species composition: cattle only, sheep only and both cattle and sheep. This reflects the different types of farm found in UK and allows the impact of the different host species on the basic reproductive number to be assessed in greater detail.

RESULTS
For all farm types, values of R 0 O1 were associated with higher values of the probability of transmission from host to vector (b), the vector to host ratios (m C or m S ) and temperature (T; see electronic supplementary Taylor (1986) and Melville et al. (1996) duration of viraemia (sheep) mean 1/r S 16.4 duration of viraemia based on experimental infection and virus isolation in embryonated chicken eggs; parameters estimated by fitting a gamma distribution to data presented in papers referred to in right-hand column Goldsmit et al. (1975) and Veronesi et al. (2005) no. of stages Gerry & Mullens (2000), Wittmann et al. (2002) and Mullens et al. (2004) no. stages Birley & Boorman (1982), Braverman et al. (1985),  and Wittmann et al. (2002) material, appendix B for the full results). However, the effect of the vector to host ratios was weaker for farms with both cattle and sheep. For the remaining parameters, there were few differences between the distributions of the parameters yielding R 0 !1 or R 0 O1 for any of the farm types. Cattle-only farms yielded the highest values for R 0 (assessed in terms of the median and maximum values, and the proportion of samples yielding R 0 O1), followed by sheep-only farms, then mixed cattle and sheep farms (table 2). This effect of farm type is shown by the differences in empirical cumulative distribution functions (CDFs) for R 0 (figure 1a,c,e) and the distribution of R 0 at different sampled temperatures (figure 1b,d, f ). Although the magnitude of R 0 differed according to farm type, its dependence on temperature did not (figure 1b,d, f ). Below 108C, R 0 was zero because the extrinsic incubation rate is zero (table 1). Above 108C, R 0 increased to a peak value in the range 20-258C, after which it decreased.
PRCCs were calculated between each input parameter and the corresponding value of the basic reproduction number (figure 2). For all farm types, the strongest correlation was between temperature, which influences biting rate, EIP and vector mortality rate, and R 0 , with a PRCCO0.7. Similarly, the probability of transmission from host to vector (b) was significantly, positively correlated with R 0 (PRCCZ0.3) for all farm types. For farms which had only a single host species, there was a significant positive correlation between the vector to host ratio and R 0 (PRCCZ0.3). However, the strength of correlation diminished when two host species were present on a farm, only that for cattle remained significantly different from zero ( p!0.05). For the remaining parameters, any correlations were weak and typically did not differ significantly from zero ( pO0.05).

DISCUSSION
The basic reproduction number, R 0 , provides a quantitative framework in which to address the question of level of risk posed by a disease to a population. It allows the identification, first, of what needs to be known in order to assess risk and, second, which of these factors are most important in determining the magnitude of risk. In the case of the transmission model for BTV developed in this paper, (2.1)-(2.6), 12 parameters were needed to compute R 0 . For three of these parameters (biting rate, EIP and vector mortality rate), it was possible to include explicit temperature dependence (table 1). Uncertainty and sensitivity analyses identified temperature as the most important factor in determining the magnitude of R 0 , though the probability of transmission from host to vector and the ratio of vectors to both cattle and sheep were also important (figures 1 and 2). The importance of temperature in determining R 0 reflects the number of temperature-dependent processes involved in the transmission of BTV: the biting rate (Mullens & Holbrook 1991;Mullens et al. 2004), the EIP Wittmann et al. 2002;Mullens et al. 2004) and the vector mortality rate Wittmann et al. 2002). Moreover, these three temperaturedependent parameters were identified as being important if they were included directly in the uncertainty and sensitivity analyses instead of temperature (results not shown).
As with all uncertainty and sensitivity analyses, these conclusions are valid only over the parameter ranges considered. Wherever possible, these ranges (table 1) were derived to reflect the situation in UK, but in many cases they were drawn from data for a range of countries and vector species. Parameters relating to infection in the host are likely to be robust, but those relating to the vector are likely to vary with species and climate. For some parameters (time interval between blood meals, probability of transmission from host to vector and vector mortality rate), estimates were obtained for UK Culicoides species (Birley & Boorman 1982;Carpenter et al. 2006), but in most cases they were derived from outbreak areas, notably the USA (principal vector species, Culicoides sonorensis) and the Mediterranean basin (C. imicola; see references in table 1 for details). Although most effort should focus on estimating the five key parameters (biting rate, EIP, vector mortality rate, probability of transmission from host to vector and the ratio of vectors to cattle and sheep) and, in particular, their dependence on temperature and other environmental factors, it should be noted that all parameters require more robust estimates and, in particular, estimates that are directly applicable to UK.
The biting rate has often been cited as a critical parameter for the transmission of vector-borne diseases (McDonald 1957;Anderson & May 1991;Dye 1992;Lord et al. 1996), largely because it influences transmission from both host to vector and vector to host. The biting rate for Culicoides midges is related to the time interval between blood meals (table 1), and increases with temperature (Mullens & Holbrook 1991;Mullens et al. 2004). Furthermore, Gerry et al. (2001) showed that the host biting rate (defined as the product of the biting rate and the vector to host ratio) varied seasonally, most probably reflecting the influence of temperature . In addition, a threshold value for the host biting rate was identified below which no seroconversions were observed in sentinel calves , which reflects the dependence of R 0 on this parameter (see, for example, equation (3.1)).
The EIP for BTV exhibits strong temperature dependence Wittmann et al. 2002;Mullens et al. 2004), varying from 26 days at 158C to 4 days at 308C . However, these values for the EIP are based on laboratory experiments in which Culicoides were maintained at a constant temperature; yet temperatures will clearly fluctuate in the field. As BTV is not able to replicate to transmissible levels within the vector at low temperatures , it will be important to extend the model so that the temperature dependence of the EIP is incorporated in a more realistic way. Relatively small changes in the survivorship of a vector can, in principle, greatly affect the transmission of a disease agent, such as BTV, which requires a period after infection of a vector before transmission is possible. Seasonal variation in survivorship has been demonstrated for several Culicoides species (Birley & Boorman 1982;Braverman et al. 1985; and, in particular, the vector mortality rate has been shown to increase with temperature Wittmann et al. 2002). The vector mortality rate will determine the probability that a vector will survive the EIP. Similarly, the biting rate is related to the time interval between blood meals, and the probability of a vector surviving this period will depend on survivorship. Consequently, the importance of the vector mortality rate arises from its impact on these processes. Moreover, the temperature dependence of all three of these parameters means that temperature has a strong influence on the magnitude of R 0 (figures 1 and 2).
Estimates for the probability of transmission from host to vector varied from 0.1 to 15% (table 1). The lower estimates were for field-caught C. sonorensis midges in the USA (Nunamaker et al. 1997;Gerry et al. 2001), while the higher estimates were for Culicoides obsoletus group midges caught in UK and subsequently fed on infected blood in the laboratory (Carpenter et al. 2006). The difference in the probability of transmission may reflect species or population differences in oral susceptibility to BTV, or may result from the methods used to determine competence. Moreover, viral titres in the host vary during the course of an infection (Goldsmit et al. 1975;Veronesi et al. 2005), which may result in changes in the probability of transmission to the vector (Bonneau et al. 2002). In the model, however, host viral titre is assumed to be constant throughout viraemia, which may result in the overestimation of the basic reproduction number.
For farms with only a single host species (either cattle or sheep), the vector to host ratio was an important determinant of the magnitude of the basic reproduction number, with higher ratios associated with higher values of R 0 (figure 2a,b; cf. equation (3.1)). For farms with both host species, the positive association remained, but was weaker (figure 2c). Moreover, R 0 tended to be highest on cattle-only farms, lower on sheep-only farms and lowest on mixed farms (table 2 and figure 1). This represents a dilution effect (Lord et al. 1996) because sheep are assumed to be a poorer host compared with cattle in terms of the duration of viraemia. There are, however, aspects to the impact of the number of cattle and sheep on the magnitude of R 0 , which are not apparent from the uncertainty and sensitivity analyses. In particular, adding cattle to a sheep farm can result in an increase in R 0 over a limited range for the number of cattle introduced, beyond which R 0 decreases as the vector to host ratio decreases. This effect depends on vector preference and the relative duration of viraemia in the two host species (see, for example, equation (3.1); see also Lord et al. (1996)).
The transmission model, (2.1)-(2.6), considers only a single vector. However, there are likely to be multiple species involved, which may differ in their ability to transmit BTV. Both the C. obsoletus and C. pulicaris groups of midges are found in UK (Carpenter et al. 2006) and have been implicated in the transmission of BTV ( Torina et al. 2004;Purse et al. 2006). Variation in oral susceptibility to BTV infection has been shown in C. obsoletus group midges caught in different regions of UK (Carpenter et al. 2006), which could represent different proportions of those species within the group possessing different levels of competence. Furthermore, Torina et al. (2004) presented tentative evidence based on light-trap catches that there may be differences between the C. obsoletus and C. pulicaris groups in terms of vector abundance required for transmission to occur. These features could be incorporated in the model by extending it to include multiple vector species or, alternatively, to use location-specific parameter values based on local species composition. There are 24 known serotypes of BTV (Mertens 1999) and, in this study, data for several serotypes have been combined when estimating certain parameters, for example, the duration of host viraemia or the EIP (see references in table 1 for details). However, there may be differences among serotypes and, in particular, interactions between host, virus and vector, which influence the ability of BTV to invade a host population. For example, BTV serotypes 1, 2, 4, 9 and 16 have all been introduced to Europe, but have not extended much beyond the northern rim of the Mediterranean basin (Mellor & Wittman 2002;Purse et al. 2005), while BTV serotype 8 has appeared and spread in northern Europe (Albers et al. 2007).
The principal motivation for computing the basic reproduction number, R 0 , was to provide a quantitative framework in which to assess the risk of BTV to UK livestock. Uncertainty analysis (figure 1) provides evidence that there is a risk (i.e. there were parameter combinations which yielded R 0 O1), which is perhaps unsurprising in the light of the recent BT outbreak in northern Europe (Albers et al. 2007). Furthermore, the analysis suggested that the risk is greatest when the temperature is between 15 and 258C (figure 1b,d, f ). At lower temperatures, BTV is unable to replicate to transmissible levels, while at higher temperatures the vector is less likely to survive for long enough to complete the EIP. However, the transmission model ignores a large number of heterogeneities which will influence both spatial and temporal variation in R 0 , for example, the distribution of ruminant and vector species, and regional differences in environmental factors, in particular, temperature. A full risk assessment will need to take these factors into account and produce spatio-temporal maps for R 0 .
In addition to its use in predicting risk, the basic reproduction number can also be used to inform control strategies. In particular, those measures which affect parameters identified as most important by the sensitivity analysis are most likely to have the greatest effect. Any measures which reduce the biting rate, for example, could potentially reduce transmission and, hence, R 0 . This could possibly be achieved through the use of dips or pour-ons Doherty et al. 2004) and, consequently, it is essential to estimate the reduction in biting rate that can be achieved when using these treatments to assess whether or not they are likely to be effective control measures during a BTV outbreak.
The focus of the current study has been the risk of an invasion of BTV occurring in British livestock and, as such, it does not address the question of whether or not BTV could become endemic in UK. The transmission model, (2.1)-(2.6), could be used to assess the risk of persistence, but to do so will require a greater understanding of the dynamics of the vector population. Moreover, consideration would need to be given to mechanisms by which BTV could overwinter, such as covert persistence in immune cells in the vertebrate host ( Takamatsu et al. 2003) or persistence in the larval stages of the vector, following the transovarial transmission (White et al. 2005).