key: cord-326631-7gd3hjc3
authors: Ma, Junling; Earn, David J. D.
title: Generality of the Final Size Formula for an Epidemic of a Newly Invading Infectious Disease
date: 2006-04-08
journal: Bull Math Biol
DOI: 10.1007/s11538-005-9047-7
sha: 
doc_id: 326631
cord_uid: 7gd3hjc3

The well-known formula for the final size of an epidemic was published by Kermack and McKendrick in 1927. Their analysis was based on a simple susceptible-infected-recovered (SIR) model that assumes exponentially distributed infectious periods. More recent analyses have established that the standard final size formula is valid regardless of the distribution of infectious periods, but that it fails to be correct in the presence of certain kinds of heterogeneous mixing (e.g., if there is a core group, as for sexually transmitted diseases). We review previous work and establish more general conditions under which Kermack and McKendrick's formula is valid. We show that the final size formula is unchanged if there is a latent stage, any number of distinct infectious stages and/or a stage during which infectives are isolated (the durations of each stage can be drawn from any integrable distribution). We also consider the possibility that the transmission rates of infectious individuals are arbitrarily distributed—allowing, in particular, for the existence of super-spreaders—and prove that this potential complexity has no impact on the final size formula. Finally, we show that the final size formula is unchanged even for a general class of spatial contact structures. We conclude that whenever a new respiratory pathogen emerges, an estimate of the expected magnitude of the epidemic can be made as soon the basic reproduction number ℝ(0) can be approximated, and this estimate is likely to be improved only by more accurate estimates of ℝ(0), not by knowledge of any other epidemiological details.

Whenever a serious infectious disease emerges or re-emerges in a human population, a matter of immediate interest is the likely magnitude of the outbreak. This is often called the expected final size of the epidemic (Bailey, 1975; Anderson and Watsons, 1980; Anderson and May, 1991; Andersson and Britton, 2000; Diekmann and Heesterbeek, 2000) , which we denote Z. Of course, if the emergent pathogen becomes endemic in the population then there will never really be a "final" size, but in that case the matter at issue is the likely magnitude of the first major wave of cases.

To our knowledge, a formal, mathematical argument leading to an estimate of an epidemic's expected final size Z was first published in the landmark 1927 paper by Kermack and McKendrick (1927) . The formula for Z that Kermack and McKendrick obtained [Eq. (5) below] depends only on the basic reproduction number, R 0 (the expected number of secondary cases caused by a typical primary case in a fully susceptible population). However, without further analysis, it is unclear to what extent this formula depends on the particular assumptions that Kermack and McKendrick made in constructing their model, and hence whether it bears a strong relationship to the final size in realistic situations.

One implicit assumption of the Kermack and McKendrick (1927) analysis is that infectious periods are exponentially distributed. Anderson and Watsons (1980) showed that the final size formula remains valid if infectious periods follow a Gamma distribution. In Section 1.3 of their book, Diekmann and Heesterbeek (2000) generalized this statement to cover an arbitrary distribution of infectious periods.

Another key assumption of the Kermack and McKendrick (1927) analysis is that the host population is homogeneously mixed. It is well known that if this assumption is dropped then the final size will not necessarily be given by the standard formula. In particular, the existence of a core group, as for sexually transmitted diseases, gives rise to a different formula for the final size (Anderson and May, 1991; Diekmann and Heesterbeek, 2000) .

In this paper, we explore the generality of the standard Kermack-McKendrick final size formula (5). We begin with a pedagogical review of the main results that have been obtained previously, adding model structure in steps. We then proceed to generalize these results in three new directions, showing that the standard formula remains valid (i) regardless of the number of distinct infectious stages, (ii) if the mean contact rate is itself arbitrarily distributed and (iii) for a large class of spatially heterogeneous contact structures. We conclude that the Kermack-McKendrick formula (5) applies in great generality, making it an extremely useful relationship in practice.

The standard SIR model (Anderson and May, 1991; Diekmann and Heesterbeek, 2000; Brauer and Castillo-Chavez, 2001 ) is represented by a system of three ordinary differential equations,

Here, S, I, and R denote the proportions of the population that are susceptible, infectious and recovered, respectively. Recovered individuals are assumed to be immune to reinfection. Since S + I + R = 1, one of the three equations is redundant. The two parameters in the model are the transmission rate (β) and the recovery rate (γ ). The mean infectious period is T = 1/γ . The basic reproduction number is R 0 = βT. The proportion of the population that is infected can increase, and hence an epidemic can occur, if and only if R 0 > 1. In Eqs.

