The Interaction between Vector Life History and Short Vector Life in Vector-Borne Disease Transmission and Control

Epidemiological modelling has a vital role to play in policy planning and prediction for the control of vectors, and hence the subsequent control of vector-borne diseases. To decide between competing policies requires models that can generate accurate predictions, which in turn requires accurate knowledge of vector natural histories. Here we highlight the importance of the distribution of times between life-history events, using short-lived midge species as an example. In particular we focus on the distribution of the extrinsic incubation period (EIP) which determines the time between infection and becoming infectious, and the distribution of the length of the gonotrophic cycle which determines the time between successful bites. We show how different assumptions for these periods can radically change the basic reproductive ratio (R0) of an infection and additionally the impact of vector control on the infection. These findings highlight the need for detailed entomological data, based on laboratory experiments and field data, to correctly construct the next-generation of policy-informing models.

to modelling vector life histories in the main text including references to equations in the main text. Wherever possible we keep the same notation as the main text even when discussing the general theory rather than using 'standard' notation for renewal process theory (e.g. as used in Grimmett and Stirzaker [15; main text]). Our aim is to immediately make connection between the general theoretical ideas and the specific application on vector modelling.

Theory: Laplace transforms, moment generating functions and renewal processes
We describe the underlying theory used in this paper. The results on renewal processes and moment generating functions can be referenced in many standard texts on probability theory (e.g. [15; main text]), the results on Laplace transforms are also standard (e.g. [1]).

Laplace transforms and moment generating functions
We will require some results for Laplace transforms of functions. We define the Laplace transform of a function g:ĝ and the Laplace convolution of two functions g, h, (g * h)(t) = t 0 g(s)h(a − s) ds. (2) We use three standard results for Laplace transforms (e.g. [1]): g * h =ĝĥ (Convolution to product), (4) d dt g(t) = θĝ(θ) − g(0) (Differential rule).
In this work we will only need to consider positive random variables X > 0 which we use to model uncertain durations of time. Each random variable has a distribution specified by their distribution function F X (x) = P(X ≤ x). If F X is differentiable then there exists a density function for X which is is the derivative of the distribution function f X = F X . The density function can be thought of as the infinitesimal probability density around values the random variable can take, f X (x) dt ≈ P(X ∈ [x, x + dx]). For each random variable X we can define a moment generating function (MGF) φ X (θ), The MGF of many standard random variables can be calculated analytically and the various moments of X can be calculated from the derivatives of φ X . If the positive random variable has a density function f X we note that the MGF of X is the Laplace transform of its density function, (7)

Gamma distributed durations and survival probabilities.
In this work we develop theory that applies to any positive duration X > 0 however for concreteness we apply the theory to gamma distributions. We parameterise the gamma distribution in a slightly non-standard way, fixing the mean E[X] and dispersion d X = VAR[X]/(E[X] 2 ), giving a gamma density function, Where the normalisation constant Γ is the gamma function. The parametrisation we use for the gamma distribution relates to standard parametrisations: the gamma shape parameter is 1/d X and the gamma scale parameter is E[X]d X .
In this work we fix mean durations to known entomological parameters, but allow the dispersion value to vary since this quantity is generally under-reported or uncertain in the literature. We consider a range of dispersions d X ∈ [0, 1], this scales between the exponentially distributed duration with mean E[X] (d X = 1) and the fixed duration of length E[X] (d X = 0). To see this: 1) which is the density function of an exponential random variable with mean E[X]. 2) The MGF for the gamma distribution is: The limit as d X → 0 can be found by relating (9) to a standard result, We note that in the limit d X → 0 the MGF of the gamma distribution duration X becomes the MGF of the constant 'random' duration X = E[X] i.e. X is fixed at its average with no variation.
In the main text we demonstrate that the probability of a vector surviving for some random duration X > 0 whilst undergoing constant mortality at rate µ is the MGF of X evaluated at the vector mortality rate [equation (14; main text)], P(Vector surviving for random duration X) = φ X (µ).

(11)
This is a general result which connects MGFs to survival probabilities.

