Epidemiology and control of maedi-visna virus: Curing the flock

Maedi-visna (MV) is a complex lentiviral disease syndrome characterised by long immunological and clinical latencies and chronic progressive inflammatory pathology. Incurable at the individual level, it is widespread in most sheep-keeping countries, and is a cause of lost production and poor animal welfare. Culling seropositive animals is the main means of control, but it might be possible to manage virus transmission effectively if its epidemiology was better quantified. We derive a mathematical epidemiological model of the temporal distributions of seroconversion probabilities and estimate susceptibility, transmission rate and latencies in three serological datasets. We demonstrate the existence of epidemiological latency, which has not explicitly been recognised in the SRLV literaure. This time delay between infection and infectiousness apparently exceeds the delay between infection and seroconversion. Poor body condition was associated with more rapid seroconversion, but not with a higher probability of infection. We estimate transmission rates amongst housed sheep to be at about 1,000 times faster than when sheep were at grass, when transmission was negligible. Maternal transmission has only a small role in transmission, because lambs from infected ewes have a low probability of being infected directly by them, and only a small proportion of lambs need be retained to maintain flock size. Our results show that MV is overwhelmingly a disease of housing, where sheep are kept in close proximity. Prevalence of MV is likely to double each year from an initial low incidence in housed flocks penned in typically-sized groups of sheep (c. 50) for even a few days per year. Ewes kept entirely at grass are unlikely to experience transmission frequently enough for MV to persist, and pre-existing infection should die out as older ewes are replaced, thereby essentially curing the flock.


d dt
Pr(T > t) = −λ(t)Pr(T > t) The solution of this differential equation is where ξ is a dummy integration variable. The cumulative distribution function (CDF) of the susceptible ewe's infection time, T , is F T (t; t 0 ) = Pr(T ≤ t). By definition, Pr(T ≤ t) = 1 − Pr(T > t), therefore, the CDF of T is F T (t; t 0 ) = 1 − Pr(T > t) Let f T (t; t 0 ) be the probability density function of the ewe's infection time. It is given by the derivative of its cumulative distribution function, We assume a Gamma distributed seroconversion period. However, we also assume that the seroconversion rate in winter 1982, when ewes were in poor condition, may be different than at all other times when ewes were healthy. The CDF of a Gamma distributed random variable X is F X (x) = Pr(X ≤ x) = γ(a,bx) Γ(a) where γ(a, bx) = bx 0 ξ a−1 e −ξ dξ is the lower incomplete gamma function, Γ(a) is the gamma function, a is the shape of the Gamma distribution and b is the rate. To model a seroconversion rate, b(t), that varies in time t through the experiment, we replace the term bx in the CDF with the function This is the cumulative seroconversion rate between times t and t + x. When b(t) = b is constant we recover Z(x; t) = bx which is independent of t as expected. Thus is the CDF of the seroconversion period, X, of a ewe infected at time t. With a time varying seroconversion rate seroconversion period is not Gamma distributed. However, except in winter 1982, seroconversion rate is constant and so seroconversion period will be Gamma distributed. The PDF of X is given by Substituting in Equation 9 gives Applying the chain rule gives Substituting in Equation 8 gives from which we obtain the final form of the PDF Let f S (s; t 0 ) be the distribution of seroconversion time S, which is the sum of infection time and seroconversion period, i.e., S = T + X. This distribution is, therefore, the convolution of the distribution of infection time T , and the distribution of seroconversion period X If a ewe is known to have seroconverted sometime between times s − and s + , the probability of it doing so is the integral of its distribution of seroconversion time between these two times, i.e., If a susceptible ewe never seroconverted between introduction into the flock and its last serological test at time t r , its probability of doing so is It now remains to solve the above integrals and to specify the number of infectious ewes I(t).

