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 reduction in high-dimensional, stochastic systems is an important and fundamental problem in nonlinear dynamical systems. In this article, we present a general theory of stochastic model reduction, which is based on a normal form coordinate transform. This nonlinear, stochastic projection allows for the deterministic and stochastic dynamics to interact correctly on the lower-dimensional manifold so that the dynamics predicted by the reduced, stochastic system agrees well with the dynamics predicted by the original, high-dimensional stochastic system. Although the method may be applied to any physical or biological system with well-separated time scales, here we apply the method to an epidemiological model. We show that when compared with the original stochastic epidemic model, the reduced model properly captures the initial and recurrent disease outbreaks, both in amplitude and phase. This long-term accuracy of the reduced model allows for the application of effective disease control where phase differences between outbreak times and vaccine controls are important. Additionally, in practice, one can only measure the number of infectious individuals in a population. Our method allows one to predict the number of unobserved exposed individuals based on the observed number of infectious individuals.

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 dynamics3 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–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 manifold) was performed in Refs. 19–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 low-dimensional models, and an analysis of the solutions found using the high and low-dimensional systems is performed.

FIG. 1.

Schematic demonstrating the injection of a common noise into both the high-dimensional system and its associated low-dimensional system.

FIG. 1.

Schematic demonstrating the injection of a common noise into both the high-dimensional system and its associated low-dimensional system.

Close modal

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

Ṡ(t)=μβI(t)S(t)μS(t)+σ1ϕ1(t),
(1a)
Ė(t)=βI(t)S(t)(α+μ)E(t)+σ2ϕ2(t),
(1b)
İ(t)=αE(t)(γ+μ)I(t)+σ3ϕ3(t),
(1c)
Ṙ(t)=γI(t)μR(t)+σ4ϕ4(t),
(1d)
where σi is the standard deviation of the noise intensity Di=σi2/2. Each of the noise terms ϕi describes a stochastic, Gaussian white force that is characterized by the correlation functions
ϕi(t)=0,
(2a)
ϕi(t)ϕj(t)=δ(tt)δij.
(2b)

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 R0=αβ/[(α+μ)(γ+μ)]>1. The interpretation of R0 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

(Se,Ee,Ie)=(1,0,0)
(3)

and corresponds to a disease free or extinct equilibrium state. The second fixed point corresponds to an endemic state and is given by

(S0,E0,I0)=((γ+μ)(α+μ)βα,μα+μμ(γ+μ)αβ,μα(γ+μ)(α+μ)μβ).
(4)

To ease the analysis, we define a new set of variables, S¯, E¯, and I¯, as S¯(t)=S(t)S0, E¯(t)=E(t)E0, and I¯(t)=I(t)I0. 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

dS¯dτ=βμI¯S¯(α0+μ2)(γ0+μ2)α0I¯α0μ3β(α0+μ2)(γ0+μ2)S¯,
(5a)
dE¯dτ=βμI¯S¯+(α0+μ2)(γ0+μ2)α0I¯+μ2[α0βμ(α0+μ2)(γ0+μ2)](α0+μ2)(γ0+μ2)S¯(α0+μ2)E¯,
(5b)
dI¯dτ=α0E¯(γ0+μ2)I¯,
(5c)
dμdτ=0,
(5d)
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

P=[110α0+γ0γ000α0+γ0γ001].
(6)

Using the fact that (S¯,E¯,I¯)T=P(U,V,W)T, then the transformation matrix leads to the following definition of new variables: U, V, and W,

U=γ0α0+γ0E¯,
(7a)
V=S¯+γ0α0+γ0E¯,
(7b)
W=I¯+E¯.
(7c)

The application of the transformation matrix to Eqs. (5a)–(5c) leads to the transformed evolution equations given by

dUdτ=α0U+μ2(γ0Vα0U)α0+γ0(γ0+μ2)(α0+μ2)[(α0+γ0)U+γ0W]α0(α0+γ0)μβα0+γ0(γ0W+(α0+γ0)U+μ2α0γ0(γ0+μ2)(α0+μ2))(U+V),
(8a)
dVdτ=α0Uμ2(γ0Vα0U)α0+γ0(γ0+μ2)(α0+μ2)[(α0+γ0)U+γ0W]γ0(α0+γ0)μβα0γ0(α0+γ0)(γ0W+(α0+γ0)U+μ2α0γ0(γ0+μ2)(α0+μ2))(U+V),
(8b)
dWdτ=α0U(γ0+μ2)(U+W)+(γ0+μ2)(α0+μ2)[(α0+γ0)U+γ0W]α0γ0μ2V+μβγ0(γ0W+(α0+γ0)U+μ2α0γ0(γ0+μ2)(α0+μ2))(U+V),
(8c)
dμdτ=0.
(8d)

