Model overview

The adult lives of the vectors are modelled as a series of independent and identically distributed (i.i.d) random waiting periods—the gonotrophic cycles—which we label G_n (n ≥ 0). These define the time between each successful bite until death, and death is assumed to occur at constant rate μ (Fig 1). This means that the lifetime biting process of the vector follows a renewal process [15] stopped at the end of the vector’s life. The starting point of the gonotrophic cycles for each vector is taken to be its emergence as an adult. If a vector is inoculated with the pathogen, then only bites after a random duration E (the EIP) are potentially infectious. The product of the probabilities of transmission from infectious host to vector and from infectious vector to host (the vector competence) is given as V. The classical assumption is that the G_n are exponentially distributed (G ∼ exp(α)) and E = 1/σ is constant. Here we drop these assumptions and generate a generalised reproductive ratio () for a single midge species spreading a pathogen amongst a single host species (Fig 1). Assuming that the average host infectious duration is 1/γ, expressed in explicitly probabilistic terms: (3) The factorisation in eq (3) is appropriate because the future biting of a vector after its inoculation can be taken as independent of its past life-history. The probability of a vector surviving its EIP duration and the average biting rate over the vector population can be solved in terms of the moment generating functions (MGFs) of the random life-cycle distributions for the vector. The MGF for a random variable X is defined as, (4) From an entomological modelling perspective the MGF evaluated at θ > 0 is the probability of a vector surviving a random period defined by X whilst undergoing constant mortality rate θ. The two MGFs we will use in our analysis are ϕ_E, the MGF of the EIP duration, and ϕ_G, the MGF of the gonotrophic cycle duration. Therefore the probability of a vector surviving its EIP, (5) The classical assumption that the EIP is fixed at σ⁻¹ for each inoculated vector in fact leads to the smallest estimate of P_E over any choice of distribution for the EIP with average duration σ⁻¹ [eq (15)]; that is the classical estimate is almost certainly an over-estimate for the proportion of inoculated vectors that die before becoming actively infectious.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Reproductive ratio and vector life history examples.

The generalised reproductive ratio, , is derived by computing the expected number of vertebrate hosts infected through one complete generation of the parasite due to a single infected host. The average quantities underlying depend upon assumptions on the chance life events of individual vectors. Three example realisations of vector lives are given for three extremal modelling assumptions (Right). 1) The classical model. The vector waits for an exponentially distributed period between bites (circles), once inoculated with the pathogen the vector has a fixed latency period (green line) until the end of its EIP after which it makes infectious bites. 2) The constant model. Both the EIP duration and the time between bites are fixed. 3) The Markov model. Both the time between bites and the EIP duration are exponentially distributed. In each case the vector’s life ends with a fixed per-capita mortality rate.

https://doi.org/10.1371/journal.pcbi.1004837.g001

The average number of bites, or equivalently survived gonotrophic cycles, in the lifespan of a vector is, (6) A traditional modelling assumption is that each vector bites according to a Poisson process at rate α(each gonotrophic cycle is exponentially distributed with a mean cycle duration of 1/α) giving a population biting rate per host of Mα where M is the vector-to-host population ratio. A convenient feature of modelling vector bites as a Poisson process is that the chance of the vector biting at any moment in time is completely independent of the vector age, a, since emerging as an adult. This is not true in general; it is biologically plausible that having bitten successfully a vector is less likely to bite shortly afterwards, even if the average time between bites remains 1/α. This could be experimentally confirmed by observing that gonotrophic cycles are under-dispersed compared to an exponential distribution. As a consequence vectors of different ages will have different chances of biting: vector age dependence in the model emerges from considering non-exponentially distributed gonotrophic cycles without any assumption on vector senescence. We will demonstrate that for a general choice of gonotrophic cycle distribution the equilibrium biting rate of the vector population upon hosts is the expected total number of bites made in per vector lifetime divided by the vector life expectancy. This has a concise mathematical expression: (7) Eq (7) is consistent with the classical prediction of biting rating from the vector population in that if the biting process is Poisson (waiting times G are exponentially distributed) then, (8) Although eq (7) holds in general, eq (8) is specific to the Poisson assumption.