(1), recruitment of new susceptibles through birth or immigration is ignored, as is loss of individuals through mortality or emmigration. This approximation is reasonable provided the timescale of the epidemic is much shorter than the timescale over which demographic turnover occurs.

No recruitment ensures that the disease will eventually burn out, i.e., I(∞) = 0. To see formally that we must have I(∞) = 0, note that the positive orthant is invariant so all solutions of Eqs. (1) lie in the non-negative, bounded set defined by S, I, R ≥ 0 and S + I + R = 1. Observing that

we see that S + I is decreasing whenever I > 0. However, S + I is bounded below by 0; hence it has a limit. Moreover, Eq.

(2) implies that d dt (S + I) is bounded because I is bounded. Hence lim t→∞ d dt (S + I) = 0, so Eq. (2) yields I(∞) = 0. This model (1) is sufficiently simple that an exact solution can be obtained for the phase portrait. Forming dI/dS and integrating yields

An epidemic ends when no infectives remain. Consequently, we can find the final size Z by setting I = 0 in Eq. (3) and solving for Z = S(0) − S. This yields the implicit relation

(an explicit form is given in Appendix A). The general formula (4) applies regardless of the proportions of the population that are susceptible and infective initially. A less than fully susceptible population must be taken into account if some individuals have been vaccinated or retain immunity from previous exposures. However, in the important special case where a new pathogen enters a fully susceptible population, we have both I(0) 1 and S(0) ∼ 1. In the limit I(0) → 0, S(0) → 1,

which is the usual final size formula.

The practical utility of the final size formula would appear to be limited by the fact that it is derived from the overly simplistic standard SIR model. In the remainder of this paper we show that, in fact, the formula is valid under very general circumstances. We generalize the model in several steps, in order to explain the ideas clearly and to avoid exposing the reader to a notational burden at the outset.

The standard SIR model ignores any delay between the time an individual is infected and the onset of infectiousness. This latent period is often substantial compared with the infectious period (Anderson and May, 1991) . In addition, some pathogens (notably HIV) give rise to a sequence of distinct stages of infection, each yielding a different transmission rate (Redfield et al., 1986; Seligmann et al., 1987) .

It is common to incorporate latency by including an "exposed" class (E), yielding an SEIR model (Schwartz and Smith, 1983; Earn et al., 2000a) . However, noting that the exposed stage can be regarded as an infectious stage during which the transmission rate happens to be zero, we lose no generality by restricting attention to multiple infectious stage models, which we refer to as SI n R models if there are n infectious stages. Generalizing the SIR model (1) to include multiple infectious stages, the equations for an SI n R arė

where I i is the density of individuals in the ith infectious stage, β i is the transmission rate for contact with individuals in this stage, and 1 γi is the mean duration of this stage.

The sum of Eqs. (6) is zero, implying that S + n i=1 I i + R is invariant (and hence that Eq. (6d) is superfluous, as in the basic SIR model). To see that the final size is well-defined in this model, we generalize the argument given for the simple SIR model in Section 2. Let J k = S + k i=1 I i . Then d dt J k = −γ k I k , so I k > 0 implies J k is decreasing and lim t→∞ d dt J k = 0. Hence I k (∞) = 0 for all k.

Before a new disease is introduced into a population, everyone is susceptible, i.e., S = 1, I 1 = · · · = I n = 0 and R = 0, which is known as the disease free equilibrium (DFE). The stability of the DFE determines whether the new disease can cause an epidemic. Hyman et al. (1999) proved that the quantity

determines the stability of the DFE, i.e., the DFE is locally stable if R 0 < 1 and unstable if R 0 > 1. Note R 0 as defined by (7) is the sum of the products of the transmission rate and mean duration in each infectious stage, i.e., it is the basic reproduction number for this model. We now proceed to show that Eqs. (6) yield the same final size formula as (5). Anderson and Watsons (1980) established this for the special case in which each stage has the same transmission rate and duration (β i = β and γ i = γ for all i). Their method can easily be generalized to arbitrary β i and γ i as follows:

Let G k = n i=k+1 I i + R, k = 1, . . . , n − 1 and G n = R. Then Eqs. (6c) and (6d) imply d dt G k = γ k I k for each k = 1, . . . , n. Consequently, Eq. (6a) implies that the function

is a constant of the motion. For a newly invading disease, S(0) → 1, I k (0) → 0 and

Recalling that Z = R(∞) = 1 − S(∞), and Eq. (7), the final size formula Eq. (5) immediately follows.