The Jacobian of Eqs. (8a)–(8d) to zeroth order in μ and evaluated at the origin is

[(α0+γ0)0γ02(α0+γ0)000α0γ0(α0+γ0)000000000],
(9)

which shows that Eqs. (8a)–(8d) may be rewritten in the form

dxdτ=Ax+f(x,y,μ),
(10)
dydτ=By+g(x,y,μ),
(11)
dμdτ=0,
(12)

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,

A=[(α0+γ0)],B=[0α0γ0(α0+γ0)00].
(13)

Therefore, the system will rapidly collapse onto a lower-dimensional manifold given by center manifold theory.32 Furthermore, we know that the center manifold is given by

U=h(V,W,μ),
(14)

where h is an unknown function.

Substitution of Eq. (14) into Eq. (8a) leads to the following center manifold condition:

h(V,W,μ)VdVdτ+h(V,W,μ)WdWdτ=α0h(V,W,μ)+μ2[γ0Vα0h(V,W,μ)]α0+γ0(γ0+μ2)(α0+μ2)[(α0+γ0)h(V,W,μ)+γ0W]α0(α0+γ0)μβα0+γ0(γ0W+(α0+γ0)h(V,W,μ)+μ2α0γ0(γ0+μ2)(α0+μ2))(h(V,W,μ)+V).
(15)

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 μ:

h(V,W,μ)=h0+h2V+h3W+hμμ+h22V2+h23VW+h33W2+hμ2μV+hμ3μW+hμμμ2+,
(16)

where h0, h2, h3, 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

U=γ02(α0+γ0)2W+O(ϵ3),
(17)

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

dVdτ=μ2γ02α0W(α0+γ0)3μ4α0W(α0+γ0)2γ0μ2Vα0+γ0(γ0+μ2)α0Wα0+γ0βα02μ(α0+γ0)2(W+μ2(α0+γ0)(γ0+μ2)(α0+μ2))(Vγ02W(α0+γ0)2),
(18a)
dWdτ=μ2γ02W(α0+γ0)2+μ4Wα0+γ0μ2V+βμα0α0+γ0(W+μ2(α0+γ0)(γ0+μ2)(α0+μ2))(Vγ02W(α0+γ0)2).
(18b)

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

S¯̇(t)=βI¯S¯(α+μ)(γ+μ)αI¯αμβ(α+μ)(γ+μ)S¯+σ1ϕ1(t),
(19a)
E¯̇(t)=βI¯S¯+(α+μ)(γ+μ)αI¯+μ[αβ(α+μ)(γ+μ)](α+μ)(γ+μ)S¯(α+μ)E¯+σ2ϕ2(t),
(19b)
I¯̇(t)=αE¯(γ+μ)I¯+σ3ϕ3(t).
(19c)

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

dUdτ=α0U+μ2(γ0Vα0U)α0+γ0(γ0+μ2)(α0+μ2)[(α0+γ0)U+γ0W]α0(α0+γ0)μβα0+γ0(γ0W+(α0+γ0)U+μ2α0γ0(γ0+μ2)(α0+μ2))(U+V)+σ4ϕ4,
(20a)
dVdτ=α0Uμ2(γ0Vα0U)α0+γ0(γ0+μ2)(α0+μ2)[(α0+γ0)U+γ0W]γ0(α0+γ0)μβα0γ0(α0+γ0)(γ0W+(α0+γ0)U+μ2α0γ0(γ0+μ2)(α0+μ2))(U+V)+σ5ϕ5,
(20b)
dWdτ=α0U(γ0+μ2)(U+W)+(γ0+μ2)(α0+μ2)[(α0+γ0)U+γ0W]α0γ0μ2V+μβγ0(γ0W+(α0+γ0)U+μ2α0γ0(γ0+μ2)(α0+μ2))(U+V)+σ6ϕ6,
(20c)
where
σ4ϕ4=μγ0α0+γ0σ2ϕ2,
(21a)
σ5ϕ5=μσ1ϕ1+μγ0α0+γ0σ2ϕ2,
(21b)
σ6ϕ6=μσ2ϕ2+μσ3ϕ3.
(21c)
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¯, E¯, and I¯ using the transformation given by Eqs. (7a)–(7c). We shift S¯, E¯, and I¯, respectively, by S0, E0, and I0 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.

