A genetic epidemiological model to describe resistance to an endemic bacterial disease in livestock: application to footrot in sheep

Selection for resistance to an infectious disease not only improves resistance of animals, but also has the potential to reduce the pathogen challenge to contemporaries, especially when the population under selection is the main reservoir of pathogens. A model was developed to describe the epidemiological cycle that animals in affected populations typically go through; viz. susceptible, latently infected, diseased and infectious, recovered and reverting back to susceptible through loss of immunity, and the rates at which animals move from one state to the next, along with effects on the pathogen population. The equilibrium prevalence was estimated as a function of these rates. The likely response to selection for increased resistance was predicted using a quantitative genetic threshold model and also by using epidemiological models with and without reduced pathogen burden. Models were standardised to achieve the same genetic response to one round of selection. The model was then applied to footrot in sheep. The only epidemiological parameters with major impacts for prediction of genetic progress were the rate at which animals recover from infection and the notional reproductive rate of the pathogen. There are few published estimates for these parameters, but plausible values for the rate of recovery would result in a response to selection, in terms of changes in the observed prevalence, double that predicted by purely genetic models in the medium term (e.g. 2–5 generations).


Introduction
Preventive measures and lost production due to endemic disease form an important component of the costs of production in many livestock production systems [1], and they also affect animal welfare and marketability of breeding stock. It is well known that, for many diseases, resistance has a genetic component [2] and selection for disease resistance has long been considered a promising way to reduce disease prevalence e.g. [3].
Selection for resistance to an infectious disease has the added benefit that it may reduce the pathogen burden, especially when the population under selection is the main reservoir of pathogens. This will lead to an addi-tional reduction in prevalence, in addition to the direct genetic effect, as a result of reduced contamination from infectious animals, e.g. [4].
The phenotype used in selection for disease resistance is often a score, which includes a class of healthy animals, and one or more classes of affected animals. In the case of endemic diseases the vast majority of animals at any one time may be classified as healthy, and this limits the opportunity for intense selection. In a threshold model, that is appropriate for this type of data, prevalences that are much lower than 50% also lead to low heritabilities on the observed scale. A successful selection programme can therefore be expected to decrease the subsequent response to selection through increased resistance and decreased pathogen burden.
Anderson and May [5] describe the spread of a microparasitic (viral or bacterial) infection through a population of animals using a so-called SIR model, based on the rates at which susceptible (S) animals are infected (I) and then recover or are removed (R). A key parameter is R 0 , which is the number of secondary infections caused by the first infected animal. One or more of these rates can be under genetic control and hence affect R 0 . This model can be extended in various ways; for instance Bishop and Mac-Kenzie [6] have described how a disease that is spread from animal to animal may or may not lead to an epidemic, and Nath et al. [7] have explored the consequences of selection to alter different model parameters. To make these models more applicable to typical livestock bacterial infections, Bishop et al. [8] have considered a disease in which the pathogen survives for some time in the environment (E) from where it can infect susceptible animals. This was termed a SEIR model.
In the above models, the assumption is that recovered animals are no longer susceptible to the disease, and in a closed population without re-infection these models will therefore always lead to zero prevalence, either because there are no susceptible animals left or because the disease has died out. This outcome is inappropriate for typical endemic diseases, which often have a more or less stable prevalence over time. This stable prevalence, the endemic equilibrium, may be due to recovered animals losing their resistance and becoming susceptible again, or it may simply be a consequence of a continued introduction of new susceptible animals, e.g. offspring.
Building on the SEIR model [8], we introduce two new aspects to the model: (i) a period of latency (L) in which animals are infected but not yet infectious and (ii) loss of immunity so that recovered (R) animals can revert back to susceptible. This creates a SELIRS model. Further, we consider the SELDCRS model in which animals can be dis-eased and infectious (D) or a carrier, which is infectious but no longer clinically diseased (C). This model can be used to distinguish direct effects related to the number of diseased animals from indirect effects related to pathogen burden, by manipulating relative rates associated with D and C.
Footrot is an infectious disease of sheep caused by bacteria (Dichelobacter nodosus) that survive in soil for a limited time. The prevalence in adult sheep in Britain is around 6% [9]. Resistance to footrot has been shown to be heritable [10][11][12], and selection for increased resistance is feasible.
The aim of this study was to develop SELIRS and SELD-CRS epidemic models, in which pathogens survive in the environment for a limited time. The models were applied to footrot in sheep and used to predict the changes in prevalence of footrot over time if selection is for resistance to the disease, accounting for the disease dynamics. The predicted progress in terms of reduction in prevalence was compared with a model that ignores epidemiological effects.