Discretising time
It is highly desirable that the above integrals have analytical solutions to allow rapid (order of hours) sampling of the posterior distribution. However, given the time-varying transmission rate β(t), and the time-varying seroconversion rate b(t), the above integrals are not analytically solvable. In order to obtain analytical solutions we can discretise time into 1 month intervals (for Houwers) or 1 week intervals (for Peterson) which roughly match the sampling intervals. (To not interrupt the flow of our argument we only refer to monthly intervals, but the same arguments hold for weekly intervals.) This means we can treat the transmission and seroconversion rates as piecewise-constant such that in month i, transmission rate β i , and seroconversion rate b i , are constant. We also assume that the number of infectious ewes I i , in month i is constant. This implies that the force of infection λ i = β i I i , in month i is also constant. We first discretise the distribution of infection time f T (t; m 0 ). The probability of a susceptible ewe being uninfected at time t (Eqt. 2) can be partitioned into monthly intervals. We set the units of time to be months with 1 month equalling 365 /12 = 30.4 days. Let m 0 be the month a susceptible ewe was introduced into the flock and let u n (t; m 0 ) be the probability of it being uninfected at time t in month n (i.e., t ∈ [n, n + 1)). This probability is given by Define u m0,n = exp − n−1 i=m0 λ i , which is the probability of a ewe being uninfected before month n, and substituting into Eqt. 19 we obtain u n (t; m 0 ) = u m0,n e −λn [t−n] (20) Let f T,n (t; m 0 ) be the piecewise distribution of infection time T , in month n. This is obtained by substituting Eqt. 20 into Eqt. 7: f T,n (t; m 0 ) = λ n u n (t; m 0 ) = λ n u m0,n e −λn[t−n] Second, we discretise the distribution of seroconversion period f X (x; t). Let Z n,m (x; t) be the cumulative seroconversion rate between times t and t + x such that t is in month n (i.e., t ∈ [n, n + 1)) and t + x is in month m (i.e, t + x ∈ [m, m + 1), n ≤ m). It is given by which, with some rearranging, is Let f X,n,m (x; t) be the piecewise distribution of seroconversion period X, in month m given infection at time t in month n. This is obtained by substituting Eqt. 25 into Eqt. 14: Finally we discretise the the distribution of seroconversion time S. Let f S,m (s; m 0 ) be the piecewise distribution of seroconversion time S in month m (s ∈ [m, m + 1)). This is given by Eqt. 15 partitioning the integral into 1 month intervals we obtain For convenience we define substituting in Eqt. 25 and defining B n,m = m−1 Substituting Eqts. 22 and 26 into Eqt. 28 and moving constants outside the integrals we obtain and substitute into Eqt. 31 to obtain where the constants of integration have cancelled. Consider the indefinite integral from Eqt. 32 Applying a change of integration variable dz dt = −b n , using the fact from Eqt. 30 that t − n = − z bn +

Bn,m+bm[s−m] bn
and moving constant factors out of the integral, we obtain is the lower incomplete gamma star function. By using this form of the incomplete gamma function we remove the denominator b n − λ n from Eqt. 36 which will cause numerical instability when b n ≈ λ n . Substituting into Eqt. 36 gives where the integration constant has been removed as it cancels in Eqt. 33. Also note that v m,m (s, s) = 0 in Eqt. 33 because z m,m (s, s) = 0. Let P m0,m = m+1 m f S,m (s; m 0 ) ds be the probability that a sheep introduced in month m 0 seroconverts in month m. Substituting in Eqt. 33 gives Moving the integrals into the summation and moving constant factors out of the integrals gives We define the indefinite integral and substitute into Eqt. 40 to obtain where the constants of integration have cancelled each other. Note that w m,m (m, m) = 0. Substituting Eqt. 38 into Eqt. 41 we obtain Applying a change of integration variable dz and moving constant factors out of the integral, we obtain The solution of which is If a susceptible ewe, introduced in month m 0 , tested negative in month m − and then positive in month m + the probability of it doing so is If a susceptible ewe, introduced in month m 0 , never seroconverted until its last test in month m r , the probability of it doing so is The equations for p + and p − are found by applying Eqts. 30, 45 and 42. These equations are functions of the force of infection λ i for each month i which, in turn, depend on the number of infectious ewes I i in each month i. In the next section we describe how we find I i . Amparo et al. (2016) provides an algorithm to numerically solve the lower incomplete gamma star function.