FIG. 2.

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 [Eqs. (1a)–(1c)] (red solid line), and computed using the transformed stochastic system of equations of the SEIR model [Eqs. (20a)–(20c)] (blue dashed line). The standard deviation of the noise intensity used in the simulation is σi=0.0005, i=1,,6.

FIG. 2.

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 [Eqs. (1a)–(1c)] (red solid line), and computed using the transformed stochastic system of equations of the SEIR model [Eqs. (20a)–(20c)] (blue dashed line). The standard deviation of the noise intensity used in the simulation is σi=0.0005, i=1,,6.

Close modal

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 σ6ϕ6 of the transformed system are new, independent noise processes with different variance than the σ1ϕ1, σ2ϕ2, and σ3ϕ3 noise processes found in the original system. If we carefully take the noise realization from the original system, transform this noise using Eqs. (21a)–(21c), and use this realization to solve the transformed system, then the two solutions would be identical.

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):

dVdτ=μ2γ02α0W(α0+γ0)3μ4α0W(α0+γ0)2γ0μ2Vα0+γ0(γ0+μ2)Wα0α0+γ0βα02μ(α0+γ0)2(W+μ2(α0+γ0)(γ0+μ2)(α0+μ2))(Vγ02W(α0+γ0)2)+σ5ϕ5,
(22a)
dWdτ=μ2γ02W(α0+γ0)2+μ4Wα0+γ0μ2V+βμα0α0+γ0(W+μ2(α0+γ0)(γ0+μ2)(α0+μ2))(Vγ02W(α0+γ0)2)+σ6ϕ6.
(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¯, E¯, and I¯ using the transformation given by Eqs. (7a)–(7c). We shift S¯, E¯, and I¯, respectively, by S0, E0, and I0 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 continues 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.

FIG. 3.

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 based on the deterministic center manifold with a replacement of the noise terms [Eqs. (22a) and (22b)] (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.

FIG. 3.

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 based on the deterministic center manifold with a replacement of the noise terms [Eqs. (22a) and (22b)] (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.

Close modal

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 pSI 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 103.

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

dKL=i,jPi,j|log(Pi,jQi,j)|,
(23)

where Pi,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 Qi,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 1010 to each Pij and Qij. This eliminates the possibility of having a Qij=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 dKL=0. The two histograms shown in Figs. 4(a) and 4(b) have a relative entropy of dKL=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.

FIG. 4.

Histogram of probability density pSI of the S and I values found using (a) the complete stochastic system of transformed equations for the SEIR model [Eqs. (20a)–(20c)] and (b) 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)]. 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 103.

FIG. 4.

Histogram of probability density pSI of the S and I values found using (a) the complete stochastic system of transformed equations for the SEIR model [Eqs. (20a)–(20c)] and (b) 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)]. 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 103.

Close modal

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 publications19–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,X1,X2)T coordinate system using a near-identity stochastic coordinate transform given as

U=Y+ξ(Y,X1,X2,τ),
(24a)
V=X1+η(Y,X1,X2,τ),
(24b)
W=X2+ρ(Y,X1,X2,τ),
(24c)
where the specific form of ξ(Y,X1,X2,τ), η(Y,X1,X2,τ), and ρ(Y,X1,X2,τ) 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, X1-dynamics, and X2-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