Consideration of a special case of the model in the previous section allows us to infer that the distribution of stage durations in any infectious stage need not be exponentially distributed. The final size formula (5) holds true if the stage durations follow a Gamma distribution, which has probability density

Here, the shape parameter k is a positive integer and the scale parameter φ > 0. The mean is k/φ and the variance is k/φ 2 . This family of distributions includes the exponential distribution (k = 1), nearly normal distributions (k large) and the Delta distribution, which yields a fixed duration (k → ∞).

The term −γ I in Eq. (1b) for the simple SIR model implies that the durations of infectious periods are exponentially distributed with mean 1/γ (Brauer and Castillo-Chavez, 2001) . Since the sum of exponentially distributed random variables is Gamma distributed, a standard trick (Anderson and Watsons, 1980; Lloyd, 2001a,b) enables us to replace an exponential distribution with mean 1/γ by a Gamma distribution g k,kγ (which also has mean 1/γ ). The trick is to replace a single infectious stage with k identical exponentially distributed substages of mean duration 1/(kγ ). The model then has precisely the form the SI n R model (6), but with the sequence of infectious stages being mathematical artefacts rather than biologically meaningful.

Since this substage trick can be applied equally well to any infectious stage, Anderson and Watson's (1980) conclusion that the final size in an SIR model with Gamma distributed infectious periods is given by the usual formula (5) now generalizes to an arbitrary number of stages, each with Gamma distributed durations.

Although expressed differently, in their original paper Kermack and McKendrick (1927) studied a particular limit of the Gamma distributed multiple stage model described in the previous section.

Suppose we have a sequence of n infectious stages, each with the same, fixed stage duration τ . Thus, this is the limit k → ∞ in the Gamma distribution (10), with φ = 1/τ . During stage i, the transmission rate is β i (i = 1, . . . , n). If we now let the number of infectious stages increase (n → ∞) and the length of each stage decrease (τ → 0), keeping the total infectious period (T = nτ ) constant, then we arrive at a model in which the transmission rate (β) varies continuously through the infectious period. In such a model, an individual's infectivity depends on his/her stage-age, i.e., the amount of time since initial infection.

Using a stage-age SIR formulation, Kermack and McKendrick (1927) derived the final size formula (5) and showed that its form is independent of how the transmission rate depends on stage-age.

We have seen that the final size formula (5) is the same if stage durations are distributed according to any member of the family of the Gamma distributions (10). This suggests that any distribution of stage durations will yield the same final size.

Unlike the situation for Gamma distributed infectious periods, with an arbitrary distribution the time-evolution of the infectious class can no longer be expressed using ordinary differential equations (ODEs). Instead, we must construct a system of integro-differential equations (Feng and Thieme, 2000) . Let U(t) be the rate at which individuals become infectious (enter class I) at time t ≥ 0, i.e.,

The equation for the susceptible class is the same as in the standard SIR model,

We also define U(t) for t < 0 since this will allow us to specify the initial conditions. For t < 0, U(t) is the rate at which individuals who were still infectious at the initial time (time 0) became infectious at time t. Hence, the number of infectious individuals at time t = 0 is

Let f (t) be the probability density of the infectious period (so, for the special case of a Gamma distribution, f (t) is given by g k,φ (t) in Eq. 10). For an individual who becomes infectious at time τ ≥ 0, the probability density that s/he recovers (leaves the I class) at time t > τ is

However, for individuals who became infectious at time τ < 0, the probability density to leave class I at time t must be conditioned by still being infectious at time 0,

Hence, the rate of change of the number of infectious individuals at time t is

Those who leave class I enter class R, hence

In Appendix B, we show that if the infectious period is exponentially distributed then Eqs. (11) reduce to the the standard SIR model (1). As in the standard SIR model (1), the population size is invariant (S + I + R = 1), S ≥ 0, and S(∞) exists. Establishing that I(∞) = 0 (so the disease will burn out regardless of the infectious period distribution) requires a little more work. Since f (t) and l(τ, t) are probability densities, ∞ 0 f (t) dt = ∞ 0 l(τ, t) dt = 1; consequently, switching the order of integrations, we find

Next, integrating from 0 to ∞ on both sides of Eq. (11f), we have I(∞) = 0.

To relate the final size in the present model to the usual formula (5), we need an expression for the basic reproduction number R 0 . Feng and Thieme (2000) have

We now divide Eq. (11b) by S and integrate, yielding