A second quantity that is key to calculating the basic reproductive ratio is the expected number of infectious vector bites. Infectious bites are those that occur after the latent EIP duration—we denote this quantity B_I (see eq (20) below). In order to make at least one infectious bite an inoculated vector must survive its EIP and any remaining gonotrophic period after its EIP finishes. This implies three key considerations:

1. A significant proportion of inoculated vectors (1 − P_E; eq (5)) will die before the end of their EIP and contribute no infectious bites.
2. Inoculated vectors that survive their EIP might still die before they complete their remaining gonotrophic period between the end of the EIP and the vector’s first infectious bite.
3. The remaining period will be typically shorter than a complete gonotrophic cycle. Therefore, an infectious vector that survives its EIP will, on average, achieve more infectious bites than a typical emerging adult vector.

Because a vector’s inoculating bite also initiates its EIP it is possibly to calculate B_I exactly using an integral formula [eq (20)]. We can also approximate B_I explicitly, using the asymptotic distribution of remaining time until next bite, (9) The approximation is exact if the gonotrophic cycle is exponentially distributed.

Putting Eqs (7) and (9) together gives the generalised reproductive ratio: (10) and the approximate generalised reproductive ratio (): (11) Eq (11) explains the discrepancy in the introduction. If the EIP duration is fixed at E = 1/σ then ϕ_E(μ) = exp{−μ/σ} which recovers the classical [eq (1)] whereas assuming an exponentially distributed EIP (ϕ_E(μ) = σ/(σ + μ)) recovers the alternative R₀ [eq (2)]. A key point is that the generalised predictions for the reproductive number depends on the MGFs of the EIP and gonotrophic cycle. MGFs depend upon the full distribution, therefore can be estimated using all the information from an entomological study rather than just summary statistics such as mean (or variance) of gonotrophic cycle duration. In Results we will concentrate on gamma distributed durations which are determined by their mean and dispersion around the mean, however there is no theoretical necessity to use these distributions in general.

Vector life duration and survival probabilities

Each vector has an age a, which defines the time since its first gonotrophic cycle was initiated after the end of its juvenile period. After reaching adulthood the vector has a random exponentially distributed lifetime L with mean described by the probability density function, (12) Exponentially distributed lifetimes, which derive from a constant hazard, have the advantage of a memoryless property. Having survived to some age a the vector’s expected remaining life is still 1/μ; that is the age of vector doesn’t effect its hazard of dying. Since we acknowledge that the age of the vector effects its biting rate we necessarily consider the age distribution of the vector population. In principle the vector age distribution could be seasonally varying, however we focus on the age distribution at population equilibrium. The age of a vector selected at random from its population at equilibrium is a random variable A. The equilibrium age probability density (f_A(a)) can be calculated using the microcosm principle for population processes [20] which states that the proportion of a population in some state is proportional to the mean lifetime an individual would spend in that state. In this case, the microcosm principle gives, (13) For exponentially distributed life-times, the equilibrium population density of vector age is equal to the life duration density i.e. f_L = f_A, but this is not generally true.

After successful inoculation with the pathogen the vector will become infectious after a random EIP duration E with density function f_E. The memoryless property of the vector lifetime distribution implies that the probability of surviving its EIP (P_E) is given by the MGF of the EIP duration evaluated at the vector mortality rate: (14) Where (.) is an indicator function, taking the value 1 if the statement is true or 0 otherwise. Since exp{−μx} is a convex function in x then Jensen’s inequality gives that for any distribution of E (c.f. inequality (2.23) in Klebaner [21]): (15) Eq (15) shows that the probability of an inoculated vector surviving a fixed EIP is the lowest possible estimate for any EIP distribution.

Vector biting process

Vector biting follows a natural repeating cycle: the vector finds a host to bite at the end of each of its gonotrophic cycles and, once successfully fed, initiates a new gonotrophic cycle which includes locating a suitable oviposition site, oviposition, oogenesis and then subsequent host seeking. We model the gonotrophic cycle durations as a collection of independent and identically distributed (i.i.d) random variables labelled G_n, each with a common density function, f_G. The (random) number of bites the vector has managed by age a, denoted B(a), is the same as the number of gonotrophic cycles completed, (16) Processes such as this that count the number of i.i.d. waiting periods completed up to a given time are known as renewal processes, and their probabilistic properties have been extensively analysed [15]. For vectors dying at constant rate μ, the chance of surviving each gonotrophic cycle is independent of the number of cycles survived in the past. Eq (14) gives the probability of a vector surviving its EIP, by exactly the same argument the vector survives each gonotrophic cycle with probability ϕ_G(μ) up until its death therefore the lifetime number of bites by a vector is distributed geometric, with mean ϕ_G(μ)(1 − ϕ_G(μ))⁻¹.

