key: cord-0011560-xqvpvebr
authors: Parra-Rojas, César; House, Thomas; McKane, Alan J.
title: Stochastic epidemic dynamics on extremely heterogeneous networks
date: 2016-12-19
journal: Phys Rev E
DOI: 10.1103/physreve.94.062408
sha: de53bb2b749d677e458221edaa8f17912a691d1d
doc_id: 11560
cord_uid: xqvpvebr

Networks of contacts capable of spreading infectious diseases are often observed to be highly heterogeneous, with the majority of individuals having fewer contacts than the mean, and a significant minority having relatively very many contacts. We derive a two-dimensional diffusion model for the full temporal behavior of the stochastic susceptible-infectious-recovered (SIR) model on such a network, by making use of a time-scale separation in the deterministic limit of the dynamics. This low-dimensional process is an accurate approximation to the full model in the limit of large populations, even for cases when the time-scale separation is not too pronounced, provided the maximum degree is not of the order of the population size.

Mathematical models are used throughout infectious disease epidemiology to understand the dynamical processes that shape patterns of disease [1, 2] . While early models did not include complex population structure, modeling approaches now frequently let the epidemic spread take place on a network, which enables greater realism than a model in which all individuals mix homogeneously. It does, however, pose many technical problems for model analysis, particularly the question of how heterogeneity in the number of links each individual participates in-their degree-influences the epidemic [3] [4] [5] .

Similarly, even though some of the earliest mathematical models of infectious diseases were stochastic, accounting for the chance nature of transmission of infection [6] , much of the applied modeling that followed was deterministic and based on nonlinear differential equations [7] . More recent applied work has, however, recognized the importance of using stochastic epidemic models [8] and also of the development of associated methodology [9] .

The difficulty in mathematically analyzing models that include both stochastic elements and network structure can be a reason for not including these factors, but we prefer to include them, and subsequently to systematically reduce the complexity of the resulting model. This is the approach we adopt in this paper; the reduction process being made possible because of the existence of a separation of timescales: many variables "decaying" away rapidly, leaving a few slow variables, which govern the medium-to long-term dynamics.

Compartmental models of epidemics typically assume that the majority of the population starts susceptible to infection, with a small number infectious, who then spread infection to others before recovering. A key distinction is between susceptible-infectious-susceptible (SIS) dynamics in which recovery returns an individual to the susceptible compartment and susceptible-infectious-recovered (SIR) dynamics in which recovery removes an individual from a further role in the next epidemic [4] , with the former being used to model, e.g., sexually transmitted infections other than HIV and the latter, e.g., pandemic influenza [2, 7, 10] .

In the theoretical analysis of epidemic models, a crucial quantity corresponds to the basic reproductive ratio R 0verbally described as the expected number of secondary cases infected per primary case early in the epidemic [11] . Depending on the value of R 0 , either O(N ) individuals will experience infection in a population of size N -the model is then said to be supercritical-or O(1) individuals will experience infection in a population of size N -the model is then subcritical-thus defining an epidemic threshold.

In reality, the contact network on which the disease spreads will often change over the course of the epidemic; while approaches exist in which the dynamics of both the disease and the network are considered [12] , it is more common to consider two limiting cases. The first of these is a static approach (also called "quenched") in which the network is assumed to evolve much more slowly than the disease and which is typically approached analytically through the use of pair approximation and related techniques [12, 13] . The second is a dynamic limit (also called "annealed" or "discrete heterogeneous") in which the network is assumed to evolve much more quickly than the epidemic, and can be described by an effective network characterized by its degree distribution, in which all individuals sharing the same degree are considered equivalent. This case can be analyzed through use of a set of degree-indexed differential equations (often called the "heterogeneous mean field" approach) provided the maximum degree is not too large [14] [15] [16] .

When the distribution of degrees in the population is highly variable-a situation that appears to be supported empirically [17, 18] -it was recognized that the epidemic may not exhibit straightforward critical behavior. This happens because as the population size N becomes large, extremely small levels of infectiousness can lead to large epidemic sizes [19] or more accurately speaking the critical level of infectiousness can depend very sensitively on the largest degree, K, in the network [20, 21] . The behavior of highly heterogeneous network epidemics near criticality continues to generate interest in both the physics and mathematics literature [4, 15, 22, 23] .

In this paper we investigate the stochastic behavior of heterogeneous network epidemics over time. We study a network in the dynamic limit, characterized by a power-law degree distribution such that the probability of an individual having degree k is given by d k ∝ k α , with −3 < α < −2, although the method of analysis is applicable to other distributions. This type of network has the property that the basic reproductive ratio, which is found to be proportional to the second moment of the degree distribution k 2 , diverges in the limit of infinite populations, leading to the absence of an epidemic threshold. This is evidently not the case when the population-and thus the degree cutoff-is finite, but the second moment of the distribution can still be extremely large for sufficiently large N , and heterogeneity can play an important role in the dynamics of the system.

For the case of large but finite population size, we derive a two-dimensional stochastic differential equation (SDE) approximation to full SIR dynamics, which reduces to an analytically solvable one-dimensional system early in the epidemic. We perform simulations using a power-law degree distribution with a maximum degree cutoff K, which show that our approach provides a good approximation provided K is not too close to the population size. There have, to our knowledge, been no directly comparable studies of this kind. Perhaps those closest have been, first the study of finite-size scaling of the SIS model near the critical point to produce a phenomenological one-dimensional SDE [24] , second, the derivation of a one-dimensional SDE for the SIS model from first principles, but using an adiabatic approximation that would only be expected to hold near the critical point [25] , and finally the study of a four-dimensional SDE for the SIR model on a static network that is much more complex than ours and less amenable to analysis [26] .