Calculating the expected number of infectious ewes
In Peterson the number of infectious ewes equals the number of donor ewes which are assumed to be infectious from the start of the experiment until they are removed from their flocks in weeks 50 or 51. There are 10 donors in the Texels and 9 in the Blessumers (see Supplementary spreadsheet Peterson data).
In Houwers the situation is more complicated because infected recipient ewes become infectious after a latent period. We use a discrete-time, infectious disease model to calculate the expected number of infectious ewes I i , and the force of infection λ i , in all months i = 0, . . . , M where M is the last test month.
We do not know the latent period L of an infected ewe. We therefore treat L as a Gamma distributed random variable with mean l and standard deviation σ. Let F L (y) = Pr(L ≤ y) = P (y; l 2 /σ 2 , σ 2 /l) be the cumulative distribution function of latent period L, where P (y; a, k) is the lower incomplete gamma function with shape a and scale k.
The contribution I i,n+j , a susceptible ewe i, introduced into the flock in month m 0,i , infected in month n (n = m 0,i , . . . , m r,i − 2) and removed in month m r,i , makes to the expected number of infectious ewes in Parameter Prior β field Exp(0.02) β housed Gamma(2, 0.1) β unhealthy Gamma(2, 1) τ Uniform(1, 4) l Uniform(0, 60) σ Gamma(2, 1.5) α Uniform(0.001, 100) µ healthy Exponential(24) µ unhealthy Exponential(24) month n + j (j = 1, . . . , m r,i − 1 − n) is the product of the probability the ewe is infected at time T in month n and the probability the ewe's latent period is less than or equal to j − 1 months, i.e., Substituting in Eqt. 20 gives where λ n = β n I n is the force of infection in month n and I n = i∈ewes I i,n is the expected number of infectious ewes in month n. Note that infectious ewes are assumed to have been removed in the last month they are tested and so do not contribute to the force of infection in that month.
We assume that donor ewe with ID 72 is infectious from the start of the experiment and so contributes 1 to the expected number of infectious ewes for each month until it is removed in January 1981. We test whether the other donor (ID 22) is infectious or latent at the start of the experiment.