The Poisson process, which is associated with constant rates and hence ODE models, is an important special case of a renewal process. The implicit assumption that the average biting rate per vector is some fixed rate α, independent of the age and time since last bite of the vector, is equivalent to assuming that each gonotrophic cycle is exponentially distribution with average duration . However, we do not need to restrict our attention to this special case. Over a long life duration the elementary renewal theorem [15] gives that the expected bites per unit time will converge on the inverse of the average duration α, (17) Eq (17) potentially explains why it is uncommon in the mathematical epidemiology literature to model epidemic contact processes as general renewal processes, as opposed to Poisson processes, since if the life time of an individual is long compared to its inter-contact periods a Poisson process is a reasonable approximation. However, for short lived vectors this is not true; the average biting rate per vector will crucially depend upon the age distribution amongst the population.

We can now write the average biting rate, denoted α(a), for vectors at age a in terms of probability densities: (18) Where is the density function for the sum of gonotrophic cycles (see supporting information S1 Text). The equilibrium population average biting rate per vector is therefore, (19) Here we have used standard results on the MGF of sums of independent random variables and geometric sums. The equilibrium biting rate differs from the long term biting rate α because of the non-exponentially distributed gonotrophic cycles and the typically short lives of the vector, such that many vectors may die before they seek their first blood meal (Fig 2). Note that Eq (19) implies that in the long life limit irrespective of the gonotrophic cycle distribution, and hence results are independent of the assumed distribution.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The effect of dispersion on the important epidemiological quantities in VBDs.

The equilibrium average biting rate per midge (left), the probability of an inoculated midge surviving its EIP (right) for Culicoides genus midges. The mortality rate is μ = 0.22 per day (Culicoides genus midge mortality rate at 20°C [13]), the asymptotic biting rate α = 0.18 per day (Culicoides genus midge inverse mean gonotrophic cycle duration at 20°C [22]), mean EIP σ⁻¹ = 8.4 days (the posterior median for Culicoides genus midge mean EIP duration [19]). The red dots indicate the results derived from the classical Ross-Macdonald assumptions, other colour dots represent EIP dispersion estimates in the literature: BTV-1 spread by C. bolitinos midges [19] and Trypanosoma brucei brucei spread by tsetse flies [23]. Each quantity is found to be increasing with increasing dispersion.

https://doi.org/10.1371/journal.pcbi.1004837.g002

Vectors that become inoculated via a bloodmeal from an infectious host have bitten at least once, therefore the population equilibrium biting rate is not appropriate to model their subsequent biting process. The expected number of future bites by an inoculated vector after its EIP, denoted B_I, can be calculated directly as, (20) That is the expected infectious bites from an inoculated vector are due to the aggregate of the biting rate at all times t after the inoculating bite, weighted by both the probability that the vector is still alive at that time and has completed its EIP (Fig 3). Eq (20) is numerically solvable using that α(•) is the rate of change of the solution to the renewal equation with waiting time distribution G (see supporting information S1 Text). Probabilistically, the exact value for B_I is difficult to compute because the probability of the vector surviving from the end of its EIP until its next bite depends upon the distribution of E in a complex fashion. However, it is clear that only the small proportion of vectors that are by chance comparatively long-lived contribute to B_I when the EIP is typically longer than life expectancy. Therefore it can be argued that the asymptotic biting rate α (eq 17) is appropriate for such vectors conditional on having survived their EIP, i.e. (21) Eq (21) is exact when the gonotrophic cycle durations are exponentially distributed. In the supporting information we make the argument for Eq (21) as a general approximation to Eq (20) more rigorous by considering the long-time distribution of remaining time after the EIP until next bite. We also consider the special case of each gonotrophic cycle being of fixed, rather than random, duration (S1 Text).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Relative reproductive ratio and expected biting from infected vectors.