Therefore, the usual final size formula (5) will follow if ∞ 0 I(t) dt = T Z. This is plausible, intuitively, since the integral ∞ 0 I(t) dt sums all infectious individuals, weighted by the time they are infectious; this must also equal the mean infectious period T times the number of individuals infected over the course of the entire epidemic (the final size Z). More rigorously, in Appendix B we prove Lemma 6.1. For the model specified by Eqs. (11), in the limit (11). If the initial number of infectious individuals is small (I(0) → 0), then the final size of the epidemic (Z) is given by the classical formula (5). In this formula, the basic reproduction number is R 0 = βT, where T is the mean infectious period.

With modest additional effort, we can generalize Theorem 6.1 to the multiple stage SI n R model, i.e., the usual final size formula (5) is still valid if there are multiple infectious stages and the durations of each stage are arbitrarily distributed. Recall that latent stages can be considered infectious stages with zero transmission rate, so we do not make any explicit reference to latency. As in all the situations considered above, we assume that vital dynamics (births and deaths) can be ignored and that recovery entails lifelong immunity.

Suppose the ith infectious stage has duration distribution with probability density f i (t) and mean T i . Thus, in the special case that durations in each stage follow Gamma distributions (10), we have f i (t) = g ki φi (t) and

The the rate of change of the susceptible class is given by Eq. (6a), as when all stage durations are exponentially distributed. Similar to Section 6, the equation for each infectious stage can be writteṅ

where U i (t) is the number of individuals entering class I i at time t, and

The initial conditions are specified by

Individuals who leave the last infectious stage recover (into class R), hencė

When the infectious periods and the latent period are exponentially distributed, i.e., f i (t) = γ i e −γi t , the technique introduced in Appendix B can be applied to show that the model (13) is indeed the standard multiple-stage SIR model. In their analysis of the SI n R model with arbitrarily distributed stage durations (13), Feng and Thieme (2000) have shown that the basic reproduction number is

and that the disease free equilibrium is unstable if and only if R 0 > 1.

Similar to previous sections, we find that S + i I i + R = 1 is invariant, S ≥ 0, and S(∞) exists. Integrating Eq. (13a) from 0 to t and and switching the order of integrations, we have I(t) ≥ 0; similarly, integrating from 0 to ∞, we have I(∞) = 0. Thus, the disease will also eventually burn out.

Theorem 7.1. Consider the multiple stage SI n R epidemic model with arbitrarily distributed stage durations, specified by Eqs. (13) . If the initial number of infectious individuals is suffciently small, i.e., in the limit i I i (0) → 0, the final size Z is given by the unique solution of Eq. (5), or explicitly by Eq. (A.2) , in which R 0 is given by Eq. (14) .

Proof. To prove this, we divide by S on both sides of (6a) and integrate with respect to time,

Using Lemma 6.1, Eq. (15) becomes

Suppose I i (0) = i 1. Switching the order of integration,

Also, integrating Eq. (13d) and switching the order of integrations, we have

i.e., summed over all time, the number of individuals who enter the (i + 1)th stage is the same as the number who enter the ith stage. We then integrate Eq. (6a). From the definition of U 1 (t) in Eq. (13c), we have

for all i = 1, . . . , n. Notice that the final size Z = S(0) − S(∞), so Eq. (16) becomes

During the epidemic of severe acute respiratory syndrome (SARS) in 2003 (Poutanen et al., 2003; Low and McGeer, 2003) , one aspect of observed transmission that received a great deal of attention was the existence of super-spreaders (WHO SARS Update 27, 2003), i.e., individuals who infect many more people than the average (e.g., at least an order of magnitude greater than R 0 ). In this section, we prove that the existence of the type of super-spreaders that occurred during the SARS outbreak in 2003 has no effect on the final size formula.

One mechanism for the generation of super-spreaders is that the infectious period distribution could be bimodal, such that a small proportion of individuals are infectious much longer than average. This situation is already covered by Theorem 7.1, which deals with arbitrary distributions of stage durations.

In this section, we consider a different type of generalization of the SIR model. In all of the models we have considered thus far, it has been assumed implicitly that all infectious individuals have the same transmission rate. If some individuals actually have a much higher than average transmission rate then they will be super-spreaders. This could arise, for example, because of individual variation in viral load (reflecting variation in immune response) or as the result of chance occurence of environmental circumstances that promote disease spread (as apparently occurred during the SARS epidemic (WHO SARS Update 33, 2003)). Whatever the origin, we focus here on variation in transmission rate that is manifested only after an individual is infected. In the following sections we consider the effects of intrinsic variation in contact rates among individuals, which results from inhomogeneous contact network structure in the host population.

We now proceed to develop and analyze an SIR model with a distribution of transmission rates and show that the usual final size formula (5) applies. Extension of our results to a multiple stage SI n R model is straightfoward and similar to the way we generalized our results for the SIR with arbitrary infectious period distribution to an SI n R model with arbitrary stage duration distributions.