Renewal processes
A renewal process B = (B(t), t ≥ 0) is a stochastic model for representing the number of randomly occurring events that have arrived within each time period ([0, t], t ≥ 0). Events occur instantaneously with a definite arrival time, we denote the arrival time for the n th event T n . The waiting time between time 0 and the first event arrival is denoted G 1 and subsequent waiting durations between the (n − 1) th and n th event arrivals are denoted G n . Therefore, and the renewal process is defined at each time in terms of the waiting durations, The key underlying assumption that differentiates renewal processes from other event counting processes is that the set of waiting durations {G n } n≥1 are treated as positive identically and independently (i.i.d.) random variables. Assuming i.i.d. waiting durations implies three properties of renewal processes: 1. Each waiting duration G n has the same distribution function as a common random variable G with distribution function F G (t) = P(G ≤ t).
2. Positive waiting durations implies that F G (0) = 0, f G (t) = 0 for t < 0 and the density function of Where f * n G is the n-fold Laplace convolution of the density f G with itself (a simple corollary of theorem 4.8.1 [15; main text]). The distribution function F Tn (t) = P(T n ≤ t) can be found from integration, 3. The independence of the waiting durations {G n } n≥1 implies that if we condition on an event arriving at some time T then the duture event arrivals are probabilistically equivalent to a renewal process started at time T . The process is said to renew at each event arrival.

The renewal function.
The expected number of events that have arrived by time t is called the renewal function, The elementary renewal equation [15; main text, Eq. 3 Section 10.2] gives that m(t)/t → 1/E[G] as t → ∞ but we will need finite time results. We will use two different, but equivalent, representations for the renewal function: • • The two solutions are identical representations but will be useful in analysing biting from two different sections of the vector population: 1) biting from the susceptible population of vectors and 2) biting from the inoculated vectors.

Remaining time until next event.
At any time t we can define the remaining time until next event as a random variable: There is an important limit theorem for the distribution function of remaining time until next event [15; main text, Eq. 5 section 10.3]: , as t → ∞. (20) The limit theorem (20) states that if we pick t a long time after the start of a renewal process (t 0) then the random time remaining until the next event after t will be approximately distributed according to the random variable R(∞). This approximation becomes exact as t → ∞.

Delayed renewal processes.
A common extension to the renewal process model is to assume that the first waiting period has a different distribution to the other waiting durations; that is G 1 is independent of, but not identically distributed to, (G n ) n≥2 which share the common distribution function F G . Traditionally, this modelling approach has been used when the first waiting time is expected to be on average longer than subsequent waiting durations and so the extended model is called a delayed renewal process. The renewal function for the delayed renewal process is denoted m d and has analogous solutions to a standard renewal process : • m d can be solved as an infinite sum [15; main text, Eq.14 section 10.4], Where the delayed distribution functions for the arrival times are: • The delayed renewal function can also be represented as the implicit solution to the integral equation (23)

Application: Modelling vector life histories
In this section we apply the standard theory presented above to modelling the epidemiologically relevant events in the life of a disease vector: the frequency and timing of biting, the duration of the EIP for inoculated vectors and the life duration. We use the model to calculate two key quantities for R 0 prediction: • The expected rate of bites per host from the susceptible vector population per host (M E[α(A)]).
• The average number of potentially infectious bites from inoculated vectors on the host population (B I ).
Other variables that effect R 0 , namely vector competence, vector-to-host population ratio and mean viraemic duration for hosts have a direct multiplicative effect (Eq. 3; Main text). We make three specific modelling decisions for the vector life histories: • The bites of a vector after its emergence as an adult are assumed to occur at the arrival times of a renewal process B(t), with the renewal waiting durations being the gonotrophic cycle lengths (G n ) n≥0 .
Each gonotrophic cycle duration is modelled as an independent gamma distribution with known mean duration E[G] = 1/α, i.e. the inverse of the asymptotic biting rate, and dispersion d G . The density function f G is given by equation (8).
• If a vector is inoculated then only bites after an EIP duration E are potentially infectious. Variation between vectors in EIP duration is modelled by treating E > 0 as a gamma distributed random variable with known mean duration E[E] = 1/σ, i.e. the inverse of the pathogen development rate, and dispersion d E , identically distributed for each vector with density function f E given by equation (8).
• Each vector has a density-and disease-independent mortality rate µ which is equivalent to an exponentially distributed lifespan L ∼ exp(µ).