The generalised reproductive ratio for bluetongue virus relative to its classical estimate for a range of dispersions (d_G, d_E) (top left) and the relative error of the analytic approximation ) (top right) calculated using the same parameters as in Fig 2 (i.e. for midge life histories assuming a constant 20°C temperature). The lower parameter sets (middle/bottom, left and right) give illustrations of the expected life histories of inoculated midges at indicated dispersion parameters. The grey filled curve gives the expected biting rate of midges at each time after the inoculating bite conditional on their survival; the total grey area up to that time is the expected number of bites by the surviving midge. The expected number of actively infectious bites per inoculated midge (B_I) is calculated by integrating over the biting rate re-weighted at each time by the survival probability (blue curve; right hand axis) and the incubation probability (upper bars; green indicating a low probability and red indicating a high probability). The biting rate predicted by modelling bites as a Poisson process (d_G = 1) is shown (red curve). Low dispersion in gonotrophic cycles (parameter sets A,C) predicts that inoculated midge biting is expected to peak around multiples of the average gonotrophic cycle length (α⁻¹ = 5.5 days) whereas higher dispersion (parameter sets B,D) predicts rapid convergence to the Poisson process biting rate. Low dispersion in EIP predicts a sharp transition at which virtually all inoculated midges will become infectious (parameter sets C,D). In general, higher dispersion favours transmission however the details can be important. The approximation is accurate for all values of d_E only when d_G > 0.5, otherwise can be either an underestimate or overestimate of .

https://doi.org/10.1371/journal.pcbi.1004837.g003

Dispersion and R₀ for VBDs spread by midges

Using a gamma distribution to interpolate between exponentially distributed (when the gamma shape parameter is 1) and constant periods (when the gamma shape parameter diverges to infinity with fixed gamma mean) is a common modelling choice due to its comparative simplicity and flexibility with a pedigree in epidemic modelling [28]. The dispersion of a gamma distributed random variable X, d_X, is defined as, (23) The dispersion of the random period X interpolates between fixed (d_X = 0) and exponential (d_X = 1), it is also possible for a duration to be over-dispersed (d_X > 1) but we restrict attention to under-dispersal in this work. The dispersion is also usefully incorporated into the MGF of a gamma distributed random variable X: (24) Modelling E and G as gamma distributed with dispersions, d_E and d_G respectively allows us to calculate the population average biting rate per host and probability of surviving its EIP as (25) and (26) respectively. The mortality rate, mean gonotrophic cycle duration and mean EIP duration for BTV are all strongly temperature dependent for C. sonorensis with the reproductive ratio being maximal at around 20–25°C [9]. At 20°C, using realistic parameters, only a minority of midges are expected to survive the average EIP duration of BTV; the chance long-lived tail of the midge population are the sole contributors to the outbreak. However, the actual survival probability for the EIP varies from the 15.6% survival (d_E = 0; the classical assumption that all EIP durations are identical) to 35.0% survival (when d_E = 1). The population averaged biting rate per midge is similarly sensitive, and varies from 0.18 bites per midge per day (d_G = 1; the classical assumption that the gonotrophic cycle is exponentially distributed) to 0.09 bites per midge per day (d_G = 0) (Fig 2).

Turning our attention to the relative reproductive ratio (R_rel), and using the numerically exact form of B_I (eq (20) and Fig 3), we find that this increases with increasing dispersion of either distribution. As such, the reproductive ratio is maximised by both distributions being exponential (R_rel = 2.24) and is lowest when both are fixed (R_rel = 0.50) (Fig 3). The classical for BTV can therefore be either a significant overestimate or a significant underestimate compared to the true prediction that accounts for dispersion.

The approximate is exact when the midges are expected to bite according to a Poisson process (d_G = 1). Whenever the gonotrophic cycle duration has high dispersal (d_G > 0.5) the approximate reproductive ratio is a good estimate across the entire EIP dispersion range (Fig 3). However can be both an underestimate or an overestimate of when d_G < 0.5 depending on an complex interplay between the distributions of E and G. In particular, when the EIP duration is nearly exponential (d_E ≈1) and the gonotrophic cycle is under-dispersed (d_G < 0.5), with the worst relative error being when (d_E = 1, d_G = 0). This can be understood through the failure of the assumption that only chance long-lived midges with biting rates well approximated by the asymptotic rate α survive their EIP.

