Response of a deterministic epidemiological system to a stochastically varying environment
See allHide authors and affiliations

Edited by Simon A. Levin, Princeton University, Princeton, NJ (received for review October 16, 2002)
Abstract
Fluctuations in the natural environment introduce variability into the biological systems that exist within them. In this paper, we develop a model for the influence of random fluctuations in the environment on a simple epidemiological system. The model describes the infection of a dynamic host population by an environmentally sensitive pathogen and is based on the infection of sugar beet plants by the endoparasitic slimemold vector Polymyxa betae. The infection process is switched on only when the temperature is above a critical value. We discuss some of the problems inherent in modeling such a system and analyze the resulting model by using asymptotic techniques to generate closedform solutions for the mean and variance of the net amount of new inoculum produced within a season. In this way, the variance of temperature profile can be linked with that of the inoculum produced in a season and hence the risk of disease. We also examine the connection between the model developed in this paper and discrete Markovchain models for weather.
1. Introduction
Most analytical work on stochasticity in epidemics has focused on demographic variability (1, 2). However, in many instances, environmental variation can be a critical influence on the development of an epidemic. A wide range of economically and environmentally important fungi and invertebrates have a strong sensitivity to environmental factors. Many species of fungi, including the dampingoff fungi (Rhizoctonia solani), mildews (Blumeria graminis), and rusts (Puccinia spp.), have threshold temperatures and humidity levels for germination to occur. Some pathogenic nematodes that cause severe disease in staple crops exhibit a critical sensitivity to soil moisture content, becoming inactive at low levels (3). A variety of insect pests and parasites are strongly influenced by environmental switching, most notably by diapause, a suspension of development often triggered by changes in temperature, light levels, or humidity (4). Organisms with the sensitivities described above may be acting directly as pathogens on the hosts or as biocontrol agents, predating pests on a host (5, 6). Environmental variables such as temperature and moisture levels have both a predictable mean trend over time and a short timescale random component. Both aspects of this variation can be transmitted to the disease process through the sensitivity of the organisms involved. Hence the variability of the environment is fed through to the state of the epidemic.
In this paper, we derive solutions for the evolution of the mean and variance for newly generated infectious material for a broad class of epidemics under the influence of stochastic switching. Our model is based on a dynamically changing host population as susceptibles are produced or become infected or resistant; inoculum germinates in response to an environmental threshold and amplifies after infection of the host plant. This analysis depends, in part, on the effective correlation time of the stochastically switching system. This correlation time is the mean time over which the switching process remains in a single state. We explicitly link our results with Markovchain weathergenerating models, enabling us to connect the statistics of environmental variation with those of resulting epidemic. The results presented here have broad applicability to a range of pests and diseases.
Our work on this problem arises from modeling the amplification, transmission, and invasion of rhizomania disease in the United Kingdom (7–10). Rhizomania is a highly persistent soilborne disease of sugar beet with a serious economic impact for growers. The disease is caused by the beet necrotic yellow vein virus, which is transmitted via the soilborne endoparasitic slimemold vector, Polymyxa betae. A critical feature of the amplification of the inoculum in the host plant is the ability of its fungal vector to attack the host plant's roots. This process is characterized by a cutoff temperature of ≈12°C, below which the fungus is unable to enter the roots (11, 12). As a result, the extent of disease in the host is very sensitive to climate and sowing date. A stochastically varying temperature can cause the root system to move back and forth through the critical temperature, thus switching the infectivity of the vector on and off. In this way, variability in temperature is transmitted directly into the hostvector system and hence to the quantity of new inoculum produced by the host. This mechanism is discussed in more detail in Section 2.2.
In Section 2.1, we introduce the deterministic system in detail and discuss the effect of varying temperature on its evolution. We then derive a model of the full system as a stochastic differential equation (SDE) and discuss some of its assumptions and limitations in Section 2.3. We analyze the resulting SDE in Section 3 by using asymptotic techniques to get closedform solutions for the mean and variance of the solutions in terms of the properties of the driving environmental variable. In Section 4, we show how these results can be linked with a specific Markovchain weathergeneration model. In Section 5, we compare the results of the analysis with numerical simulations of the system and examine how the analysis can lead to testable hypotheses about the behavior of the system.
2. The Model
2.1. The Deterministic System. The biological system addressed is a susceptibleinfectious epidemic driven by primary infection from a reservoir of inoculum (X) in which the force of infection is switched on and off above and below a critical temperature. The model presented below is motivated by work done on a sugar beet–P. betae model described and analyzed in refs. 7–10. The model comprises the following system of ordinary differential equations: [1][2][3] with initial conditions, n(0) = s(0) = n_{0}, i(0) = 0. The variables n, s, and i represent the total population, the susceptible population at risk from the disease, and the inoculum generated by the infected population, respectively. The function, f(n, t), represents the rate of introduction of new hosts at any time, t, and encompasses very many commonly used growth functions. For numerical work, we use a monomolecular function, f(n, t) = r(1  n), which is both representative and convenient for calculation. Clearly, Eq. 1 can be solved independently of the other equations, and hence we can replace the function, f, in Eq. 2 with n_{t}(t), the derivative of the total population. We assume that newly generated inoculum is released as fully active for the next season. In the context of rhizomania, the variables n, s, and i represent the fraction of the total root length in a particular state. The model contains only four parameters: m, the rate at which susceptibles become resistant to infection; X, the initial inoculum loading; Λ_{m}, the force of infection, and Q, the amplification factor for new inoculum. Temperature dependence enters through the parameter, Λ_{m}, as follows: [4] That is, for temperature less than the critical value, T_{c}, infection is effectively turned off, whereas for temperature above this value, the force of infection is constant.
2.2. The Fluctuating Environment. In the present context, we focus on temperature as the critical aspect of the environment. However, more generally, a range of possible factors could introduce stochastic effects into the system, such as rain, sunlight, and animal vector movement. These are discussed in more detail in Sections 4 and 6. During the growing season, the mean temperature rises considerably and passes through the critical temperature of the parasite. The obvious deterministic approach to this changing environment is to consider the parasite to be “switched off” until the mean temperature reaches T_{c} and “switched on” afterwards. However, if we consider temperature to be a stochastically varying quantity, as experienced by the parasite, then the system will by switched on and off randomly while the variation in temperature spans T_{c}. Given that temperature can be taken to vary about the mean between definable extremes, the rising temperature profile can be divided into three temporal phases (Fig. 1) with respect to the parasite.