Definition of epidemic models
Case 1: Infectious and diseased animals are equivalent An overview of all symbols and abbreviations used and their definitions is given in Table 1. Consider a population of N individuals which, at time j, consists of S susceptible, L latently infected, I infected and R recovered animals. It is assumed that it is only category I animals that show clinical signs of disease and are infectious. Environmental contamination is quantified by the concept of an infectious dose; therefore at time j there are E infectious doses of the pathogen in the environment. Following Bishop et al., [8], and approximating the discrete process (e.g. daily steps) by a continuous one, the SELIRS model is defined by the following five differential equations: ( 1 ) dL/dt = ES-L.
where  is the rate at which susceptible animals become infected (latent) per infectious dose,  is the rate of loss of immunity,  is the rate at which latent animals develop clinical signs and become infectious,  is the rate at which infectious doses are shed by infected animals,  is the rate at which infectious doses die (other than by host animals),  is the rate at which a host animal physically removes infectious doses (e.g. by ingestion, adherence to the animal or squashing them) and  is the rate at which infectious animals become immune. All rates are nonnegative.
Two properties of this model are of importance, the notional reproductive rate (R') and the equilibrium prevalence. In standard epidemiological models the basic reproductive ratio, R 0 , is the number of infections caused by a single infected animal in a wholly susceptible population, during the course of its infectious period. The equivalent in the SELIRS model is the number of secondary infections due to a single infectious animal, for the time period over which the infectious material remains in the environment. Bishop et al. [8] have derived an expression for the notional reproductive rate in a SEIR model as: where  =  + N, i.e. the total rate at which the infection is removed from the environment. The same expression may be used as an approximation to R' in an SELIRS model (see Appendix 1), however it is not exact as the loss of immunity potentially increases the number of secondary infections in situations where the environmental contamination is long-lived and the period of immunity is short.
An equilibrium state will be reached when the number of animals in each state is the same from one day to the next. In other words, the rates of change for the numbers of each category of animal (equations 1, 3, 4, 5) are all zero, i.e. dS/dt = dL/dt = dI/dt = dR/dt = 0. It can be shown (see Appendix 1) that the corresponding equilibrium number of diseased and infectious animals (I*) is: Combined with (6), it follows that the equilibrium prevalence (p*), i.e. I*/N, is Hence, after rearrangement, R' may be defined as a function of the equilibrium prevalence as follows: Note that in equation (7   The prevalence of the disease as observed from clinical signs R' The notional reproductive rate of the infectious disease the number of units required for an animal to become infected. This means that there is no need for an exact definition of the infective dose.

Case 2: Infectious animals with and without clinical signs
In many instances animals may be infectious, even when no clinical signs of disease are apparent, a phenomenon in some circumstances referred to as 'carrier status'. Defining clinical disease and infectious status as two separate but overlapping categories also allows greater flexibility in the exploration of the disease dynamics. This is achieved in the SELDCRS model, in which animals can be diseased and infectious (D) or infectious carriers who no longer show clinical signs (C), with N = S+L+D+C+R. Equations (1) and (3) describing change in S and L remain the same, with the remaining equations being: dD/dt = L-D. (10) dC/dt = D-C. (11) dR/dt = C-R.
Where the new symbols are:  is the rate at which diseased animals no longer show clinical signs and move to the carrier state, and  is the rate at which carrier state animals cease to be infectious. The total time an animal is infectious is 1/+1/, rather than simply 1/ as in the SELIRS model. This accounts for the supposition that infected animals may stop showing clinical signs yet continue to be infectious, rather than assuming that only animals with clinical signs are infectious. For the SELIRS and SELDCRS models to be equivalent in terms of transmission of infection then the total infectious period must be the same in both models, i.e. 1/ = (1/+1/). Note that this is proposed solely as a theoretical tool to distinguish between the effects of recovery on the animal it self (clinical signs) and on the population (infectious), which are confounded in the SELDIRS model. It is not proposed as a selection strategy. Table 2 shows a contrast of similar parameters in the SELIRS and SELDCRS models.
For purposes of comparing the SELDCRS model with quantitative genetic models that ignore the transmission of infection (see below), the pathogen burden (E) in the population can be made independent of changes in the rate at which diseased animals apparently recover () if any reductions in the number of diseased animals (D) are compensated by an increase in the number of carriers (C), i.e. the total time that an animal is infectious (1/+1/) is kept constant. Properties of the SELDCRS model are derived in Appendix 1, with R', the reproductive rate of the disease, being defined as: and the equilibrium prevalence p*, being defined as: Equation (14) can alternatively be written as: There are two differences between the equations defining equilibrium prevalence in the SELIRS and SELDCRS models (equations 8 and 15). First, if the assumption is invoked that the total infectious period is kept constant, then with decreasing 1/ the term R' is constant in SELD-CRS but will decrease with decreasing 1/ in SELIRS. Secondly, under the same assumption of constant infectious period, inspection of equations 8 and 15 reveals that the SELDCRS model will lead to the lower equilibrium prevalence, given the same values for all other parameters. The equilibrium prevalence will be relatively insensitive to change in R' if it has a high value (i.e. not close to 1).

Predicting responses to selection
Improvement of resistance to disease by genetic selection in its simplest form, i.e. disregarding information from relatives, consists of selection of those animals that do not show clinical signs at the time point at which the observations are made. The expected response to selection on such a binary trait depends on the heritability of resistance and the prevalence. This can be calculated assuming a threshold ('all or none') model with an underlying normally distributed liability with heritability h L 2 , which depends on the heritability of the binary trait (h 01 2 ) as: where p is the prevalence of the binary trait and z is the ordinate of the standardised normal distribution corresponding to p [13].
For the case where the number of healthy animals available for breeding exceeds the number required for selection, and again disregarding information from relatives, the response to selection is calculated as if the selected proportion is equal to 1-p, and on the underlying scale the response is R L = ih L 2 with i the selection intensity corresponding to 1-p. Because of reductions in prevalence, the response to selection is not linear but decreases with increasing values of 1-p.
In the context of the disease dynamics, selection for 'healthy' animals may be thought of as selection on one or more of the following: a. The recovery rate  or  (i.e. infectious or diseased animals recover quicker, e.g. due to an acquired immune response) b. The susceptibility of animals  (i.e. susceptible animals are more resistant to infection per se, requiring exposure to a greater number of infectious doses before becoming infected) c. The rate of shedding of bacteria  (i.e. infectious animals spread fewer bacteria) -this is comparable to selection for reduced egg counts with nematode infections d. The rate at which immunity is lost  (i.e. longer lasting immunity) e. The rate at which latently infected animals become infectious  (i.e. longer latency, as seen in scrapie).
Note that  and  occur only as products in R' and p*, therefore their effects can be considered jointly in this model.
In the simplest analogy to the threshold model, one can think of the model parameters, particularly  and  (i.e. recovery rate) or  (susceptibility) as being analogous to the liability, with heritable between-animal variation. This being the case, we can then assume that the liability and one or more of these parameters have the same distribution and heritability, and responses to selection can be calibrated between the threshold model and the SELIRS or SELDCRS models. Note that because of (6) selection on  (or ) is equivalent to selection on R', given a constant .
The response to selection for a binary disease trait can now be estimated in two ways; based on the threshold model and on the SELIRS or SELDCRS model. The threshold model is fully described by the prevalence and heritability of resistance, but the epidemiological models require a greater number of parameters. While estimates for most are available, prediction of the response to selection also requires knowledge of variation in the parameter under selection. One way to approach this is to parameterise the threshold model and the epidemiological models so that (for example) they give the same predicted response to one round of selection, and then study the longer-term predicted selection responses from these models.
In the threshold model, the response to selection on the observed scale is equal to the change in prevalence p* 1p* 0 , with p* 1 being the prevalence predicted in the threshold model after one round of selection. For the epidemiological model to predict the same response to selection we use the standard selection response equation ih 2  x = p 1 *p 0 *, so that  x = (p 1 *-p 0 *)/ih 2 , ( 1 6 ) with  x being the implied phenotypic standard deviation of the trait under selection. This standard deviation can now be used for parameters in the SELIRS and SELDCRS models, in order to calibrate selection responses. Some properties of responses to selection for these models and their dependence on various parameters are derived in Appendix 1.
If selection is on 1/R', or one of its animal-components other than , i.e.  or , and assuming that 1/R' is normally distributed, the phenotypic standard deviation  x that leads to the same response to one round of selection in the SELIRS model as the threshold model (p* 1 -p* 0 ) can be calculated (see Appendix 1) as with R' = 20,  = 0.025 day -1 and p 0 = 0.5.

Application to footrot
To investigate the predicted response to selection in the various models and parameter combinations, the equilibrium prevalence of footrot was calculated over 20 rounds of selection on 1/ and 1/R' in the SELDCRS and SELIRS models, using discrete generations and a heritability of 0.3, assuming a normal distribution for 1/, 1/ and 1/R' with constant underlying variances and ignoring the Bulmer effect (whereby genetic variation is reduced as a consequence of selection). With the initial prevalence of footrot set at 0.08, the response to one round of selection using a threshold model was calculated to estimate  x according to (16) and (17). Then (14) and (8) were used to estimate R' given the other parameters. Changes in the parameter under selection (1/, 1/ or 1/R') that would result from selection of a random sample of healthy animals were calculated for each generation and the new value inserted in (14) or (8) as appropriate to obtain the equilibrium prevalence in the next generation.