Vector survival probabilities
In order to bite each vector must survive its gonotrophic cycle and potentially infectious bites can only occur after an inoculated vector survives its EIP. The survival probabilities implied by the model can be calculated directly from equations (9) and (11): P(Vector survives a gonotrophic cycle) : P(Vector survives its EIP) : This derives P E [equation (26; main text)] used in the main text.

The biting renewal process
We now apply the results from section 1.2 for both biting from the susceptible vector population per host and for biting from inoculated vectors.

Biting from susceptible vectors.
Susceptible adults vectors may or may not have bitten, the only knowledge we have of their biting pattern is that they emerged at some time. Therefore the relevant time for the susceptible vectors biting process is their age a since emergence. The expected rate of biting for a vector at age a is the time derivative of the renewal function (16) which we denote, We use solution 1 for the renewal equation (17) Equation (27) justifies (equation 18; main text). The expected biting rate from the entire vector population will depend on the age distribution of population; that is the age of a vector chosen at random from the adult vector population. We denote this randomly chosen age A with density f A . At equilibrium the age distribution is f A (a) = µ exp(−µa) (equation 13; main text). For vector to host ratio M the expected biting rate from the susceptible population per host is then: Using the argument in equation (19; main text). Applying our specific result for gamma distributed waiting times from section 2.1 to (28) we get: This derives the biting rate per host used in the main text [equation (25; main text)].

Biting from inoculated vectors.
Vectors are inoculated by successful transmission during a bite on an infectious host. Therefore, inoculated vectors must have survived for long enough to make at least one successful bite. The biting process renews at this time (property 3 of renewal processes). Exponentially distributed lifespans have the memoryless property, the future life expectancy for the vector after the inoculating bite is still exponentially distributed exp(µ). Therefore the relevant time for the future biting of the vector is time since the inoculating bite t, the age of the vector at the inoculating bite is irrelevant to the future chance of biting. We use a numerical technique to solve for the the expected biting rate based on the Laplace transform of solution 2 of the biting renewal equation,m The expected biting rate at time t after the inoculating bite is defined as α(t) = d/ dtm(t) and m(0) = 0.
We use the Laplace differential rule [Equation (5)] to construct the Laplace solution for the expected biting rate for inoculated vectors,α We solve for the age varying biting rate α(a) along a pre-defined mesh of a values by using the Talbot numerical Laplace inversion algorithm 1 [2]. The random number of infectious bites over the remaining lifetime of the vector are those which occur after a random latency period E but before the death of the vector L. This leads directly to a stochastic integral expression,

Choose a mesh of t values of size
2. Use the Talbot Laplace numerical inversion algorithm on equation (31) to solve: 3. Use a quadrature rule on the mesh points to numerically evaluate (33).