Phase 1, in which temperature is always below T_{c} and the infection process is always quiescent (t_{0} to t_{1} in Fig. 1).

Phase 2, in which temperature varies across the critical temperature and the infection process switches between quiescent and active (t_{1} to t_{2} in Fig. 1).

Phase 3, in which the temperature is constantly above T_{c} and the infection process is always active (t_{2} to t_{3} in Fig. 1).
Clearly the dynamics of Phase 1 are straightforward, assuming deterministic initial conditions. Phase 3 is once again deterministic but must describe the evolution of the mean and variance of the variables as they stand at the end of Phase 2. Phase 2, however, requires a new approach to integrate the stochasticity into the dynamics.
2.3. Stochastic Differential Equation (SDE) Model. An SDE is used to model the continuoustime stochastic process. The SDE comprises the deterministic mean behavior to which a noise term is added in the form of an infinitesimal Wiener process. The deterministic solution of Eqs. 1–3 is easily calculated in terms of the function, f(n, t). It is now necessary to formulate the infection switching in terms of the Wiener process. Problems arise from the fact that the switching process is binary; i.e., it is either “on” or “off.” In trying to formulate a continuoustime representation of switching based on a Wiener process, one must take a limit as dt → 0. As this limit is approached, the highfrequency elements in the Wiener process spectrum dominate, causing the variance and hence stochastic terms in the equation to vanish. Viewed from a physical or biological standpoint, the nature of the problem is clear. In the context of the sugarbeet rhizomania system, switching at high frequencies would require soil to heat up and the behavior of the parasite to change almost instantaneously. A more accurate description of the process would require either that the sharp step from “off” to “on” at T_{c} be made into a smooth function or that the physical system be given some inertia in changing state. Of these two, the second can be approximated by assuming that there exists an effective correlation period over which the state of the system will not change. The system samples only the random variable between these periods. In Section 4, this approach is compared directly to discrete firstorder Markovchain models for generating stochastic data, as used in weather modeling (13, 14). The strong correspondence between the SDE formulation and such models allows us to directly include such environmental influences into SDE descriptions.
To derive a continuoustime stochastic differential equation representation of the system described above, we first calculate the mean and variance of the process. Consider the system described by Eqs. 1–3 at time, t,instate(n, s, i). Let p_{+}(t) be the probability that the temperature is above T_{c} at time t, and let ΔT, the effective correlation time, define the discrete time interval at which the system checks the temperature. We can interpret ΔT as a property of either the deterministic system or the driving stochastic variable. Within the deterministic system, it can represent the time for the force of infection to respond to a change in temperature. In the case of the vector for rhizomania, this could be the characteristic time scale for soil heating. From the viewpoint of temperature, it can represent a correlation time over which the temperature remains unchanged. This possibility is discussed in detail in Section 4 with regard to discrete Markovchain models.
The effective correlation time is assumed to be small compared to the time scale of evolution of the deterministic system. Considering just the changes over ΔT of the stochastically varying elements, s and i, with probability, 1  p_{+}(t), and with probability, p_{+}(t). This gives a mean of [5] and a variance of [6] We want these variances to be matched by the SDE, where dW(t) represents the infinitesimal Wiener deviate at time, t. Only one infinitesimal Wiener process is required, because there is only a single stochastic process operating. Considering again the period, ΔT, and assuming a slow evolution for the deterministic processes, we arrive at a mean of and using the identity, (15). By comparison with Eqs. 5 and 6, we have Therefore, the SDE describing the dynamics of Phase 2 is [7] [8] [9] where . We have taken advantage of the degree of freedom in the signs of B_{s} and B_{i} to allow for material removed from s to be added to i.
3. Analysis
By using Eqs. 1–3 for Phases 1 and 3 and Eqs. 7–9 for Phase 2, it is possible to describe the evolution of a single realization of the system through all three phases in the form of a stochastic integral (15). This solution is, in general, too complex to yield any practically useful information about the statistics of the system, such as mean and variance. To simplify the solution, we consider an asymptotic limit in which the stochasticity is small, thereby allowing a solution through asymptotic techniques. Nondimensionalizing with respect to time and grouping parameters gives [10] [11] [12] for the SDE and [13] [14] [15] for the deterministic system, where and using the identity, . The parameter, λ′ = 0 in phase 1 and λ′ = ε in Phase 3.
Because the stochastic effects of Phase 2 enter through the force of the infection process and the parameter, λ_{m}, the smallnoise limit is analyzed by taking the limit in which the parameter, ε → 0. We look for solutions of Eqs. 10–12 as a series in ε, i.e., Substituting these series into the SDE gives for the first two terms [16] [17] where n_{1}, s_{1} represent values at the start of Phase 2, and τ is measured from the start of Phase 2. We are interested only in the first term of ĩ, which is ĩ_{b} [18] It is easy to calculate the mean and variance of these quantities by using the identity, (15), giving [19] [20] [21] The O(ε^{2}) for these terms reflects the fact that the infection process and hence its products are O(ε). As yet, we have made no assumptions about the form of p_{+}(τ) and hence ϕ(τ′). The simplest assumption is that p_{+}(τ) is constant through Phase 2, in which case ϕ can be removed from the integrals in Eqs. 19–21. The simplest approximation to the timedependent behavior illustrated in Fig. 1 is a linear rate of change for p_{+}, which captures the smooth transition in p_{+} from 0 at the start of Phase 2 to Phase 1 at the end. [22] and [23] where Δτ_{2} is the duration of Phase 2.
We can make a further simplification to the above expressions for variance when Phase 2 is short in comparison to the rate of change of s_{a}. We can approximate [24] where τ is measured from the start of Phase 2. The derivative of s_{a} at the start of the phase can be expressed in terms of n_{τ} and s_{a} through Eq. 14. With this approximation, the integrals in Eqs. 19–21 become analytically tractable. For example, the expression for the variance of ι̃ (Eq. 20) becomes which is simply the integral of a polynomial in τ.
The evolution of the mean and variance of s and i in the deterministic third phase is addressed in the Supporting Appendix, which is published as supporting information on the PNAS web site, www.pnas.org. Essentially, the development of these quantities in Phase 3 can be expressed as functions of the mean, variance, and covariance of s and i at the end of Phase 2.
4. Comparison of Results with Discrete FirstOrder MarkovChain Models
Discrete firstorder Markovchain models are widely used in model weathergenerating systems as a simple and accurate method of generating daily rainfall patterns (13). The seminal example is a model of rainfall in Tel Aviv (14). We will use this model to illustrate the connection between firstorder Markovchain models and that derived in Section 3. Using the relationships between SDE and Markovchain models derived in this section, we can correctly parameterize an SDE to simulate the stochastic influence of a discrete Markovchain environmental variable.
The rainfall model assumes the weather to be in one of two states on any day: raining or dry. Hence the day is the fundamental unit of time, t_{0}, over which the state of the system remains unchanged. The state on any day depends only on that of the previous day through the conditional probabilities,
It can be shown that over a period of n days the mean and variance of the cumulative length of time of rainy days, r, are given by [25] [26] as n gets large, where π is the absolute probability of rain on any day,
These expressions can be directly compared to those for the mean and variance of i in Section 3. The quantity, i, represents the accumulation of inoculum generated by the stochastically switching infection process and is entirely analogous to the accumulation of rainy days in the Markovchain weather model. If we consider the dynamics of the underlying model to be static through Phase 2, we can express the mean and variance of i as [27] [28] in dimensional form. The parameter grouping λ_{m}Xs_{a} is a rate appearing in both expressions, converting time into the units of i. Hence we can make a direct comparison between the time spent in an infectious state in Phase 2 and the time spent raining in the Markovchain model, and we can associate the parameters of the two models:
SDE model ⇔ Markovchain model [29] [30] Both p_{+} and π are absolute probabilities, whereas the quantities in Eq. 30 are both functions of correlation parameters.
We can use this basis, therefore, to construct SDE models for systems governed by discrete Markovchain stochastic processes over relatively long periods of time. Simulations show that the relative error in the variance as predicted by the SDE model decreases with increasing time. This decrease would of course be accompanied by a corresponding increase in error from the asymptotic approximation.
5. Numerical Results
The theory developed in Section 3 was tested against numerical simulations by using a fourthorder Runge–Kutta algorithm to integrate Eqs. 1–3 over the three phases of the system's evolution. During the stochastic phase, the temperature of the system for any correlation period was chosen randomly according to the probability, p_{+}. For all the following simulations, the temperature phases were defined as The parameters used were those presented in Table 1 unless stated otherwise. A monomolecular growth function was used for simulations for algebraic tractability.
Table 2 below compares values for variance of inoculum production at the end of Phase 3 between the numerical solution and the theoretical predictions from Section 3 for different correlation periods, ΔT. Predictions for the mean were found to be very accurate for all values of the period. For variance, the accuracy increases as the correlation period gets smaller. These results also illustrate the problem discussed in Section 2.3: as the correlation period is decreased toward zero, the variance in the system also vanishes.
Fig. 2 shows realizations of the stochastic process produced by the simulator. The thicker lines represent the mean and two standard deviations on either side as calculated from theory (Eq. 20). The region between these lines represents approximately a 98% confidence interval for the realizations, assuming a normal distribution. This assumption can be seen to be approximate from the fact that the lower curve goes negative at the start of the second phase. Nevertheless, they function well as bounds on the behavior of the realizations. Fig. 3 shows more clearly the variance of generated inoculum over the three phases of evolution. In particular, it shows the continuing change in variance in the third deterministic phase.
Simulation and theory are in close agreement through the entire range of p_{+}. The shape of the curve reflects the dependence of the variance of inoculum on ϕ^{2} (Eq. 20), where , giving a characteristic parabolic shape. A linear rate of change in probability (Eq. 22) is a more realistic form for p_{+} to take, as discussed in Section 3. Fig. 4 was generated by using this linear form for p_{+} and the linear approximation for s_{a} (Eq. 24). It shows the error in the variance of the infected roots against the length of Phase 2. Clearly, the variance error increases with the length of Phase 2, showing that the improvement in the accuracy of the SDE described in Section 3 is outweighed by the accumulating errors from the asymptotic approximation and the assumption of linearity in s_{a}.
The inclusion of environmental variability allows us to make useful statistical statements about how disease will develop in the host. Such predictions can be tested against experiments in the field. As an example, we consider the probability, P(i > i_{c}), that the inoculum generated in a season exceeds some critical level. i_{c}, such as that resulting in symptoms in the host or detection in the environment. In Fig. 5, these are calculated as a function of growth rate and force of infection, respectively, from the analysis in Section 3. The probability of exceeding the threshold as a function of force of infection shows an error function shape, reaching P = 0.5 as the mean of i = i_{c}. For a deterministic representation, the probability will switch from 0 to 1 at this point. Fig. 5B shows P(i > i_{c}) as a function of the growth rate, r. For low values of r, the slow generation of new roots limits the quantity of inoculum generated. For high values, the majority of roots have already entered the resistant class before T_{c} is reached and the disease process can develop strongly. Hence P(i > i_{c}) has a local maximum as a function of r. For these particular parameter values, the deterministic value for i never exceeds i_{c}. The stochastic model, however, indicates that there is a quite high probability (P > 0.4) of breaching the threshold for certain values of r.
6. Conclusion
The development of a model for this system falls into two distinct parts: first, the development of a conceptual model for the real system, and second, the development of a mathematical description of the conceptual model.
Fokker–Planck equations and stochastic differential equations are two different ways of representing the same underlying continuoustime stochastic process (15, 16). Stochastic differential equations describe the evolution of the random variables representing the state of the system, whereas Fokker–Planck equations describe the evolution of the probability distribution for the state of the system. From an analytic standpoint, they have a wellbehaved deterministic limit as the effect of stochasticity is reduced. For systems such as the present one, in which the deterministic solutions are available, there is the possibility of generating asymptotic series around the small noise limit. In the analogous limit for Fokker–Planck equations, no such convenient form is obtained.
Results from Section 3 show that it is possible to construct an SDE model for the stochasticswitching model developed in Section 2.3. Moreover, a closedform solution to this system of equations can be found in the asymptotic limit and hence an analytical approximation for the stochastic process. This approach has broad applicability to a range of epidemics and growth processes subject to environmentally controlled stochastic switching. In the present example, analytic tractability rests on the simplicity of the underlying epidemic model. For more complex systems, numerical techniques would be required.
Numerical results in Section 5 show that the solutions for the statistics of the inoculum, i, calculated in Section 3 and the Supporting Appendix describe those of the simulation very well. Predictions for the mean value of the inoculum are very accurate, whereas that for the variance are typically <5% out for the parameter values used. This accuracy is maintained across the whole range of p_{+}. There is an implicit assumption in the construction of the SDE that the Wiener process can be used to describe what is the essentially binomial process of switching. One problem can be seen in Eq. 12, where dW(t) can be negative and large enough to allow di to be negative. Clearly, in the real system, di ≥ 0. Hence, although the statistics predicted by the theory are reliable, individual realizations may be unrepresentative. As the stochastic period lengthens, the accuracy of the estimate of variance declines (Fig. 4), because of the accumulation of errors from the approximation.
In Section 4, we explicitly compare the behavior of the SDE model to that of a discrete firstorder Markovchain model. The concept of Markov chains underlies many stochastic processes, and discrete firstorder chains are widely used to model weather phenomena. Strong parallels can be drawn between such models and the SDE systems in terms of their statistics, in particular, the correlation period and the conditional probabilities of the Markov chain. That is, where t_{0} is the length of the discrete time unit. Using this relationship and the approach outlined in this paper, we can include the stochastic influence of such Markovchain environmental phenomena on evolving deterministic systems subject to stochastic switching.
For the specific model treated in this paper, we have considered stochastic forcing through fluctuating temperature in detail. This forcing is mediated through an extreme sensitivity to temperature in the parasite. A similar sensitivity to temperature would be found in plants vulnerable to frost damage and infection. In this case, the cutoff temperature is 0°C. Periods of rainfall or high humidity can also be triggers for infection and, as has been illustrated in Section 4, rainfall patterns are accurately simulated by Markovchain models and hence SDEs. Other candidates include insolation and the movement of animal vectors.
By explicitly including stochastic effects into a deterministic model, we address properties of the system that are not recoverable from a purely deterministic description. The probability of exceeding a threshold shown in Fig. 5 A and B is an example of such a property. Most biological experiments yield results of a statistical result of this nature based on many replicates, which can be tested against the probability distributions predicted by the model. Hence we can examine directly the nature and strength of the interaction between environmental phenomena and the biological systems they control.
The work presented here illustrates some of the difficulties inherent in representing the influence of stochastic processes on timecontinuous systems. Simply taking the limit as dt → 0 to geta continuous representation can fail to capture the variability of the system. It can “disappear” in the limit. This is the result of the fact that the variance of the Wiener process is proportional to time. Table 2 illustrates this point well. As the correlation period, ΔT, is made smaller, and variance shrinks proportionally. From a practical point of view, the problem lies with the infinitesimal properties of the Wiener distribution combined with the stepfunction nature of the switch. It could be said that the Weiner process is capable of finite changes in infinitesimally small periods of time. Because of inertia, no natural process behaves in this way. Hence for a switching process, it is necessary to include inertia artificially. In the present case, we included the effective correlation time to represent thermal inertia. An alternative approach is to replace the switching function with a smooth alternative. This would better represent the response of the system to an environmental influence. What form these functions should take and what effect their form will have on the transmission of variance into the system are the subjects of further work.
Acknowledgments
We gratefully acknowledge funding from the Biotechnology and Biological Sciences Research Council (J.E.T.) and the Royal Society and Leverhulme Trust (C.A.G.).
Footnotes

