Modeling Marek’s disease virus transmission: A framework for evaluating the impact of farming practices and evolution

Marek’s disease virus (MDV) is a pathogen of chickens whose control has twice been undermined by pathogen evolution. Disease ecology is believed to be the main driver of this evolution, yet mathematical models of MDV disease ecology have never been confronted with data to test their reliability. Here, we develop a suite of MDV models that differ in the ecological mechanisms they include. We fit these models with maximum likelihood using iterated filtering in ‘pomp’ to data on MDV concentration in dust collected from two commercial broiler farms. We find that virus dynamics are influenced by between-flock variation in host susceptibility to virus, shedding rate from infectious birds, and cleanout efficiency. We also find evidence that virus is reintroduced to farms approximately once per month, but we do not find evidence that virus sanitization rates vary between flocks. Of the models that survive model selection, we find agreement between parameter estimates and previous experimental data, as well as agreement between field data and the predictions of these models. Using the set of surviving models, we explore how changes to farming practices are predicted to influence MDV-associated condemnation risk (production losses at slaughter). By quantitatively capturing the mechanisms of disease ecology, we have laid the groundwork to explore the future trajectory of virus evolution.


Introduction
Marek's disease virus (MDV), the causative agent of Marek's disease (MD), imposes a substantial economic burden on chicken meat and egg production, costing the worldwide poultry industry in excess of 1 billion USD per year (Morrow and Fehler, 2004). Historical control measures have at least twice been undermined by virus evolution, leading to speculation that future evolution could undermine current control (Nair, 2005). The ecology This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/). of the disease appears to be the driving force behind past evolution, with explanations invoking vaccination (Witter, 1997;Atkins et al., 2013a;Read et al., 2015), rearing period duration (Atkins et al., 2013a;Rozins and Day, 2017), and virus persistence during downtime between bird flocks (Rozins and Day, 2017). Understanding the ecology of the virus is thus a key component in predicting whether and when control efforts will lose efficacy. Such an understanding is also crucial in developing immediate responses should the efficacy of current control measures wane. Yet the ecology of MDV is poorly understood. This is perhaps most clearly exemplified by the conventional wisdom that the virus is ubiquitously found on industrialized poultry farms (Office International des-Epizooties, 2010;Dunn, 2013), despite recent surveillance data suggesting that the virus may not be present on a large fraction of farms (Groves et al., 2008;Wajid et al., 2013;Walkden-Brown et al., 2013;Bettridge et al., 2014;Kennedy et al., 2015bKennedy et al., , 2017Ralapanawe et al., 2015).
Mathematical models of disease ecology can provide valuable insight into infectious disease dynamics. Such models quantitatively relate changes in ecology to changes in disease dynamics, which is particularly useful when experimental manipulation is unethical or, as with commercial-scale chicken rearing, financially costly. Models provide cheap and safe opportunities to explore the impact of system manipulation on pathogen control, and this approach has been applied to MDV (Atkins et al., 2013a,b;Day, 2016, 2017). The reliability of a model, however, can only be assessed by challenging it with data, and this has never been done for MDV. Here we develop a suite of models to describe MDV dynamics on commercial broiler farms, and we use model selection methods to identify the ecological mechanisms that are most important to explaining MDV dynamics in the field.
Poultry intended for consumption are inspected and condemned for a condition called "leukosis" at the time of processing. This condition can be caused by various diseases, but in chickens reared for meat, leukosis is almost exclusively caused by MD (Sharma, 1985). Current rates of condemnation due to leukosis are extremely low (Kennedy et al., 2015b), but future reductions in MD vaccine efficacy due to virus evolution might cause leukosis rates to increase, as was documented with the erosion of vaccine efficacy in the past (Witter, 1996). A method to relate changes in farming practices to changes in risk of condemnation would therefore be a useful tool should virus evolution continue along the trajectory of the past.

