Bayesian analysis of botanical epidemics using stochastic compartmental models
See allHide authors and affiliations

Edited by R. James Cook, Washington State University, Pullman, WA (received for review February 5, 2004)
Abstract
A stochastic model for an epidemic, incorporating susceptible, latent, and infectious states, is developed. The model represents primary and secondary infection rates and a timevarying host susceptibility with applications to a wide range of epidemiological systems. A Markov chain Monte Carlo algorithm is presented that allows the model to be fitted to experimental observations within a Bayesian framework. The approach allows the uncertainty in unobserved aspects of the process to be represented in the parameter posterior densities. The methods are applied to experimental observations of dampingoff of radish (Raphanus sativus) caused by the fungal pathogen Rhizoctonia solani, in the presence and absence of the antagonistic fungus Trichoderma viride, a biological control agent that has previously been shown to affect the rate of primary infection by using a maximumlikelihood estimate for a simpler model with no allowance for a latent period. Using the Bayesian analysis, we are able to estimate the latent period from population data, even when there is uncertainty in discriminating infectious from latently infected individuals in data collection. We also show that the inference that T. viride can control primary, but not secondary, infection is robust to inclusion of the latent period in the model, although the absolute values of the parameters change. Some refinements and potential difficulties with the Bayesian approach in this context, when prior information on parameters is lacking, are discussed along with broader applications of the methods to a wide range of epidemiological systems.
Compartmental models, in which individuals pass through a series of classes from susceptible through exposed to infected and removed, are fundamental to the study of the dynamics, invasion, persistence, and control of plant and animal diseases (1, 2). More recently there has been renewed interest in the incorporation of stochasticity into these models to capture the inherent variability of epidemics (3–7). By taking account of this variability, stochastic models provide a theoretical framework by which to predict the probability of disease invasion after introduction of inoculum or to compare the risks of failure of different disease control strategies. This demands knowledge of the values for the model parameters. Major new challenges therefore now lie in testing stochastic epidemiological models against experimental data to identify and estimate critical parameters associated with the control of disease. Parameter estimation, however, is not easy, and although it is an essential bridge between theory and application, it has tended to lag behind model development. Estimation is complicated by restrictions in epidemic data that often comprise short time series or have unobserved compartments (for example when individuals are infected but not infectious and cannot be distinguished from susceptibles). Here, however, we introduce and test a method based on a Markov chain Monte Carlo (MCMC) algorithm that allows a stochastic compartmental epidemic model to be fitted to experimental observations within a Bayesian framework. The approach allows the uncertainty in unobserved compartments to be represented in the parameter posterior densities and has broad applicability to many other epidemic systems. It is tested here for the control of a botanical epidemic that combines well regulated, replicated, and repeatable epidemics along with a parsimonious and biologically plausible model that includes uncertainty in observed compartments and quenching of susceptibles as they become resistant to infection (4, 8).
Although still in its infancy, rigorous parameter estimation for stochastic models of biological processes is now the subject of intense research activity, much of it exploiting the modern widespread availability of computing power via computational methods such as MCMC (9, 10). New statistical methodology for model fitting has already made an impact on the study of botanical epidemics for diseases of agricultural and plantation crops (e.g., refs. 4 and 11–13). These models all focused on a simple transition from susceptible to infected plants. We now address parameter estimation when allowance is made for an exposed class, whereby there is a latent period between infection of a plant and that plant becoming infectious. This builds on previous work by Gibson et al. (4) in which we fitted a stochastic epidemic model to experimental data from replicated microcosm experiments on epidemics of dampingoff caused by Rhizoctonia solani (Kühn) in radish (Raphanus sativus L.), in the presence or absence of the biological control agent Trichoderma viride (Pers ex Gray). R. solani is a ubiquitous soilborne fungal parasite that causes a wide range of economically important diseases on many different crop species. The Rhizoctoniaradish system is also an ideal experimental model allowing measurement of fast, replicated epidemics under controlled conditions in which the effects of demographic stochasticity on epidemics can be quantified (4).
The model considered by Gibson et al. (4) was a simple SI form for transitions between susceptible (S) and infected (I) individuals, given by: It represents a natural stochastic reformulation of a deterministic version previously considered by Kleczkowski et al. (14) with dual sources of infection in which susceptibility is quenched over time. Thus plants may be infected from a reservoir of inoculum in soil (primary infection, with rate r_{p} ) or by planttoplant transmission (secondary infection, with rate r_{s} ). Host susceptibility, s(t), decays with time, to capture the increasing resistance of plants characteristic to infection by R. solani and many other seedling diseases.
Thanks to the simplicity of the model, given the experimental data it was possible to derive an exact form for the parameter likelihood in terms of generalized Erlang distributions (15). Moreover, since the model represented only a single type of transition (from susceptible to infectious), the timevarying susceptibility could be easily incorporated in the likelihood calculation by a nonlinear transformation of the axis of time. Individual parameters were then estimated by computing profile likelihoods, and the results were compared between treatments with and without addition of the biocontrol agent (4). The study provided strong evidence that the presence of T. viride affects the rate of primary infection, whereas estimates of other parameters varied little between the treatments. The inference rests, however, on the simple twocompartmental model in which no allowance is made for a latent period between infection and subsequent transmission of infection. Here, we introduce the exposed or latent class (E) into the estimation of parameters so that we now estimate four parameters. These include the rate of primary infection, the rate of secondary infection, the rate of change in host susceptibility, and the latent period. This requires the implementation of a method for parameter estimation. Specifically we ask: is there evidence for a latent period? How is this affected by uncertainty in the observed classes? Is the inference about biocontrol robust to inclusion of a latent period?
The inclusion of the exposed class renders the computation of a parameter likelihood difficult. The timevarying susceptibility cannot now simply be accommodated via a transformation of time (4), since the transitions from E to I would no longer respect the new time scale. Moreover, it is not possible to derive an analytic solution for the parameter likelihood. This means that it is not so convenient to use profile likelihoods (16) to make inferences about individual parameters. However, several authors (e.g., refs. 10, 17, and 18) have successfully used Markov chain methods within a Bayesian framework to handle the problem of fitting nonlinear stochastic models to partial data, when the resulting likelihoods are intractable. This is the approach adopted here.
Inclusion of the exposed class complicates the interpretation of the observations. With the SI model, a natural assumption is made that observations record the number of plants in I at each time. However, in the presence of the additional E class the question arises as to whether the observed infections (determined by visual inspection) represent plants in I, or plants that are in either E or I. This uncertainty is common to many plant diseases. It could be accommodated by extending the SEI model by partitioning the exposed class, E, into two subclasses E_{u} and E_{o} , denoting, respectively, unobservable and observable exposed plants. Newly exposed plants would enter E_{u} before passing to E_{o} at the time when their symptoms became observable. This approach would add complexity to the model and, instead, we prefer to conduct and compare analyses under the two distinct assumptions that the infections observed at time t record I(t) and E(t) + I(t), respectively. These two scenarios can be considered to represent the extremes of the possible interpretation of the data.
Experimental Methods
The experimental data were derived from microcosm experiments on epidemics of Rhizoctonia dampingoff of radish, in both the presence and absence of the biological control agent T. viride. These data, shown in Fig. 1, record the observed development of dampingoff of radish seedlings, caused by R. solani, which attacks hosts in the early stages of development, plants rapidly becoming resistant to parasitism as they mature (19). Experimental data were obtained for five replicate microcosms with (Fig. 1b ) and without (Fig. 1a ) addition of the antagonistic fungus T. viride. We denote these two treatments as T+ and T, respectively. For more detailed descriptions of the microcosm experiments see refs. 4 and 14. In summary, each replicate curve shows the evolution of numbers of diseased seedlings in a population of 50 seedlings grown in clear plastic boxes. R. solani (AG 2.1, R5, IMI 385769) was added in the form of 10 mycelial discs placed randomly in each box. In the case of the curves of Fig. 1b , T. viride was added next to each host plant, in the form of a single colonized poppy seed. Host genotype and density, water availability, light, and temperature were strictly and identically controlled in each replicate. The boxes were opened at various times, and cumulative numbers of dampedoff seedlings were recorded.
Stochastic SEI Model with TimeDecaying Susceptibility
The model is an extension of the model of Gibson et al. (4). We denote the fixed population size by N and let S(t), E(t), and I(t) denote the number of susceptible, exposed, and infected individuals, respectively, at time t, so that N = S(t) + E(t) + I(t). At time t, the occurrence of a new exposure (i.e., the transition of an individual from S to E) is governed by the probabilistic rule, Here s(t) denotes the timevarying susceptibility of the host population as in Eq. 1 . The inclusion of this term is necessary biologically since as radish seedlings grow they become increasingly resistant to infection by R. solani. Such quenching is typical of many host–parasite systems. A key feature of the model is the presence of both primary and secondary infection processes in the model. These are controlled by the rate parameters r_{p} and r_{s} , respectively. The mechanism for secondary infection is based on a simple massaction principle of contacts between the I(t) infectious plants and the remaining S(t) susceptible plants. We assume that the timevarying susceptibility takes the form of an exponential decay s(t) = e^{at} . Here, we consider a fixed latent period of length μ, so that a plant that is exposed at time t will become infectious at time t + μ. The model, therefore, has four free parameters: the primary and secondary infection rates, r_{p} and r_{s} ; the latent period, μ; and the exponential decay rate of susceptibility, a. On setting μ = 0, we recover the SI model with decay of susceptibility considered by Gibson et al. (4).
Suppose now that we obtain data from a single replicate y = {(t_{i} , n_{i} ), 0 ≤ i ≤ k}, where n_{i} = I(t_{i} ) and (t _{0}, n _{0}) = (0, 0). This represents the situation where the plants scored as positives are those that are infectious. Now, given the model parameters r_{p} , r_{s} , a, and μ, the only unobserved aspects of the process are the precise times of the E → I transitions, denoted by z = (z_{1} , z_{2} ,... z_{nk}), z_{1} < z_{2} <... z_{nk} , that occur in the observation period [0, t_{k} ], the exposure times then being determined from z after subtraction of the fixed latent period μ. Working in a Bayesian statistical framework, for any parameter prior distribution π(r_{p} , r_{s} , a, μ), given the observations y , we can estimate (r_{p} , r_{s} , a, μ) by investigating the joint posterior density of the parameters and missing data components The term L(r_{p} , r_{s} , a, μ z ) is the parameter likelihood (given z ) and is proportional to the probability density function of z specified by parameters (r_{p} , r_{s} , a, μ). Here z ranges over those sets of infection times that are consistent with y , i.e., for which t_{j} ≤ z_{nj+1} ≤ z_{nj+2}... ≤ z_{nj+1} ≤ t_{j+1}, 0 ≤ j ≤ k  1. Inferences on individual parameters are made via their respective posterior marginal distributions (defined by appropriate integrals of the joint distribution specified by 3). When y represents several replicates considered jointly, the corresponding likelihood is obtained as the product of the singlereplicate likelihoods. Note that there is an implicit assumption of independence between replicates (so that the joint probability density function of the corresponding z s is the product of the densities for the individual replicates).
Given z , the corresponding unseen exposure times are given by x = (x _{1}, x _{2},... x_{nk} ), where x_{i} = z_{i}  μ and the parameter likelihood L(r_{p} , r_{s} , a, μ z ) can be derived as Eq. 3 cannot be integrated directly to obtain marginal parameter densities. However, (see e.g., refs. 17 and 18) we can simulate samples from the joint density π(r_{p} , r_{s} , a, μ, zy ) by using MCMC methods. Specifically, we investigate π(r_{p} , r_{s} , a, μ, zy ) by using a standard Metropolis sampler that applies updates to the parameters and the missing infection times sequentially. A single iteration of the chain comprises modifications to each parameter in turn, followed by modifications to the unknown infection times. To update r_{p} from its current value , we first propose a value drawn uniformly from . The window size, Δ _{p} , is a parameter that influences how efficiently the Markov chain explores its stationary distribution and is chosen by experimentation. The new value is then accepted with probability Otherwise it is rejected . Updates to r_{s} , a, and μ are made in an analogous fashion, with window sizes for the proposal distributions given by algorithm parameters Δ _{s} , Δ _{a} , and D_{m} , respectively.
Updates to the set of infection times z ^{i} are achieved by proposing and accepting/rejecting changes to the times of infections that occur between a particular pair of successive observation times t_{j} and t_{j+1} , for a given replicate. New times for infections occurring in (t_{j} , t_{j+1} ) are proposed uniformly from the set to give a candidate set of infection times for the replicate, z ′. The new set of times is accepted with probability Since replicates are assumed to be independent, it is necessary only to consider the contribution of the replicate to which the proposed change applies when computing the acceptance ratio. This procedure is applied to each pair of successive observation times for every replicate, to generate a new set of infection times for the set of replicates z ^{i} ^{+1}.
When the data record the values of E(t) + I(t), then the MCMC procedure used is similar to the above except that the missing events to be updated by the Markov chain are the precise exposure times x_{i} . The latter must satisfy t_{j} ≤ x_{nj+1} ≤ x_{nj+2}... ≤ x_{nj+1} ≤ t_{j+1} . Furthermore, a likelihood that differs slightly from Eq. 6 must be used, since the model assumes that all exposures occurring before t_{k} are observed [instead of t_{k}  μ for the case where data record I(t)]. Therefore for a given set of exposure times, x , the likelihood is given by Throughout this article, we work with noninformative prior densities for the parameters. Specifically, we consider that a priori the model parameters are independently distributed so that π(r_{p} , r_{s} , a, μ) = π(r_{p} )π(r_{s} )π(a)π(μ), where the prior for each parameter is uniform over bounded regions (cf. ref. 17). An alternative approach would be to use vague exponential or gamma priors (e.g., ref. 18). For the case where the data record E(t) + I(t) the likelihood fails to integrate over the space {(r_{p} , r_{s} , a, μ)  r_{p} , r_{s} , a, μ ≥ 0} so that we cannot use an improper uniform prior over this unbounded parameter space.
The Markov chain described above mixes well, and the acceptance rates for proposed changes to infection/exposure times are typically of the order of 90% or more, with acceptance rates of the order of 50% for changes to parameters. We found that simulation runs consisting of 501,000 iterations of the chain, with the first 1,000 iterations discarded to allow burnin produced consistent estimates of posterior densities of parameters over replicate simulation runs. The results illustrate that even simple MCMC approaches are powerful tools in this context and give acceptable performance, with comparatively little finetuning.
Results
We first compare parameter estimates and inferences about biological control yielded by the Bayesian treatment of the SI model with the results of an analysis using profile likelihoods by Gibson et al. (4). Fig. 2 shows histogram estimates of posterior densities π(r_{p}  y ), π(r_{s}  y ), and π(a y ), as estimated by the MCMC algorithm, where y represents the replicate data shown in Fig. 1, conditioned on the latent period μ = 0 (corresponding to the SI model). Upper bounds for r_{p} , r_{s} , and a were set at 0.08, 0.04, and 0.4 respectively. It is clear that the posterior marginal modes for the three parameters are close to their maximumlikelihood values (r_{p} = 0.2265, 0.00074; r_{s} = 0.0118, 0.0102, and a = 0.167, 0.127, for T and T+, respectively) given by Gibson et al. (4). Moreover, the greatest disparity between the treatments in the posterior densities occurs for the rate of primary infection rate, r_{p} . This confirms the conclusion of Gibson et al. (4) that the biological control agent affects epidemics by reducing the rate of primary infection.
We repeat the analysis, but allow our prior density for μ to be uniform rather than concentrated on μ = 0 and assume that the observations record I(t) (so that individuals in E and S are indistinguishable in the experiment). This implies that μ ≤ 1, since some replicates show nonzero incidence of infection after 1 day for both treatments, so that our uniform prior for μ is restricted to the interval 0 ≤ μ ≤ 1. Fig. 3 shows the posterior marginal densities π(r_{p}  y ), π(r_{s}  y ), π(a y ), and π(μ y ) for the two treatments. From these we see that, although the densities for the rate parameters are generally concentrated on higher values than in Fig. 2, the same qualitative conclusion emerges. There is clearly the most divergence between the posterior densities for r_{p} , supporting the notion that primary infection is influenced by the presence of T. viride.
Finally, we carry out the analysis for the situation where the observations represent the sum of individuals in I(t) + E(t). In this case, because plants that are infected but not infectious are scored as positive with respect to R. solani, the latent period need not be bounded above by 1. Inspection of the posterior marginal densities in Fig. 4 shows that, again, the largest divergence in the marginal densities between the two treatments is seen for the parameter r_{p} . Compared with the analysis for observations of I(t) alone, π(r_{s}  y ) is suggestive of higher values of the secondary infection parameter r_{s} . Moreover, π(μ y ) for both treatments is indicative of a nonzero latent period, with the posterior marginal modal values of μ exceeding 2 days for both treatments. There is also evidence of bimodality in π(μ y ) for the case where T. viride is added. Bimodality in a marginal posterior density can arise for several reasons and can be influenced by the choice of prior distributions. The interpretation of the phenomenon is not straightforward and should be the subject of further investigation.
Discussion
Several qualitative conclusions emerge from the analyses. A consistent finding across all of the models considered is that there is an obvious difference between the treatments in posterior marginal density of r_{p} , which is concentrated around lower values when T. viride is added (T+). Differences between marginal densities for other parameters are much less marked. This confirms the main conclusion of Gibson et al. (4) that T. viride operates by inhibiting the primary infection of plants by inoculum in the soil.
Here, for simplicity and ease of analysis, we assumed a fixed latent period. Clearly other choices for the distribution of latent periods may be made. The exponential distribution, commonly used in epidemic modeling, may be criticized on the grounds of biological realism since it is restrictive in the shape of distribution it can represent. A more flexible and realistic alternative is the gamma distribution used by O'Neill and Becker (20). Our principal aim here was to assess the robustness of the main conclusion of Gibson et al. (4) (namely that the biological control agent influences primary infection) to the inclusion of the latent state in the model, and we believe that the simplicity afforded by the fixed period makes its choice sensible.
It is important to consider the parameter likelihood for the model in Fig. 4 in more detail to avoid overinterpreting the densities of Fig. 4. Intuitively, we can anticipate a bias toward larger values of μ in π(μ y ); as μ increases, the integral of S(t)I(t) over the observation period decreases. Now, this quantity represents the total contact between infectious and susceptible individuals and determines the amount of secondary infection. Therefore, as μ increases, the dependence of L(r_{p} , r_{s} , a, μ y ) on r_{s} will decrease until μ exceeds the length of the observation window, so that I(t) = 0 throughout. At this point L(r_{p} , r_{s} , a, μ y ) no longer depends on r_{s} and μ. An immediate consequence of this is that unbounded uniform priors for r_{s} and μ cannot be used since the likelihood would not be integrable. Moreover, it suggests that a high value of π(μ y ) (obtained by integration of the likelihood) could be caused by the local flatness of the likelihood with respect to the other parameters and does not necessarily imply the likelihood itself taking high values. Consequently, it is possible for the properties of the marginal posterior distributions π(μ y ) and the profile likelihood L_{p} (μ y ) (obtained by maximizing the likelihood with respect to the other three parameters) to be very different. We suggest, therefore, that, although the analysis of Fig. 4 supports the earlier qualitative conclusions regarding the effect of T. viride on r_{p} , it would not be reasonable to draw conclusions regarding the value of the latent period μ without further investigation of the sampling properties of L(r_{p} , r_{s} , a, μ y ) and without further experimental data. However, independent empirical observations based on measuring the time of transfer between pairs of donor and recipient plants (21) suggest that a period of this order (i.e., 1–2 days) is realistic.
These analyses serve to illustrate how stochastic modeling, allied with powerful computational methods of inference, provide a powerful tool to analyze experimental data that typically offer only an incomplete account of the processes studied. Here, the regularity and completeness of the sampling procedures mean that the problems presented by the missing data components are mild, and that effective MCMC samplers can be designed in a very simple manner by using nothing more complex than Metropolis–Hastings methods. The approach used to construct the chain is not the most sophisticated one available, and the algorithm parameters have not been selected to optimize its convergence properties. Indeed, there are many enhancements that could be made to the algorithms presented here that would improve their efficiency and enable them to be applied in more testing scenarios.
We do not claim to be testing the limits of the computational approach in any way. Rather, the constraints that must be placed on the strength of the conclusions from this study derive from the nature of the data and the availability of prior knowledge of the process. For example, our ability to demonstrate the existence of a latent state, or characterize its duration, is undermined to an extent by the uncertainty in how the observations should be interpreted. In the light of this, we believe that significant further understanding of host–pathogen systems will come about because of advances in experimental methodology, supported by, but not replaced by, mathematical and statistical analyses of the kind reported here.
Acknowledgments
We gratefully acknowledge funding from the Biotechnology and Biological Sciences Research Council.
Footnotes

