key: cord-0000311-0pe8691t authors: Forgoston, Eric; Billings, Lora; Schwartz, Ira B. title: Accurate noise projection for reduced stochastic epidemic models date: 2009-10-29 journal: Chaos DOI: 10.1063/1.3247350 sha: 7ab1a7cdf14acc362661a984688a378665040313 doc_id: 311 cord_uid: 0pe8691t We consider a stochastic susceptible-exposed-infected-recovered (SEIR) epidemiological model. Through the use of a normal form coordinate transform, we are able to analytically derive the stochastic center manifold along with the associated, reduced set of stochastic evolution equations. The transformation correctly projects both the dynamics and the noise onto the center manifold. Therefore, the solution of this reduced stochastic dynamical system yields excellent agreement, both in amplitude and phase, with the solution of the original stochastic system for a temporal scale that is orders of magnitude longer than the typical relaxation time. This new method allows for improved time series prediction of the number of infectious cases when modeling the spread of disease in a population. Numerical solutions of the fluctuations of the SEIR model are considered in the infinite population limit using a Langevin equation approach, as well as in a finite population simulated as a Markov process. The interaction between deterministic and stochastic effects in population dynamics has played, and continues to play, an important role in the modeling of infectious diseases. The mechanistic modeling side of population dynamics is well known and established. 1, 2 These models typically are assumed to be useful for infinitely large, homogeneous populations, and arise from the mean field analysis of probabilistic models. On the other hand, when one considers finite populations, random interactions give rise to internal noise effects, which may introduce new dynamics. Stochastic effects are quite prominent in finite populations, which can range from ecological dynamics 3 to childhood epidemics in cities. 4, 5 For homogeneous populations with seasonal forcing, noise also comes into play in the prediction of large outbreaks. [6] [7] [8] Specifically, external random perturbations change the probabilistic prediction of epidemic outbreaks as well as its control. 9 When geometric structure is applied to the population, the interactions are modeled as a network. 10, 11 Many types of static networks which support epidemics have been considered. Some examples include small world networks, 12 hierarchical networks, 13 and transportation networks of patch models. 14 In addition, the fluctuation of epidemics on adaptive networks, where the wiring between nodes changes in response to the node information, has been examined. 15 In adaptive network models, even the mean field can be high dimensional, since nodes and links evolve in time and must be approximated as an additional set of ordinary differential equations. Another aspect of epidemic models, which is often of interest, involves the inclusion of a time delay. The delay term makes the analysis significantly more complicated. However, it is possible to approximate the delay by creating a cascade consisting of a large number of compartments. 16 For example, one could simulate the delay associated with a disease exposure time with several hundred "exposed" compartments. These model examples are just a few of the types of very high-dimensional models that are currently of interest. As a result of the high dimensionality, there is much computation involved, and the analysis is quite difficult. In particular, real-time computation is not currently possible. However, there are usually many time scales that are well separated ͑due typically to a large range in order of magnitude of the parameters͒ when considering such high-dimensional problems. In the presence of well-separated time scales, a model reduction method is needed to examine the dynamics on a lower-dimensional space. It is known that deterministic model reduction methods may not work well in the stochastic realm, which includes epidemic models. 17 The purpose of this article is to examine a method of nonlinear, stochastic projection so that the deterministic and stochastic dynamics interact correctly on the lower-dimensional manifold and predict correctly the dynamics when compared with the full system. Because the noise affects the timing of outbreaks, it is essential to produce a low-dimensional system that captures the correct timing of the outbreaks as well as the amplitude and phase of any recurrent behavior. We will demonstrate that our stochastic model reduction method properly captures the initial disease outbreak and continues to accurately predict the outbreaks for time scales, which are orders of magnitude longer than the typical relaxation time. Furthermore, in practice, real disease data include only the number of infectious individuals. Our method allows us to predict the number of unobserved exposed individuals based on the observed number of infectious individuals. For stochastic model reduction, there exist several potential methods for general problems. For a system with certain spectral requirements, the existence of a stochastic center manifold was proven in Ref. 18 . Nonrigorous stochastic normal form analysis ͑which leads to the stochastic center mani-fold͒ was performed in Refs. [19] [20] [21] [22] . Rigorous theoretical analysis of normal form coordinate transformations for stochastic center manifold reduction was developed in Refs. 23 and 24. Later, another method of stochastic normal form reduction was developed, 25 in which any anticipatory convolutions ͑integrals into the future of the noise processes͒ that appeared in the slow modes were removed. Since this latter analysis makes the construction of the stochastic normal form coordinate transform more transparent, we use this method to derive the reduced stochastic center manifold equation. Figure 1 shows a schematic demonstrating our approach to the problem. We consider a high-dimensional system along with its corresponding reduced low-dimensional system. In this article, two types of low-dimensional system are discussed: a reduced system based on deterministic center manifold analysis and a reduced system based on a stochastic normal form coordinate transform. Regardless of the type of low-dimensional system being considered, a common noise is injected into both the high-dimensional and lowdimensional models, and an analysis of the solutions found using the high and low-dimensional systems is performed. In this article, as a first study of a high-dimensional system, we consider the susceptible-exposed-infected-recovered ͑SEIR͒ epidemiological model with stochastic forcing. As previously mentioned, we could easily consider a SEIR-type model where the exposed class was modeled using hundreds of compartments. Since the analysis is similar, we consider the simpler standard SEIR model to demonstrate the power of our method. Section II provides a complete description of this model. Section III describes how to transform the deterministic SEIR system to a new system that satisfies the spectral requirements needed to apply the center manifold theory. After the theory is used to find the evolution equations that describe the dynamics on the center manifold, we show in Sec. IV how the reduced model that is found using the deterministic result incorrectly projects the noise onto the center manifold. Section V demonstrates the use of a stochastic normal form coordinate transform to correctly project the noise onto the stochastic center manifold. A discussion and the conclusions are contained, respectively, in Secs. VI and VII. We begin by describing the stochastic version of the SEIR model found in Ref. 26 . We assume that a given population may be divided into the following four classes that evolve in time: ͑1͒ Susceptible class s͑t͒ consists of those individuals who may contract the disease. ͑2͒ Exposed class e͑t͒ consists of those individuals who have been infected by the disease but are not yet infectious. ͑3͒ Infectious class i͑t͒ consists of those individuals who are capable of transmitting the disease to susceptible individuals. ͑4͒ Recovered class r͑t͒ consists of those individuals who are immune to the disease. Furthermore, we assume that the total population size, denoted as N, is constant and can be normalized to S͑t͒ + E͑t͒ + I͑t͒ + R͑t͒ = 1, where S͑t͒ = s͑t͒ / N, E͑t͒ = e͑t͒ / N, I͑t͒ = i͑t͒ / N, and R͑t͒ = r͑t͒ / N. Therefore, the population class variables S, E, I, and R represent fractions of the total population. The governing equations for the stochastic SEIR model are Schematic demonstrating the injection of a common noise into both the high-dimensional system and its associated low-dimensional system. where i is the standard deviation of the noise intensity D i = i 2 / 2. Each of the noise terms i describes a stochastic, Gaussian white force that is characterized by the correlation functions Additionally, represents a constant birth and death rate, ␤ is the contact rate, ␣ is the rate of infection, so that 1 / ␣ is the mean latency period, and ␥ is the rate of recovery, so that 1 / ␥ is the mean infectious period. Although the contact rate ␤ could be given by a time-dependent function ͑e.g., due to seasonal fluctuations͒, for simplicity, we assume ␤ to be constant. Throughout this article, we use the following parameter values: = 0.02͑yr͒ −1 , ␤ = 1575.0͑yr͒ −1 , ␣ =1/ 0.0279͑yr͒ −1 , and ␥ =1/ 0.01͑yr͒ −1 . Disease parameters correspond to typical measles values. 26, 27 Note that any other biologically meaningful parameters may be used as long as the basic reproductive rate R 0 = ␣␤ / ͓͑␣ + ͒͑␥ + ͔͒ Ͼ 1. The interpretation of R 0 is the number of secondary cases produced by a single infectious individual in a population of susceptibles in one infectious period. As a first approximation of stochastic effects, we have considered additive noise. This type of noise may result from migration into and away from the population being considered. 28 Since it is difficult to estimate fluctuating migration rates, 29 it is appropriate to treat migration as an arbitrary external noise source. Also, fluctuations in the birth rate manifest itself as additive noise. Furthermore, as we are not interested in extinction events in this article, it is not necessary to use multiplicative noise. In general, for the problem considered here, it is possible that a rare event in the tail of the noise distribution may cause one or more of the S, E, and I components of the solution to become negative. In this article, we will always assume that the noise is sufficiently small so that a solution remains positive for a long enough time to gather sufficient statistics. Even though it is difficult to accurately estimate the appropriate noise level from real data, our choices of noise intensity lie within the huge confidence intervals computed in Ref. 29 . The case for multiplicative noise will be considered in a separate paper. Although S + E + I + R = 1 in the deterministic system, one should note that the dynamics of the stochastic SEIR system will not necessarily have all of the components sum to unity. However, since the noise has zero mean, the total population will remain close to unity on average. Therefore, we assume that the dynamics is sufficiently described by Eqs. ͑1a͒-͑1c͒. It should be noted that even if E͑t͒ + I͑t͒ = 0 for some t, the noise allows for the reemergence of the epidemic. One way to reduce the dimension of a system of equations is through the use of deterministic center manifold theory. In general, a nonlinear vector field can be transformed so that the linear part ͑Jacobian͒ of the vector field has a block diagonal form where the first matrix block has eigenvalues with positive real part, the second matrix block has eigenvalues with negative real part, and the third matrix block has eigenvalues with zero real part. 30, 31 These blocks are associated, respectively, with the unstable eigenspace, the stable eigenspace, and the center eigenspace. If we suppose that there are no eigenvalues with positive real part, then orbits will rapidly decay to the center eigenspace. In order to make use of the center manifold theory, we must transform Eqs. ͑1a͒-͑1c͒ to a new system of equations that has the necessary spectral structure. The theory will allow us to find an invariant center manifold passing through the fixed point to which we can restrict the transformed system. Details regarding the transformation can be found in Sec. III A, and the computation of the center manifold can be found in Sec. III B. Our analysis begins by considering the governing equations for the stochastic SEIR model given by Eqs. ͑1a͒-͑1c͒. We neglect the i i ͑t͒ terms in Eqs. ͑1a͒-͑1c͒ so that we are considering the deterministic SEIR system. This deterministic system has two fixed points. The first fixed point is and corresponds to a disease free or extinct equilibrium state. The second fixed point corresponds to an endemic state and is given by To ease the analysis, we define a new set of variables, S, Ē , and Ī, as S͑t͒ = S͑t͒ − S 0 , Ē ͑t͒ = E͑t͒ − E 0 , and Ī͑t͒ = I͑t͒ − I 0 . These new variables are substituted into Eqs. ͑1a͒-͑1c͒. Then, treating as a small parameter, we rescale time by letting t = . We may then introduce the following rescaled parameters: ␣ = ␣ 0 / and ␥ = ␥ 0 / , where ␣ 0 and ␥ 0 are O͑1͒. The inclusion of the parameter as a new state variable means that the terms in our rescaled system which contain are now nonlinear terms. Furthermore, the system is augmented with the auxiliary equation d / d = 0. The addition of this auxiliary equation contributes an extra simple zero eigenvalue to the system and adds one new center direction that has trivial dynamics. The shifted and rescaled augmented system of equations is where the endemic fixed point is now located at the origin. The Jacobian of Eqs. ͑5a͒-͑5d͒ is computed to zeroth order in and is evaluated at the origin. Ignoring the components, the Jacobian has only two linearly independent eigenvectors. Therefore, the Jacobian is not diagonalizable. However, it is possible to transform Eqs. ͑5a͒-͑5c͒ to a block diagonal form with the eigenvalue structure that is needed to use the center manifold theory. We use a transformation matrix P consisting of the two linearly independent eigenvectors of the Jacobian along with a third vector chosen to be linearly independent. There are many choices for this third vector; our choice is predicated on keeping the vector as simple as possible. This transformation matrix is given as Using the fact that ͑S , Ē , Ī͒ T = P · ͑U , V , W͒ T , then the transformation matrix leads to the following definition of new variables: U, V, and W, The application of the transformation matrix to Eqs. ͑5a͒-͑5c͒ leads to the transformed evolution equations given by The Jacobian of Eqs. ͑8a͒-͑8d͒ to zeroth order in and evaluated at the origin is which shows that Eqs. ͑8a͒-͑8d͒ may be rewritten in the form where x = ͑U͒, y = ͑V , W͒, A is a constant matrix with eigenvalues that have negative real parts, B is a constant matrix with eigenvalues that have zero real parts, and f and g are nonlinear functions in x, y, and . In particular, Therefore, the system will rapidly collapse onto a lowerdimensional manifold given by center manifold theory. 32 Furthermore, we know that the center manifold is given by where h is an unknown function. Substitution of Eq. ͑14͒ into Eq. ͑8a͒ leads to the following center manifold condition: In general, it is not possible to solve the center manifold condition for the unknown function h͑V , W , ͒. Therefore, we perform the following Taylor series expansion of h͑V , W , ͒ in V, W, and : where h 0 , h 2 , h 3 , h ,¯are unknown coefficients that are found by substituting the Taylor series expansion into the center manifold condition and equating terms of the same order. By carrying out this procedure using a second-order Taylor series expansion of h, the center manifold equation is where ⑀ = ͉͑V , W , ͉͒ so that ⑀ provides a count of the number of V, W, and factors in any one term. Substitution of Eq. ͑17͒ into Eqs. ͑8b͒ and ͑8c͒ leads to the following reduced system of evolution equations, which describe the dynamics on the center manifold We now return to the stochastic SEIR system given by Eqs. ͑1a͒-͑1c͒. The shift in the fixed point to the origin will not have any effect on the noise terms, so that the stochastic version of the shifted equations is As Eqs. ͑19a͒-͑19c͒ are transformed using Eqs. ͑7a͒-͑7c͒, the ␣ = ␣ 0 / scaling, the ␥ = ␥ 0 / scaling, and the t = time scaling, the noise also is scaled so that the stochastic transformed equations are given by where 043110-5 Noise projection for epidemic models Chaos 19, 043110 ͑2009͒ The stochastic terms 4 , 5 , and 6 in Eqs. ͑20a͒-͑20c͒ are still additive Gaussian noise processes. However, Eqs. ͑21a͒-͑21c͒ show how the transformation has acted on the original stochastic terms 1 , 2 , and 3 to create new noise processes, which have a variance different from that of the original noise processes. Also note that we have suppressed the argument of 4 , 5 , and 6 in Eqs. ͑20a͒-͑20c͒. The time scaling means that these noise terms should be evaluated at . The system of equations given by Eqs. ͑20a͒ and ͑21c͒ is an exact transformation of the system of equations given by Eqs. ͑1a͒-͑1c͒. We numerically integrate the original stochastic system of the SEIR model ͓Eqs. ͑1a͒-͑1c͔͒ along with the transformed stochastic system ͓Eqs. ͑20a͒-͑20c͔͒ using a stochastic fourth-order Runge-Kutta scheme with a constant time step size. The original system is solved for S, E, and I, while the transformed system is solved for U, V, and W. In the latter case, once the values of U, V, and W are known, we compute the values of S, Ē , and Ī using the transformation given by Eqs. ͑7a͒-͑7c͒. We shift S, Ē , and Ī, respectively, by S 0 , E 0 , and I 0 to find the values of S, E, and I. Figure 2 compares the time series of the fraction of the population that is infected with a disease I, computed using the original stochastic system of equations of the SEIR model with the time series of I computed using the transformed stochastic system of equations of the SEIR model. Although the two time series shown in Fig. 2 generally agree very well, there is some discrepancy. This discrepancy is due to the fact that the noise processes 4 4 , 5 5 , and It is tempting to consider the reduced stochastic model found by substitution of Eq. ͑17͒ into Eqs. ͑20b͒ and ͑20c͒, so that one has the following stochastic evolution equations ͑that hopefully describe the dynamics on the stochastic center manifold͒: ͑22b͒ One should note that Eqs. ͑22a͒ and ͑22b͒ also can be found by naïvely adding the stochastic terms to the reduced system of evolution equations for the deterministic problem ͓Eqs. ͑18a͒ and ͑18b͔͒. This type of stochastic center manifold reduction has been done for the case of discrete noise. 27 Additionally, in many other fields ͑e.g., oceanography, solid mechanics, fluid mechanics͒, researchers have performed stochastic model reduction using a Karhunen-Loève expansion ͑principal component analysis, proper orthogonal decomposition͒. 33, 34 However, this linear projection does not properly capture the nonlinear effects. Furthermore, one must subjectively choose the number of modes needed for the expansion. Therefore, even though the solution to the reduced model found using this technique may have the correct statistics, individual solution realizations will not agree with the original, complete solution. We will show that Eqs. ͑22a͒ and ͑22b͒ do not contain the correct projection of the noise onto the center manifold. Therefore, when solving the reduced system, one does not obtain the correct solution. Such errors in stochastic epidemic modeling impact the prediction of disease outbreak when modeling the spread of a disease in a population. Using the same numerical scheme previously described, we numerically integrate the complete stochastic system of transformed equations of the SEIR model ͓Eqs. ͑20a͒-͑20c͔͒ along with the reduced system of equations that is based on the deterministic center manifold with a replacement of the noise terms ͓Eqs. ͑22a͒ and ͑22b͔͒. The complete system is solved for U, V, and W, while the reduced system is solved for V and W. In the latter case, U is computed using the center manifold equation given by Eq. ͑17͒. Once the values of U, V, and W are known, we compute the values of S, Ē , and Ī using the transformation given by Eqs. ͑7a͒-͑7c͒. We shift S, Ē , and Ī, respectively, by S 0 , E 0 , and I 0 to find the values of S, E, and I. Figures 3͑a͒ and 3͑b͒ compare the time series of the fraction of the population that is infected with a disease I, computed using the complete stochastic system of transformed equations of the SEIR model ͓Eqs. ͑20a͒-͑20c͔͒ with the time series of I computed using the reduced system of equations of the SEIR model that is based on the deterministic center manifold with a replacement of the noise terms ͓Eqs. ͑22a͒ and ͑22b͔͒. Figure 3͑a͒ shows the initial transients, while Fig. 3͑b͒ shows the time series after the transients have decayed. One can see that the solution computed using the reduced system quickly becomes out of phase with the solution of the complete system. Use of this reduced system would lead to an incorrect prediction of the initial disease outbreak. Additionally, the predicted amplitude of the initial outbreak would be incorrect. The poor agreement, both in phase and amplitude, between the two solutions con-tinues for long periods of time, as seen in Fig. 3͑b͒ . We also have computed the cross correlation of the two time series shown in Figs. 3͑a͒ and 3͑b͒ to be approximately 0.34. Since the cross correlation measures the similarity between the two time series, this low value quantitatively suggests a poor agreement between the two solutions. Using the same systems of transformed equations, we compute 140 years worth of time series for 500 realizations. Ignoring the first 40 years of transient solution, the data are used to create a histogram representing the probability density p SI of the S and I values. Figure 4͑a͒ shows the histogram associated with the complete stochastic system of transformed equations, while Fig. 4͑b͒ shows the histogram associated with the reduced system of equations with a replacement of the noise terms. The color-bar values in Figs. 4͑a͒ and 4͑b͒ have been normalized by 10 −3 . One can see by comparing Fig. 4͑a͒ with Fig. 4͑b͒ that the two probability distributions qualitatively look the same. It is also possible to compare the two distributions using a quantitative measure. The Kullback-Leibler divergence or relative entropy measures the difference between the two probability distributions as where P i,j refers to the ͑i , j͒th component of the probability density found using the complete stochastic system of transformed equations ͓Fig. 4͑a͔͒ and Q i,j refers to the ͑i , j͒th component of the probability density found using the reduced system of equations ͓Fig. 4͑b͔͒. In our numerical computation of the relative entropy, we have added 10 −10 to each P ij and Q ij . This eliminates the possibility of having a Q ij =0 in the denominator of Eq. ͑23͒ and does not have much of an effect on the relative entropy. If the two histograms were identical, then the relative entropy given by Eq. ͑23͒ would be d KL = 0. The two histograms shown in Figs. 4͑a͒ and 4͑b͒ have a relative entropy of d KL = 0.0391, which means that the two histograms, while not identical, are quantitatively very similar. However, one cannot rely entirely on the histograms alone to say that the solutions of the complete system and the reduced system agree. As we have seen in Figs. 3͑a͒ and 3͑b͒, the two solutions have differing amplitudes and are out of phase with one another. It is important to note that these features are not picked up by the histograms of Fig. 4 . To project the noise correctly onto the center manifold, we will derive a normal form coordinate transform for the complete stochastic system of transformed equations of the SEIR model given by Eqs. ͑20a͒-͑20c͒. The particular method we use to construct the normal form coordinate transform not only reduces the dimension of the dynamics, but also separates all of the fast processes from all of the slow processes. 25 This technique has been modified and applied to the large fluctuations of multiscale problems. 17 Many publications [19] [20] [21] [22] discuss the simplification of a stochastic dynamical system using a stochastic normal form transformation. In some of this work, 19, 22 anticipative noise processes appeared in the normal form transformations, but these integrals of the noise process into the future were not dealt with rigorously. Later, the rigorous theoretical analysis needed to support normal form coordinate transforms was developed in Refs. 23 and 24. The technical problem of the anticipative noise integrals also was dealt with rigorously in this work. Even later, another stochastic normal form transformation was developed. 25 This new method allows for the "͓removal of͔ anticipation… from the slow modes with the result that no anticipation is required after the fast transients decay" ͑Ref. 25, p. 13͒. The removal of anticipation leads to a simplification of the normal form. Nonetheless, this simpler normal form retains its accuracy with the original stochastic system. 25 We shall use the method of Ref. 25 to simplify our stochastic dynamical system to one that emulates the long-term dynamics of the original system. The method involves five principles, which we recapitulate here for completeness: ͑1͒ Avoid unbounded secular terms in both the transformation and the evolution equations to ensure a uniform asymptotic approximation. ͑2͒ Decouple all of the slow processes from the fast processes to ensure a valid long-term model. ͑3͒ Insist that the stochastic slow manifold is precisely the transformed fast processes coordinate being equal to zero. ͑4͒ To simplify matters, eliminate as many as possible of the terms in the evolution equations. ͑5͒ Try to remove all fast processes from the slow processes by avoiding as much as possible the fast time memory integrals in the evolution equations. In practice, the original stochastic system of equations ͑which satisfies the necessary spectral requirements͒ in ͑U , V , W͒ T coordinates is transformed to a new ͑Y , X 1 , X 2 ͒ T coordinate system using a near-identity stochastic coordinate transform given as where the specific form of ͑Y , X 1 , X 2 , ͒, ͑Y , X 1 , X 2 , ͒, and ͑Y , X 1 , X 2 , ͒ is chosen to simplify the original system according to the five principles listed previously, and is found using an iterative procedure. To outline the procedure, we provide details for a simple example in Appendix A. Several iterations lead to coordinate transforms for U, V, and W along with evolution equations describing the Y-dynamics, X 1 -dynamics, and X 2 -dynamics in the new coordinate system. The Y-dynamics have exponential decay to the Y = 0 slow manifold. Substitution of Y = 0 leads to the coordinate transforms All of the stochastic terms in Eqs. ͑25a͒-͑25c͒ consist of integrals of the noise process into the past ͑convolutions͒, as given by Eqs. ͑26͒ and ͑27͒. These memory integrals are fast-time processes. Since we are interested in the long term slow processes and since the expectation of G equals e −Ꭽ ‫ء‬ E͓͔, where E͓͔ = 0, we neglect the memory integrals and the higher-order multiplicative terms found in Eqs. ͑25a͒-͑25c͒ so that Note that Eq. ͑28a͒ is the deterministic center manifold equation, and at first order, matches the center manifold equation that was found previously ͓Eq. ͑17͔͒. Substitution of Y = 0 and neglecting all multiplicative noise terms and memory integrals using the argument from above ͑so that we consider only first-order noise terms͒ leads to the following reduced system of evolution equations on the center manifold: The specific form of F and G in Eqs. ͑29a͒ and ͑29b͒ is complicated and is therefore presented in Appendix B. We numerically integrate the complete stochastic system of transformed equations of the SEIR model ͓Eqs. ͑20a͒-͑20c͔͒ along with the reduced system of equations that is found using the stochastic normal form coordinate trans form ͓Eqs. ͑29a͒, ͑29b͒, ͑B1a͒, and ͑B1b͔͒. The complete system is solved for U, V, and W, while the reduced system is solved for X 1 = V and X 2 = W. In the latter case, U is computed using the center manifold equation given by Eq. ͑28a͒. As before, once the values of U, V, and W are known, we compute the values of S, Ē , and Ī using the transformation given by Eqs. ͑7a͒-͑7c͒. We shift S, Ē , and Ī, respectively, by S 0 , E 0 , and I 0 to find the values of S, E, and I. 5. ͑Color online͒ Time series of the fraction of the population that is infected with a disease I, computed using the complete stochastic system of transformed equations of the SEIR model ͓Eqs. ͑20a͒-͑20c͔͒ ͑red solid line͒, and computed using the reduced system of equations of the SEIR model that is found using the stochastic normal form coordinate transform ͓Eqs. ͑29a͒, ͑29b͒, ͑B1a͒, and ͑B1b͔͒ ͑blue dashed line͒. The standard deviation of the noise intensity used in the simulation is i = 0.0005, i =4,5,6. The time series is shown for ͑a͒ t =0 to t = 40 and for ͑b͒ t =40 to t = 100. computed using the complete stochastic system of transformed equations of the SEIR model ͓Eqs. ͑20a͒-͑20c͔͒ with the time series of I computed using the reduced system of equations of the SEIR model that is found using the stochastic normal form coordinate transform ͓Eqs. ͑29a͒, ͑29b͒, ͑B1a͒, and ͑B1b͒ ͔. Figure 5͑a͒ shows the initial transients, while Fig. 5͑b͒ shows the time series after the transients have decayed. One can see that there is an excellent agreement between the two solutions. The initial outbreak is successfully captured by the reduced system. Furthermore, Fig. 5͑b͒ shows that the reduced system accurately predicts recurrent outbreaks for a time scale that is orders of magnitude longer than the relaxation time. This is not surprising since the solution decays exponentially throughout the transient and then remains close to the lower-dimensional center manifold. Since we are not looking at periodic orbits, there are no secular terms in the asymptotic expansion, and the result is valid for all time. Additionally, any noise drift on the center manifold results in bounded solutions due to sufficient dissipation transverse to the manifold. The cross correlation of the two time series shown in Fig. 5 is approximately 0.98, which quantitatively suggests there is an excellent agreement between the two solutions. Using the same systems of transformed equations, we compute 140 years worth of time series for 500 realizations. As before, we ignore the first 40 years worth of transient solution, and the data are used to create a histogram representing the probability density p SI of the S and I values. Figure 6͑a͒ shows the histogram associated with the complete stochastic system of transformed equations, while Fig. 6͑b͒ shows the histogram associated with the reduced system of equations found using the normal form coordinate transform. The color-bar values in Figs. 6͑a͒ and 6͑b͒ have been normalized by 10 −3 . As we saw with Figs. 4͑a͒ and 4͑b͒, the probability distribution shown in Fig. 6͑a͒ looks qualitatively the same as the probability distribution shown in Fig. 6͑b͒ . Using the Kullback-Leibler divergence given by Eq. ͑23͒, we have found that the two histograms shown in Figs. 6͑a͒ and 6͑b͒ have a relative entropy of d KL = 0.0953. Since this value is close to zero, the two histograms are quantitatively very similar. In addition to computing the cross correlation between the solution of the original system and the solutions of the two reduced systems for i = 0.0005, we have computed the cross correlation for the case of zero noise as well as for noise intensities ranging from = 5.0ϫ 10 −10 to = 5.0ϫ 10 −5 . These cross correlations were computed using time series from t = 800 to t = 1000. For the deterministic case ͑zero noise͒, the cross correlation between the time series which were computed using the original system and the reduced system based on the deterministic center manifold is 1.0, since the agreement is perfect. The cross correlation between the original system and the reduced system found using the stochastic normal form is also 1.0. Figure 7 shows the cross correlation between the original system and the two reduced systems for various values of . One can see in Fig. 7 that the solutions found using the reduced system based on the deterministic center manifold compare poorly with the original system at very low noise values. Furthermore, as the noise increases, the agreement between the two solutions gets worse. On the other hand, Fig. 7 shows that the solutions computed using the reduced system found using the normal form coordinate transform compare very well with the solutions to the original system across a wide range of small noise intensities. We have demonstrated that the normal form coordinate transform method reduces the Langevin system so that both the noise and dynamics are accurately projected onto the lower-dimensional manifold. It is natural to consider ͑a͒ the replacement of the stochastic term by a deterministic period drive of small amplitude and ͑b͒ the extension to finite populations. These cases are discussed, respectively, in Secs. VI A and VI B. . 6 . ͑Color online͒ Histogram of probability density p SI of the S and I values found using ͑a͒ the complete stochastic system of transformed equations for the SEIR model with mortality ͓Eqs. ͑20a͒-͑20c͔͒ and ͑b͒ the reduced system of equations of the SEIR model with mortality that is found using the stochastic normal form coordinate transform ͓Eqs. ͑29a͒, ͑29b͒, ͑B1a͒, and ͑B1b͔͒. The histograms are created using 100 years worth of time series ͑starting with year 40͒ for 500 realizations, and the color-bar values have been normalized by 10 −3 . Cross correlation between time series found using the original stochastic system of transformed equations and the reduced system of equations based on the deterministic center manifold ͑"circle" markers͒ and cross correlation between time series found using the original stochastic system of transformed equations and the reduced system of equations based on the stochastic normal form coordinate transform ͑"square" markers͒. The cross correlation is computed using time series from t = 800 to t = 1000. A single time series realization of the noise might be thought of as a deterministic function of small amplitude driving the system. One could rederive the normal form for such a deterministic function. However, since our derived normal form holds specifically for the case of white noise, we show that a simple replacement of the stochastic realization with a deterministic realization does not work. As an example, one could consider the following sinusoidal functions: 1 1 = cos͑10t͒/8000, ͑30a͒ 2 2 = sin͑4t͒/8000, ͑30b͒ where 4 4 , 5 5 , and 6 6 are given by Eqs. ͑21a͒-͑21c͒. Using Eqs. ͑30a͒-͑30c͒ or some other similar deterministic drive, the solution computed using the reduced system based on the deterministic center manifold analysis will agree perfectly with the solution computed using the complete system of equations. On the other hand, since the reduced system based on the normal form analysis was derived specifically for white noise, the transient solution found using this reduced system will not agree with the solution found using the complete system. It is possible to find a normal form coordinate transform for periodic forcing, but the normal form will be different than the one derived in this article for white noise. Figures 8͑a͒ and 8͑b͒ compare the time series of the fraction of the population that is infected with a disease I, computed using the complete system of transformed equations of the SEIR model ͓Eqs. ͑20a͒-͑20c͔͒ with the time series of I computed using the reduced system of equations of the SEIR model that is found using the stochastic normal form coordinate transform ͓Eqs. ͑29a͒, ͑29b͒, ͑B1a͒, and ͑B1b͔͒, but where the stochastic terms of both systems have been replaced by the deterministic terms given by Eqs. ͑30a͒-͑30c͒. Figure 8͑a͒ shows the initial transients, while Fig. 8͑b͒ shows a piece of the time series after the transients have decayed. One can see in Figs. 8͑a͒ and 8͑b͒ that although the two solutions eventually become relatively synchronized with one another, there is a poor agreement, both in phase and amplitude, throughout the transient. The solutions to the original system and both reduced systems are continuous solutions based on an infinite population assumption and are found using Langevin equations having Gaussian noise. It is interesting to examine the effects of general noise by using a Markov simulation to compare solutions of the original and reduced systems. The complete system in the original variables ͑see p. 2͒ will evolve in time t in the following way: transition rate ͑s − 1,e + 1,i͒ ␤si/N ͑s,e − 1,i + 1͒ ␣e ͑s,e,i − 1͒ ␥i ͑s + 1,e,i͒ N ͑s − 1,e,i͒ s ͑s,e − 1,i͒ e ͑s,e,i − 1͒ i. Using a total population size of N =10ϫ 10 6 , we have performed a Markov simulation of the system. After completing the Markov simulation, we divided s, e, and i by N to find S, E, and I. Figure 9͑a͒ shows a time series, after the transients have decayed, of the fraction of the population that is infected with a disease I. The results reflect both the mean and the frequency of the deterministic system. Performing the simulation for 500 realizations allows us to create a histogram representing the probability density p SI of the S and I values. This histogram is shown in Fig. 9͑b͒ , and one can see that the probability density reflects the amplitude, which var- ies with the population size of S and I. The color-bar values in Fig. 9͑b͒ have been normalized by 10 −4 The complete system in the transformed variables has the stable endemic equilibrium at the origin. To bound the dynamics to the first octant, we use the fact that s Ն 0, e Ն 0, and i Ն 0 to derive the appropriate inequalities for the transformed discrete variables u, v, and w. These inequalities can be found in Appendix C as Eq. ͑C1͒. These inequalities enable us to define new discrete variables Y 1 , Y 2 , and Y 3 given by Eqs. ͑C2a͒-͑C2c͒ in Appendix C. In the Y i variables, we define evolution relationships similar to those found in Eq. ͑31͒. The complete transformed system will evolve in time according to the transition and rates given by Eq. ͑C3͒ in Appendix C. After performing a Markov simulation of Eq. ͑C3͒ with a population size of N =10ϫ 10 6 , we can compare the dynamics of the transformed system with the dynamics of the original system by transforming the Y i variables in the time series back to the original s, e, and i variables. Dividing by N yields S, E, and I. Figure 10͑a͒ shows a time series, after the transients have decayed, of the fraction of the population that is infected with a disease I. The mean and the frequency agree with those found from the Markov simulation of the original system. We have performed the simulation for 500 realizations, and a histogram representing the probability density p SI is shown in Fig. 10͑b͒ . The color-bar values in Fig. 10͑b͒ have been normalized by 10 −4 . One can see in Fig. 10͑a͒ that the relative fluctuations of the I component have nearly doubled. While the fluctuation size was 0.152 for the original system, it is 0.310 for the transformed system. Additionally, the two histograms shown in Figs. 9͑b͒ and 10͑b͒ have a relative entropy of d KL = 0.9519, which means they are not in agreement. Because the simulation of the stochastic dynamics in the complete system of transformed variables does not qualitatively ͑or quantitatively͒ resemble the original stochastic system, we cannot expect that the reduced system will agree with either the original or the transformed systems. Therefore, much care should be exercised when extending the model reduction results ͑which show outstanding agreement͒ derived for a specific type of noise in the limit of infinite population to finite populations with a more general type of noise. We have considered the dynamics of a SEIR epidemiological model with stochastic forcing in the form of additive Gaussian noise. We have presented two methods of model reduction, whereby the goal is to project both the noise and the dynamics onto the stochastic center manifold. The first method uses the deterministic center manifold found by neglecting the stochastic terms in the governing equations, while the second method uses a stochastic normal form coordinate transform. Since the original system of governing equations does not have the necessary spectral structure to employ either deterministic or stochastic center manifold theory, the system of equations has been transformed using an appropriate linear transformation coupled with appropriate parameter scaling. At this stage, the first method of model reduction can be performed by computing the deterministic center manifold equation. Substitution of this equation into the complete stochastic system of transformed equations leads to a reduced system of stochastic evolution equations. The solutions of the complete stochastic system of transformed equations as well as the reduced system of equations were computed numerically. We have shown that the individual time series does not agree because the noise has not been correctly projected onto the stochastic center manifold. However, by comparing histograms of the probability density p SI of the S and I values, we saw that there was a very good agreement. This is caused by the fact that although the two solutions are out of phase with one another, their range of amplitude values is similar. The phase difference is not represented in the two histograms. This is a real drawback when trying to predict the timing of outbreaks and leads to potential problems when considering epidemic control, such as the enhancement of disease extinction through random vaccine control. 35 To accurately project the noise onto the manifold, we derived a stochastic normal form coordinate transform for the complete stochastic system of transformed equations. The numerical solution to this reduced system was compared with the solution to the original system, and we showed that there was an excellent agreement both qualitatively and quantitatively. As with the first method, the histograms of the probability density p SI of the S and I values agree very well. It should be noted that the use of these two reduction methods is not constrained to problems in epidemiology, but rather may be used for many types of physical problems. For some generic systems, such as the singularly perturbed, damped Duffing oscillator, either reduction method can be used since the terms in the normal form coordinate transform which lead to the average stochastic center manifold being different from the deterministic center manifold occur at very high order. 17 In other words, the average stochastic center manifold and deterministic center manifold are virtually identical. For the SEIR model considered in this article, there are terms at low order in the normal form transform, which cause a significant difference between the average stochastic center manifold and the deterministic manifold. Therefore, as we have demonstrated, when working with the SEIR model, one must use the normal form coordinate transform method to correctly project the noise onto the center manifold. In summary, we have presented a new method of stochastic model reduction that allows for impressive improvement in time series prediction. The reduced model captures both the amplitude and phase accurately for a temporal scale that is many orders of magnitude longer than the typical relaxation time. Since sufficient statistics of disease data are limited due to short time series collection, the results presented here provide a potential method to properly model real, stochastic disease data in the time domain. Such longterm accuracy of the reduced model will allow for the application of effective control of a disease where phase differences between outbreak times and vaccine controls are important. Additionally, since our method is general, it may be applied to very high-dimensional epidemic models, such as those involving adaptive networks. From a dynamical systems viewpoint, the reduction method has the potential to accurately capture new, emergent dynamics as we increase the size of the random fluctuations. This could be a means to identify new noise-induced phenomena in generic stochastic systems. Infectious Diseases of Humans ͑Oxford University Press The Mathematical Theory of Infectious Diseases ͑Griffin Random Dynamical Systems ͑Springer-Verlag Differential Equations, Dynamical Systems, and Linear Algebra ͑Academic Introduction to Applied Nonlinear Dynamical Systems and Chaos ͑Springer-Verlag Applications of Centre Manifold Theory ͑Springer-Verlag The authors benefited from the comments and suggestions of anonymous reviewers. We gratefully acknowledge the support from the Office of Naval Research and the Air Force Office of Scientific Research. E. We consider the system given byThe iterative procedure begins by lettingx Ϸ X, ͑A2a͒and by finding a change to the y coordinate ͑fast process͒ with the form y = Y + ͑,X,Y͒ +¯,where and G are small corrections to the coordinate transform and the corresponding evolution equation. Substitution of Eqs. ͑A2a͒ and ͑A3b͒ into Eq. ͑A1b͒ gives the equation Replacing YЈ = ‫ץ‬Y / ‫ץ‬ with −Y + G ͓Eq. ͑A3b͔͒, noting that ‫ץ‬X / ‫ץ‬ =0 ͓Eq. ͑A2b͔͒, and ignoring the term ‫ץ‬ / ‫ץ‬Y · G since it is a product of small corrections leads to Equation ͑A5͒ must now be solved for G and . In order to keep the evolution equation ͓Eq. ͑A3b͔͒ as simple as possible ͓principle ͑4͒ of Sec. V͔, we let G = 0, which means that the coordinate transform ͓Eq. ͑A3a͔͒ is modified by = X − X 3 . Therefore, the new approximation of the coordinate transform and its dynamics are given bywhere = ͉͑X , Y , , ͉͒ so that provides a count of the number of X, Y, , and factors in any one term.For the second iteration, we seek a correction to the x coordinate ͑slow process͒ with the formwhere and F are small corrections. Substitution of Eqs. ͑A6a͒ and ͑A7b͒ into Eq. ͑A1a͒ leads toReplacing XЈ = ‫ץ‬X / ‫ץ‬ with F ͓Eq. ͑A7b͔͒, replacing ‫ץ‬Y / ‫ץ‬ with −Y ͓Eq. ͑A6b͔͒, and ignoring the term ‫ץ‬ / ‫ץ‬X · F since it is a product of small corrections gives the equationmust now be solved for F and . As in the first step, we employ principle ͑4͒ and keep the evolution equation ͓Eq. ͑A7b͔͒ as simple as possible. However, since the terms ͑X − X 3 ͒ located on the right-hand side of Eq. ͑A9͒ do not contain or Y, these terms must be included in F. Therefore, one piece of F will be F = ͑X − X 3 ͒.The remaining deterministic term on the right-hand side of Eq. ͑A9͒ contains Y. This term can therefore be integrated into . The equation to be solved iswhose solution is given as =−Y.To abide by principle ͑4͒, we would like to integrate the stochastic piece on the right-hand side of Eq. ͑A9͒ into by solving the equation ‫ץ/ץ‬ = . ͑A11͒However, the solution of Eq. ͑A11͒ is given bywhich has secular growth like a Wiener process. Since this would violate principle ͑1͒, we must let F = .Putting the three pieces together yields =−Y and F = ͑X − X 3 ͒ + . Therefore, the new approximation of the coordinate transform and its dynamics are given byThe construction of the normal form continues by seeking corrections and F to the x coordinate transform and the X evolution using the updated residual of the x equation ͓Eq. ͑A1a͔͒, and by seeking corrections and G to the y coordinate transform and the Y evolution equation using the updated residual of the y equation ͓Eq. ͑A1b͔͒. The specific form of F and G in Eqs. ͑29a͒ and ͑29b͒ is given as The complete system in the transformed variables has the stable endemic equilibrium at the origin. To bound the dynamics to the first octant, we transform the new variables by using the original properties of s Ն 0, e Ն 0, and i Ն 0, so thatTherefore, we define the following new variables:In these variables, we define evolution relationships similar to Eq. ͑31͒. The complete transformed system will evolve in the following way:transition rate