↵* To whom correspondence should be sent at the present address: Department of Infectious Disease Epidemiology, Faculty of Medicine, Imperial College, London W2 1PG, United Kingdom. Email: j.truscott{at}ic.ac.uk.

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

Abbreviation: SDE, stochastic differential equation.
 Received October 16, 2002.
 Accepted May 20, 2003.
 Copyright © 2003, The National Academy of Sciences
References
 ↵
Andersson, H. & Britton, T. (2000) Stochastic Epidemic Models and Their Statistical Analysis (Springer, Berlin).
 ↵
Renshaw, E. (1991) Modelling Biological Populations in Space and Time (Cambridge Univ. Press, Cambridge, U.K.).
 ↵
Grant, J. A. & Villani, M. G. (2003) Environ. Entomol. 32, 8087.
 ↵
Tauber, M. J., Tauber, C. A. & Masaki, S. (1986) Seasonal Adaptations of Insects (Oxford Univ. Press, Oxford).
 ↵
Fransen, J. J. (1990) in Whiteflies: Their Bionomics, Pest Status and Management, ed. Gerling, D. (Intercept, Andover, MA).
 ↵
 ↵

Truscott, J. E., Webb, C. R. & Gilligan, C. A. (2000) Bull. Math. Biol. 62, 366393.
 ↵
 ↵
Asher, M. J. C. (1993) in The Sugar Beet Crop: Science into Practice, eds. Cooke, D. & Scott, R. (Chapman & Hall, London), Chapter 9, pp. 313346.
 ↵
Blunt, S. J., Asher, M. J. C. & Gilligan, C. A. (1992) Plant Pathol. 41, 257267.
 ↵
Wilks, D. S. & Wilby, R. L. (1999) Prog. Phys. Geogr. 23, 329357.
 ↵
 ↵
Gardiner, C. W. (1985) Handbook of Stochastic Methods (Springer, Berlin).
 ↵
Risken, H. (1989) The Fokker–Planck Equation: Methods of Solution and Applications (Springer, Berlin).