↵ † To whom correspondence should be addressed. Email: g.j.gibson{at}ma.hw.ac.uk.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviation: MCMC, Markov chain Monte Carlo.
 Copyright © 2004, The National Academy of Sciences
References

↵
Anderson, R. M. & May, R. M. (1991) Infectious Diseases of Humans: Dynamics and Control (Oxford Univ. Press, Oxford).

↵
Gilligan, C. A. (2002) Adv. Bot. Res. 38 , 164.

↵
Andersson, H. & Britton, T. (2000) Stochastic Epidemic Models and Their Statistical Analysis (Springer, New York).
 ↵

Renshaw, E. (1991) Modelling Biological Populations in Space and Time (Cambridge Univ. Press, Cambridge, U.K.).

↵
Truscott, J. E. & Gilligan, C. A. (2003) Proc. Natl. Acad. Sci. USA 100 , 90679072. pmid:12857951
 ↵
 ↵
 ↵

↵
Gibson, G. J. (1997) J. R. Stat. Soc. C Appl. Stat. 46 , 215233.

Gibson, G. J. (1997) Phytopathology 87 , 13911346.
 ↵

↵
Kleczkowski, A., Bailey, D. J. & Gilligan, C. A. (1996) Proc. R. Soc. London Ser. B 263 , 777783.

↵
Johnson, N. L. & Kotz, S. (1970) Distributions in Statistics: Continuous Univariate Distributions 2 (Wiley, New York).

↵
BarndorffNielsen, O. E. & Cox, D. R. (1994) Inference and Asymptotics (Chapman & Hall, London).

↵
Gibson, G. J. & Renshaw, E. (1998) J. IMA Math. Appl. Med. Biol. 15 , 1940.
 ↵

↵
Deacon, J. W. (1997) Modern Mycology (Blackwell, Oxford).
 ↵

↵
Kleczkowski, A., Gilligan, C. A. & Bailey, D. J. (1997) Proc. R. Soc. London Ser. B 264 , 979984.