Model construction
We model the transmission and persistence of MDV within and between flocks of broiler chickens on commercial poultry farms. Our models are constructed assuming standard rearing practices in Pennsylvania, United States. These practices are fairly standard for commercial poultry rearing across much of the developed world.
Industrial-scale rearing of broiler chickens on farms tend to follow an "all-in, all-out policy," meaning that all chickens within a house are reared as a single-aged cohort of birds. We refer to a cohort of birds that occupy a single house on a farm as a flock. Birds are placed on litter that consists of wood chips and sawdust at one-day-old, and the birds are reared in this environment until they are ready for processing. Birds in houses are provided ad libitum food and water. Temperature, humidity, and air quality are maintained by a combination of active ventilation through fans or wind tunnels and heating. Flocks are collected for processing when sufficient time has elapsed for birds to reach a particular target weight.
While chickens are being reared, houses accumulate "chicken dust," a by-product of farming that consists of bits of food, epithelial cells, dander, bacteria, and feces (Collins and Algers, 1986;Pandey et al., 2016). The amount of dust produced by birds increases as birds grow (Islam and Walkden-Brown, 2007;Atkins et al., 2013a). Infectious MDV can be contained in this dust (Carrozza et al., 1973), being shed with the epithelial cells of infectious chickens and transmitted through the inhalation of virus-contaminated dust (Colwell and Schmittle, 1968). The concentration of MDV in dust can be measured through quantitative polymerase chain reaction (qPCR) (Baigent et al., 2005(Baigent et al., , 2016Islam et al., 2006). Our model is constructed with this type of data in mind. We thus track the infection status of birds as well as total dust and total virus quantities.
Coming from extremely hygienic hatcheries, chickens are unexposed to MDV when first placed in a house. Shedding of virus from a bird can begin as early as one week post exposure to virus and tends to reach maximal levels two to three weeks post exposure (Islam and Walkden-Brown, 2007;Read et al., 2015). Once reached, virus shedding stabilizes at peak levels for the duration of a broiler chicken's life (Islam and Walkden-Brown, 2007;Read et al., 2015). Shed virus can infect other chickens, causing the pathogen to spread to other hosts in the flock. Even if virus were absent on a farm, it is possible that it may be introduced from outside sources, for example through dispersal in the air from nearby farms, on feed trucks, by service technicians, or by other farm visitors.
Typical commercial broiler farms vaccinate against MD by using bivalent vaccination (Morrow and Fehler, 2004). Although vaccinated birds can still be infected with MDV and can still shed MDV (Witter et al., 1971;Islam et al., 2008;Ralapanawe et al., 2016), vaccination greatly reduces clinical signs of disease (Witter et al., 1971). This, along with other measures to ensure bird health, means that total mortality from hatch to processing is typically minimal (≈3% and ≈8% in the two farms used for model inference below -in line with the national average of 4.8%, National Chicken Council, 2016). We therefore assume that bird mortality is negligible in our model. A schematic representation of the infection dynamics are shown in Fig. 1, corresponding to the following set of mathematical equations: where a c μ a LN − σ a 2 2 , σ a .
In the above equations, S(t) is the number of birds susceptible to virus infection at time t. E n (t) is the number of exposed birds in class n at time t where there are 5 total exposed classes (see below), I(t) is the number of infectious birds at time t, D(t) is the total milligrams of dust in the house at time t, and V(t) is the total number of virus copies in the house at time t. α c is the log normally distributed transmission rate of the virus in flock c determined by mean μ α and adjusted scale parameter σ α , β is the transition rate between exposed classes, γ is the rate at which dust and virus are removed from the house through ventilation, and δ is the rate at which virus decays in a house. All birds in the house produce dust at a rate of d(τ), which is a function of flock age τ. Keeping in mind that bird mortality is assumed negligible, the total number of birds in the house at any time is equivalent to the initial number placed S 0 . Infectious birds from flock c also produce virus at a concentration of a c per mg of dust produced, which is log normally distributed with mean μ a and adjusted scale parameter σ a . Additionally, to allow for the possible reintroduction of virus on farms, a quantity of M μ infectious virus copies are introduced at rate M r .
We model the exposed class as a series of compartments to create an incubation period of the virus that follows a gamma distribution (Wearing et al., 2005). In practice, we used N = 5 incubation classes with a transition rate β chosen to produce virus shed rates over time that were consistent with average viral copy number in feather tips of experimentally infected birds (Read et al., 2015), which is strongly correlated with viral shedding intensity (Baigent et al., 2013).
The above model describes disease dynamics within a flock of chickens, but a complete model of MDV dynamics must also incorporate the persistence of dust and virus between flocks. We follow the practice of nesting within-flock disease dynamics into a discrete-time model to capture the multiple timescales of dynamics (Dwyer et al., 2000;Elderd et al., 2008).
MDV is highly stable in the environment, and ventilation fans are typically off between flocks of chickens. As a consequence, the duration of downtime between flocks of birds is unlikely to strongly influence virus dynamics. Rather, the reduction in virus and dust between flocks is likely the result of cleaning. This cleaning process can be highly variable between flocks and between different farms. On some farms, bedding material, or "litter", is reused for several flocks, up to a year or even more. Between flocks on these farms, only caked or wet areas of litter are removed and top dressed with new litter. On other farms, litter is completely changed between each flock. Full removal of litter involves removing all bedding material down to the clay floor or cement floor. Also between flocks, farmers may or may not use forced air blowers to blow down chicken dust from fans, walls, and other above ground surfaces. Blow downs are sometimes followed by a "thermal fog", comprised of a vaporized disinfectant and an insecticide. Typically once per year, farmers will do a "wet clean and disinfect", which is a wash down by spraying all surfaces with high pressure, low volume water. These wet cleans are typically followed by spraying with disinfectant. The outcome of the cleanout process is that dust and virus might be removed concurrently through mechanical processes, or virus may degrade without the removal of dust due to the action of disinfectants. We thus model the initial virus and dust concentration in flocks using the following model: where C c Beta(μ C , ν C ), f c Beta(μ f , ν f ) .
Here, ψ c and ω c are the respective total dust and total virus at the start time of flock c. l c−1 is the load out time of flock c −1, making D(l c−1 ) and V(l c−1 ) the respective dust and virus at the end of the previous flock. C c is the fraction of dust and virus removed between flocks c −1 and c by the physical removal of dust, and f c is the fraction of persisting virus that is removed between flocks c −1 and c due to chemical disinfectants and other environmental decay factors. Note that C c and f c are beta distributed random numbers parameterized by their means (μ C , μ f ) and sample sizes (ν C ,ν f ). The initial number of birds susceptible to virus infection at the start of flock c is equal to the total number of birds placed in the house, and the initial number of exposed and infectious birds are equal to 0. A list of all model parameters is provided in Table 1.
To account for the fact that bird population sizes are finite, and for computational convenience, the above differential equations were effectively replaced with their corresponding probabilistic transition equations (Kennedy et al., 2014(Kennedy et al., , 2015a, where time was discretized to units of one day. Bird numbers were discretized to integer values, and transitions between classes were determined by binomial probabilities. For example, the number of birds that became exposed in a time step Δt was distributed Binomial(S(t), 1 −e −α c V(t)Δt). Virus and dust continued to be treated as continuous variables. For example, the amount of virus that was removed through ventilation in a time step Δt was equal to γV(t)Δt.
The above model allows for variation between flocks of chickens in virus transmission rate, virus shed rate from infectious birds, dust and virus cleanout efficiency during downtime, virus decay rate during downtime, and virus reintroduction. These mechanisms for generating variation can however be removed from the model. Variation in virus transmission rate or in virus shed rate can be removed by respectively assuming that σ α = 0 or σ a = 0. Variation in the cleanout efficiency of dust and virus or in the decay rate of virus during downtime can be respectively removed by assuming that ν C or ν f are extremely large. Virus reintroduction can be removed by assuming that M r = 0. Note that this additionally removes the parameter M μ from the model. To test for the importance of between-flock variation in transmission rate, virus shed rate, dust and virus cleanout efficiency, and virus degradation, and for the importance of virus reintroduction, we therefore generated every version of the model that includes or excludes each of these factors (32 models in total), and we compare the fit of each of these models to data. Our null model is the model that lacks all of these factors.

Parameter estimation, model evaluation, and model comparison
To visualize the basic dynamics imposed by the model, we generated simulations of the most complex model using an illustrative parameter set (Table 1), and we examined these simulations to identify common features in the model realizations. We then used model inference to estimate parameters and determine how well the above model formulations describe data of MDV dynamics on commercial poultry farms. Several of our model parameters are already known from previous experiments. These include the rate of dust production as a function of bird age d(τ) (Atkins et al., 2013a), the incubation period of the virus (Islam and Walkden-Brown, 2007;Read et al., 2015), which is determined by β, and the mean transmission rate of the virus μ α (Atkins et al., 2011(Atkins et al., , 2013b. Placement dates, load out dates, and flock sizes were fixed at known values from the respective datasets being modeled. Point estimates of the remaining parameters were determined from the data, as described below. The data used in this study are qPCR data reflecting the virus copy number (VCN) per mg of dust. These data were originally collected as part of a surveillance study quantifying the spatial and temporal variation of MDV across farms in Pennsylvania (Kennedy et al., 2017). We use two representative datasets from that study: the data from Farm A House 1, and the data from Farm E House 4 (Fig. 2). These datasets were chosen because they were among the most exhaustively sampled farms and houses, and because they include highly disparate patterns of MDV dynamics. We therefore fit these datasets separately. The To construct a likelihood function, we first determined the marginal likelihood of virus being detectable by qPCR as a function of mean virus concentration. We treated 100 virus copies per mg of dust as the limit of detection by qPCR, because below this concentration, virus detection is unreliable (Kennedy et al., 2017). Using the longitudinal data from Kennedy et al. (2017), we found that the probability of detection in biological replicates can be well described by a probit regression of the mean log 10 VCN/mg of dust plus 1. Using a generalized linear model with binomial data and a probit link we found an intercept of −2.206 and a slope of 1.555 (Fig. 3). For data that did not exceed our limit of detection, the likelihood was one minus the value resulting from this regression assuming the mean given by the model. For data that did exceed the limit of detection, the likelihood was the probability of detection multiplied by the probability of observing the particular virus concentration seen in the data (data were log 10 plus 1 transformed), where the standard deviation was determined by a linear regression of the standard deviation with respect to the mean. We thus modeled this component of likelihood as a normal distribution with mean μ given by the model and standard deviation 1.127 -0.151μ (Fig. 3).
Parameter estimation was performed using maximum likelihood approximated by iterated filtering. This was implemented using the 'mif2' function of the 'pomp' package in the R statistical computing language (King et al., 2016(King et al., , 2017R Core Team, 2017). Fitting arguments used in 'mif2' included 'Np = 1000', 'Nmif = 1500', 'cool-ing.fraction.50 = 0.5', and 'rw.sd = 0.02' for all regular parameters or 'rw.sd = 0.2' for initial value parameters. Arguments 'Np', 'Nmif', and 'cooling.fraction.50' were modified for any model that repeatedly had trouble finding high likelihood parameter space. This was performed until at least 20 parameter sets were generated that had log likelihoods within 50 points of the current observed maximum likelihood, to ensure the our fitting routine had multiple opportunities to explore reasonable parameter space. The likelihood of each final parameter set was determined using particle filtering. In practice this was implemented by averaging the likelihood scores from 10 iterations of 'pfilter' with 'Np = 10,000'. All code and data are publicly available (https://github.com/dkenned1/KennedyDunnRead).
Models were compared using Akaike's Information Criterion (AIC) (Burnham and Anderson, 2002). To aid in comparison, we present ΔAIC and model weights (Burnham and Anderson, 2002). We further explored the ability of the models to explain the data by simulating 5000 realizations from the weighted set of reasonable models, defined as any model with weight greater than 0.01. Each model was simulated using its respective maximum likelihood parameter estimates. We present envelopes that encompass central quantiles from these simulations.
MD takes four to twelve weeks post exposure to develop, and so the risk of leukosis condemnation is likely to positively correlate with the fraction of birds exposed to virus at least 30 days prior to load out. Hereafter, we refer to this fraction of birds exposed to virus 30 or more days before loadout as the "condemnation risk". We present envelopes that encompass the fraction of birds exposed to MDV over time for the maximum likelihood estimated parameters from the set of reasonable models. We explore the impact of altering features of poultry rearing by exploring how two fold increases or decreases in model parameters relative to maximum likelihood estimates alter condemnation risk (rearing duration is altered by plus or minus 5 days).
To explore whether the parameter estimates that best explain the data from Farm A House 1 also provide a reasonable estimate for the data from Farm E House 4, and to explore the impact of rearing duration on MDV dynamics, we repeat these methods, applying the fitted models from Farm A House 1 to the rearing parameters (i.e. placement dates and load out dates) from Farm E House 4. We also do the reverse to see whether the best fit models and parameters from Farm E House 4 provide a reasonable explanation for the data from Farm A House 1.

Results
The structure of the above model comes from a qualitative understanding of MDV natural history. Relating this qualitative understanding to quantitative dynamics, however, at a minimum requires simulation of the model with parameter values within their respective plausible ranges. In Fig. 4, we show simulated realizations of the model for conventional duration (40 days) and long duration (80 days) rearing periods using the illustrative parameter set from Table 1. In Fig. 5 we show the envelopes produced by simulating the model 5000 times. These two figures demonstrate at least five features of MDV dynamics that emerge from our mechanistic model. First, the sampled virus dynamics are inherently stochastic. This stochasticity can be visualized as the jaggedness in the curves, and it becomes less pronounced as virus concentrations increase (Fig. 4). Second, different realizations of the model generate variable virus concentrations even when all simulations are performed with a single set of parameters (Fig. 4). Third, virus densities within flocks tend to form "U" or "J" shaped trajectories, with virus concentrations starting relatively high, decreasing early on during the rearing period and increasing back to high concentrations later (Fig. 5). Fourth, the sharpness of the "U" shaped trajectories vary between different flocks and realizations, sometimes generating less pronounced troughs than other times (Fig. 4). Fifth, the rearing duration can strongly alter virus dynamics (Fig.  5).
By fitting our suite of models to the data from Farm A House 1, we found that three models have ΔAIC < 2, and are thus essentially indistinguishable ( Table 2). Each of these three best models include variation between flocks in the transmission rate of the virus σ α , in the shedding intensity of virus from infectious birds σ a , and in the cleanout efficiency of dust and virus between bird flocks ν C , suggesting that these features of the model are important to explaining virus dynamics. The decay of virus between flocks ν f and the reintroduction rate of virus M r on the other hand are each excluded from at least one of the best three models, and therefore might not need to be included in models of virus dynamics on this farm.
We find somewhat different results when we examine the Farm E House 4 data. While we again find that several models provide comparable fits to the data according to AIC (Table 3), we find that for this dataset, every model that includes reintroduction of virus (M μ and M r ) has a better AIC score than every model that excludes this mechanism. This result suggests that virus reintroduction is important to explaining dynamics on this farm. In addition, we find that virus reintroduction alone is unable to explain the breadth of variation present on the farm, as noted by the fact that model M1 which includes only stochastic virus introduction is not in our set of models with weights greater than 0.01. Rather, we find that each of the models that provide reasonable explanations for the data also include variation between flocks in either virus transmission rate or virus shed rate from infectious birds.
In Fig. 6 we show the parameter estimates that emerged as the maximum likelihood estimates from the set of reasonable models as defined by having model weights greater than 0.01. Using these parameter estimates in simulating our model we generate envelopes that shows the variation in the data expected given the model structure and maximum likelihood parameter sets. These simulations incorporate both process error (variation in virus dynamics between simulated realizations) and observation error (variation in the data due to measurement error). Note the differences in the envelopes that are generated from these different datasets (Fig. 7). On Farm A House 1, virus concentrations generally remain at high levels across many different realizations. On Farm E House 4 on the other hand, the virus concentrations quickly fall to low levels, but with reasonable chances for outbreaks to reemerge in future flocks due to virus reintroduction. These patterns are similarly reflected in the data. In Fig. 8, we show how the fraction of birds exposed to virus changes over time for the two farms using maximum likelihood estimates from the set of reasonable models, revealing that substantially more birds are exposed to virus on Farm A House 1 than on Farm E House 4. Fig. 9 shows how changes in farming practices (i.e., changes in parameter values of the model) are predicted to impact condemnation risk. Note that despite large differences between the two farms in overall condemnation risk, altering most parameters has similar effects on both farms. The transmission rate scale σ α , the mean virus shed rate μ a , the virus shed rate scale σ a , the mean between-flock cleanout efficiency μ C , and the number of birds placed S 0 have relatively large effects on condemnation risk for both farms. The ventilation rate γ, the sample size of cleanout efficiency ν C , and the mean decay rate of virus between flocks μ f have relatively large impacts on only one of the farms. The decay rate of virus during the rearing period δ, the sample size of virus decay between flocks ν f , the rate of virus reintroduction M r , the quantity of virus reintroduced M μ , and the initial virus population size V 0 have relatively little effect on the condemnation risk for either farm.
Applying our maximum likelihood parameter estimates from Farm A House 1 to the rearing conditions of Farm E House 4, and vice versa, allows us to visualize the relative impacts of differences in model parameters and differences in rearing duration on virus dynamics. Fig. 10 shows that the maximum likelihood parameter estimates inferred from virus concentration data on one of these farms do not accurately predict virus concentrations on the other farm, which can be seen in the mismatch between the prediction envelopes and the data. Comparisons between Figs. 7 and 10 show that if Farm A House 1 were to shorten rearing durations to those of Farm E House 4, it would lead to an approximate 10 fold reduction in median virus concentration. If Farm E House 4 were to extend rearing durations to those of Farm A House 1, it would lead to an increase of median virus concentrations from undetectable levels to between 10 2 and 10 3 virus copies per mg of dust.

Discussion
Here, we constructed mechanistic models to describe MDV dynamics across multiple flocks of chickens reared on a commercial chicken farm. We used data on virus concentration to identify the mechanisms that drive variation in virus dynamics between farms (Tables  2 and 3), and to provide the first data-inferred estimates of the parameters that define these mechanisms (Fig. 6). Our models were able to capture the dynamics of two different datasets, one in which virus intensities tended to be much higher than in the other (Figs.  2 and 7). The different dynamics in these two datasets arise from slight differences in model structure, parameter estimates, and rearing practices. We then showed that altering the parameters of these models influences condemnation risk (Fig. 9), giving insight into how MDV losses might be managed should virus evolution undermine existing vaccines.
Two previous sets of models have been developed to describe MDV dynamics. The first set by Atkins et al. (2013a,b) uses an individual-based approach that captures many of the biological details particular to MDV dynamics. The second set by Day (2016, 2017) uses an impulsive-differential-equation based approach that allows the models to be completely communicated through a set of equations. Our modeling approach is a combination of these two approaches. Like the Rozins and Day approach, our models have an underlying impulsive differential equation structure, making it easy to describe them. Like the models of Atkins and colleagues, our models allow for stochasticity and biological details that are not captured by simpler models. The key difference between the approach presented here and previous approaches, however, is that the parameters of our models have been determined by statistical integration with field data. Our best models are therefore those which have withstood confrontation with data, and these models are therefore more likely to produce accurate predictions than untested models.
Our previous surveillance on MDV dynamics has revealed that virus concentrations in dust tend to form U-shaped trajectories (Kennedy et al., 2017). We attributed this to dilution of virus early in the rearing period when few birds are shedding virus-contaminated dust followed by increasing concentration of virus later in the rearing period when many birds are shedding virus-contaminated dust (Kennedy et al., 2017). Here, we have shown that this qualitative shape naturally arises from the structure of our models (Fig. 5), and that our models quantitatively explain such patterns in field data (Fig. 7).
We used model comparison methods to identify the features of poultry farming that lead to variation in MDV dynamics between flocks and between farms. For the Farm A House 1 data, we found that all reasonable models included flock-to-flock variation in virus transmission rates σ α , virus shedding rates σ a , and cleanout efficiency between flocks ν C , thus highlighting the importance of these mechanisms for explaining dynamics on this farm (Table 2). These results suggest that flocks of birds differ in their susceptibility to virus and in their shedding rate should they become infectious. High susceptibility and shedding rates might be rapidly assessed in the field by determining the virus concentration in dust samples or feather tips (Baigent et al., 2016). Our results also suggest that cleanout efficiency varies on this farm and that this variation has had detectable effects on virus dynamics. For the Farm E House 4 data, we found that all reasonable models included flock-to-flock variation in either virus transmission rates σ α or virus shedding rates σ a , and they also all included the stochastic reintroduction of virus M μ and M r (Table 3). This again suggests that flock-to-flock variation in either disease susceptibility or virus shedding is affecting virus dynamics. Our finding that virus is stochastically reintroduced on this farm, also suggests that eradication of virus may be quite challenging. Even if a control program were capable of clearing infection from a farm, there appears to be substantial risk of reintroduction, meaning that control measures might need to be sustained indefinitely to maintain an infection-free farm.
Although our two datasets identified different mechanisms as key to driving MDV dynamics, this should not be surprising. Such an outcome might arise from true differences in the mechanisms that generate variation between farms, but it also might arise due to other differences in the data. For example, the importance of stochastic reintroduction of virus is likely realized on Farm E House 4, because virus concentrations are low in that dataset. Similar levels of virus reintroduction might also occur when virus is common such as on Farm A House 1, but because virus is common, this mechanism would have only a negligible impact on the dynamics of the virus. This outcome highlights that the absence of a mechanism in a reasonable model does not mean that the mechanism is not acting, but only that it is not important to explaining those particular data. Of the mechanisms tested by our approach, variable decay of virus between flocks ν f is the only mechanism that was not necessary to explain either of the datasets studied.
Examination of the maximum likelihood parameter estimates that arose from the sets of reasonable models reveals that there is substantial overlap of the range of these maximum likelihood parameter estimates across the datasets (Fig. 6). There are only three parameters where the ranges of the maximum likelihood estimates do not clearly overlap between the two farms. These are the virus shed rate scale σ a , the virus reintroduction rate parameter M r and the virus reintroduction quantity parameter M μ . We note that variation in virus shed rate σ a might be hard to distinguish from virus reintroduction on Farm E House 4, leading to an overestimate of this parameter. Likewise, the rate and quantity of virus reintroduction might be hard to quantify on Farm A House 1 because these parameters are likely obscured by the high virus concentrations on this farm. Both lines of reasoning are supported by the observations that 4 out of 12 of the reasonable models for Farm E House 4 do not include the parameter σ a , and 1 out of the 3 reasonable models for Farm A House 1 do not include the parameters M r and M μ , implying a great deal of uncertainty regarding the point estimates of these parameters.
The parameter estimates that arose from model fitting give further insight into model reliability, virus management, and perhaps even virus evolution. Regarding model reliability we can compare our parameter estimates to those measured in previous studies. To our knowledge, the only parameter in our model that was directly measured through lab experiments is the mean virus shed rate μ a . Two studies have provided estimates of this value from bivalent vaccinated chickens. Islam and Walkden-Brown (2007) found that birds shed 6.36 log 10 virus copies per mg of dust. Atkins et al. (2011) tested three virus strains of different virulence rank and found respective peak shedding rates of 6.04, 7.06, and 7.12 log 10 virus copies per mg of dust. Our estimate (mean = 6.30, s.d. = 0.61 log 10 virus copies per mg of dust) is in close agreement with these previous studies. The decay rate of virus δ has never been measured, but two studies have looked at the infectiousness of dust after storage at room temperature. Jurajda and Klimes (1970) found no change in virus infectivity after 44 days and Witter et al. (1968) found that 8 of 14 infectious dust samples maintained infectivity after 112 days. Our median estimate of δ implies a half life of viral persistence of 66 days (range 20-425 days), consistent with those previous experiments, further suggesting that our model is reasonably capturing the biology of the system.
Our estimates of virus reintroduction rate M r on Farm E House 4 imply a median reintroduction probability of 3.2% per day (range 1.7-5.1%). This equates to approximately one virus reintroduction event per month. Based on these estimates, a typical flock of chickens is likely to be exposed to virus even if the virus is not present at the time of placement. However, the timing of virus reintroduction is an important factor in determining whether virus amplification will take place. Introductions late in the rearing cycle, for example, are unlikely to amplify to high levels before birds are removed for processing, highlighting the importance of good biosecurity early in the rearing cycle. Nevertheless, the frequent reintroduction rate suggests that maintaining the absence of virus is likely to require rearing conditions that limit the long term amplification of virus.
Our estimates of the parameters dictating between-flock removal of virus and dust μ C and decay of virus μ f suggest that the cleanout and sanitization processes between flocks are fairly efficient, with the median estimate for cleanout efficiency at 95%, and the median estimate of virus decay at 69% (note that the median for μ C includes only the estimates from Farm A House 1 because the estimates of this parameter from Farm E House 4 are highly uncertain across models). This relatively high efficiency implies that current practices are good at reducing dust and virus between flocks. This efficiency, however, may be a cause for concern given that Rozins and Day (2017) found efficient cleaning between flocks can promote the evolution of increased MDV virulence by decreasing the lifespan of the pathogen and thus decreasing the cost of virulence.
We found that condemnation risk is highly sensitive to the mean cleanout efficiency μ C and flock size S 0 when compared against the impact of the mean virus decay between flocks μ f , the ventilation rate γ, or the within flock virus decay rate δ. This suggests that, should vaccine efficacy be reduced by pathogen evolution, condemnation rates might best be managed by focusing efforts on blowing down dust, changing litter between flocks, and reducing flock sizes, as opposed to improving chemical treatments between flocks or increasing ventilation rates within flocks. Reducing the frequency M r or intensity M μ of virus reintroduction from other farms also appears to have little influence on condemnation risk. Nevertheless, M r would become a key parameter should the goal change from reducing to eliminating MD related condemnation.
We also found that condemnation risk is highly sensitive to the transmission rate scale parameter σ α , and the virus shed rate scale parameter σ a . Because the scale parameters σ α and σ a are directly related to variability between flocks in transmission rate and virus shed rate, our analysis demonstrates that condemnation risk is highest when flocks have little variability between them in transmission rate and virus shed rate. This observation highlights the value of intervening when a farm is having problems, even if the intervention is economically unsustainable in the long term, because a reduction in the transmission and virus shedding rates within a single flock can have large impacts on condemnation risk in future flocks. Currently, such a strategy is practiced through remedial actions taken in response to MD breaks such as changing the vaccination program to include Rispens vaccine. Cycling flocks of chickens between standard broilers and MD-resistant broiler might have similar benefits (Hunt and Dunn, 2013).
Our results highlight the importance of rearing duration on virus dynamics. Long rearing durations give the virus ample time to spread through a flock, and reach high concentrations in dust (Figs. 4 and 10). The historical decline in MDV-associated condemnation is strongly correlated with declines in rearing durations, making it difficult to tease apart the relative impact that rearing duration and other changes in biosecurity and vaccine administration have had on disease control (Kennedy et al., 2015b). Nevertheless, our model shows that rearing duration is one of several factor that can potentially be manipulated to control problems with disease should they reappear.
We note that the data used in this study were collected from farms in the United States, where universal bivalent vaccination of broiler chickens is standard. In other parts of the world, chickens are often given either different versions of the MD vaccine, or even no vaccine (Dunn and Gimeno, 2013). These differences would likely alter the ecology of the virus, and in turn, our parameter estimates. For example, vaccination is known to suppress virus shedding rates (Islam et al., 2008;Atkins et al., 2011;Read et al., 2015). But of course, there are numerous differences between poultry rearing in the United States and other places in the world, suggesting that many potential factors in addition to vaccination, may lead to differences in disease dynamics. Even in settings similar to those analyzed here, estimates of reintroduction rates and other parameters may be useful to operation managers, for example, by identifying farms that would benefit from improved biosecurity such as those located in areas of high chicken farm densities. We thus encourage application of our models to novel datasets, and to facilitate this, all code is available online (https:// github.com/dkenned1/KennedyDunnRead).
Given the history of MDV evolution undermining vaccine efficacy, any long term control strategy must carefully consider the ecology of disease transmission, and its impact on pathogen evolution. Our hope is that the data-tested models that we present here will lay the ground work for this approach. Future work might use these models to ask how altering rearing practices, approximated by changing model parameters, affects the profitability of poultry farming as well as how this might alter competition between virus strains and in turn the evolution of the virus. Schematic of the model. States and parameters are as described in the main text and Table 1.
Solid lines indicate transitions between model classes. Dashed lines indicate that producing dust and virus does not cause birds to leave their current model class. Dotted lines indicate the between flock persistence of dust and virus. Note that without altering the model, we depict the exposed class as a single group, where the time until an exposed host becomes infectious is gamma distributed with shape equal to 5 and rate equal to β. Data from Kennedy et al. (2017). Under the assumption that noise in the data is homoscedastic and normally distributed, bars show 95% confidence intervals around maximum likelihood estimates of virus copy number per mg of dust (Kennedy et al., 2017). Note that in rare circumstances, the error bars do not overlap the data points because the log mean virus concentration differs slightly from the maximum likelihood virus concentration. This discrepancy is partially due to the assumption of homoscedasticity in the data, an assumption that we relax to fit our models. The empirically derived likelihood function using data from Kennedy et al. (2017). Shown on the left is the probability of virus detection across biological replicates as a function of measured log mean virus concentration. Points are the raw data, and the line is the best fit probit regression line. Shown on the right is the observed standard deviation of biological replicates from the same study, discarding when log mean virus concentrations were below 10 3 . Points are again the data, and the solid line is the best fit linear regression line. Ten simulated realizations of the model shown in Fig. 1 using an illustrative parameter set ( Table 1). As in Fig. 2, red intervals show periods of down time, when birds are absent from houses. Both panels show realizations using the same model parameters, but where rearing duration is 80 days (top) or 40 days (bottom). Note the differences in dynamics between these panels, as well as the differences in dynamics across realizations and within realizations between flocks. These differences are the combined result of observation error, stochasticity introduced by the probabilistic spread of virus through a flock, and differences between flocks in transmission rates and virus shedding rates. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.) Envelope of model realizations generated from 5000 simulations with the illustrative parameter set (Table 1). Figure layout is identical that of Fig. 4. Black solid lines show the median model prediction at each time point. Grey shaded regions show the respective 50%, 75%, 95% and 99% central quantiles. Maximum likelihood parameter estimates for each of the models with AIC weights greater than 0.01. Shown in black are the models fit to the Farm A House 1 data. Shown in red are the models fit to the Farm E House 4 data. Each shape shows a different model. For parameters that are missing from a particular model, no point is shown. Fit of the model to the data. Red intervals as in Fig. 4   Fraction of birds present in a house that are currently in the exposed or infectious class. Note that we arbitrarily set the value to zero when birds were absent from the house. See Fig.  7 for interpretation of colors and lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.) Simulated condemnation risk for Farm A House 1 (black) and Farm E House 4 (red) when altering the model parameters (solid horizontal lines mark condemnation risk at maximum likelihood estimates). "0.5x" on the x-axis denotes halving the parameter value, and "2x" denotes doubling the parameter value, with the exception of the mean cleanout efficiency μ C and the mean virus degradation between flocks μ f . For these latter parameters "0.5x" is the case where twice as much virus or dust is held over between flocks, bounded to be non-negative (i.e. argmax(2μ −1, 0)), and "2x" is where hold over virus or dust is cut in half (i.e. μ + 1/2). Respectively, "−5" and "+5" denote decreasing or increasing the rearing duration τ max by 5 days. Note that the direction of the effect for each parameter is always the same between Farm A House 1 and Farm E House 4, but some parameters have disproportionately more effect on the condemnation risk of one farm than the other. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.) Fit of the model to the data when model weights and parameter estimates from Farm A House 1 are applied to Farm E House 4 (top), and the reverse (bottom). Note that the panels are swapped relative to Fig. 7 to facilitate comparison of the impact of rearing duration on virus dynamics. The lack of agreement between model predictions and data suggest that the parameter values governing virus transmission differ at least slightly between these farms.  Table 2 AIC table for Farm A House 1 dataset.   Table 3 AIC table for Farm E House 4 dataset.