Let I(β, t) be the distribution (at time t) of infectious individuals that have transmission rate β. Note that we assume implicitly that a given individual retains the same transmission rate throughout his/her infectious period, and hence that this transmission rate is determined as soon as infection occurs. Let q(β) be the probability that a newly infectious individual has transmission rate β. Thus, the mean transmission rate is

The rate of change of the proportion of the population that is susceptible iṡ

Upon infection, individuals leave the susceptible class S and enter the infectious class I, a proportion q(β) of them having a transmission rate β. Upon recovery, they leave the infectious class at rate γ . Thus,

In order to ensure that the model is well posed, we must establish that if I(β, 0) > 0 then I(β, t) > 0 for all time. Suppose this is not true. Then there existsβ and T such that ∞ 0 I(β, T) dβ > 0, I(β, t) > 0 for t < T, but I(β, T) = 0. Substituting these conditions into Eq. (21b), we obtainİ(β, T) > 0. But this contradicts the hypothesis that I(β, t) decreases from I(β, 0) > 0 to I(β, T) = 0. Hence I(β, t) > 0 for all t ≥ 0.

We now check that the disease will always burn out in this model, as expected.

which is strictly negative. Hence, S + I tot is decreasing and bounded below by zero. Hence lim t→∞ d dt (S + I tot ) = 0, i.e., I tot (∞) = 0. For the model specified by Eqs. (21), the basic reproduction number is

In Appendix B,we prove that the disease free equilibrium (DFE) is unstable if and only if R 0 > 1. To find the final size, we divide by S on both sides of (21a) and integrate,

Substituting Eq. (21a) in Eq. (21b) we can writė

Integrating this we obtain

Since we have established that I(β, ∞) = 0, we can write

Substituting (25) into (24) and using (23), we obtain

Hence, since the final size Z = S(0) − S(∞), we have It is possible to generalize this theorem for an SI n R model with arbitrarily distributed stage durations and arbitrarily distributed transmission rates in each stage. We do not go through the details because no new ideas are required. With Theorems 7.1 and 8.1 in hand, the extra effort required for the generalization is merely one of keeping track of notation.

In all of the models that we have discussed so far, we have assumed that the population is homogeneously mixed. In the remainder of this paper, we explore the significance of heterogeneous mixing for the final size formula. In the present section we identify an important class of spatially structured models for which the standard formula (5) remains valid.

A simple but important example of heterogenous mixing occurs if the population is divided into a number of spatially isolated patches (e.g., cities), often called a metapopulation (Hanski and Gilpin, 1997) . In ecological models, coupling among patches in a metapopulation usually occurs as a result of migration (Earn et al., 2000b) . In the present context, inter-patch coupling occurs because individuals travel temporarily from their home patch to other patches. On such journeys, susceptible individuals might contact infectious individuals and become infected, and infectious individuals might contact susceptibles and infect them. If we restrict attention to a single infectious stage with exponentially distributed infectious period then the model is specified as follows.

Here, S i , I i , and R i (i = 1, . . . , n) are the proportions of the population in patch i that are susceptible, infectious recovered, respectively. Summing Eqs. (26), we see that S i (t) + I i (t) + R i (t) = 1 is invariant in all patches. We denote by N i the number of individuals in patch i, so the total population size is N = i N i . It is also convenient to use B = (β i j ) to denote the n × n transmission matrix. Similar to our analysis of other models in previous sections, it is straightforward to prove that S i (∞), R i (∞), I i (∞) exist, that S i (∞) + R i (∞) = 1 and that I i (∞) = 0. Since β i j ≥ 0, the dominant eigenvalue λ of B is real and positive (Horn and Johnson, 1985) . If we linearize the model (26) at the disease free equilibrium (DFE) (S i = 1, I i = 0, R i = 0), we see that the stability of the DFE is determined by λ. If λ/γ > 1, then the DFE is unstable, otherwise it is stable. Thus, the threshold for exponential growth of cases (an epidemic) is determined by the value of R 0 = λ/γ , which can be interpreted as the basic reproduction number (van den Driessche and Watmough, 2002) .

The final size in patch i is Z i = S i (0) − S i (∞). We can obtain a system of n coupled algebraic equations for the Z i if we divide by S i on both sides of Eq. (26a) and integrate from 0 to ∞. These algebraic equations contain ∞ 0 I i dt, which can be eliminated by integrating Eq. (26c) and noting that Z i = R i (∞) − R i (0). In the limit I i (0) → 0, S i (0) → 1, we obtain