U=σ4G(ϕ4)+γ02(σ6G(ϕ6)X2)(α0+γ0)2+μ[γ0β[σ5X2G(ϕ5)X1X2+σ6X1G(ϕ6)](α0+γ0)2σ4α0βγ0G2(ϕ4)[(α0+γ0)X1+γ0X2](α0+γ0)(γ0+μ2)(α0+μ2)]+μ2[γ0[X1α02X2+2σ6G(ϕ6)]α0(α0+γ0)σ4G2(ϕ4)(γ0+μ2)(α0+μ2)(2γ03+α03α0+γ0+α02γ02(α0+γ0)(γ0+μ2)(α0+μ2))σ5γ0G(ϕ5)(α0+γ0)2]+μ3[β[X1+σ5G(ϕ5)](α0+γ0)2σ4βG2(ϕ4)(γ0+μ2)(α0+μ2)(α0γ0α0+γ0+γ0X2+X1(α0+γ0))+β(σ6X1G(ϕ6)+σ5X2G(ϕ5)X1X2α0)α0(α0+γ0)]+O(μ4),
(25a)
V=X1+μ[σ4α0βX1G(ϕ4)α0+γ0+σ4α0βX2G(ϕ4)(α0+γ0)2]+μ2[σ4G(ϕ4)(α02+α0γ0+γ02)γ0(α0+γ0)2]+μ3[σ4α0βG(ϕ4)γ0(α0+γ0)2+σ4βX2G(ϕ4)γ0(α0+γ0)+σ4βX1G(ϕ4)γ02]+O(μ4),
(25b)
W=X2+μ[σ4βX1G(ϕ4)γ0σ4βX2G(ϕ4)(α0+γ0)]+μ2[σ4G(ϕ4)(α02+α0γ0+γ02)α0γ0(α0+γ0)]+μ3[σ4βG(ϕ4)γ0(α0+γ0)σ4(α0+γ0)βX1G(ϕ4)α0γ02σ4βX2G(ϕ4)α0γ0]+O(μ4),
(25c)
where
(26)
G(ϕ)=eτϕ=τexp[(τs)]ϕ(s)ds,
=α0γ0(α0+γ0)(α0+μ2)(γ0+μ2)
and

G2(ϕ)=eτeτϕ.
(27)

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

U=γ02X2(α0+γ0)2μβX1(α0+γ0)(μ2(α0+γ0)+γ0X2(α0+γ0)+μ2X2)+μ2γ0(α0+γ0)(X12X2α0),
(28a)
V=X1,
(28b)
W=X2.
(28c)
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:

dX1dτ=F(X1(τ),X2(τ)),
(29a)
dX2dτ=G(X1(τ),X2(τ)).
(29b)
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 X1=V and X2=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¯, E¯, and I¯ using the transformation given by Eqs. (7a)–(7c). We shift S¯, E¯, and I¯, respectively, by S0, E0, and I0 to find the values of S, E, and I.

Figures 5(a) and 5(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 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.

FIG. 5.

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.

FIG. 5.

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.

Close modal

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 pSI 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 103.

FIG. 6.

Histogram of probability density pSI 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 103.

FIG. 6.

Histogram of probability density pSI 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 103.

Close modal

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 dKL=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×1010 to σ=5.0×105. 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 σ.

FIG. 7.

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.

FIG. 7.

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.

Close modal

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.

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(10πμt)/8000,
(30a)
σ2ϕ2=sin(4πμt)/8000,
(30b)
σ3ϕ3=cos(10πμt)/8000,
(30c)
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.

FIG. 8.

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)] (red solid line), and computed using the reduced system of equations of the SEIR model that is found using the normal form coordinate transform [Eqs. (29a), (29b), (B1a), and (B1b)] (blue dashed line). The stochastic terms in both systems have been replaced by the deterministic terms given by Eqs. (30a) and (30b). The time series is shown from (a) t=0 to t=25 and from (b) t=65 to t=70.

FIG. 8.

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)] (red solid line), and computed using the reduced system of equations of the SEIR model that is found using the normal form coordinate transform [Eqs. (29a), (29b), (B1a), and (B1b)] (blue dashed line). The stochastic terms in both systems have been replaced by the deterministic terms given by Eqs. (30a) and (30b). The time series is shown from (a) t=0 to t=25 and from (b) t=65 to t=70.

Close modal

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:

transitionrate(s1,e+1,i)βsi/N(s,e1,i+1)αe(s,e,i1)γi(s+1,e,i)μN(s1,e,i)μs(s,e1,i)μe(s,e,i1)μi.
(31)

Using a total population size of N=10×106, 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 pSI 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 varies with the population size of S and I. The color-bar values in Fig. 9(b) have been normalized by 104

