Two-Host, Two-Vector Basic Reproduction Ratio (R 0) for Bluetongue

Mathematical formulations for the basic reproduction ratio (R 0) exist for several vector-borne diseases. Generally, these are based on models of one-host, one-vector systems or two-host, one-vector systems. For many vector borne diseases, however, two or more vector species often co-occur and, therefore, there is a need for more complex formulations. Here we derive a two-host, two-vector formulation for the R 0 of bluetongue, a vector-borne infection of ruminants that can have serious economic consequences; since 1998 for example, it has led to the deaths of well over 1 million sheep in Europe alone. We illustrate our results by considering the situation in South Africa, where there are two major hosts (sheep, cattle) and two vector species with differing ecologies and competencies as vectors, for which good data exist. We investigate the effects on R 0 of differences in vector abundance, vector competence and vector host preference between vector species. Our results indicate that R 0 can be underestimated if we assume that there is only one vector transmitting the infection (when there are in fact two or more) and/or vector host preferences are overlooked (unless the preferred host is less beneficial or more abundant). The two-host, one-vector formula provides a good approximation when the level of cross-infection between vector species is very small. As this approaches the level of intraspecies infection, a combination of the two-host, one-vector R 0 for each vector species becomes a better estimate. Otherwise, particularly when the level of cross-infection is high, the two-host, two-vector formula is required for accurate estimation of R 0. Our results are equally relevant to Europe, where at least two vector species, which co-occur in parts of the south, have been implicated in the recent epizootic of bluetongue.


Introduction
Mathematical formulations for the basic reproduction ratio (R 0 ) -defined as the average number of secondary infections produced by a typical primary infection in an otherwise totally susceptible population [1] -exist for several vector-borne diseases including those with one host and one vector, such as malaria [2] and those with two hosts and one vector, such as zoonotic trypanosomiasis [3], African horse sickness [4] and bluetongue [5,6]. To date, with the exception of Lopez et al. [7], almost no attention has been paid to developing mathematical formulations of R 0 where there are both multiple hosts and multiple vectors. However, this is a common situation: for trypanosomiasis in Africa for example, two or more species of tsetse fly vector often co-exist; while for both African horse sickness and bluetongue in southern Africa, two competent vectors (Culicoides imicola and C. bolitinos) are frequently trapped together. Other diseases transmitted by multiple vectors include dengue [8], Japanese encephalitis [9] and malaria [10].
This situation may also apply to the recent European outbreak of bluetongue, which has caused the deaths of well over a million sheep. The outbreak began in 1998 in regions of southern Europe where the Afrotropical midge, C. imicola, occurs. Starting in 1999, it was also detected in Balkan countries where C. imicola was not known, thereby implicating local Culicoides species, such as the obsoletus and pulicaris groups, as vectors. Since these co-occur with C. imicola over the latter's European range [11], it was reasonable to suspect that they may transmit the virus alongside C. imicola in some places; and epidemiological evidence for this was later provided in Sicily [12]. Subsequently, both BTV 1 and BTV 8 have been transmitted in regions with indigenous vectors, both with and without C. imicola. It is therefore quite likely that two or more vector species have co-transmitted BT virus in several parts of Europe in recent years.
Given the likely widespread existence of multivector, multihost disease systems, we derive and analyse R 0 for the simplest of these: a two-host, two-vector system, using bluetongue as an example. We particularly wish to investigate the effect on R 0 of measurable parameter values relating to vector abundance, vector competence and vector host preference. In order to do this we extend the work of Lopez et al. [7] by first including vector host preference and temperature-dependent transmission parameters and then studying the effect of specific parameters. We consider the effect of the following: (1) vector to host ratio, which is linked to vector abundance and varies with species and temperature, as evidenced by C. imicola and C. bolitinos in RSA [13,14]; (2) probability of transmission from host to vector, which is linked to vector competence and also varies with species and temperature [15]; (3) vector host preference, which differs between species [14,16,17]. Importantly, we also work directly with R 0 , rather than the threshold T proposed by Lopez et al. [7], which is not valid in all regions of feasible parameter space. The notation we use allows us to make direct comparisons with a previously published two-host, one-vector formula for bluetongue. We illustrate our results by parameterising the model for a specific disease system, namely  [16,17], but prefers cattle [23]. C. bolitinos feeds on cattle and horses [16,17] and breeds in cattle dung [14].
m ij N j /H i ratio of vectors of type j to hosts of type i Many areas: m C1 = m S1 = 500 (0-5000) m C2 = m S2 = 50 (0-500) Colder high-lying areas: m C1 = m S1 = 50 (0-100) m C2 = m S2 = 500 (0-5000) In general, C. imicola is approx. 10 times more abundant than C. bolitinos [15]. In colder, high-lying areas, C. imicola is approx. 10 times less abundant than C. bolitinos [14]. bluetongue in South Africa. We use the situation in South Africa because of the availability of extensive distribution data, together with detailed experimental results on the relative vector competencies of the two main vector species [15]. Similar data for different European bluetongue vectors do not exist. It is known that several European vector species transmit bluetongue virus and that there are differences in host preference between these species. For example, Garros et al. [18] show that C. chiopterus prefers to feed on cattle while C. obsoletus is more of a generalist. However, nothing is known of their respective vector competencies. This highlights the need for a two-host, two-vector formula for R 0 as well as experimental work to establish the vector competence of each species. Although we have focussed on the situation in South Africa, the framework and general results presented here are equally relevant to Europe.