These equations form a special case of a more general final size formula discussed by Diekmann and Heesterbeek (2000) (exercises 6.17 and 6.19). In their solution to exercise 6.19, Diekmann and Heesterbeek (2000) sketch an argument that implies that if R 0 > 1 then there is a unique non-trivial (i.e., non-zero) solution to Eqs. (27). The final size for the whole metapopulation, Z, is not in general the simple average of the Z i 's in Eq. (27). Rather, each Z i must be weighted by the population size N i of patch i, i.e.,

Consequently, there will always be particular values of the N i that make Z in Eq. (28) equal to the solution of the standard formula (5). However, since population sizes are given a priori, the problem of interest is to find necessary and sufficient conditions on the transmission matrix (B) that ensure that the sum in Eq. (28) is equal to the final size obtained from the standard formula (5), regardless of the values of the N i . In Eq. (28), Z will have no dependence on the patch population sizes if and and only if

Given this, Eq. (27) implies

i.e., the transmission rate is the same in all patches. If, on the other hand, we start with Eq. (30) as an assumption, then inserting Eq. (29) in Eq. (27) Restricting attention to situations in which the transmission rate is the same in every patch may seem to be an extremely stringent constraint. In fact, this is often the situation of interest. For example, individuals who normally reside in Boston, Philadelphia or New York can expect to contact similar numbers of people each day, even though the New York metropolitan area has a population that is several times larger than the other two. More generally, if patches represent large cities in a given country, then it is reasonable to expect that the number of contacts an individual has will not depend strongly on the size of the city.

Given Eq. (30), we can write

where P i j is the probability, for an individual in patch i, that a given contact is with an individual in patch j, and β is the transmission rate (the product of the number of contacts per unit time and the probability that a contact between a susceptible and an infectious individual leads to transmission). Note that since the rows of the matrix P are probability distributions, it follows that n j=1 P i j = 1, i = 1, . . . , n,

i.e., P is a stochastic matrix (Horn and Johnson, 1985) . For convenience, we make the mild assumption that P i j > 0 for all i and j, i.e., for any pair of patches, there is a non-zero probability that residents of the two patches will come into contact. It then follows that if the final size Z is positive then the final size Z i is strictly positive in every patch. To see this, note that if Z i = 0 for some i then Eq. (27) 

Since P i j > 0 implies β i j > 0, Eq. (33) implies Z j = 0 for all j.

Since P is a positive, stochastic matrix, its largest eigenvalue is 1 (Horn and Johnson, 1985) . Consequently, the largest eigenvalue of B = β P is β and the basic reproduction number is R 0 = β γ . Therefore, we have Theorem 9.1. Consider the multi-patch SIR epidemic model specified by Eqs. (26) . Suppose the transmission rate in each patch is the same, i.e., n j=1 β i j = β, for each i. Then, in the limit that the initial proportion of individuals that are infectious is small, i.e., I i (0) → 0 for each i, the final size Z of an epidemic is given implicitly by Eq. (5), or explicitly by Eq. (A.2) , where R 0 = β γ .

At this point, it will come as no surprise that with some effort in keeping track of notation, Theorem 9.1 can be generalized to a multi-patch SI n R model with arbitrarily distributed stage durations and arbitrarily distributed transmission rates. 

Equations (26), which specify the spatially heterogeneous model we analyzed in the previous section, can be interpretted as representing other types of social heterogeneities, which are known to yield different final size formulae (Gart, 1968; Dwyer et al., 2000; Diekmann and Heesterbeek, 2000; Andreasen, 2003) . In this section, we briefly discuss two other interpretations, which lead to a transmision matrix that is not proportional to a stochastic matrix and will not yield the usual final size formula (5).

Rather than specifying heterogeneities of transmission resulting from spatial structure, as in the previous section, suppose that heterogeneities arise from age-structured mixing patterns. One possible formulation is to categorize the population by discrete age cohorts, which is typically motivated by mixing patterns of children in schools. If we consider a time-scale short enough that transfer from one cohort to the next can be ignored, then the standard "realistic agestructured model" (Schenzle, 1984) reduces to Eqs. (26) , where the transmission matrix B now refers to contact patterns among different age cohorts. As in the case of the spatial patch interpretation, it might be reasonable to assume that there is a stochastic matrix P such that

where the interpretation is that individuals in cohort i have a proportion P i j of their contacts with individuals in cohort j, and β i is the transmission rate for individuals in cohort i. Unlike the spatial interpretation, however, it would be hard to justify taking β i to be the same for each i. In particular, young children tend to come into much closer contact with their classmates than teenagers or adults do with members of their cohorts. As a result, Theorem 9.1 will not apply and the final size will not be given by the usual formula (5).