Predicted progress in the threshold, SELDCRS and SELIRS models
Predicted responses to selection according to the SELD-CRS model with selection on 1/ (the time an animal is diseased and infectious) and the SELIRS model with selection on 1/R' are identical for values of  = 0.2 and R' = 2.91 (Fig. 1). Both responses are larger than the threshold model prediction after the first round (when they were fixed at the same value). This is the result of differences in the relationship between the trait under selection and prevalence which is close to linear for SELIRS and SELD-CRS but not for the threshold model. As expected, the SELIRS model predicts a significant additional response from selection on 1/ compared to the SELDCRS model, in which the total infection period is held constant in this parameterisation. All graphs have a similar shape, showing diminishing returns at lower prevalences.

Sensitivity of the predicted response to selection to  and R'
A value of R' of just greater than 1 leads to predicted prevalence of footrot quickly going to 0 (Fig. 2). The reason is that R' drops below 1, so that the infection is not expected  [11] to be maintained in the population. In contrast, very high values of R' (such as 18) do not lead to a noticeable additional predicted response. With typical values for R' for footrot of about 1.5 an additional response in prevalence of over 2% of animals can be expected after a few generations, but this difference decreases again in later generations.
For more extreme cases, with R' and p* large and  small, a different picture emerges with the SELIRS model actually predicting a lower response to selection for larger  in early generations than SELDCRS or selection on R' (Fig.  3). This is because the variation in 1/ is relatively small and changes have little effect on R'.
Sensitivity to  and  Additional calculations (results not shown) have confirmed that the predicted response to selection on  or R' does not depend on values of the rate at which latentlyaffected animals become infectious () or the rate at which recovered animals lose immunity (), given , R' and p*.