Model Equations
Equations describing the dynamics of a two-host, two-vector system are given below, whilst the variables and parameters of the model are defined and described in Table 1. For clarity, we have adopted a similar notation to that used by Gubbins et al. [6]. In short, hosts can be either susceptible, infectious or recovered (and in this case immune), whilst vectors can be either susceptible, latent or infectious. The proportions of susceptible, infectious and recovered hosts are denoted by x i , y i and z i respectively, whilst the numbers of susceptible, latent and infectious vectors are denoted by S j , L j and I j respectively (N j in total). Susceptible hosts of type i [where i can be either C (cattle) or S (sheep)] become infectious at rate l Hi , which is the sum over vector types indicated by j [where j can be either 1 (C. imicola) or 2 (C. bolitinos)] of b j a j (w ij m ij I j =N j ). The third term is composed of m ij the ratio of vectors of type j to hosts of type i, I j N j the proportion of vectors of type j that are infectious and w ij the proportion of vectors of type j attracted to hosts of type i (i.e. reflecting the preference of vector type j for host type i). So, the third term gives the average number of infectious vectors of type j attracted to a host of type i (after taking into account vector type j's preference for host type i). This is multiplied by a j , the (temperature-dependent) biting rate of vectors of type j, and b j , the probability of transmission from a vector of type j to a host given an effective contact. Similarly, susceptible vectors of type j become latent at rate l Vj , which is the sum over host types (indicated by i) of b j a j (w ij y i ). The third term is the probability of a vector of type j being attracted to an infectious host of type i. This is multiplied by a j , the (temperature-dependent) biting rate of vectors of type j, and b j , the probability of transmission from a host to a vector of type j given an effective contact. An infectious host remains infectious until it either recovers (at rate r i ) or dies (at rate d i ). After a short extrinsic incubation period (on average 1/n j ), latent vectors become infectious. They remain infectious until they die, which occurs at rate m j . Susceptible vectors are added to the system at rate r j . The model assumes that there is no seasonal aspect to vector recruitment or population size and that there is no latent period in hosts, recovered animals are immune and the host population remains constant except for losses due to diseaseinduced mortality. Hosts. Vectors.

Basic Reproduction Ratio
The ability of the pathogen to spread can be expressed in terms of the basic reproduction ratio R 0 . Mathematically, R 0 is the dominant eigenvalue of the next-generation matrix K. For vectorborne transmission models like the one described above, where matrix A describes vector to host transmission and matrix B describes host to vector transmission (see Appendix 1 in File S1). We could work directly with the characteristic equation DK{lID~0. However, there are significant advantages in using a result shown in Appendix 2 in File S1, namely that R 0 is the square root of the dominant eigenvalue of BA (a 464 submatrix of K 2 ). Not only is BA smaller than K but also its elements have an obvious biological interpretation in terms of R ij , the average number of infectious vectors of type i produced by one infectious vector of type j (necessarily in two generations). It is such biological interpretation that we seek. The utility of working with BA is doubtless associated with the argument that, in contrast to directlytransmitted infections, for vector-borne infections R 2 0 makes more sense biologically [2] and is in fact what is measured in the field (i.e. two-generation 'like' to 'like' transmission).
Following the above procedure we find that where specifically Two-Host, Two-Vector Basic Reproduction Ratio Note that the proportion of vectors of type j attracted to hosts of type i where v ij is a measure of vector type j's preference for host type i and H i is the total number of hosts of type i. For two host species (cattle C and sheep S), this can be rewritten as where vector host preference is now denoted by s j . In terms of vector to host ratios, where m Cj~Nj =H C and m Sj~Nj =H S , the first of these is As Gubbins et al. [6], for clarity we replace w Cj and w Sj with w j and 1{w j respectively. From Gubbins et al. [6], we know that R 0~ffi ffiffiffiffiffiffi ffi R 11 p and R 0~ffi ffiffiffiffiffiffi ffi R 22 p for two-host, one-vector systems involving vector type 1 and vector type 2, respectively. For the two-host, two-vector system we find from equations (2) that when w 1~w2 , . If vectors 1 and 2 are identical in terms of parameter values (i.e. all parameter values for vector species 1 equal those for vector species 2) and equal in number, then R 11~R22~R12~R21 and so R 0~ffi ffiffiffiffiffiffiffiffiffi 2R 11 p . In other words, the two-host, two-vector R 0 is greater than the two-host, one-vector R 0 by a factor of ffiffi ffi 2 p . So, in this case, if it were assumed that there was only one vector species transmitting the infection, the basic reproduction ratio would be underestimated (because the average number of relevant vectors per host would be underestimated). When the vectors are not identical, as the level of cross-infection R 12 R 21 increases, R 0 increases from ffiffiffiffiffiffiffi ffi R 11 p (when R 12 R 21~0 and R 11 wR 22 ), through ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Results
For the two-host, one-vector system, Gubbins et al. [6] found that the parameters with a significant effect on R 0 were temperature (T) [via biting rate, extrinsic incubation period (EIP), vector mortality rate], the probability of transmission from host to vector (b) [which was not temperature-dependent in Gubbins et al. [6]] and the ratios of vectors to hosts (m C and m S ).
For the two-host, two-vector system, we propose to focus on the effects on R 0 of varying the ratios of vectors to hosts (m C1 , m C2 , m S1 , m S2 ) [linked to vector abundance], the probabilities of transmission from host to vector (b 1 , b 2 ) [linked to vector competence and temperature-dependent in our model] and the vector host preferences (s 1 , s 2 ). Our aim is to provide a general framework for two-host, two-vector approaches to bluetongue; however there is a paucity of data. There is one situation, South Africa, where C. imicola and C. bolitinos coexist across most of the country and for which we do have data. We undertake the analysis with reference to this.
As shown by Meiswinkel et al. [13], there are many areas (e.g. Western Cape, western part of the Eastern Cape, Mpumalanga, Gauteng and Limpopo Province) where C. imicola is 10 to 100 times more abundant than C. bolitinos. However, there are areas, in particular in the cooler high-lying areas of the Free State, where C. bolitinos is approximately 10 times more abundant than C. imicola [14] and Venter et al. [19] suggest that C. bolitinos may play an important role in the transmission of BTV in these areas. Paweska et al. [15] demonstrate that, regardless of incubation temperature (10, 15, 18, 23.5 or 30uC), the mean virus titre/midge, infection rate and proportion of infected females with transmission potential (i.e. virus titre/midge $10 3 TCID 50 , where TCID 50 (tissue culture infectious dose 50) is the amount of virus that will infect 50% of midges inoculated with it) are significantly higher in C. bolitinos than in C. imicola and suggest that, because of its significantly higher vector competence, C. bolitinos could be the primary vector in areas where it occurs in lower numbers than C. imicola, as well as in these cooler regions. Here, abundance is expressed through the ratios of vectors to hosts (m C1 , m C2 , m S1 , m S2 ), while vector competence is expressed through the probabilities of transmission from host to vector (b 1 , b 2 ). Regarding vector host preferences, there is evidence [14,16] that many Culicoides species prefer to feed on cattle and some suggestion that C. bolitinos may not feed on sheep at all [16,17].

Estimating b 1 and b 2
Fu et al. [20] show that only midges containing $10 3 TCID 50 release detectable amounts of virus in their saliva. So, first we define 'infectious' as 'having a virus titre $10 3 TCID 50 '. Next we obtain from Table 2 of Paweska et al. [15], for several different temperatures, the proportion of vectors remaining that are infectious [i.e. (number of infectious vectors)/(number of initially susceptible vectors known to have fed on infected blood and still be alive after the incubation period)]. Each of these data points is equal to b j at a given temperature. By fitting curves to the data, we can find temperature-dependent functions for b 1 (C. imicola) and b 2 (C. bolitinos). Exponential curves of the form b j~pj exp(q j T) were fitted using a nonlinear least-squares method with bisquare weighting of the residuals. The coefficients and goodness of fit statistics are given in Table 2. The curves (shown in Appendix 3 in File S1) adequately describe the relationships between b 1 , b 2 and temperature over this range of temperatures. Note that as temperature varies from 10 to 30uC, the ratio b 2 /b 1 varies from 9.85 to 12.91 (i.e. the probability of transmission of C. bolitinos is always about 10 times greater than that of C. imicola).

Other parameter estimates
The estimates for m C1 , m C2 , m S1 and m S2 were based on catch sizes and species composition reported in Venter & Meiswinkel [14] and Venter et al. [17]. They are rough estimates designed to reflect the relative orders of magnitude of each vector species. As in Guis et al. [21], we assume that catch size approximates ratio of vectors to hosts. Estimates for the remaining parameters were taken from Gubbins et al. [6]. Details are given in Table 1. Six of these estimates a j , n j and m j (i.e. three for each vector species) depend on temperature. They are positive and increase monotonically for temperatures between 10.4uC and 35.5uC.

Effect of differences in m ij and b j
In order to focus on the effects on R 0 of differences in vector competence, vector abundance and vector host preference, the two species of vectors are assumed to be the same in every way except b j , m ij and s j . We first consider the effect of differences in m ij for fixed b j and s j .
First note that, as shown in equation (3), w Cj varies with m Cj , m Sj and s j . However, when m Cj equals m Sj , w Cj (and hence w j ) depends on s j alone. In Figure 1 (and Figure 2 below), m S1 and m S2 equal m C1 and m C2 respectively and s 1 = s 2 = 0.5. Consequently w 1 and w 2 are fixed at 0.67. The transmission probabilities b 1 and b 2 are determined (as described in Table 1) by temperature, which is 25uC in Figure 1A and 15uC in Figure 1B. The parameters m C1 and m C2 vary independently along the x and y axes respectively. In Figure 1A, we can clearly see that R 0 is greater when vector 2 (C. bolitinos) is more abundant than vector 1 (C. imicola) [i.e. when m C2 is greater than m C1 ] than when the reverse is true. For example, when m C1 = 50 and m C2 = 500 (i.e. in the top left-hand corner), R 0 is 7.0, whereas, when m C1 = 500 and m C2 = 50 (i.e. in the bottom right-hand corner), R 0 is only 3.1. This large difference is due to the fact that b 2 is approximately 10 times greater than b 1 and illustrates the balance between vector abundance and vector competence for C. imicola and C. bolitinos. The same relationship is observed when the temperature is 15uC (i.e. in Figure 1B), but R 0 is much smaller at this temperature.
It is also clear from Figure 1 that omitting one vector species (i.e. being constrained to one axis) leads to underestimation of R 0 and when that species has a significantly higher vector competence (as does vector species 2) the degree of underestimation can be dramatic.
In Figure 2, b 1 and b 2 vary with temperature, as described in Table 1. The vector to host ratios m S1 and m S2 equal m C1 and m C2 respectively, which vary simultaneously with x as described below: for 0ƒxƒ1: Relative abundance (m C1 /m C2 = N 1 /N 2 ) is therefore described by the hyperbola m C1 m C2~5

50z450x
{1: Figure 2A shows how R 0 varies with relative abundance. Curves 1 (green) and 2 (blue) show the relationship when temperature is fixed at 15uC and 25uC, respectively. Curve 3 (red) is produced when temperature (T~25{10x for 0ƒxƒ1) and relative abundance vary simultaneously with x. Hence, curve 3 (red) shows the relationship when temperature varies from 15uC to 25uC as relative abundance varies from 0.1 (when C. bolitinos is 10 times more abundant than C. imicola) to 10 (when C. imicola is 10 times more abundant than C. bolitinos). By definition, curve 3 is constrained to start at the same point as curve 1 and end at the same point as curve 2. Curve 3 can be thought of as a path across the landscape, moving from the cooler high-lying regions where C. bolitinos dominates to the warmer low-lying regions where C. imicola dominates. Along this path, temperature (and hence b j ) and relative abundance vary simultaneously. Figure 2B shows how R 0 varies with relative abundance (y axis) and temperature (x axis) separately. It clearly shows that, for a fixed temperature, R 0 always decreases as relative abundance of C. imicola increases. It also shows that, for fixed relative abundance, R 0 initially increases with temperature, but starts to decrease again beyond 31uC. Curve 3 in Figure 2A corresponds to moving from A to B across the surface in Figure 2B. Along this path, the highest R 0 corresponds to a relative abundance of approximately 1.4 and a temperature of approximately 21.1uC. In this case, for temperatures greater than 21.1uC, R 0 drops with rising temperature because, at the same time, the less competent vector is replacing the more competent one.
In summary, we find that high vector competence can compensate for low vector abundance and that temperature, which determines the transmission probabilities b 1 and b 2 and also influences the abundance and composition of vector species, has a marked effect on R 0 .

Effect of differences in s j
We now consider the important effect that vector host preference (s j ) has on R 0 . In order to focus on the effect of s j , we assume that the vector species differ only in s j . Also, to ensure that there is no advantage to choosing cattle over sheep (or vice versa), we set r C = r S , d C = d S and m C1 = m S1 = m C2 = m S2 . When s j = 0, the proportion of vectors of type j attracted to cattle (w j ) equals 1. When s j = 1 (i.e. no preference), w j just depends on the relative numbers of each host species, with a greater number resulting in a greater share of the vectors. As s j R', w j R0 where all vectors are attracted to sheep. To prevent loss of detail as s j R', in Figure 3 we use a j rather than s j , where a j is an alternative measure of vector host preference such that s j~aj =(1{a j ). From this formula we can see that as a j varies from 0 to 1, s j varies from 0 to ' and that a j = 0.5 indicates no preference. Figure 3 clearly shows that the minimum value of R 0 lies on a straight line running from (a 1 = 0, a 2 = 1), where w 1~1 and w 2~0 (i.e. vector type 1 feeds exclusively on cattle while vector type 2 feeds exclusively on sheep), through (a 1 = 0.5, a 2 = 0.5), where w 1~0 :5 and w 2~0 :5 (i.e. neither vector has a preference and so both vector species are equally distributed between both host species), to (a 1 = 1, a 2 = 0), where w 1~0 and w 2~1 (i.e. vector 1 feeds exclusively on sheep while vector 2 feeds exclusively on Table 2. Coefficients and goodness of fit statistics for exponential curves of the form b j~pj exp(q j T), where j can be 1 (C. imicola) or 2 (C. bolitinos), fitted to data extracted from Paweska et al. [15]. cattle). Any deviation from this line results in a higher R 0 . This figure clearly shows two things: firstly, that different combinations of vector host preferences can result in the same R 0 ; second, that when both vectors prefer the same host species, R 0 is greater. This result is important because it shows that, when the vector species differ only in host preference and the host species are equally good as hosts (in this case, they share the same infectious period and the same pathogen-induced mortality rate) and equally abundant, overlooking vector host preference can result in an underestimation of R 0 . In Figures 1 and 2 and Table 3, we used s 1 = s 2 = 0.5 (which corresponds to a 1 = a 2 = 1/3) as many Culicoides species prefer to feed on cattle [14,16]. However, while it is clear that C. imicola also feeds on sheep (and even horses and pigs too), there is some evidence that C. bolitinos does not -instead feeding exclusively on cattle and horses [16,17]. A strong association between C. bolitinos and cattle is further suggested by the fact that C. bolitinos breeds in cattle dung [14], rather than soil like C. imicola. In terms of R 0 , if C. bolitinos does not feed on sheep, then s 2 = 0 (i.e. a 2 = 0) and the true value of R 0 will be higher than our estimates based on s 2 = 0.5.

R 0 approximations
The need for the two-host, two-vector formula is further emphasised when we consider several approximations based on the two-host, one-vector formula for R 0 , which is Suppose we have information on two vectors that are circulating in the same area and feeding on the same host populations. We might think it reasonable to assume that the vectors are acting independently, merely feeding on the same hosts. In which case, one option would be to calculate R 0 for each vector separately and add them together. We refer to this approximation as R 0,sum (i.e. R 0,sum~ffi ffiffiffiffiffiffi ffi R 11 p z ffiffiffiffiffiffiffi ffi R 22 p ). It incorporates the idea that there is no cross-infection. Table 3 contains the true value of R 0 (i.e. calculated using the two host, two vector formula) and the value of R 0,sum under different scenarios. In examples a and b the temperature is 25uC and C. imicola is 10 times more abundant than C. bolitinos (representing warmer low-lying areas), whereas in examples c and d the temperature is 15uC and C. bolitinos is 10 times more abundant than C. imicola (representing cooler highlying areas). In a and c, the vectors are assumed to have the same preference for cattle (i.e. s 1 = s 2 = 0.5) so w 1~w2 . In b and d, C. imicola is assumed to have no preference for a particular host, while C. bolitinos is assumed to feed exclusively on cattle (i.e. s 1 = 1,  . Effect on R 0 of relative abundance and temperature. In (A) R 0 is plotted against relative abundance (m C1 /m C2 = N 1 /N 2 ), which varies from 0.1 (when C. bolitinos is 10 times more abundant than C. imicola) to 10 (when C. imicola is 10 times more abundant than C. bolitinos). Temperature is either fixed at 25uC or 15uC or varies from 15uC to 25uC as relative abundance varies from 0.1 to 10. In (B) R 0 is plotted against ln(relative abundance) and temperature. Parameter values (1 = C. imicola, 2 = C. bolitinos): b 1 = 0.9, b 2 = 0.9, s 1 = 0.5, s 2 = 0.5, r C = 1/20.6, r S = 1/16.4, d C = 0, d S = 0.005, m S1 = m C1 , m S2 = m C2 , a 1 , a 2 , m 1 , m 2 , n 1 , n 2 , b 1 and b 2 are determined by temperature. doi:10.1371/journal.pone.0053128.g002 In these examples, R 0,sum consistently overestimates R 0 by between 5% and 45%.
An alternative approach would be to pool the vectors (e.g. m C~mC1 zm C2 and m S~mS1 zm S2 ) and use average values (e.g. b~(b 1 zb 2 )=2 and w~(w 1 zw 2 )=2). In the examples in Table 3, we have assumed that the vectors are very similar and so share many parameter values. In fact, we have assumed that they differ only in vector to host ratio (m ij ), host to vector transmission probability (b j ) and vector host preference (s j ). So, R 0,ave can be obtained by substituting m C , m S , b and w into equation (4). Surprisingly, the examples reveal that R 0,ave can sometimes overestimate and sometimes underestimate the true value by a significant amount.
Another possible approximation is obtained by first calculating the lower and upper bounds given by R 0,minb and R 0,maxb and then taking the weighted average (R 0,wtsum ). R 0,minb is calculated in the same way as R 0,ave except that the host to vector transmission probability (b) takes the minimum value (b 1 , where b 1 vb 2 ), rather than the average, and the proportion of vectors attracted to cattle (w) takes the minimum value (w 1 , where w 1 vw 2 ), rather than the average. R 0,maxb is the equivalent calculation using the maximum host to vector transmission probability (in this case b 2 ) and the maximum proportion of vectors attracted to cattle (in this case w 2 ) . R 0,wtsum is then the weighted sum of R 0,minb and R 0,maxb (i.e. ½m 1 =(m 1 zm 2 )R 0,min b z½m 2 =(m 1 zm 2 )R 0,max b , where m 1~mC1 zm S1 and m 2~mC2 zm S2 ). We can see from Table 3 that R 0,wtsum can provide a fairly good approximation to R 0 . In our examples, it consistently underestimates R 0 , but never by more than 19% and sometimes by as little as 2%. Alternative formulations in terms of R 11 and R 22 are given in Table 3 for comparison with R 0,sum . Note however that, for R 0,ave , R 0,minb , R 0,maxb and R 0,wtsum , the formula given is for w 1~w2 only. There is insufficient space to give the more general expression.
These examples suggest that, even when the vectors are very similar and share many parameter values, simply summing the contribution from each vector species (R 0,sum ) will lead to overestimation of R 0 and that the degree of overestimation can be large. R 0,wtsum appears to provide a more consistent estimate. Table 3 also shows that intuitive approximations like R 0,ave can be very misleading, sometimes underestimating and sometimes overestimating the true value of R 0 .

Discussion
We have presented an expression for R 0 for a two-host, twovector system and demonstrated its sensitivity to parameters relating to vector abundance, vector competence and vector host preference. We have shown that high vector competence can offset low vector abundance and that, where high vector competence and high vector abundance coincide, R 0 can reach high values. We have also shown that the highest value of R 0 does not always coincide with the highest b j values. Earlier work using a one-host, one-vector formulation showed that when a, m and n vary with temperature, R 0 at first increases with temperature then decreases [22]. We observed the same behaviour when using the slightly different temperature-dependent functions described in Table 1. Figure 2B shows that this relationship is maintained when b 1 and b 2 also increase with temperature.
As shown in Figure 3, vector host preference has an interesting effect on R 0 . When the vector species differ only in host preference and the host species are equally good as hosts and equally   { For the approximations R 0,minb , R 0,wtsum , R 0,ave and R 0,maxb , the formula given is for w 1~w2 only. There is insufficient space to give the more general expressions. doi:10.1371/journal.pone.0053128.t003 abundant, a preference for one host species can increase R 0 if the total feeding rate is maintained. When both vectors prefer the same host species, R 0 will increase. When the preferred host benefits transmission (e.g. by having a longer infectious period, like cattle with bluetongue), then R 0 will increase further. However, if the preferred host is less beneficial or more abundant, then R 0 will decrease. In this model, the vector species do not directly interact. They merely feed upon the same pool of susceptible hosts. So, we might expect a simpler formulation expressed in terms of the two-host, one-vector R 0 for each species to provide a good approximation to R 0 . We considered several possibilities and found that simply summing the contribution from each vector species (R 0,sum ) leads to overestimation of R 0 , while using average values (R 0,ave ) can lead to under or overestimation. A more consistent estimate was provided by R 0,wtsum . However, this approximation relies on the fact that the vectors differ in m ij , b j and s j only. When the vectors differ in many ways, we can see from equation (1) that the twohost, one-vector formula will provide a good approximation when the level of cross-infection between vector species is very small. As this approaches the level of intraspecies infection, a combination of the two-host, one-vector R 0 for each vector species (i.e. ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 11 zR 22 p ) becomes a better estimate. Otherwise, particularly when the level of cross-infection is high, the two-host, two-vector formula is required for accurate estimation of R 0 . The results of this work demonstrate the need for a two-host, two-vector formula for R 0 in areas that support two significant vectors, particularly where those vectors differ in many ways. Further extensions of this model would be required for areas where there were more than two important vectors. Northern Europe could be one such area. Both C. pulicaris and C. obsoletus transmit bluetongue is this region. However, both of these vectors are in fact vector species groups containing multiple vector species (e.g. the C. obsoletus group contains four distinct vector species). At the moment, there is insufficient information about differences in vector competence between these species to be able to use this R 0 formula (or an extension of it) in this region.

Supporting Information
File S1 Supporting information.