The outline of the paper is as follows. In Sec. II we introduce the model, and formulate it first at the level of individuals (the microscale) and then at the level of populations, but retaining a stochastic element (the mesoscale). We begin to explore the reduction of the mesoscopic-level SDEs in Sec. III and complete the process in Sec. IV. In Sec. V we perform a further reduction, which allows us to make contact with a model that has already been discussed in the literature. In Sec. VI these reduced models are used to find the distribution of epidemic sizes, and we conclude with a discussion of our results in Sec. VII. There are two appendices where technical details relevant to Secs. III and IV have been relegated.

In this section we formulate the model first at the microscale, that is, in terms of individuals, and from this derive a mesoscopic version of the model where the variables are continuous and represent the fractions of types of individuals who are infected or who are susceptible to infection. The general procedure to do this is reviewed in Ref. [27] and has been previously applied to models of epidemics on networks [28, 29] , albeit in situations where the nodes of the network contained many individuals, rather than just one, as in the present context.

As discussed in the Introduction, individuals are located at the nodes of a network through which they are exposed to infection from individuals at neighboring nodes. Since we work in the dynamic limit, individuals do not alter their number of contacts when acquiring the infection or recovering from it, and the network plays no role other than encoding the information about how many connections to other individuals a given individual possesses. In practice, this means that the network is not explicitly generated, and it is only characterized through its degree distribution.

An individual is labeled according to whether it is (i) susceptible, infected, or recovered, and (ii) the degree of the node where it is located. Thus, the variables at the microscale are S k , I k , and R k , k = 1, . . . ,K. We will assume that the population is closed, i.e., that S k (t) + I k (t) + R k (t) = N k , with total population N = K k=1 N k , at all times, so that we can remove one of the variables:

Two types of events can occur: the infection of a susceptible individual, whenever an individual of type I comes in contact with one of type S-contacts between individuals i and j having degrees k i and k j , respectively, are governed by a Poisson process with rate βk i k j ; and the recovery of an infectious individual with rate γ . In this paper we will be interested in the properties of an epidemic that takes place over a much shorter time frame than demographic processes and hence will ignore the birth and death of individuals, leaving the only interactions between the individuals to be:

(1) Infection of an individual at a node of degree k. The transition rate for this process is given by

where the summation is over infected individuals labeled according to the degree of the node on which they are located. The only arguments of T 1 , which are shown explicitly, are those involving individuals on nodes of degree k, since only these change in this process; those on the right represent the initial state, and those on the left represent the final state.

(2) Recovery of an infected individual at a node of degree k. This proceeds at a rate given by γ , and so the transition rate for this process is given by

The dynamics of the process can be described by writing down an equation for the probability that the system is in a state {S l ,I l |l = 1, . . . ,K} at time t, which will be denoted by P ({S l ,I l },t). This is the master equation [30] , which equates dP /dt to the net rate of increase of P due to the processes given by Eqs. (1) and (2) , and is given by

[T 1 (S k ,I k |S k + 1,I k − 1)P (S k + 1,I k − 1,t) + T 2 (S k ,I k |S k ,I k + 1)P (S k ,I k + 1,t)]

[T 1 (S k − 1,I k + 1|S k ,I k ) + T 2 (S k ,I k − 1|S k ,I k )]P (S k ,I k ,t).