In Section 8 where we discussed super-spreaders, we specifically excluded the types of super-spreaders that occur for sexually transmitted diseases, namely individuals who always have a higher rate of contact with others, regardless of whether they happen to be infected. When such "core groups" (Yorke and Hethcote, 1984; Anderson et al., 1986) exist, Eqs. (26) can still be used, with the interpretation that the stratification of the population is by social group rather than spatial region or age cohort. Like the age-structured situation, we can safely assume that Eq. (34) holds but not that β i is the same for every i as in Eq. (30). Indeed, a large difference in β i is precisely what defines the core groups. Again, we have a situation where the usual final size formula will not apply. Core groups are well known to have a critical role in the transmission dynamics of sexually transmitted diseases, so it is not surprising that they will affect the final size of epidemics.

The well-known formula (5) for the expected final size of an epidemic is valid in remarkably general circumstances. Previous work (e.g., Kermack and McKendrick (1927) ; Anderson and Watsons (1980) ; Diekmann and Heesterbeek (2000) ) established that the formula is valid in an SIR model with an arbitrarily distributed infectious period (Theorem 6.1). Here we have shown, in addition, that the standard formula (5) is invariant to the number of latent and infectious stages of disease (Theorem 7.1), the distributions of transmission rates within stages (Theorem 8.1), and even to common spatial contact heterogeneities (Theorem 9.1).

The invariance of the final size formula has important practical implications. Typically, the time at which one wishes to estimate the expected final size of an epidemic is long before enough information has been gathered to estimate the distributions of latent and infectious periods or other epidemiological details. Our theorems provide rigorous support for estimates of the expected magnitude of epidemics based solely on estimates of the basic reproduction number R 0 , which is the only parameter that appears in the final size formula (5).

It should be noted that the final size formula refers only to the ensemble average size of an epidemic for a disease with a given R 0 . Different stochastic realizations of the same process will lead to different final sizes and our analysis says nothing about the variance or higher moments of the final size distribution. Kurtz (1980) proved that in the limit that the population size goes to infinity, the ensemble mean of the stochastic SIR model with arbitrarily distributed infectious period converges to the solutions of the integro-differential equation model specified in our Eqs. (11). Hence, the standard final size formula gives the mean final size of the stochastic SIR model. However, it is difficult to deduce the distribution of the final size, and this remains an open problem for stochastic SIR models with arbitrarily distributed infectious periods. For stochastic SEIR models with Gamma distributed latent and infectious periods, Anderson and Watsons (1980) derived a normal approximation to the final size distribution that is valid in the limit of large population size.

If there is little variation in the length of the latent period, and the infectious period is short, the stochastic SEIR model can be approximated by a chain binomial model (Bailey, 1975) , which is a discrete time Markov chain. The time step for this chain is equal to the fixed-length latent period, since the infectious period is presumed infinitesimal (all contacts occur instantaneously at the end of the latent period). For such models, the only parameter is the probability q that a contact between an infectious individual and a susceptible individual leads to infection. With the assumption that q = 1 − e −R0 , Von Bahr and Martin-Lof (1980) and Scalia-Tomba (1985) showed that final size distribution for the traditional chain binomial model is asymptotically normal in the limit of large population size. Scalia-Tomba (1986) and Andersson (1999) studied more general chain binomial models (with heterogeneous contact structures equivalent to those we discussed in Sections 9 and 10) and found asymptotic final size distributions. Their results imply that if q = 1 − e −R0 then the ensemble mean final size in these chain binomial models is the solution of our Eqs. (27).

We emphasized two circumstances under which the usual final size formula will not apply. One was the case of age-structured heterogeneity of transmission; while this would influence the final size, it is unlikely that we would be able to parameterize the age-structured contact patterns sufficiently well to improve on the final size estimate generated with the assumption of homogeneous mixing. The other case we discussed was the existence of social core groups, which are extremely important in the epidemiological dynamics of sexually transmitted diseases; in this case, some estimate of the difference in transmission rates in different social groups would likely be needed to obtain a useful estimate of the final size.