Controlling VBDs spread by midges

The most important feature of the reproductive ratio as a measure of transmission intensity is that it yields a threshold quantity for disease persistence, and hence provides a quantitative insight into control. If an agency tasked with disease control can cause the reproductive ratio to fall below unity for an extended period of time then eradication can be achieved locally. One method for reducing the reproductive ratio is by vaccinating the host population; in the case of BTV there was a mass vaccination of northern European cattle during 2008 which was ultimately effective at eliminating BTV incidence in northern Europe. However, a vaccine for a novel or unexpected strain might not be immediately available, as occurred in northern Europe in 2006 which was invaded by BTV-8 serotype rather than the BTV-2 or BTV-4 serotypes which circulate in southern Europe [17], in which case other control measures need to be considered.

Reducing the midge population by using insecticides would seem to be an obvious solution, if this could be done whilst respecting standards of user and environmental safety. Wide-spread application of insecticide around a farm has not been recorded as successfully causing significant local midge population reduction [29], but it has been suggested that targeting spraying of an insecticide with long residual life in proximity to cattle might be more successful [30]. To model the effect of such insecticide spraying proximate to cattle we assume that the midge mortality is increased by an excess mortality μ_e but the overall midge to host ratio (M) remains unchanged, due to the inability to target the entire midge population. If we make the classical assumptions (d_E = 0, d_G = 1), then the critical excess mortality () required to achieve eradication can be determined as the solution to: (27) This is derived from eq (1), by including the excess mortality and insisting that the resultant reproductive ratio is one.

For disease control using vaccination there is a simple relationship between the critical vaccination coverage of hosts and the reproductive ratio; such that the critical vaccination coverage amongst hosts implied by a classical estimate is . For eradication via increasing midge mortality the relationship between the critical excess mortality and the reproductive ratio estimate derived using the classical assumptions (e.g. [26]) solves, (28) Note that depends on the reproductive ratio estimate , the midge mortality rate μ and the incubation rate σ in contrast to the critical vaccination coverage of hosts which depends only on . A natural consideration is the effect that generalising the distribution of the EIP and the gonotrophic cycle durations has upon the predicted critical excess mortality or critical vaccination coverage needed for eradication.

We begin by considering excess mortality and assume that we have access to a classical estimate of the pre-control reproductive ratio . In the generalised setting, the midge population biting rate depends on the precise nature of the excess mortality so we consider two separate scenarios:

• Scenario 1: The insecticidal spray increases the mortality rate of midges that have successfully fed on cattle, which are assumed to remain nearby these hosts, and it also reduces biting from the susceptible midge population by increasing their per capita mortality. This implies that the critical excess mortality () is the solution to: (29)
• Scenario 2: The insecticidal spray increases the mortality rate of midges that have successfully fed on cattle, which remain nearby these hosts, but it does not reduce biting from the susceptible midge population which have yet to be affected by the insecticide. This implies that the critical excess mortality () solves the following equation: (30)

Where is the expected number of bites after its EIP by a midge undergoing excess mortality . It should be noted that Eqs (29) and (30) are defined in terms of an estimate of the reproductive ratio based on classical assumptions, which can be related to a generalised reproductive ratio via Eq (22). For either control scenario calculating an approximate critical excess mortality () by using the approximate expression for B_I Eq (21) is significantly easier numerically. The quality of the analytic approximation to the critical excess mortality has similar dependency on the dispersion parameters as the approximate reproductive ratio : the approximation is broadly good when d_E < 0.5 or when the gonotrophic cycle is more dispersed than the EIP (S1 Fig).