Discussion
A model was developed to predict the response to selection for resistance to an endemic bacterial disease, and this model was then applied to footrot in sheep. While the exact description of the epidemiological process requires a great number of parameters, many of which are poorly known, the prediction of relative responses to selection depends on only a few parameters. By using the threshold model to predict the response to one round of selection, and setting this as the standard, the only epidemiological parameters required are the rate at which animals recover () and the reproductive rate (R'), alongside the heritability of resistance to the disease (or the parameter under selection).
It was shown that if the notional reproductive rate R' is under selection, changes in the prevalence are proportional to changes in 1/R'. If selection is for larger  (i.e. quicker recovery) the SELIRS model predicts a response that consists of the direct effect of animals recovering more quickly, plus an additional component arising from the resulting lower pathogen burden.
Predicted response to selection for resistance to footrot depending on the model and the trait under selection The expectation might be that if fewer animals get infected, there should also be a reduction in pathogens. Formula (7) shows that the equilibrium prevalence only depends on , i.e. susceptibility, through R'.
An explanation is that a decrease in susceptibility only means that it takes longer for animals to get infected when facing the same challenge (1/ days), but eventually they still become infected and shed the same number of bacteria. The advantage of a lower susceptibility, i.e. longer time until infection, is that it takes animals longer to complete the full SELIRS cycle, thereby reducing the number of animals that are diseased at any point in time prior to the equilibrium being reached.
Raadsma et al. [11] and Nieuwhof et al. [10] have estimated the heritability for resistance to footrot based on the genetic variation within a population, rather than the response to selection. Therefore, these estimates are independent of possible effects of reduced pathogen burden in selected populations.
Little information is available on the value of the recovery rate , and one reason is that prompt treatment of affected animals means it is not fully expressed (i.e. 1/ is censored). The prevalence without treatment could then be considerably higher and R' larger. This scenario was investigated with an initial prevalence of 50%, and it was found that in early generations the expected additional effects are in fact negative, but the effect becomes positive later on. The reason for this negative effect, best understood by comparing (8) and (14) is that with increasing recovery rate in the SELDCRS model animals spend increasing times in the C phase, i.e. they are no longer considered diseased but continue to spread bacteria. This slows down the whole cycle, with there being fewer susceptible animals compared to the SELIRS model, where recovered animals move directly to the phase of immunity. For our application to footrot this effect may be considered an artefact of the model, occurring only under extreme assumptions, rather than of biological importance, but it may be relevant for other diseases. Previously, MacKenzie and Bishop [14] have shown that in the SIR model applied to viral diseases, if R 0 is high then it may take many generations of selection before the expected number of animals infected during an epidemic is expected to decrease. This, also, is a scenario in which the Predicted response to selection for resistance to footrot depending on the model used, the notional reproductive rate R' and the initial recovery rate  Knowledge of the traits that show genetic variation is clearly important. While the current results suggest that this can be done by comparing the variation within a population with the response to selection, in practice this will be very difficult because it would require an unselected control or population with the same environmental conditions, but not affected by decreased shedding of bacteria by the selected animals. An alternative would be to estimate parameters directly from the length of time various epidemiological stages last, in selected and unselected populations following (deliberate) infection, comparable to the figures in Table 3.
In the absence of any estimates for the genetic variance of resistance to footrot, this study used the threshold model to standardise response to selection. Following the concept of an underlying normally distributed trait in the threshold model, normal distribution were assumed for 1/R' and 1/. The inverse of the recovery rate 1/ is a length of time, and it seems plausible that it has a positive skewness, with no negative values, and some animals taking extremely long to recover. A positive skewness looks likely for R' as well, especially for scenarios where mean R' is in the critical range just above 1, with some animals potentially being extremely infectious. Positively skewed distributions for the inverse 1/ and 1/R' would for instance occur if  and R' were normally distributed. Under these scenarios, relative responses to selection can be recalculated with appropriately altered selection intensities. However, it should be remembered that a normally distributed liability in the threshold model is also an assumption that can be challenged.
It was shown that, under the prevailing assumptions, the reduction in prevalence at a given R' does not depend on the rate of loss of immunity  and the rate of conversion of latently infected animals to the infectious state . This does not mean that these parameters are not important for the potential genetic progress; it implies that once the response to one round of selection is known it is possible In the current study, a constant environment and a homogeneous population have been assumed. In practice, environmental conditions will vary and this may affect survival of bacteria in the environment or the animals' phenotypes, while there may also be different classes of animals, e.g. adults and offspring with different phenotypes with regard to the disease. All these variations can be investigated based on the SELIRS equations, but may require running of a dynamic algorithm that calculates daily prevalence, rather than relying on equilibria. In rapidly changing environments the time to reach equilibrium will become an important factor, with potentially no equilibrium being attained by the time prevalence is measured or selection decisions are made. In a stable environment, changes in parameters as a result of selection will lead to only small changes in the expected equilibria so that a new equilibrium can quickly be established.
Selection in this study was based on own performance and one observation per animal with disease resistance as the only breeding goal. In practice, information from relatives and repeated measurements will increase the response to selection. On the one hand, assuming that resistance to footrot is not genetically correlated to any other traits under selection, selection on an index of traits will decrease the expected response to selection for resistance. On the other hand, disease information on relatives will greatly improve the potential selection response rates for improved resistance. While all these considerations affect the magnitude of the response to selection, essentially by changing the 'h 2 ' term in the response equations, they do not alter the nature of the epidemiological effects of selection. Therefore, simple extrapolation is appropriate.
The models developed in this study are used to consider an endemic bacterial disease with bacteria being transmitted through the environment, where they can only survive for a limited period of time. The models can be applied to a variety of diseases and host species, where these conditions apply. The general trend of results is in fact similar to that seen for a different disease, ruminant gastro-intestinal parasitism, as shown by Bishop and Stear [15]. One difference is that these authors had better estimates of some traits, especially the rate at which animals spread infection, as this is captured in the faecal egg count trait.
Based on the current study it can be expected that selection for resistance to footrot in sheep will be more considerably effective, especially in the medium term, than purely genetic models predict. There are, however, many other important issues to consider in a practical breeding programme, such as obtaining consistent disease scores across a population of sufficient size and the simultaneous selection for other traits, which may be correlated, on a phenotypic or genetic level, with resistance to footrot.
In summary, this paper presents a novel epidemic model, applied to footrot in an attempt to explore likely responses to selection. A key parameter for the model, and also from a biological perspective, is the recovery rate. Given the long time that it takes animals to recover from the disease without human intervention, low values for the rate of recovery () seem likely. If this is indeed the trait under selection when selecting for increased resistance, then the response to selection in terms of observed prevalence, including effects of reduced pathogen burden, could in the medium term be double that predicted by purely genetic models.

Derivation of R' for the SELIRS model
Assuming that N is large, so that S is approximately equal to N, in the SEIR model an infected animal sheds  infectious doses over 1/ days, these doses survive for 1/ days infecting N/ daily so that R' = N/. A more formal derivation is given in [8].
The extra L step in the SELIRS model does not affect this, as all latently infected animals will (sooner or later depending on ) become diseased. For most parameter values, the loss of immunity (R animals reverting to S) does not affect the number of secondary infections, but extreme parameter values (long-lived environmental contamination combined with a short period of immunity) may lead to more secondary infections.

Derivation of numbers of animals in various categories at the equilibrium in the SELIRS model
At the equilibrium (denoted by *) all derivatives, dI/dt etc are equal to 0, so that from (4) it follows that: I* = L* = (N-S*-I*-R*).

Summary parameters in the SELDCRS model
Following Appendix A, R' for the SELDCRS model can be calculated from the product of number of infectious doses produced at rate  by an infectious animal over the duration of the infection (1/+1/), the time these infections survive in the environment (1/) and the number of animals that will be infected at a rate of  N, so that: As with the SELIRS model in Appendix B, the equilibrium prevalence p* can be calculated by setting all derivatives to 0. Combining (3) and (10) then gives S* = D*/E*, and from (9) and (11) it follows that so that: From (10), (11) and (12), L* = D*/, C* = D*/ and R* = D*/. Substituting these values in L+D+C+R = N-S and rearranging yields: Therefore, from (14) the equilibrium prevalence is Derivation of the predicted response to selection on  or  and demonstration that the rate of response in p* depends on  and , given R' and p* in the SELDCRS and SELIRS models SELDCRS model From (15) , with (K set to be constant).
Therefore . Now consider the derivative d(1/p*)/d. Solving this equation we obtain: Therefore: Thus, in the SLEDCRS model the predicted response in p* to changes in  is a linear function of (1/ + 1/), and it also depends on R' and p* but under the model assumptions it is independent of .

SELIRS model
Define a new parameter Q' so that and Q' is independent of . Since Q' = R' we have: Equation (D1) can be rearranged as Solving the derivative d(1/p*)/d we obtain: Rearranging terms yields: There are two important differences from the derivative for p* in the SELDCRS model. Firstly, there is the extra element , which is more important at low values of Q' (i.e. low values of R'). Also, the second part of the equation includes , so that higher values of  lead to higher rates of change.