In computing the final size, we have always assumed that there are no temporal changes in the transmission rate β. Temporal variations in β could result from intrinsic seasonality of contact rates, imposition of control measures, or behavioral change in response to epidemic alerts. Small seasonal variations in β cause only small perturbations to the final size formula, and substantial seasonality is generally associated with school-age children (London and Yorke, 1973) who form a small fraction of the susceptible pool when a new infection enters a population. Infection control measures, such as those adopted during the 2003 SARS epidemic (Lipsitch et al., 2003; Wallinga and Teunis, 2004) , could have a dramatic effect on the final size; in this case, rather than R 0 , the reproduction number that is relevant for the final size formula must be calculated after sustainable precautionary measures have been put into place. Individual behavioral change could alter the transmission rate through time-dependence (β = β(t), as in the case of imposed control measures) or through density-dependence (β = β(I), with β decreasing as I increases, which would arise if people tend to be more careful when they are aware of more cases). In either case, the standard final size formula might yield a poor approximation of the true final size.

All of our analysis is predicated on the assumption that vital dynamics (births and deaths), and more generally any source of new susceptible individuals, can be ignored. This approximation is not valid in circumstances where infectious periods are substantial compared with life expectancy or where immunity decays rapidly. It will, however, typically be relevant whenever a new respiratory pathogen emerges, as in the case of the SARS outbreak in 2003 (Donnelly et al., 2003) or the emergence of a new subtype of influenza (Earn et al., 2002) . We are interested in growth or decay of I(β, ), which is determined only by Eq. (D.1b). Defining X(t) = ∞ 0 β I(β, t) dβ, multiplying both sides of Eq. (D.1b) by β, and integrating with respect to β, we find

On the other hand, recall that I(t) = ∞ 0 I(β, t) dt (the total number of infectious individuals at time t). Therefore, integrating Eq. (B.2b) with respect to β, noting that q(β) is a probability density function so The origin of this system of two ordinary differential equations is locally stable if and only if ∞ 0 βq(β) dβ − γ < 0 or, equivalently, R 0 < 1. Thus, the DFE is stable if R 0 < 1 and unstable if R 0 > 1. We can expand model (21) to include multiple stages. Using the same techniques, we find that the conclusion that the epidemic threshold is given by R 0 = 1 remains valid.

On the spread of a disease with gamma distributed latent and infectious periods

Infectious Diseases of Humans: Dynamics and Control

A preliminary study of the transmission dynamics of the human immunodeficiency virus (HIV), the causative agent of AIDS

Stochastic epidemic models and their statistical analysis

The asymptotic final size distribution of multitype chain-binomial epidemic processes

Dynamics of annual influenza a epidemics with immuno-selection

The Mathematical Theory of Infectious Diseases and its Application

Mathematical models in population biology and epidemiology

Mathematical epidemiology of infectious diseases: Model building, analysis and interpretation

Epidemiological determinants of spread of causal agent of severe acute respiratory syndrome in Hong Kong

Pathogen-driven outbreaks in forest defoliators revisited: Building models from exprimental data

Ecology and evolution of the flu

Coherence and conservation

A simple model for complex dynamical transitions in epidemics

Endemic models with arbitrarily distributed periods of infection I: Fundamental properties of the model

The mathematical analysis of an epidemic with two kinds of susceptibles

Metapopulation Biology: Ecology, Genetics, and Evolution

Matrix Analysis

The differential infectivity and staged progression models for the transmission of HIV

A contribution to the mathematical theory of epidemics

Relationships between stochastic and deterministic population models

Transmission dynamics and control of severe acute respiratory syndrome

Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods

Realistic distributions of infectious periods in epidemic models: Changing patterns of persistence and dynamics. Theor

Recurrent outbreaks of measles, chickenpox and mumps. I Seasonal variation in contact rates

SARS-one year later

Identification of severe acute respiratory syndrome in Canada

The Walter Reed staging classification for HTLV-III/LAV infection

Asymptotic final size distribution for some chain-binomial processes

Asymptotic final size distribution of the multitype Reed-Frost process

An age-structured model of pre-and post-vaccination measles transmission

Infinite subharmonic bifurcation in an SEIR model

Reproduction numbers and subthreadold endemic equilibria for compartmental models fo disease transmission

Threshold limit theorems for some epidemic processes

Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures

From MathWorld-A Wolfram Web Resource. http:\math world.wolfram.com/LambertW-Function.html

Severe acute respiratory syndrome (SARS), multi-country outbreak

Severe acute respiratory syndrome (SARS), multi-country outbreak

Gonorrhea: Transmission dynamics and control

An explicit solution for Z expressed in terms of elementary functions is not possible, but Eq. (4) can be solved explicitly for Z in terms of the Lambert W function (Weisstein) , 

When the infectious period is exponentially distributed, i.e., f (t) = γ e −γ t , we have l(τ, t) = f (t). Thus Eq. (11f) can be writtenLet