Once again, the dependence of the probability distribution, P , on the state variables {S l ,I l |l = k} has been suppressed in the sum over k on the right-hand side of this equation. We will refer to the discrete model defined by Eqs. (1)-(3) as the full microscopic model, since we shortly derive a mesoscopic model from it, and subsequently further reduce this model. Once an initial state has been specified, the model is completely determined, and in principle P ({S l ,I l },t) can be found for any state {S l ,I l |l = 1, . . . ,K} at any time t. This behavior can be explored numerically via Monte Carlo simulations (e.g., using Gillespie's method [31, 32] or other approaches [16] ); however, our main interest in this paper is in approximating the master equation to obtain models that are more amenable to analysis.

A mesoscopic description of the process involves going over to new variables s k = S k /N and i k = I k /N, k = 1, . . . ,K, which are assumed continuous for large N , but retaining the stochastic nature of the system by keeping N large but finite. A macroscopic description, on the other hand, would involve taking the limit N → ∞, which renders the process deterministic, and which eliminates all effects of the original discrete nature of the system.

The construction of the mesoscopic model involves expanding the master equation, Eq. (3), in powers of 1/N and truncating this expansion at order O(N −2 ). This is discussed in many places in the literature, but here we follow Ref. [27] , which gives explicit forms for the functions A (k) μ and B (k) μν , μ,ν = 1,2, which define the mesoscopic model. Carrying out this procedure, the master equation becomes a Fokker-Planck equation [33] of the form

where

We have also introduced x = (x (1) , . . . ,x (K) ). The explicit form of the functions A (k) μ and B (k) μν are [27] A (k)

where

are the original transition rates (scaled by N) given in Eqs. (1) and (2) , but written in terms of the continuous state variables s k and i k . The final stage of the construction of the mesoscopic form of the model is to adopt a coarser time scale by introducing a new time variable τ = t/N. From Eq. (4) it can be seen that without this change of variable, the macroscopic (N → ∞) limit eliminates the right-hand side of the equation, giving a trivial macroscopic limit. With the change of variable, the first term on the right-hand side survives (but without the factor N −1 ), giving a Liouville-like equation for P , which implies a nontrivial deterministic dynamics.

A clearer way to see this, and in fact a more intuitive formulation of the mesoscopic dynamics to that of the Fokker-Planck equation, is to use the SDE, which is equivalent to Eq. (4). We utilize the general result that a Fokker-Planck equation of the form given by Eq. (4) is equivalent to the SDE dx (k)

where the η (k) μ are Gaussian random variables with zero mean and a correlation matrix given by the matrix B (k) . In terms of the variables s k and i k , the SDEs take the form

where

Here the noise is to be interpreted in the sense of Itō [33] . It is clear from Eqs. (7) and (8) that the limit N → ∞ eliminates the noise terms leaving deterministic equations, and we recover Eq. (2.6) from Ref. [16] . It is convenient to make the noise in the SDEs explicitly multiplicative. Since the rates

, μ = 1, 2, and then transform to new noise variables,

It is straightforward to check that ρ (k) μ (τ )ρ (l) ν (τ ) = δ μν δ kl δ(τ − τ ), and therefore that all the state dependence has been transferred from the correlator to the SDEs. Equations (7) and (8) correspond to the full mesoscopic model.

In order to describe the evolution of the entire population, we need to consider the SDEs for all k = 1, . . . ,K; that is, we have to work with 2K equations. In this paper, we will be interested in networks where individuals can pick degrees from a truncated Zipf distribution of the form

and we let 1 K N − 1. Since we are interested in the case when K is large, the number of equations making up the full mesoscopic model are so large that it can make any direct treatment impractical. Therefore, in the following sections we will be concerned with reducing the number of equations, and we start by looking at the deterministic limit of the system. This will also have a crucial role to play in the reduction of the stochastic model, described in Sec. IV.

The deterministic limit that controls the dynamics of the macroscopic version of the model is found by taking the N → ∞ limit of Eqs. (7) and (8), and so eliminating the noise terms. For notational convenience we reorder the entries of the state vector, so that x(τ ) = (s 1 (τ ), . . . ,s K (τ ),i 1 (τ ), . . . ,i K (τ )), and write Eqs. (7) and (8) in the form d x/dτ = A(x). To find the fixed points of the dynamics we set the time derivatives on the left-hand side of these equations to zero; these fixed points will be denoted by an asterisk and obey A(x * ) = 0. It is immediately clear that there is a set of stable fixed points with i * k = 0 and s * k undetermined, for all k, corresponding to a disease-free state.

For a given initial condition, the fixed point is uniquely determined, however, it may not be stable. To investigate the stability of the fixed point (s * 1 , . . . ,s * K ,0, . . . ,0) we perform a linear stability analysis. 2 and is found to have the form

where J is a K × K matrix with entries J kl = βs * k kl and I is the K × K unit matrix. The eigenvalues and eigenvectors of J {2K} are straightforward to determine. Any vector that has the last K entries equal to zero is an eigenvector with eigenvalue zero, reflecting the degeneracies of the system. Similarly, a vector with the first K entries set equal to zero and the last K entries equal to ψ k , k = 1, . . . ,K, is an eigenvector with eigenvalue −γ if K l=1 J kl ψ l = 0, that is, if K l=1 lψ l = 0. These then yield another degenerate set of (K − 1) eigenvalues −γ . The final eigenvalue can be found by equating the trace of J {2K} to the sum of its eigenvalues. In summary, the 2K eigenvalues are

while a suitable set of corresponding eigenvectors are listed in Appendix A. From Eq. (13) we see that the condition for the fixed points to be (marginally) stable is

As mentioned in the Introduction, there are two possible outcomes for the deterministic dynamics depending on the parameter values: either a negligible or a significant number of individuals in the population will get infected during the course of the epidemic. The quantity that determines which of these cases we are in is often referred to as the basic reproductive ratio; the final size of the epidemic, in turn, can be obtained from the s coordinates of the fixed point, s * 1 , . . . ,s * K . In order to determine these in the general heterogeneous case, we look first at the homogeneous case, where K = 1. If we set an initial condition with only a few infective individuals in an almost completely susceptible population, we will have i(τ = 0) = i 0 1 and s(τ = 0) = s 0 ∼ 1, where s = s 1 and i = i 1 ; with this, the N → ∞ limit of Eq. (8) takes the form di/dτ = (β − γ )i, early in the epidemic. From this we obtain

where R 0 = β/γ corresponds to the basic reproductive ratio. This defines an epidemic threshold, such that for R 0 > 1 the number of infectives grows exponentially in the early stages of the epidemic, while for R 0 < 1 the epidemic does not take off and its final state corresponds to (s * ,i * ) = (s 0 ,i 0 ) ∼ (1,0). For the case when R 0 > 1, we can find the s coordinate of the fixed point by dividing Eq. (8) by Eq. (7) with K = 1; this yields

If we set again an initial condition i 0 1, with s 0 ∼ 1, we can integrate the equation above to find

where we have used that the final state of the epidemic corresponds to i * = 0. Therefore, to a very good approximation, we can find s coordinate of the fixed point as the smallest solution to

note that s * = 1 is always a solution, and the only one for R 0 < 1. The total epidemic size will be given by r * = 1 − s * , which corresponds to the fraction of recovered individuals at the end of the outbreak. This is due to the fact that the initial number of recovered individuals is zero and that every individual that gets infected during the course of the epidemic will eventually recover and remain permanently recovered. Now, a significant simplification can be made in the deterministic limit of the heterogeneous case, which reduces the 2K dynamical equations to two equations (see Refs. [34] [35] [36] , which consider the static case). This is achieved by use of the Ansatz s k (τ ) = d k θ (τ ) k ; we also define λ(τ ) such that λ(τ ) ≡ k ki k (τ ). With these two variables, Eqs. (7) and (8), in the case when N → ∞, can be rewritten independently of k as

where we have introduced φ(x) = K k=1 k 2 d k x k . While the new variables θ and λ are justified mathematically through the fact that they simplify the equations, they also have a physical interpretation. For the case of θ , this is related to the recently introduced concept of a test node [12, 37] , which plays an analogous role in epidemic dynamics to the concept of a test charge in electrostatics. If we introduce a single initially susceptible individual with one link into the system, and assume that the population is so large that this individual has negligible influence on the epidemic dynamics, then θ (τ ) is the probability that such a "test node" has avoided infection at time τ . The quantity λ is related to the force of infection that is defined in epidemiology as the per-capita rate at which susceptible individuals become infective [10] . Our λ(τ ) is the network generalization of this, being instead the per-link rate at which susceptible individuals become infective at time τ . We therefore argue that where a nonnetwork SIR epidemic can have its dynamics represented in phase space with coordinates (s(τ ),i(τ )), for the network model the physically relevant coordinates, designed to account for link heterogeneity, are (θ (τ ),λ(τ )).

If we also introduce the probability generating function of the degree distribution G(

. In this reduced model, there is a line of stable fixed points, uniquely determined by the initial conditions, which satisfy λ * = 0. The stability matrix for these fixed points is the 2 × 2 matrix

which has a similar form to the stability matrix for the full system given by Eq. (12) . Once again, an eigenvector exists that has the last entry equal to zero. This corresponds to an eigenvalue zero, reflecting the existence of a line of fixed points. The trace of J {2} then gives the second eigenvalue, and the two eigenvalues are thus

with the eigenvectors given in Appendix A. The stability condition is now

Now we can find analogous expressions to those in Eqs. (15) and (18) for the heterogeneous case. We start by noting that, in terms of θ , the total number of susceptible individuals will be given by

This means that a completely susceptible population will correspond to θ being equal to unity. A small initial number of infectives will, in turn, correspond to λ 1. In addition from Eq. (20) , dλ/dτ = γ (ρ − 1)λ, where ρ = βφ(1)/γ = β k 2 /γ , at early times. Therefore,

From Eq. (25) it is clear that there is a critical value of ρ equal to 1; if ρ is larger than 1 the infection grows, and if it is below 1 it does not. Thus, ρ is the heterogeneous basic reproductive ratio, generalizing R 0 = β/γ in the homogeneous case. Central to our investigation is the fact that this quantity does not have a deterministic threshold for networks characterized by power-law degree distributions in the case of interest, i.e., when −3 < α < −2.

In order to find the heterogeneous final size, we proceed again as in the homogeneous case; dividing Eq. (20) by Eq. (19), we find

From this, the heterogeneous analog of Eq. (17) is given by

where we again set an initial condition with a very small number of infectives, λ 0 1, in an almost completely susceptible population, θ 0 ∼ 1, and the final state is given by λ * = 0. The θ component of the fixed point will then correspond to the smallest solution to

with ρ the heterogeneous basic reproductive ratio defined above. The epidemic size, by definition, will be given by r * = 1 − G(θ * ). Furthermore, the s k components of the fixed point correspond to s * k = d k θ * k . Our interest in this paper is in the (finite N ) stochastic dynamics of the model, where a complete reduction of the kind that led to Eqs. (19) and (20) is not possible. In order to see this, we proceed in the same manner as in the deterministic case: we start from Eq. (7), this time without taking the N → ∞ limit, and use the Ansatz s k (τ ) = d k θ (τ ) k ; summing over k and dividing by G (θ ), we arrive at

where we have introduced ξ 1 

. Since the noises ρ (k) 1 (τ ) have zero mean and are independent of each other, we can rewrite ξ 1 (τ ) as

where ζ 1 (τ ) has zero mean and unit variance,

With this, we arrive at a reduced equation for θ given by

As in the deterministic case, then, we can express the K equations for the s k variables, described by Eq. (7), entirely in terms of θ and λ via Eq. (31). However, we will see that the same does not hold for the dynamics of the i k variables, as we have anticipated.

In a similar fashion as with the s k variables, we take the equation for the i k variables, Eq. (8), multiply it by k and sum over k, to obtain

Let us now denote the noise terms appearing in the equation for λ above as ξ 2 (τ ) = − k kσ (k) 1 ρ (k) 1 (τ ) and ξ 3 (τ ) = k kσ (k) 2 ρ (k) 2 (τ ). These two noise terms, like ξ 1 (τ ), have zero mean; also, it can be readily verified that the correlations between all three of them are given by ξ μ (τ )ξ ν (τ ) = δ(τ − τ )B μν , where now μ,ν = 1, 2, 3, with

We can now, as in Sec. II, express ξ 1 (τ ),ξ 2 (τ ) and ξ 3 (τ ) in terms of three independent normally distributed random variables, ζ 1 (τ ), ζ 2 (τ ), and ζ 3 (τ ), with zero mean and ζ μ (τ )ζ ν (τ ) = δ(τ − τ )δ μν , μ, ν = 1, 2, 3. The definition of ζ 1 (τ ) in terms of ξ 1 (τ ) has already been given by Eq. (30) ; the other two are 

The main feature we note from this is that the twodimensional system of equations cannot be closed. This is so becauseσ 3 is still dependent on i k for all values of k, and so it is not possible to achieve a complete reduction of the stochastic system. Therefore, when we take the intrinsic noise present in the system into account, a partially reduced model of (K + 1) equations-Eq. (31) for θ , with λ = k ki k , and Eq. (8) for i k -is the best we can achieve.

Due to this, the exploration of the epidemic dynamics for extremely heterogeneous cases, in which K can take very large values, is rendered impractical: even after reducing the number of equations, we will still be left with a very-high-dimensional system. However, progress can be made by observing that, under certain conditions, the deterministic limit of the model presents a separation of time scales, as we shall see in the following section.

We will end this section by carrying out this partial reduction in the deterministic case, as a preparation for the stochastic treatment. This partial reduction involves applying the above reduction to the s k variables only, by writing s k (τ ) = d k θ (τ ) k , but keeping the K i k variables. There are now (K + 1) dynamical equations, which take the form

There is again a line of fixed points corresponding to diseasefree states, i k = 0, k = 1, . . . ,K, but with θ undetermined. The (K + 1) × (K + 1) stability matrix is now

where j is a K-dimensional vector whose lth entry is βθ * l, J is a K × K matrix with entries J kl = βd k (θ * ) k kl, and I is the K × K unit matrix. The eigenvalues and eigenvectors can be found in the same way as before, the eigenvalues having the 

As we have discussed, the deterministic limit of Eqs. (7) and (8) can be greatly simplified by performing a change of variables and describing the evolution of the system in terms of θ and λ. The inclusion of demographic noise, however, hinders such a reduction of dimensionality. It is not possible to obtain a closed two-dimensional system of equations, thus enormously complicating the analysis of the stochastic epidemic, compared to its deterministic limit.

Here, we attempt to overcome this issue by exploiting the properties of the deterministic limit of the model and by finding a way to close the equations for θ and λ in the finite-N case, while still retaining the essential features of the full model. In order to do this, we start from Eqs. (36) and (37) for the partially reduced (K + 1)-dimensional deterministic system for θ and i k . Examination of the eigenvalues appearing in Eq. (39) shows that | {K+1} | < | {2,...,K} |. If the ratio ≡ | {K+1} / {2,...,K} | is small, then there is the potential to carry out a reduction of the system. In such a case, the deterministic dynamics would quickly evolve along the directions corresponding to the set of most negative eigenvalues toward a lower-dimensional region of the state space, in which the system as a whole evolves in a much slower time scale, and which contains the fixed point (θ * ,0, . . . ,0). As illustrated in Fig. 1 , showing the value of for different choices of K, we would expect this time-scale separation to occur for combinations of parameters for which the final size of the epidemic is small, i.e., for large θ * .

Since there is one most negative eigenvalue, with a degeneracy equal to (K − 1), after sufficient time has passed the evolution of the system will then be constrained to a two-dimensional surface, which we will refer to as the slow subspace. With the aim of reducing the dimensionality of the problem, then, we ignore the fast initial behavior of the system and focus the analysis on the slow surface. We do so by imposing the condition that no dynamics exists along the (K − 1) fast directions, i.e.,

where x = (θ,i 1 , . . . ,i K ), and u {j } is the j th left-eigenvector of the Jacobian J {K+1} evaluated at the fixed point (θ * ,0, . . . ,0), with θ * given by Eq. (28) . This procedure goes under the name of adiabatic elimination, or fast-variable elimination [38] , and it is a common practice in the simplification of deterministic nonlinear dynamical systems-e.g., systems of oscillating chemical reactions [39, 40] . For ease of application of the method, the left-and righteigenvectors of the Jacobian are normalized in such a way 

these are given, respectively, by Eqs. (A4) and (A6). Imposing the condition given by Eq. (40), one finds that A k+1 = Y k K l=1 lA l+1 , where k = 2, . . . ,K and Y k is given by Eq. (A5). In terms of the θ and i k variables, this reads

and this gives the behavior of i k , which limits the dynamics of the deterministic system to a slow subspace. A more convenient form can be found by multiplying Eq. (42) by k and summing from k = 2 to k = K. This allows us to relate K l=1 li l to i 1 to find the following (K − 1) expressions for i k :

It should be stressed that the (K − 1) Eqs. (43) , which we may write in the form i k (τ ) = F k (θ (τ ),i 1 (τ )), k = 2, . . . ,K, should be interpreted as the time-evolution that the i k need to follow, so that d x/dτ = 0 in the directions v {k} , k = 2, . . . ,K. This follows from the deterministic equations of motion d x/dτ = A(x) and from Eq. (40) . Thus, the only temporal evolution is in the v {1} and v {K+1} directions.

The functional forms i k (τ ) = F k (θ (τ ),i 1 (τ )), with F k defined by Eq. (43), formally only hold sufficiently near the fixed point that linearization is accurate. However, in Fig. 2 we show three-dimensional projections of the deterministic dynamics from Eqs. (36) and (37) for an epidemic on a network characterized by a power-law degree distribution with α = −2.5 and maximum degree K = 10, with heterogeneous basic reproductive ratio ρ = 1.2, which gives ≈ 0.184. It appears that the solid red line, which represents the deterministic dynamics, lies in the slow subspace for most of the time, not just at late times when the system approaches the fixed point. That this is the case can be seen from a consideration of the Jacobian of the system at an arbitrary time, and not just at the fixed point: (37); blue, yellow, green lines: stochastic differential Eqs. (31) and (8) . The initial condition corresponds to i 1,0 = 0.005, i k,0 = 0 for k = 2, . . . ,K, and G(θ 0 ) = 1 − i 1,0 . In both cases, the system is pushed back from the initial condition (lower right in both panels) toward the slow subspace described by Eq. (43), on which it seems to stay during the rest of the course of the epidemic. Note the different scales in the z axes of both panels, so the fluctuations away from the slow subspace are small in both cases; the fluctuations on the slow subspace can, however, be large, as seen from the departure of the stochastic trajectories from the deterministic dynamics.

It is straightforward to check that (0,ψ 1 , . . . ,ψ K ) is a right-eigenvector of the matrix in Eq. (44) with eigenvalue −γ , if K l=1 lψ l = 0. Therefore, remarkably, the eigenvectors v {k} of the Jacobian at the fixed point, defined in Eq. (A4), are eigenvectors of the Jacobian at all times. The existence of this negative eigenvalue with degeneracy (K − 1) means that the system remains in the plane described by the vectors v {1} and v {K+1} , since any small deviations out of the plane will tend to be pushed back again. Normally, when carrying out a fastvariable elimination, one would find that the system possesses a line of fixed points, which is quickly approached at early times, dominated by the deterministic dynamics, and along which the slow, stochastic dynamics takes place. In this case, however, the line of fixed points corresponds to an absorbing boundary of the system: the evolution of the epidemic stops altogether once it reaches the disease-free state. Therefore, it is important that the projection onto the slow variables constitutes a good approximation to the dynamics also far from the line of fixed points.

We parametrize the slow subspace using the coordinates z 1 ≡ u {1} · x and z 2 ≡ u {K+1} · x. Using Eq. (A6) and x = (θ,i 1 , . . . ,i K ), we obtain

A linear transformation of Eqs. (45) and (46) shows that an equivalent set of coordinates are θ and k ki k , previously called λ. So, in summary, in Sec. III we showed the system of 2K equations leads to a closed set of equations in θ and λ. Here, through a consideration of the slow modes of the deterministic equations, we have shown that not only can we describe the system in terms of θ and λ, but that it remains in the vicinity of a subspace that can be parametrized by these variables.

We would like to carry this reasoning over to the stochastic case, to find a reduced set of SDEs, but which are closed, i.e., which can be expressed entirely in terms of the slow variables z 1 and z 2 (or θ and λ). Before starting the actual treatment of the model, we visualize its dynamics to find whether the system evolves toward a lower-dimensional subspace, as it does in the deterministic case. Figure 2 also shows three-dimensional projections of stochastic trajectories from Eqs. (31) and (8), for a total population of N = 10 5 . The behavior is found indeed to be similar to that of the deterministic case: the dynamics seems to quickly reach, and then fluctuate around, the deterministic slow subspace, with the early behavior dominated by the deterministic dynamics. Therefore, we expect that we should be able to effectively reduce the stochastic system from (K + 1) dimensions to two dimensions, by neglecting the fast initial behavior of the epidemic.

We will do this by projecting the stochastic differential equations, Eqs. (31) and (8), onto the slow degrees of freedom of the deterministic limit of the model; that is, we will neglect the fluctuations away from the slow subspace, along the (K − 1) fast directions corresponding to the most negative eigenvalue. This can be formulated mathematically through the application of a projection operator that only picks the components of x, A(x), and η(τ ) along the slow eigenvectors of the deterministic limit, v {1} and v {K+1} [41] :

The details of the derivation are given in Appendix B, where it is shown that the SDEs given by Eqs. (31) and (35) of Sec. III are recovered, except that nowσ 3 ≡ B 33 is determined in terms of θ and λ. This closed system of stochastic differential equations for θ and λ, which we call the reduced mesoscopic model, constitutes the main result of this paper and so we reproduce it here:

We now compare the dynamics of the full microscopic model, defined by Eqs. (1)-(3), to that of the reduced mesoscopic model above. Figures 3 and 4 show the time series of the full microscopic model in terms of the variables θ and λ, and that of the reduced mesoscopic model, for N = 20 000 and different values of K having the same deterministic final size θ * , obtained from Eq. (28); Fig. 5 , in turn, compares the phase diagrams. We see that, in the case of small K, the time series of both variables in the full microscopic model are rather smooth, and are dominated by a horizontal spread due to the fact that the epidemic takes off at different moments. As heterogeneity increases, the dynamics of the system becomes noisier and noisier, as seen in the phase diagrams; however, the time series of θ continues to be smooth. This observation will be useful in the following section.

From the time series of the reduced mesoscopic model we see that, except for small K, the temporal behavior of θ and λ does not present the horizontal spread observed in the full microscopic model. However, once the epidemic has taken off, the dynamics of the full microscopic model is correctly captured when K is not so large; this can be seen from the phase diagrams. We note that the relative magnitude of the eigenvalues for this choice of parameters takes the values ≈ 0.44-0.65, and that the reduced mesoscopic model fares well despite being rather large. We find slightly worse agreement between the full and reduced models for smaller values of θ * , corresponding to larger . One further detail to add is that in the case of the reduced model there are far fewer trajectories in which the epidemic does take off, compared to the full microscopic model; that is, the reduced model results in an over-representation of very early extinctions.

Note that, for the value of N employed in these figures, K = 19 999 is the maximum possible degree cutoff. In this case, we would not expect the reduced model to yield accurate results. This is due to the fact that not even the full 2Kdimensional mesoscopic model [Eqs. (7) and (8)] provides a good approximation to the microscopic dynamics in this case, since K is of the order of N and this would need to be explicitly taken into account when performing the mesoscopic expansion leading to the Fokker-Planck Eq. (4). The implementation of the reduced mesoscopic model, however, is much more practical than that of the full mesoscopic model: while we have found that simulations of the latter take a time that, as its number of equations, grows linearly with the degree of heterogeneity of the network, for the reduced version this time, besides being much shorter, appears to scale as √ K. We end this section with some further observations regarding the reduced stochastic system. First, the existence of a (K − 1)-fold degenerate eigenvalue equal to −γ throughout the time evolution also means that stochastic trajectories, as well as deterministic ones, which venture out of the slow subspace will tend to be pushed back into the plane. Thus, we would expect the reduced set of SDEs to be a good approximation throughout the motion, and not just at late times. This may be the reason why the approximation works well even if is not particularly small. Second, in retrospect, the elimination of the fast variables in the stochastic system is not strictly necessary, since the SDEs previously obtained in Sec. III are recovered, andσ 3 could be found in terms of θ and λ using Eq. (43) obtained from the elimination of fast variables in the deterministic system. One can compare the value ofσ 3 in terms of the i k to the one in terms of λ in the deterministic limit and see that, indeed, the approximation is not very good for moderate to large values of . This does not seem to play a huge role when considered as one of many terms in the stochastic evolution of λ.

Finally, we can make a general comment about the nature of the fluctuations as the epidemic starts to grow from its initial state. In this case the relevant eigenvalues will be those of the Jacobian Eq. (44), but with i k = 0 and θ = 1, and will therefore be given by those in Eq. (39) , but with θ * set equal to 1. Therefore, the largest eigenvalue is now (ρ − 1)γ . In this case, if ρ > 1, so that an epidemic starts to grow, the dynamics is pushed away from the line (θ,0, . . . ,0) in its early stages, as can be seen from Fig. 2 , where the trajectories leave the vicinity of the absorbing state after reaching the slow subspace. 

In the previous section we derived a reduced twodimensional system with three white-noise processes, the reduced mesoscopic model Eqs. (48) and (49), which provides a good approximation to the dynamics of the full microscopic model, provided the maximum degree individuals can pick is not so large. This is true even for relatively large values of , corresponding to small separation between the eigenvalues of the system. We have also noted that the time series of θ , obtained from the full microscopic model, presents very smooth behavior over the range of parameters explored. Inspired by this observation, we attempt to further simplify the model by noting that Eqs. (48) and (49) have the form

Taking into account the fact that, in any given realization of the process, the lowest value reached by θ will be roughly around θ * , we compare the magnitudes of the noise intensities f 1 , . . . ,f 4 between θ = θ * and θ = 1 for fixed θ * and different values of K. We find that, in general, and especially during the early stages of the epidemic when θ ≈ 1, f 1 is small compared to f 2 , f 3 , and f 4 -see Fig. 6 ; the effect is more pronounced for more heterogeneous distributions with smaller deterministic final sizes, i.e., larger θ * . Therefore, the system may be described simply by a deterministic ordinary differential equation (ODE) for θ , and Eqs. (48) and (49) become

with ζ (τ ) a Gaussian noise with zero mean and ζ (τ )ζ (τ ) = δ(τ − τ ), with the noise intensity,σ , given bȳ We call this system the semideterministic mesoscopic model, which is two-dimensional like the reduced mesoscopic model, but has only one white-noise process in place of the latter's three. Here we have used that ζ μ (τ ) = 0, and ζ μ (τ )ζ ν (τ ) = δ μν δ(τ − τ ), μ,ν = 1,2,3. Now, we note that Eq. (54) corresponds to a Cox-Ingersoll-Ross (CIR) model [42] ,

where θ (τ ) is a time-dependent parameter that is the solution of the deterministic differential Eq. (53), b = 0, and

For fixed a and σ , the conditional probability distribution of λ(τ + τ ) = λ, given λ(τ ) = λ , is known in closed form, and corresponds to a noncentral chi-square distribution given by [42, 43] 

where

and where I q is the modified Bessel function of the first kind of order q. Therefore, if we assume that θ evolves discretely, staying constant from time τ to time τ + τ for reasonably small τ , we can use the result above to simulate Eqs. (53) and (54) by sampling the values of λ directly from the distribution in Eq. (59) with fixed parameters a[θ (τ )], σ [θ (τ )] at every time step. The case b = 0 corresponds to the special case of a noncentral chi-square distribution with zero degrees of freedom; thus, in order to simulate Eq. (54) we follow the method described in Refs. [44, 45] and sample from a central chi-square distribution with Poisson-distributed degrees of freedom, with the convention that the central chi-square distribution with zero degrees of freedom is identically zero [46] . Figure 7 shows the dynamics obtained in this way. While we note an additional over-representation of early extinctions, the other results are in very good agreement with the reduced mesoscopic model, Eqs. (48) and (49).

In the previous sections, we have derived a two-dimensional reduction of the 2K-dimensional SIR dynamics described by Eqs. (7) and (8), by exploiting a separation of time scales in the deterministic dynamics of the system. We have found reasonably good agreement between the temporal behaviors of the full and reduced models. Furthermore, we have also discussed an extra simplification, consisting in neglecting the noise acting on the susceptible variable θ , and describing the system in terms of one ODE and one SDE instead of two SDEs. Now we explore how accurately these reductions capture the distribution of epidemic sizes, given by the number of recovered individuals present in the system at the end of the epidemic, r ∞ , which can be obtained as r ∞ = 1 − k s k,∞ and r ∞ = 1 − G(θ ∞ ) in the full and reduced models, respectively, where s k,∞ and θ ∞ are the final values of s k and θ -not to be confused with the fixed points in the deterministic limit, s * k and θ * . Formal mathematical work has considered central limit theorems for models that generalize ours [47] ; we are interested here in explicit calculation and simulation approaches.

We obtain the distribution of r ∞ from the full model via a Sellke construction [16] , and note that it is characterized by two regions: one to the left of the interval, corresponding to early extinctions, and one to the right corresponding to the epidemic taking off. The first thing we note when comparing this to results from the two versions of the reduced model-Eqs. (48) and (49), Eqs. (53) and (54)-is that, as discussed previously, there is an over-representation of very early extinctions in the case of the latter, and therefore the relative sizes of the leftand right-most regions of the distribution are not well captured. We would not expect our diffusion approximation to capture extinction probabilities accurately, given that its derivation implicitly assumes that the infective population is large and extinctions by definition involve small numbers of infectives. We can, however, derive accurate extinction probabilities from other results we have obtained.

First, we make an argument related to "susceptibility sets" in probability [48] . For each pair of individuals in the population, we pick whether infection will be spread between them if one becomes infective and draw a link between them if it will. This forms a static graph, and if we pick a node of the graph to be initially infective, then the size of the epidemic will be the size of the network component that the initial pick is part of. This means that the size of the giant component of the graph is equal to the final size of the epidemic if it is large, and the probability of early extinction is equal to the proportion of the population in the small components of this graph. Therefore, we can conclude that we have already calculated the probability of early extinction of the epidemic starting with an individual picked uniformly at random from a very large population; it will be equal to k s * k = G(θ * ). Figure 8 shows the distribution of final sizes from the reduced and semideterministic mesoscopic models, adjusted to account for the inaccuracy in early extinctions by removing these and assuming that a fraction G(θ * ) of epidemics stop early, and a fraction 1 − G(θ * ) take the values from the mesoscopic models. These results are compared to the distribution from the full model and show that the argument we have introduced provides a simple method to adjust for the inaccuracies in early extinctions, yielding an accurate approximation overall.

Stochastic epidemics on heterogeneous networks involve two potentially large parameters. The first of these is the size of the population, N , which is proportional to the dimensionality of the microscopic description. For large populations, this motivates the derivation of a mesoscopic model based on a diffusion approximation whose errors are controlled by inverse powers of N , and which we have formulated in this paper. The second of these is the maximum degree K, which is proportional to the dimensionality of the full mesoscopic model. We have been able to use fast-mode elimination techniques to derive two-dimensional approximations to this scheme: the reduced mesoscopic model, which involves three white-noise processes; and the semideterministic mesoscopic model, which involves only one white-noise process.

Our numerical results back up our theoretical understanding that these low-dimensional systems are accurate approximations to the microscopic model, except in the case where K approaches N , although this is to be expected since our derivation of the diffusion limit implicitly assumes that K N . Nevertheless, the reduction can be applied to cases with extremely large values of K, as long as they are not too close to N , thus significantly simplifying the original model. We have observed that the elimination of fast variables in the system provides a good approximation to the full dynamics even for relatively large values of , corresponding to a small separation between time scales. We have also shown that the main inaccuracy of the reduction, relating to very early extinction probabilities of the epidemic, can be dealt with using a simple argument. In this way, we hope to have provided a useful tool in the study of both infectious diseases and spreading processes on heterogeneous networks.

[ξ 2 (τ ) + ξ 3 (τ )] (B5)

βθ * −γ + βφ(θ * ) σ 1 (z)ζ 1 (τ ) + βθ * −γ + βφ(θ * ) [σ 2 (z)ζ 2 (τ ) +σ 3 (z)ζ 3 (τ )] ,

where we obtain, from Eqs. (43) , (45) , and (46), 

χ (z)σ 1 (z)ζ 1 (τ ) +σ 2 (z)ζ 2 (τ ) +σ 3 (z)ζ 3 (τ ) .

We now formθ andλ as linear combinations ofż 1 andż 2 , using Eqs. (B1) and (B2). These are equal to the appropriate linear combinations of Eqs. (B3), (B4), (B6), and (B10). One finds that Eqs. (31) and (35) are exactly recovered, except that now σ 3 ≡ B 33 is determined, given by Eq. (B9).

The Mathematical Theory of Epidemics

Infectious Diseases of Humans

Modeling Infectious Diseases in Humans and Animals

Proc. Natl. Acad. Sci. USA

Proceedings of the 16th Annual ACM-SIAM Symposium on Discrete Algorithms

Proc. Natl. Acad. Sci. USA

Stochastic Processes in Physics and Chemistry

Markov Processes: An Introduction for Physical Scientists

Handbook of Stochastic Methods

Simulation of square-root processes

Susceptibility sets and the final outcome of stochastic SIR epidemic models

We thank Frank Ball for helpful comments on this manuscript, as well as the anonymous reviewers. C.P. 

where e i is the ith unit vector and(A2)

A set of right-eigenvectors corresponding to the eigenvalues from Eq. (22) are given by

The partially reduced system with eigenvalues given by Eq. (39) has a set of right-eigenvectors given bywhere nowIn Eqs. (A1) and (A3) we did not normalize the eigenvectors. However, in this case we do need the eigenvectors to be normalized to simplify the application of the reduction technique described in Sec. IV. This requires finding the corresponding left-eigenvectors. When finding these, we have a complication given by the fact that there is one eigenvalue, which is (K − 1)-fold degenerate, and the condition Eq. (41) will not necessarily hold in general. However, if we construct a matrix V in such a way that its ith column corresponds to v {i} , we can choose u {j } to be the j th row of the inverse matrix V −1 , if such an inverse exists. We find that indeed this is the case, and the left-eigenvectors in the partially reduced model are then given byWe note thatUsing this property, we can verify that the orthogonality condition Eq. (41) is satisfied:

In this Appendix we derive the two-dimensional reduced mesoscopic model, by projecting the dynamics of the (K + 1)dimensional partially reduced mesoscopic model onto the slow degrees of freedom of its corresponding deterministic limit.Using the projector defined in Eq. (47) on x = (θ,i 1 , . . . ,i K ) gives the slow variables defined by Eqs. (45) and (46) . Later, the inverse transformation, back to the more familiar variables θ and k ki k (= λ), will be required. It ishowever, for the moment we will use the z 1 and z 2 variables, and will denote the right-hand side of Eq. (B1) as χ (z). The deterministic slow dynamics, obtained by applying the projector to A(x), will be given bẏwhere the dot represents the time derivative. Applying now the projector to the noise terms in the finite-N case we find, for the component along v {1} ,1 ρ (k) 1 (τ ),−σ (1) 1 ρ (1) 1 (τ ) + σ (1) 2 ρ (1) 2 (τ ), . . . ,−σ (K) 1 ρ (K) 1 (τ ) + σ (K) 2 ρ (K) 2 (τ )