Parameter prior distributions
When fitting to Houwers' data, we assign parameter prior distributions that are weakly informative reflecting the lack of information about the epidemiological parameters of SRLV transmission (Table 1). We use the mean parameter estimates from the fit of Houwers' data to assign the priors for fitting Peterson's data.
There are no published estimates of field transmission rates β field . A weakly informative prior on field transmission rate can be found by assuming that all of the 23 of the 42 ewes that were infected over the 60 months of the experiment were infected by a single ewe whilst grazing. This is an overestimate because it neglects housed transmission and assumes a single infectious ewe was responsible for all infections. Assuming this gives a rate of − 1 /60 ln(1 − 23 /42) ≈ 0.01. We take a slightly larger average rate of 0.02 to make the prior less informative. We assign an exponentially distributed prior-which has its mode at zero-as it is possible that field transmission might be close to, or even, zero.
There are numerous studies that suggest that housing is an important factor in SRLV transmission, but there are no published estimates of its value. We therefore assume that housing transmission rate is about 10 times the field transmission rate of 0.01. We therefore assign a weakly informative, Gamma distributed prior with shape 2 and scale 0.1 to β housed , giving a mean transmission rate of 0.2 and zero probability of a zero transmission rate.
There are no published estimates of transmission rate in poor condition ewes, other than clinical observation that poor husbandry conditions and the presence of secondary infections are conducive to more rapid spread of MV (Prichard and McConnell, 2007). We therefore assign a weakly informative, Gamma distributed prior with shape 2 and scale 1 to β unhealthy , giving a mean transmission rate of 2.
Ewes were housed for at least one month during lambing in March and had access to housing for another three months from December to February. We therefore assign a flat uniform prior distribution from 1 to 4 to the duration of housing τ .
There are no published estimates of mean latent period l, or its variability σ. We therefore assign its mean l, a flat uniform prior between 1 and 60 months and its standard deviation σ, a weakly informative, Gamma prior with shape 2 (non-zero mode) and scale 15 months.
There is ample evidence that ewes seroconvert from a few weeks to several months after infection (e.g., Rimstad et al. (1993)). As discussed above, we have assumed that seroconversion period is Gamma distributed with shape parameter α and mean µ healthy for healthy ewes and mean µ unhealthy for poor condition ewes. As there is no prior information on the shape of this distribution we assign a weakly informative prior of Uniform(0.001, 100) to α. We assign exponentially distributed priors to µ healthy and µ unhealthy with a scale of 24 months.
The estimates from Houwers are used to assign the scale of the priors for Peterson (Table 2). In particular, the mean seroconversion period µ, is assigned a Gamma prior distribution with a mean of 8 months to reflect our inference about it from Houwers' experiment.
Parameter Prior Gamma(2, 4) 2 Calculating R 0 and the prevalence doubling time Let A be a next-generation matrix whose elements a i,j , are the expected number of ewes of age class j that are infected over the infectious lifespan of a ewe infected at age i Diekmann et al. (1990); van den Driessche and Watmough (2002). Let N g be the vector of the number of infected ewes in each age class on generation g. Then, when MV prevalence is low, the number of newly infected ewes in the next generation is The solution of this equation is N g = R g 0 v where R 0 is the dominant eigenvalue of A and v is its corresponding eigenvector whose elements are the stable proportions of infected ewes in each age class. The time between generations G, is the weighted mean ewe lifespan post infection. The weights are the stable proportions of ewes (infected plus susceptible) in each age class, and the average is taken over all age classes.
When MV prevalence is low, the growth equation of the total number of infected individuals is |N g | = R g 0 |v|. The time since initial infection is t = gG. Substituting g = t /G into this equation, taking logs and rearranging gives By setting |N g | = 2|v| and substituting into Eqt. 53 we obtain the prevalence doubling time G ln 2 ln R0 .

Markov chain convergence and mixing
In all models, three Markov chains were run for a total of 5 × 10 5 samples thinned to 5,000 samples with a burn-in of 5×10 5 samples. Markov chain convergence and mixing was assessed by examining (i) the marginal distributions of each parameter of three parallel Markov chains (ii) the autocorrelation of each chain for each parameter and (iii) the Gelman-Rubin statistics of each parameter Gelman et al. (2003). Assessments are shown for the optimal model including field transmission for Houwers (Fig. 1), Peterson's Texels (Fig. 2), Peterson's Blessumers with (Fig. 3) and without (Fig. 4)  : Houwers. Assessment of Markov chain convergence, mixing and autocorrelation for the optimal model including β field for Houwers. The lefthand panels show the trace plots of three chains as green, blue and orange dots for each parameter. Their overlap demonstrates the chains' excellent convergence, stationarity and mixing. The similar chain histograms shown as green, blue and orange bars in the second column of panels demonstrate excellent convergence between the three chains. The overlapping histograms of the first (solid lines) and second (dashed lines) halves of each chain in the third column of panels demonstrates the chains' stationarity.The autocorrelation plots in the righthand panels demonstrate the lack of autocorrelation in the three chains. Gelman Rubin statistics (Rc, Gelman and Rubin (1992)) for each parameter is given in the righthand panels and demonstrate excellent convergence. Prior distributions are represented as red lines in the second column of panels.