Approximating the number of infectious bites per inoculated vector
Equation (33) gives an exact but implicit result for the expected number of infectious bites made by an inoculated vector. However, it is often convenient to give an explicit, but approximate, result. Each vector has an independent probability φ G (µ) of surviving each gonotrophic cycle to successfully bite again, which implies that the lifetime number of bites by the vector is geometric distributed with mean φ G (µ)/(1 − φ G (µ)).
Because of the renewal property (property 3 of renewal processes), and because the vector life duration is expected to be memoryless, inoculated vectors that survive 1) their EIP E and then 2) the time remaining before its next bite R(E) will make their first bite and then are, probabilistically, equivalent to a vector starting its adult life infectious. Therefore, we can write B I (33) in probabilistic form, Where φ E+R(E) (µ) is the probability that the vector survives for the duration E + R(E) with constant mortality rate µ [equation (11)]. The random duration E +R(E) is not analytically accessible in general. However if the EIP duration is typically long compared to the convergence rate of expression (20) then we can approximate And we can make an approximation (which limits to being exact) for the probability of the vector surviving to make its first infectious bite, P(Vector survives EIP and remaining gonotrophic cycle We use equation (7) and invert the Laplace differential rule [i.e.F X (θ) =f X (θ)/θ], in order to calculate the probability of a vector surviving the limiting random remaining gonotrophic cycle duration, This leads directly to the approximation for B I used in the main text, This justifies equation (21; main text). The quality of the approximation for B I depends on how rapidly the remaining time until next bite distribution converges onto its asymptotic distribution, hence the approximation is exact when for exponentially distributed gonotrophic cycles (this can be seen by inspection of (35)). On the other hand, if the EIP duration is more dispersed than the gonotrophic cycle duration, which in turn is significantly under-dispersed compared to the exponential distribution, then the fundamental assumption that the EIP duration is long compared to the convergence rate is poor. This explains why the approximation to B I is found to be reasonable except in the parameter regime {d G < 0.5, d E > 0.5}.

Fixed inter-bite periods
The only example in the text where the above analysis is inappropriate is when the dispersion of the gonotrophic cycle is exactly zero, d G = 0, which corresponds to the time between successful bites being some fixed value T so that G n = T for all n ≥ 0. This transforms the solution to (33) into Truncating the sum (39) at some value n * gives an error which is less than exp{−n * µT }/(1 − exp{−µT }), the error can be chosen arbitrarily small by an a priori choice of n * .

Model extensions
We discuss two modelling extensions that can be easily added to the general modelling framework presented in this work. Firstly, that the first gonotrophic cycle might be intrinsically different from subsequent gonotrophic cycles and secondly, that as well as a constant mortality rate for vectors there might also be a specific risk associated with biting (both before and after an attempt).

Modified first gonotrophic cycles
It might be more realistic to model the first gonotrophic cycle separately to subsequent gonotrophic cycles. Some midge species are autogenous; that is they lay one batch of eggs before they seek a host. Essentially autogenous midge species go through two gonotrophic cycles before their first bite. The theory of delayed renewal processes (section 1.2.3) is ideal for analysing the consequences of this model modification. We assume that we can gather data on the distribution of the first gonotrophic cycle thus giving an estimate for its density function f G 1 . The modified first gonotrophic cycle alters the biting rate from the susceptible population because the probability of vectors surviving to their first bite is altered.
We repeat the analysis of section 2.2.1 using a solution to the delayed renewal equation [equation (21)]. This gives the delayed age-dependent biting rate as, The modified biting rate from the population is therefore, Inoculated vector have, by definition, succeeded in biting at least once so their future biting should still be modelled using a standard renewal process.

Vector mortality risk associated with biting
Constant mortality is an appropriate model for background vector risk such as from predation. However, biting brings its own specific additional risks which we have not modelled. For example, we could model the effect of insecticide deployment near hosts as introducing a risk of vector death both immediately before and immediately after biting. We define p 1 as the probability of the vector surviving the risk before biting and p 2 as the probability of the vector surviving the risk after biting. We denote the modified renewal equation for inoculated vectors m r . The probability that the vector survives to make its next bite after inoculation p 1 p 2 , we condition on this happening and the timing of the next biting to give, The renewal equation (42) can then be used solve for a modified B I as in section 2.2.2. If we assume that it is only the subsection of vectors that bit on the protected hosts that undergo this extra risk (as in control scenario 2 in the main text) then the biting from the susceptible population is unaffected. However, this may not be the case. Evaluating the effect of biting risk on the age-distribution of the population is possible and will be a future direction of research. MGF of random variable X φ X (θ) φ X (θ) = E[exp{−θX}] Distribution function of random variable X F X (t) F X (t) = P(X ≤ t) Distribution function of random variable X f X (t) f X (t) = ( d/ dt)F X (t) Probability vector survives random period X P X P X = φ X (µ) with mortality rate µ Biting renewal process for a vector B(t) # bites in period [0, t] n th gonotrophic cycle for a vector G n (G n ) n≥1 are i.i.d random variables G is gamma distributed with mean 1/α and dispersion d G .
Asymptotic biting rate