FIG. 9.

(a) Time series of the fraction of the population that is infected with a disease I, computed using a Markov simulation of the complete original equations of the SEIR model [Eq. (31)] and (b) a histogram of probability density pSI of the S and I values found using a Markov simulation of Eq. (31). The histogram is created using 100 years worth of data (starting with year 40) for 500 realizations, and the color-bar values have been normalized by 104.

FIG. 9.

(a) Time series of the fraction of the population that is infected with a disease I, computed using a Markov simulation of the complete original equations of the SEIR model [Eq. (31)] and (b) a histogram of probability density pSI of the S and I values found using a Markov simulation of Eq. (31). The histogram is created using 100 years worth of data (starting with year 40) for 500 realizations, and the color-bar values have been normalized by 104.

Close modal

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 s0, e0, and i0 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 Y1, Y2, and Y3 given by Eqs. (C2a)–(C2c) in Appendix  C.

In the Yi 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×106, we can compare the dynamics of the transformed system with the dynamics of the original system by transforming the Yi 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 pSI is shown in Fig. 10(b). The color-bar values in Fig. 10(b) have been normalized by 104. 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 dKL=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.

FIG. 10.

(a) Time series of the fraction of the population that is infected with a disease I, computed using a Markov simulation of the complete transformed equations of the SEIR model [Eq. (C3)] and (b) a histogram of probability density pSI of the S and I values found using a Markov simulation of Eq. (C3). The histogram is created using 100 years worth of data (starting with year 40) for 500 realizations, and the color-bar values have been normalized by 104.

FIG. 10.

(a) Time series of the fraction of the population that is infected with a disease I, computed using a Markov simulation of the complete transformed equations of the SEIR model [Eq. (C3)] and (b) a histogram of probability density pSI of the S and I values found using a Markov simulation of Eq. (C3). The histogram is created using 100 years worth of data (starting with year 40) for 500 realizations, and the color-bar values have been normalized by 104.

Close modal

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 pSI 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 pSI 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 long-term 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.

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.F. is supported by the National Research Council Research Fellowship and L.B. is supported by the National Institute of General Medical Sciences (Award No. R01GM090204). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of General Medical Sciences or the National Institutes of Health.

We consider the system given by

dxdτ=μ(y+σϕ),
(A1a)
dydτ=xx3y,
(A1b)
dμdτ=0.
(A1c)
The iterative procedure begins by letting
xX,
(A2a)
X0,
(A2b)
and by finding a change to the y coordinate (fast process) with the form
y=Y+η(τ,X,Y)+,
(A3a)
Y=Y+G(τ,X,Y)+,
(A3b)
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

Y+ητ+ηXXτ+ηYYτ=Yη+XX3.
(A4)

Replacing Y=Y/τ with Y+G [Eq. (A3b)], noting that X/τ=0 [Eq. (A2b)], and ignoring the term η/YG since it is a product of small corrections leads to

G+ητYηY+η=XX3.
(A5)

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 η=XX3. Therefore, the new approximation of the coordinate transform and its dynamics are given by