The classical formulation of the reproductive ratio Eq (1) suggests that increasing the mortality rate of vectors could be an effective control strategy since decreases faster than exponentially with increasing μ. In the generalised setting this is not necessarily true; the sensitivity of to μ can be sub-exponential, implying that increasing vector mortality is less effective at reducing the proportion of inoculated vectors expected to survive their EIP than the classical prediction would predict. On the other hand, in the generalised setting the mean vector population biting rate can decrease with increasing μ; an effect that reduces transmission and which is absent from the classical Ross-Macdonald model. In fact, the generalised critical value, , can be either greater or less than the value implied by the classical assumptions () derived from solving Eq (28), with the relative critical excess mortality () increasing for both increasing d_G and d_E (Fig 4). Unsurprisingly, the parameter region where control using spraying is expected to be harder than the classical prediction (i.e. is larger for scenario 2, which is more pessimistic about the efficacy of insecticidal spraying than scenario 1 (Fig 4). Because the critical vaccination coverage can be expressed as depending only on the reproductive ratio, the curve R_rel = 1 separates the region where it is more difficult from the region where it is less difficult to eradicate using host vaccination than the classical prediction. However, the R_rel = 1 curve is not identical to the curve defined by (Fig 4). This has an important consequence, even if (i.e. the critical host vaccination coverage is the same even when accounting for the distributions) then the critical amount of excess mortality that insecticide based control must achieve in order to eradicate may be different depending on the distributions. In fact, for culicoides biting midges the amount of necessary insecticide control predicted using classical assumptions is an underestimate everywhere along the R_rel = 1 level curve for scenario 1, and the underestimation is even more dramatic for scenario 2. In fact, whenever the generalised reproductive ratio and the classical reproductive ratio coincide, the classical model is always overly optimistic about the how little excess mortality will be required to achieve eradication. This is true for both varying (d_G, d_E) for a BTV epidemic with the classical estimate and when fixing d_E = 0.130 (the median posterior estimated dispersion for the EIP of C.bolitinos infected with BTV-1 serotype [19]) allowing to vary (Fig 4). However, it should be noted that there is an appreciable region of parameter space where eradication is easier to achieve than expected using classical modelling assumptions. If d_G is sufficiently small then in most cases we predict that it could be easier to eradicate via increasing adult vector mortality than predicted by classical modelling. Whether the findings of this paper are optimistic or pessimistic with regard to insecticidal spray efficacy will ultimately be resolved by more detailed data on biting vector gonotrophic cycle dispersion.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Critical control effort relative to classical predictions.

The relative excess mortality that would be required to control BTV using insecticides compared to its classical value given (top left). The red dot indicates the classical Ross-Macdonald prediction. The environmental temperature is assumed to be a constant 20°C giving the same parameters as in Fig 2. The critical excess mortality required to achieve eradication using insecticides as predicted by classical assumptions can be either an underestimate or an overestimate. Tracing the level curves from the classical value for critical vaccination coverage (; red curves) and critical excess mortality (; blue curves) separates the (d_G, d_E) plane into regions where control is more or less difficult compared to the classical value for both cattle vaccination and insecticide spraying. Two scenarios for the efficacy of spraying are considered: where the excess mortality reduces the population biting rate as well as reducing the life expectancy of inoculated vectors (solid blue curves) or if it solely decreases the life expectancy of inoculated vectors (dashed blue curves). For fixed (top right) and fixed d_E = 0.130 (based on C. bolitinos inoculated with BTV-1 serotype [19]; bottom) the level curve for R_rel = 1 is distinct from the level curve μ_rel = 1.

https://doi.org/10.1371/journal.pcbi.1004837.g004

Modelling mortality associated with vector biting and the first gonotrophic cycle

In the model presented here we assume that vector mortality is due to a constant background hazard rate μ and the vectors’ first gonotrophic cycle after emergence as an adult is similar to subsequent cycles. An alternative is to explicitly include heightened mortality risk associated with vector activities such as host-seeking or oviposition as has been done in modelling papers of mosquito-borne diseases (e.g. [31]). Introducing additional risk of vector mortality immediately before and after biting modifies the time-dependent biting rate α(•) by incorporating the chance that biting leads to vector death. α(•) is the solution to a modified renewal equation and therefore introducing this extra model component represents no extra difficulty compared to only having background mortality (see S1 Text).

Many midge species are autogenous: their females emerge as adults ready to produce and oviposit an initial egg batch without a bloodmeal [32]. From an midge ecology perspective this allows the midge population to persist even when hosts are frequently absent. Therefore, from an epidemiological perspective the waiting time until the first bite from a female autogenous midge on a host occurs at the end of the second gonotrophic cycle. This possibility can be included in our model since we base our analysis on renewal theory: renewal processes with an initial waiting period that is different from subsequent waiting periods (known as delayed renewal processes [15]) are theoretically standard (S1 Text).

We did not explore these biologically motivated modelling extensions further due to our focus on bluetongue diseases spread by Culicoides genus biting midges; current best estimates of midge mortality in the wild are based on estimates of the proportion of parous midges in the population [13]. This gives an estimate of the mortality rate based on the probability of surviving a gonotrophic cycle which cannot be disentangled into separate hazards. Moreover, no autogenous midge species has been implicated in the transmission of any arbovirus, including bluetongue [32].

Concluding remarks

The classical Ross-Macdonald theory of VBDs created a critically important scientific framework with a measure of generational transmission intensity (R₀) at its heart. However the specific practical insights derived from the Macdonald reproductive ratio (for example the implied efficacy of increasing the mortality rate of adult insect vectors in reducing transmission) are undermined by the growing realisation amongst epidemiologists of mosquito-borne diseases that different methods for estimating R₀ produce strongly differing results (cf. Smith et al for a discussion on the problems associated with validation of the mathematical theory of mosquito-borne pathogens [6]). The problems with theoretic cross-validation appear to have not been explicitly stated in the literature of midge-borne pathogens but the issues are essentially similar. In summary, it has become generally accepted in the modelling community that, whilst all models are simplifications of reality, the classical simplifications of Ross-Macdonald theory may produce misleading results (cf. recent reviews [5, 6]). The development of mathematical models of VBDs that deviate from Ross-Macdonald theory has generally focused on heterogeneity in transmission, whether due to heterogenous distribution of vectors across space, for example the recent northern European outbreak of BTV has inspired a number of explicitly spatial models (for example [33–35]), or because of the observation that (typically) a disproportionately greater number of bites will be distributed across a small number of vertebrates (as discussed in [26]). It seems that questions about the role of chance in generating variation in life histories between vectors, even if they are otherwise identical from a modelling perspective, has been largely overlooked despite being in principle empirically observable in the laboratory setting.

In this work, we have demonstrated that relaxing the classical assumptions on the random distributions underlying the calculation of a reproductive ratio for a VBD can have a significant impact on the prediction, even without additional modelling of heterogeneous biting or spatial vector distribution. The impact of the precise distribution of life history events is accentuated for short-lived vectors; thus the economically important vector-disease complex of BTV spread by, the typically short-lived, culicoides genus biting midges is an ideal case study for our theory. Our main results are presented in terms of moment generating functions of fairly arbitrary period distributions, so that any reasonable distribution can be used for modelling aspects of the vector life histories. However, our illustrative examples are calculated for the gamma distribution, a natural choice since it interpolates between the constant duration and the exponentially distributed duration. Gamma distributions have been used extensively to model the durations of latency (EIP) and viraemia, most commonly as a discrete sum of exponential distributions (a special case of the gamma distribution called the Erlang distribution), in both the vector-borne disease literature [9] and the directly transmitted disease literature [28]. However, the timing of epidemiological relevant contacts, (e.g. successful bites) are overwhelmingly modelled as a Poisson process throughout the mathematical modelling literature. We show that the generalisation of the contact model from a Poisson process to a renewal process, in conjunction with the typically short lives of the vector and a general EIP duration, has significant effects on the basic predictions of an otherwise standard VBD model.

The more general modelling of the vector’s biting process as a renewal process, and the EIP duration as non-constant, predicts that the classical reproductive ratio can be either a significant overestimate or underestimate compared to the generalised prediction. Moreover, the generalised modelling approach cannot be disentangled from considerations of the age distribution of the vector population. However, by assuming that the population is at equilibrium and making theoretically well-motivated approximations (as well as assuming constant mortality rates) we are still able to present the dependence of R₀ on the entomological situation in an explicit manner. For the motivating case of midges spreading BTV, the reproductive ratio is increasing with increasing variation in either the latency (EIP) duration or the time between bites (gonotrophic cycle). This can be understood by considering two consequences of the short life-span of the ‘typical’ midge: firstly, a significant proportion of midges can die whilst their expected biting rate is low (due to relaxing the Poisson biting process assumption) and secondly, that the ‘typical’ midge inoculated with BTV will not survive her EIP (which can depend sensitively on the EIP distribution). Only the ‘tail’ of the inoculated midge population actively contribute to BTV incidence amongst hosts, and therefore random variation between vector life history events becomes a crucial factor. In particular, the potential importance of the precise details about the distribution of the gonotrophic cycle duration seems to be under-recognised in both the theoretical and entomological literature. These findings echo results from at least one other model which generalises the vector biting process [36] and suggests that greater effort must be devoted to understanding and quantifying the detailed vector biology both at the individual and population levels in order to accurately parametrise these kinds of models.

The popularity of R₀ as a metric of transmission intensity is in a large part due to its uncomplicated relationship with disease eradication: if R₀ can be reduced below unity for a sufficient period of time then the disease must be eradicated. The generalised prediction of R₀ presented in this work therefore leads to a reassessment of the efficacy of standard control measures such as the mass vaccination of hosts or the insecticidal spraying against adult vectors. The predicted critical vaccination coverage amongst hosts can still be expressed as only depending on the estimate of the reproductive ratio in the standard fashion, however the critical excess mortality amongst vectors that must be achieved depends on both the estimate of the reproductive ratio and other entomological/epidemiological factors such as basic vector mortality rate and the incubation rate for the pathogen even when making the classical Ross-Macdonald assumptions about the VBD. Predicting the critical excess mortality in the more general setting is further complicated by the necessity to consider the effect that excess mortality has on the age distribution of the vector population. Trials of the ability of insecticidal spraying aimed at adult midges to reduce the local midge population have been pessimistic, leading to speculation that much more must be known about the resting places of adult midges in order to design better spraying protocols [29]. We would add to this discussion further speculation upon which part of the midge population is affected by the additional environmental hazard of a persistent insecticide. Does insecticidal spraying reduce biting from the susceptible mass of the vector population due to shifting their age profile (our scenario 1), or not (our scenario 2)? In either scenario we find that even when the reproductive value predicted by the generalised model agrees with its classical counterpart then the classical model still under-estimates how much excess mortality must be achieved in order to eradicate. None of the arguments above should be read as implying that insecticidal spraying should not be used as a control measure against a vector-borne disease in any circumstance. It might well be the case that a vaccine for hosts might be unavailable or prohibitively expensive whereas a suitable, and safe, insecticide might be readily available and cheap.

An intriguing possibility is to target the older vectors that make up the tail of the age profile using late life acting insecticides (LLAIs). This has been suggested for controlling malaria so that evolutionary pressure on the mosquito population to develop insecticidal resistance is sharply diminished [37]. Modelling the EIP as random rather than fixed always favours vector survival of EIP [eq 15]. This implies that more ‘young’ inoculated vectors will be actively infectious than expected by a model with a fixed EIP; a pessimistic finding with respect to the efficacy of LLAIs whether the LLAI is effective after a delay since exposure or simply preferentially affects older vectors. An advantage of the modelling approach taken in this work is that the action of the LLAI can be explicitly age-dependent or dependent upon real time since exposure. In Read et al [37] the action of the LLAI on the vector is modelled as occurring after a number of vector gonotrophic cycles. However, the age of a vector after a number of gonotrophic cycles is uncertain due to variation in the cycle duration. Therefore, experimental data on age-dependent (or time delayed) LLAI efficacy might be difficult to interpret in the context of a model which focuses only on the number gonotrophic cycles a vector survives. Our modelling approach doesn’t suffer from this problem and experimental data on LLAIs could be integrated into our model directly. This would certainly be an interesting direction for future research.

Mathematical models are now commonly used to underpin policy decisions concerning disease control in both developed and low- and middle-income countries. We have used BTV which affects livestock spread by Culicoides genus midges as a motivating example. However, the methodology could equally apply to Leishmania in humans spread by sandflies, which have similar entomological characteristics to midges. Similarly, while mosquitoes do not necessarily match all assumptions within our model, mosquito-borne diseases could be accommodated with only a slight increase in model complexity. Our work has highlighted the sensitivity of model predictions to entomological and epidemiological details. This points to further avenues of applied experimental or observational research: 1) the acquisition of more detailed knowledge concerning the life cycles of potential vectors, in particular going beyond the measurement of simple averages to include variability; 2) studies of the life expectancy and age profiles of vector populations in the wild; 3) a more comprehensive investigations of the response of vector populations to insecticide-based control measures. Developing highly informative policy-relevant prediction for the future control of vector-borne diseases is likely to therefore require a combination of state-of-the-art models, with meticulous quantitative studies of the insect vector.