y=Y+XX3+O(ζ2),
(A6a)
Y=Y+O(ζ2),
(A6b)
where ζ=|(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 form

x=X+ξ(τ,X,Y)+,
(A7a)
X=F(τ,X,Y)+,
(A7b)
where ξ and F are small corrections. Substitution of Eqs. (A6a) and (A7b) into Eq. (A1a) leads to

X+ξτ+ξXXτ+ξYYτ=μ(Y+XX3)+μσϕ.
(A8)

Replacing X=X/τ with F [Eq. (A7b)], replacing Y/τ with Y [Eq. (A6b)], and ignoring the term ξ/XF since it is a product of small corrections gives the equation

F+ξτYξY=μ(Y+XX3)+μσϕ.
(A9)

Equation (A9) must 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 μ(XX3) 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=μ(XX3).

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 is

YξY=μY,
(A10)

whose 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 by

ξ=μσϕdτ,
(A12)

which 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=μ(XX3)+μσϕ. Therefore, the new approximation of the coordinate transform and its dynamics are given by

x=XμY+O(ζ3),
(A13a)
X=μ(XX3)+μσϕ+O(ζ3).
(A13b)

The 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

F=[α02γ03X2+μβα02(α0+γ0)γ02(γ02α0+γ0X22+α0X1X2)+μ2(α0γ03X1+α0γ02(2α03+5α02γ0+5α0γ02+γ03)(α0+γ0)2X2α02β2γ02(α0+γ0)X12X2α02β2γ03(α0+γ0)2X1X22)+μ3(α02βγ0X1α02βγ03(α0+γ0)2X2+α02βγ02(α0+γ0)X123α0βγ03(α0+γ0)X22+α0βγ0(α03α02γ03α0γ023γ03)(α0+γ0)2X1X2)][γ0(α0+γ0)(α0+μ2)(γ0+μ2)]+σ5ϕ5μ2(α02+α0γ0+γ02)(α0+γ0)3(σ4α0ϕ4γ0+σ6γ0ϕ6α0+γ0)μ3α0β(α0+γ0)3(σ4α0ϕ4γ0+σ6γ0ϕ6(α0+γ0)),
(B1a)
G=[μ(α03βγ02(α0+γ0)X1X2α02βγ04(α0+γ0)2X22)+μ2(α02γ02X1+α02γ04(α0+γ0)2X2α02β2γ02(α0+γ0)X12X2α02β2γ03(α0+γ0)2X1X22)+μ3(α02βγ0X1α02βγ03(α0+γ0)2X2+α02βγ02(α0+μ2)(γ0+μ2)(α0+γ0)X123α0βγ03(α0+γ0)X223α0βγ02(α0+μ2)(γ0+μ2)X1X2+α02βγ0(α02+2α0γ0+3γ02)(α0+γ0)2X1X2)][α0γ0(α0+μ2)(γ0+μ2)]+σ6ϕ6+μ2σ4(α02+α0γ0+γ02)α0γ0(α0+γ0)ϕ4+μ2σ6[γ03(α0+μ2)(γ0+μ2)+α0γ0(α0+γ0)]α0(α0+γ0)3ϕ6+μ3β(α0+γ0)(σ4ϕ4γ0+σ6γ0ϕ6(α0+γ0)2).
(B1b)

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 s0, e0, and i0, so that

(C1)
uμ2Nγ0α0(α0+γ0),Nγ0(βμ3+α02+γ0α0)α0βμ(α0+γ0)v,
Nμ2(α0+γ0)γ0α0w.
Therefore, we define the following new variables:
Y1=u+Nμ2γ0α0(α0+γ0),
(C2a)
Y2=v+Nγ0(βμ3+α02+γ0α0)α0βμ(α0+γ0),
(C2b)
Y3=w+Nμ2(α0+γ0)γ0α0.
(C2c)
In these variables, we define evolution relationships similar to Eq. (31). The complete transformed system will evolve in τ the following way:

transitionrate(Y1+1,Y2,Y3)βμN(γ0α0+γ0Y2Y3+Y12)(Y11,Y2,Y3)(α0+μ2)Y1+βμN(γ0α0+γ0Y1Y3+Y1Y2)(Y1,Y2+1,Y3)μ2N+βμN(α0(α0+γ0)Y1Y3+α0γ0Y1Y2)(Y1,Y21,Y3)α0Y1+μ2Y2+βμN(α0(α0+γ0)Y2Y3+α0γ0Y12)(Y1,Y2,Y3+1)(α0+γ0)Y1+βμN(Y2Y3+(α0+γ0)γ0Y12)(Y1,Y2,Y31)(γ0+μ2)Y3+βμN(Y1Y3+(α0+γ0)γ0Y1Y2).
(C3)
1.
R. M.
Anderson
and
R. M.
May
,
Infectious Diseases of Humans
(
Oxford University Press
,
New York
,
1991
).
2.
N. T. J.
Bailey
,
The Mathematical Theory of Infectious Diseases
(
Griffin
,
London
,
1975
).
3.
G.
Marion
,
E.
Renshaw
, and
G.
Gibson
,
Theor Popul. Biol.
57
,
197
(
2000
).
4.
H. T. H.
Nguyen
and
P.
Rohani
,
J. R. Soc., Interface
5
,
403
(
2008
).
5.
P.
Rohani
,
M. J.
Keeling
, and
B. T.
Grenfell
,
Am. Nat.
159
,
469
(
2002
).
6.
D. A.
Rand
and
H. B.
Wilson
,
Proc. R. Soc. London, Ser. B
246
,
179
(
1991
).
7.
L.
Billings
,
E. M.
Bollt
, and
I. B.
Schwartz
,
Phys. Rev. Lett.
88
,
234101
(
2002
).
8.
L.
Stone
,
R.
Olinky
, and
A.
Huppert
,
Nature (London)
446
,
533
(
2007
).
9.
I. B.
Schwartz
,
L.
Billings
, and
E. M.
Bollt
,
Phys. Rev. E
70
,
046220
(
2004
).
10.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. E
63
,
066117
(
2001
).
11.
Y.
Moreno
,
R.
Pastor-Satorras
, and
A.
Vespignani
,
Eur. Phys. J. B
26
,
521
(
2002
).
12.
A.
Vazquez
,
Phys. Rev. E
74
,
056101
(
2006
).
13.
D. J.
Watts
,
R.
Muhamad
,
D. C.
Medina
, and
P. S.
Dodds
,
Proc. Natl. Acad. Sci. U.S.A.
102
,
11157
(
2005
).
14.
V.
Colizza
,
A.
Barrat
,
M.
Barthelemy
, and
A.
Vespignani
,
Bull. Math. Biol.
68
,
1893
(
2006
).
15.
L. B.
Shaw
and
I. B.
Schwartz
,
Phys. Rev. E
77
,
066101
(
2008
).
16.
W. T.
Mocek
,
R.
Rudnicki
, and
E. O.
Voit
,
Math. Biosci.
198
,
190
(
2005
).
17.
E.
Forgoston
and
I. B.
Schwartz
,
SIAM J. Appl. Dyn. Syst.
8
,
1190
(
2009
).
18.
P.
Boxler
,
Probab. Theory Relat. Fields
83
,
509
(
1989
).
19.
P. H.
Coullet
,
C.
Elphick
, and
E.
Tirapegui
,
Phys. Lett. A
111
,
277
(
1985
).
20.
E.
Knobloch
and
K. A.
Wiesenfeld
,
J. Stat. Phys.
33
,
611
(
1983
).
21.
N. S.
Namachchivaya
,
Appl. Math. Comput.
38
,
101
(
1990
).
22.
N. S.
Namachchivaya
and
Y. K.
Lin
,
Int. J. Non-Linear Mech.
26
,
931
(
1991
).
23.
L.
Arnold
,
Random Dynamical Systems
(
Springer-Verlag
,
Berlin
,
1998
).
24.
L.
Arnold
and
P.
Imkeller
,
Probab. Theory Relat. Fields
110
,
559
(
1998
).
26.
I.
Schwartz
and
H.
Smith
,
J. Math. Biol.
18
,
233
(
1983
).
27.
L.
Billings
and
I. B.
Schwartz
,
J. Math. Biol.
44
,
31
(
2002
).
28.
Mathematical Epidemiology
, edited by
F.
Brauer
,
P.
van den Driessche
, and
J.
Wu
(
Springer-Verlag
,
Berlin
,
2008
).
29.
O. N.
Bjørnstad
,
B. F.
Finkenstädt
, and
B. T.
Grenfell
,
Ecol. Monogr.
72
,
169
(
2002
).
30.
M. W.
Hirsch
and
S.
Smale
,
Differential Equations, Dynamical Systems, and Linear Algebra
(
Academic
,
New York
,
1974
).
31.
S.
Wiggins
,
Introduction to Applied Nonlinear Dynamical Systems and Chaos
(
Springer-Verlag
,
Berlin
,
1990
).
32.
J.
Carr
,
Applications of Centre Manifold Theory
(
Springer-Verlag
,
Berlin
,
1981
).
33.
A.
Doostan
,
R. G.
Ghanem
, and
J.
Red-Horse
,
Comput. Methods Appl. Mech. Eng.
196
,
3951
(
2007
).
34.
D.
Venturi
,
X.
Wan
, and
G. E.
Karniadakis
,
J. Fluid Mech.
606
,
339
(
2008
).
35.
M. I.
Dykman
,
I. B.
Schwartz
, and
A. S.
Landsman
,
Phys. Rev. Lett.
101
,
078101
(
2008
).