Using a stochastic susceptible–infected–removed meta-population model of disease transmission, we present analytical calculations and numerical simulations dissecting the interplay between stochasticity and the division of a population into mutually independent sub-populations. We show that subdivision activates two stochastic effects—extinction and desynchronization—diminishing the overall impact of the outbreak even when the total population has already left the stochastic regime and the basic reproduction number is not altered by the subdivision. Both effects are quantitatively captured by our theoretical estimates, allowing us to determine their individual contributions to the observed reduction of the peak of the epidemic.

Simple models for the spread of infectious diseases are useful for the quantitative characterization of an epidemic as well as for forecasting future infection numbers and guiding decision-making for containment. Different extensions and refined versions of these models have been created to extract various factors that may be critical for the dynamics and prevention of epidemics. Although it is well known that stochastic fluctuations can alter the dynamics as well, they are often neglected at higher infection number levels such that the contact rates and basic reproduction number become the central quantities of interest. In contrast, we investigate a situation in which stochastic effects can quantitatively change the course of an epidemic when infection numbers are large and contact rates remain unaltered. We consider an extended Susceptible–Infected–Removed (SIR) model in which a large population is subdivided into a certain number of sub-populations, each containing only a few infected individuals. For the limiting case of perfect isolation, i.e., when the epidemic evolves independently in each sub-population with no cross-infections, we derive analytical estimates for these stochastic effects that together recapitulate the results of extensive numerical simulations. Our central quantity of interest is the peak total number of simultaneously infected individuals, which we compare between the subdivided population and a single large population with an identical reproduction number. Our analysis suggests that regional isolation can resurrect certain stochastic effects and thereby contribute to effective containment, regardless of the initial distribution of infected individuals.

Generic models such as the Susceptible–Infected–Removed (SIR) model conceived by Kermack and McKendrik1 are indispensable for characterizing the bulk properties of epidemics and determining the influence of crucial parameters on the dynamics. The contact rate between individuals, which is proportional to the reproduction number R0, usually plays a crucial role, as its reduction through containment measures directly slows the spreading of the disease. On a large scale (states or countries), numbers of infections during the height of an epidemic are usually large such that deterministic mean-field descriptions are appropriate. These have been widely used to track the course of epidemics and the effect of interventions, for example, for the current spreading of COVID-19.2 

While many details about the biology and modes of infection of a specific disease are important for its dynamics in detailed models,3 even basic SIR models have been extended in various conceptual directions. Besides various general topologies of the underlying contact and mobility networks,4–6 so-called meta-population models have been used to separate the disease dynamics within local environments from its spread between them.7 It has been shown that it is possible to calculate effective quantities for the whole population, such as reproduction numbers (i.e., a threshold theorem),8 the final attack ratio,9 and criteria for persistence10 in deterministic models of such subdivided populations.

Another important deviation from simple bulk behavior arises through stochasticity (see Ref. 11 and references therein). Stochastic versions of extended SIR and related models have been used to calculate corrections to the outbreak threshold,12 consequences of stochasticity for contact tracing,13 and other control schemes,14 to only name a few. Stochastic effects are also observed in agent-based15 and meta-population models.16,17 Here, we seek to study the joint effect of subdivision and stochasticity on the overall magnitude of an epidemic for a fixed initial number of infected individuals in the total population. In general, subdivision can be expected to artificially boost fluctuations, as the infection numbers in each sub-population can be small even when the total number of infections in the entire population is large. We would like to quantify the ability of such increased stochasticity to reduce the impact of the epidemic. We deliberately refrain from applying any form of traditional containment in our model, such as further reductions in the contact rate or contact tracing.18 In particular, we design the subdivision such that the deterministic dynamics of the epidemic in the subdivided population remains unchanged compared to a single large population, as outlined in Sec. II. This allows us to compare the peak number of infected individuals in the entire population for each scenario both analytically and numerically in order to extract the specific effects of stochasticity triggered by subdivision.

We consider a population of N individuals with SIR dynamics,1 

S+Ib/NI+I,IkR,
(1)

with S, I, and R referring to susceptible, infected, and removed individuals, respectively, where removal with per capita rate k happens due to recovery, quarantine, or death. The rate b corresponds to the number of contacts per unit time an individual has with a random other individual in the population, multiplied by the probability that a contact between a susceptible and an infected individual leads to transmission. The total transition rate from S to I per unit time is, therefore, bNSI. The two rates b and k are related to the basic reproduction number R0=b/k, which is independent of population size. The deterministic epidemic threshold above which an outbreak occurs is R0=1, and we assume R0>1 throughout this study. The population is subject to the total constraint N=S(t)+I(t)+R(t), where we denote the number of individuals in each state by the same letters. For simplicity, all initial conditions assume R(0)=0 such that they are uniquely defined by N and the number of initially infected I0=I(0).

When a population of total size N is split up into Ns sub-populations, we simulate Ns separate copies of the system (1), with N, S, I and R replaced by Ni, Si, Ii, and Ri, respectively, where the index i refers to the different sub-populations and Ni=N/Ns. N=Ni, S=Si, I=Ii, and R=Ri refer to the population totals. The initial number of infected individuals is distributed either uniformly or randomly across the Ns sub-populations. All numerical results in this study are obtained from stochastic simulations of Eq. (1) using the Gillespie algorithm.19 To account for the inherent stochasticity of the system, several realizations, i.e., identical simulations with different random number generator seeds, are simulated for each parameter set. We report the number of realizations as well as distributions, averages, and standard deviations across the results as appropriate. Our main figure of interest is the peak number of infected individuals Imax or, equivalently, the peak infected fraction of the population γ=Imax/N. These could be considered a measure for the impact of the epidemic and the strain on the health care system and public health resources such as the agencies that perform contact tracing and testing.

The reaction scheme (1) results in the deterministic mean-field equations,

dSdt=bNSI,
(2a)
dIdt=bNSIkI,
(2b)
dRdt=kI,
(2c)

which give rise to two regimes in the dynamics. During the initial regime, I starts off from an initial value I(0)=I0, rises exponentially I0e(bk)t, and saturates to a peak value,

ImaxγN(1kb[1+log(b/k)])N,
(3)

where the approximation for the maximum fraction of infected individuals 0<γ<1 is valid as long as the entire population is initially susceptible; i.e., S(0)N.20 In the secondary regime when the recovery dynamics dominates, I decays to zero exponentially, as the number of susceptibles decreases below the value necessary to sustain spreading.

In this deterministic system, a subdivision into Ns smaller sub-populations of size N/Ns will have no effect since Eq. (2) remains invariant when S, I, R, and N are scaled by the common factor 1/Ns. Relative to their individual sub-population sizes Ni, the same dynamics are observed in all sub-populations and the dynamics of the population totals S=Si, I=Ii, and R=Ri are identical to those of a single large population. Therefore, the subdivision is not analogous to cutting links in a contact network but rather a redistribution of them since we assume that the contact rate b remains unchanged. This conservative assumption means that individuals in each sub-population still have the same number of contacts per unit time as they had in the large population despite the smaller number of individuals to choose from. While, in reality, the contact rate b might decrease in such a situation and deterministically reduce R0 and, therefore, Imax, we intentionally keep it constant here to extract the effects of stochasticity.

Deterministic behavior only applies if S and I are both large, particularly only after the number of infected people I has risen to appreciable levels. If I is still low, stochastic fluctuations determine whether I will “take off” and develop exponential behavior even if b>k. This effect was already considered shortly after Kermack and McKendrick introduced the original SIR model21 and is now well-known. However, in a subdivided population, it can significantly alter the course of the outbreak in the total population if the initial number of infected individuals in a single sub-population is low enough (even if the number is large in the total population). An example for populations of N=1000000 individuals split into Ns=10 sub-populations is shown in Figs. 1(a) and 1(b), along with the expected dynamics of a single large population (red curve). In one example set of ten sub-populations, only three sub-populations (blue, yellow, and green curves) experience a significant outbreak, and they are desynchronized with a broad distribution of individual peak times [Fig. 1(c)]. Spontaneous extinction and desynchronization lead to an average behavior across 100 simulations with a significantly reduced peak (turquoise curve). Note that, on average, both the undivided large population and the sum of the smaller sub-populations initially exhibit comparable exponential growth in the number of infected individuals [Fig. 1(b)]. This means that, while extinction in some sub-populations and fluctuations in timing happen early on, their effect is only seen later during the saturation phase.

FIG. 1.

Stochastic effects lower the peak in subdivided populations. (a) Time course for a population of N=1000000 with I0=10 initially infected individuals for Ns=1 large population (red) and a population split into Ns=10 sub-populations (turquoise), b=0.2, k=0.14. Shading indicates ±25% confidence intervals across 100 simulations. Sub-populations are shown from one simulation with Ns=10. (b) Enlarged plot of the initial phase for the same traces as in panel (a). (c) Distribution of peak times in sub-populations for Ns=10. The occurrence fraction indicates the fraction of sub-populations across all simulations. The dashed line indicates the analytical approximation, Eq. (8), with uniform n=I0/Ns=1.

FIG. 1.

Stochastic effects lower the peak in subdivided populations. (a) Time course for a population of N=1000000 with I0=10 initially infected individuals for Ns=1 large population (red) and a population split into Ns=10 sub-populations (turquoise), b=0.2, k=0.14. Shading indicates ±25% confidence intervals across 100 simulations. Sub-populations are shown from one simulation with Ns=10. (b) Enlarged plot of the initial phase for the same traces as in panel (a). (c) Distribution of peak times in sub-populations for Ns=10. The occurrence fraction indicates the fraction of sub-populations across all simulations. The dashed line indicates the analytical approximation, Eq. (8), with uniform n=I0/Ns=1.

Close modal

During the initial phase, we can assume that SN and that I follows a simple birth–death process with rates b for “birth” and k for “death.” We shall use this analogy for derivations throughout this study and in  Appendixes A– D. We briefly recapitulate one important result from the theory of branching processes here, namely, that an exponentially growing population that starts from an initial condition of I(0)=1 has a finite extinction probability of

P0(t)=kbe(bk)t1e(bk)tk/b,
(4)

which asymptotically approaches k/b at long times; see the derivation in  Appendix A. This means that with probability p1ext=k/b, the dynamics never enters the exponentially growing deterministic regime but decays back to zero due to number fluctuations.22 Therefore, for two independent lineages in the same population, the extinction probability is p2ext=(k/b)2, and, similarly, pnext=(k/b)n, as long as the total population is sufficiently large such that the lineages do not interfere with each other. We will use these extinction probabilities and other statistics of the birth–death process to derive analytical approximations for the effects of extinction and desynchronization on the stochastic dynamics.

1. Extinction

To obtain an estimate for the effect of extinction and the distribution of infected individuals, we add up the maximum numbers of infected individuals in the sub-populations. Each of these peaks is approximately γN/Ns but only if the infection does not stochastically become extinct during the initial stages. For large-population sizes and values of b/k that result in a significant peak, extinction usually happens well before the peak is reached in other sub-populations (see  Appendix B) such that these populations do not contribute. Therefore, on average, the contribution of each sub-population will be Is,max(n)=γ(1pnext)N/Ns, where n indicates the number of initially infected individuals in the sub-population and pnext is the probability that they go extinct without entering deterministic growth as discussed above. Therefore, the total peak number of infected individuals in all the sub-populations due to extinction is given by Imaxext=ngnIs,max(n), where gn is the number of sub-populations with n initially infected individuals. Note that Ns=ngn. Combining the above equations, we obtain

Imaxext=γNNsn(1pnext)gn=γN[1ngn(k/b)nngn].
(5)

The above result manifestly shows that

γextImaxextN=γ[1ngnNs(k/b)n]<γ
(6)

holds. Note that this reduction is exclusively due to extinction, and the simple summation of the individual maxima neglects the possible desynchronization between sub-populations, which we will consider further below.

For example, for the ideal case where each sub-population only contains at most one infected individual, we have

γ1ext=γI0Ns(1kb),
(7)

where g1=I0 is the total number of initially infected individuals in the large population (for this to make sense, NsI0 is required). Since γ corresponds to the case where the population was not split up, the peak number of infected can, therefore, be reduced by increasing the number of sub-populations Ns or by bringing b closer to k. Note that this is in addition to a potential decrease in the deterministic peak fraction γ of infected [cf. Eq. (3)] that would result if the subdivision also led to fewer contacts (i.e., a reduced rate b), which we have conservatively assumed not to be the case here.

2. Desynchronization

The independent summation of maxima in different sub-populations is a conservative estimate since fluctuations can lead to stochastic desynchronization and thus to a further reduction of the peak value. The distribution of peak times in the sub-populations from the previous example is shown in Fig. 1(c). The temporal shift between the different sub-populations can be attributed entirely to stochastic fluctuations in the initial phase of the dynamics. Assuming that this time shift accumulates while the dynamics can still be modeled as a pure birth–death process without saturation effects, we can derive the probability distribution for the deviation from the mean peak time Δtpeaktpeaktpeak as

P(Δtpeak)=k(1k/b)[1(k/b)n]×exp((bk)(τ¯+Δtpeak)kbe(bk)(τ¯+Δtpeak)),
(8)

where n is the initial number of infected individuals in the sub-population and τ¯=ln(γk/b)/(bk) with γ being the exponential of the Euler constant (see  Appendix C for details). Note that n here was only used to incorporate the extinction probability, while the shape of the distribution is based on a single initially infected individual. Nevertheless, this result is in excellent agreement with the measured distribution for randomly distributed infected individuals [see the dashed line in Fig. 1(c)].

We can then use this distribution to obtain a quantitative estimate for the additional peak reduction due to desynchronization. For this purpose, we approximate the deterministic time evolution of I in the vicinity of the peak as I(t)Nγexp(12bkγ(ttpeak)2), which is valid as long as S(t) remains of order N (see  Appendix D); i.e., b/k is not too large. In the limit of many superimposed peaks of this shape, with the variability of tpeak given by Eq. (8), the peak is reduced by an additional factor α1,

γcon=γextα,α=1+π2[R01log(R0)]6(R01)2.
(9)

The peak number of infected individuals, with both stochastic effects of the confinement taken into account, similarly becomes Imaxcon=Nγcon=Imaxext/α. It is interesting to note that this reduction factor is bounded from below by limR01α1=12/(12+π2)0.7407. The desynchronization effect is, therefore, much more limited than the extinction effect.

We consider as an example a region with a population of 8 000 000 and 500 infected individuals (I0/N6105) and assume a removal rate of k=0.14, corresponding to a realistic mean removal time of 1/k7 days for the recent epidemic23 (particularly if symptomatic individuals are quickly removed from the infectious pool through quarantining). Let us further assume that the infectious contact rate is b=0.2 (>k). This corresponds to a substantial reduction of R0 from its initial value of 2–2.524 through mild measures such as social distancing, although the epidemic would still spread exponentially, with infection numbers doubling about every 12 days.

If this population is allowed to mix homogeneously, the dynamics will evolve according to the deterministic prediction with a peak around 5% infected individuals (blue data in Fig. 2). If instead, the population is split up and the 500 infected people are distributed randomly across the sub-populations, the peak percentage of infected individuals decreases to around 3% (for 100 sub-populations of 80 000 people) or 1% (for 500 sub-populations of 16 000 people) on average (red and yellow, respectively). In all cases, the analytical estimate that only considers the extinction effect, Eq. (6), is only an upper bound for the peak percentage of infected individuals in the total population, while also considering desynchronization according to Eq. (9) yields a good estimate the typical peak values. The peak time distributions for the three different ways of splitting up the population shown in Fig. 2(c) also agree with the analytical estimate of Eq. (8). Note that these distributions are not normalized since a significant fraction of sub-populations experience extinction of the epidemic and, therefore, do not exhibit a peak. There is also a subtle, non-monotonic effect on the termination time of the epidemic [Fig. 2(d)] whose distribution is broader when the population is split up but does not change position appreciably. Note that the reduction for Ns=500 sub-populations in Fig. 2 is comparable (or even slightly lower) than the case where the 500 infected individuals are not distributed randomly across the sub-populations, but each sub-population contains exactly one infected individual. In this case (see Fig. 3), there are no sub-populations with initially zero infected individuals, implying that the reduction in the peak value compared to the large homogeneous population is strictly due to extinction and desynchronization, which are again well predicted by the analytical estimates.

FIG. 2.

Epidemics for different subdivisions of the population. N=8000000, b=0.2, k=0.14, three different values of Ns=20, 200, and 40 individual simulations for Ns=1, Ns=100, and Ns=500, respectively. (a) Time courses (solid lines) and 2.5/97.5 percentiles (shading). (b) Distributions of the observed peak percentage γcon (in the whole population). The occurrence fraction indicates the fraction of simulations. Analytical estimates for γext (dashed lines), Eq. (6), and γcon (solid lines), Eq. (9), assume gn according to a binomial distribution. (c) Distribution of peak times in the sub-populations. The inset provides an enlarged y-axis. Dashed lines indicate analytical approximation, Eq. (8), assuming a uniform n=I0/Ns for each case. (d) Distribution of termination times, defined as the time when I in the total population drops below I0.

FIG. 2.

Epidemics for different subdivisions of the population. N=8000000, b=0.2, k=0.14, three different values of Ns=20, 200, and 40 individual simulations for Ns=1, Ns=100, and Ns=500, respectively. (a) Time courses (solid lines) and 2.5/97.5 percentiles (shading). (b) Distributions of the observed peak percentage γcon (in the whole population). The occurrence fraction indicates the fraction of simulations. Analytical estimates for γext (dashed lines), Eq. (6), and γcon (solid lines), Eq. (9), assume gn according to a binomial distribution. (c) Distribution of peak times in the sub-populations. The inset provides an enlarged y-axis. Dashed lines indicate analytical approximation, Eq. (8), assuming a uniform n=I0/Ns for each case. (d) Distribution of termination times, defined as the time when I in the total population drops below I0.

Close modal
FIG. 3.

Plots analogous to Figs. 2(a) and 2(b) for the case Ns=500, but with exactly one initially infected individual in each sub-population instead of a random distribution. Analytical estimates, Eq. (6) (dashed line) and Eq. (9) (solid line), accordingly use g1=I0=500.

FIG. 3.

Plots analogous to Figs. 2(a) and 2(b) for the case Ns=500, but with exactly one initially infected individual in each sub-population instead of a random distribution. Analytical estimates, Eq. (6) (dashed line) and Eq. (9) (solid line), accordingly use g1=I0=500.

Close modal

To examine the validity of our approximations across different parameters, we varied the contact rate b and carried out numerical simulations for values of R0 ranging between 1.14 and 2. We analyzed the resulting peak magnitudes to extract the individual contributions of extinction and desynchronization, which are in excellent agreement with our predictions of Eqs. (6) and (9), as shown in Fig. 4. The contribution of extinction alone was estimated numerically by summing maxima in different sub-populations, regardless of their timing. Overall, the simulations confirm the relative importance of the extinction effect, whereas the additional reduction by desynchronization plays a smaller role. Figures 4(a) and 4(b) show the case where Ns=I0=100, i.e., number of sub-populations and initially infected individuals is the same, and exactly one infected individual is placed in each sub-population. This serves to demonstrate the maximum effect of extinction, whereas in Fig. 4(c), a large share of the peak reduction is due to sub-populations containing no infections, as I0=100<Ns=500. However, the random distribution of infected individuals for Ns=I0=500 [Fig. 4(d)] leads to a very similar result as in Fig. 4(b), although some of the reduction is due to the initial distribution (i.e., sub-population without any infections). For a high number of sub-populations Ns as in Figs. 4(c) and 4(d) (and consequently a smaller sub-population size), deviations from the theory begin to appear toward low values of b very close to k, as the timescale of the extinction process becomes comparable to that of the deterministic SIR dynamics. In this regime, the distinction between an initial stochastic phase approximated by a birth–death process and the onset of saturation effects becomes increasingly blurred, as we show analytically in  Appendix B. In particular, this affects the estimation of the extinction contribution (marked by black dots).

FIG. 4.

Peak reduction for different subdivisions and values of b. N=8000000, k=0.14. All data points represent averages across 100 independent stochastic simulations, and error bars indicate standard deviation. (a) Peak fraction of infected individuals for I0=100, Ns=100 with each sub-population initially containing exactly one infected individual. The symbol color indicates the reduction due to extinction (yellow) or both extinction and desynchronization (blue) as measured in simulations. Black symbols represent the large-population control. Yellow/blue shading and solid lines indicate analytical predictions from Eqs. (6) and (9), respectively. The black line indicates the deterministic estimate from Eq. (3). (b) Same as panel (a), plotted logarithmically and normalized by a theoretical peak fraction γ without subdivision, Eq. (3). (c) I0=100, Ns=500, each sub-population containing at most one infected individual. The yellow color now represents the reduction both to extinction and the initial distribution of infected individuals since 400 sub-populations are initialized with I=0. (d) I0=500, Ns=500, infected individuals randomly distributed across sub-populations. Black dots in (c) and (d) mark data points where the estimation of the extinction effect is affected by overlapping timescales between different processes (see the text and  Appendix B).

FIG. 4.

Peak reduction for different subdivisions and values of b. N=8000000, k=0.14. All data points represent averages across 100 independent stochastic simulations, and error bars indicate standard deviation. (a) Peak fraction of infected individuals for I0=100, Ns=100 with each sub-population initially containing exactly one infected individual. The symbol color indicates the reduction due to extinction (yellow) or both extinction and desynchronization (blue) as measured in simulations. Black symbols represent the large-population control. Yellow/blue shading and solid lines indicate analytical predictions from Eqs. (6) and (9), respectively. The black line indicates the deterministic estimate from Eq. (3). (b) Same as panel (a), plotted logarithmically and normalized by a theoretical peak fraction γ without subdivision, Eq. (3). (c) I0=100, Ns=500, each sub-population containing at most one infected individual. The yellow color now represents the reduction both to extinction and the initial distribution of infected individuals since 400 sub-populations are initialized with I=0. (d) I0=500, Ns=500, infected individuals randomly distributed across sub-populations. Black dots in (c) and (d) mark data points where the estimation of the extinction effect is affected by overlapping timescales between different processes (see the text and  Appendix B).

Close modal

Reducing the infectious contact rate b or increasing the removal rate k directly leads to a decrease of the deterministic peak fraction of infected, γ. The above analysis shows that, even without changing R0=b/k, the isolation of small sub-populations can reduce the overall peak number of infected people in the ideal case of at most one infected individual per sub-population by an additional factor of up to I0/Ns(1k/b)/α when I0/Ns<1. One contribution comes from the communities that have no infections and are now protected (I0/Ns), while another contribution comes from the possibility that an infection chain in a local community stochastically ends due to fluctuations (k/b). Stochastic desynchronization (1/α) further reduces the peak by up to about 25% according to Eq. (9). However, as shown by our estimates and confirmed by the numerical simulations, even outside this ideal scenario, a reduction can be achieved, regardless of the distribution of infected individuals across the sub-populations, and the reduction will be larger if b/k is already close to 1. It is also worth noting that, in contrast to reductions in R0=b/k, the timescale of the outbreak is not increased.

The benefits of subdivision are obvious even from a deterministic standpoint in the case where many regions initially contain no infected individuals—in this case, subdivision prevents spreading of the epidemic to disease-free communities. However, our analysis shows that this advantage persists due to stochastic extinction events and desynchronization even if the sub-populations are so large that many or all of them initially contain infections, as long as I0/Ns1. Of course, increasing Ns further is always beneficial due to the above-mentioned deterministic effect, with the trivial limiting case of one group per household (an extremely strict lockdown). In contrast, aiming at I0/Ns1 could still allow for the functioning of local socioeconomic life in fairly large sub-populations if I0 is not too large when the subdivision happens.

While extinction has been widely considered for SIR-type models11,21 and has been related to a minimum number of infections necessary to cause a “major” outbreak,14 we have shown here that, even if the dynamics in the large population is outside the stochastic regime, it is possible to resurrect these effects by artificially sub-dividing the population. Because of the strong exponential dependence of the extinction probability on n [see Eq. (5)], it is important to note that I0 denotes the true number of infections, including undetected and/or asymptomatic cases. Another aspect we have neglected here is that of cross-infections: In reality, sub-populations cannot be perfectly isolated; therefore, local extinction might only be temporary, as has been seen in studies of persistence.10,16 The calculated peak reduction would be observed in the limit of small cross-infection rates. In contrast to extinction, desynchronization does not reveal itself on the level of a single population (except as a difference in timing) and is, therefore, an emergent property of the subdivision scenario, which is likely to persist in the presence of cross-infections. In the framework presented in Sec. II A, these could be included (without changing R0) by allowing a certain fraction ξ of contacts across the entire population and only restricting the remaining fraction 1ξ to within each sub-population. We set up such a model in a separate study25 to investigate a potential realistic containment strategy.

In reality, individuals will not compensate for all avoided contacts outside the local sub-population with contacts within it, as we have conservatively assumed by keeping b constant upon subdivision. Instead, isolation will naturally lead to a reduction in b, akin to cutting links in the spreading network5 so that the effect of subdivision will be a combination of deterministic reductions in R0 and the stochastic effects presented here. Subdivision of a population can be complementary to containment measures, such as social distancing and electronic contact tracing,13,23 which still allow for the functioning of local public life. However, it also does not preclude the activation of more drastic measures in regions beginning to show deterministic exponential behavior.25 

We have benefited from discussions with J. Agudo-Canalejo, A. Bahrami, H. Bickeboeller, E. Bodenschatz, W. Brück, H. Chaté, R. Fleischmenan, T. Friede, T. Geisel, H. Grubmüller, R. Jahn, B. Mahault, V. Priesemann, T. Richter, S. Scheithauer, A. Vilfan, M. Wilczek, and R. Yahyapour. This research was supported by the Max-Planck-Gesellschaft.

The data that support the findings of this study are available within the article.

Consider a population of the infected individuals I that can undergo the following two processes:

IbI+I,Ik;
(A1)

i.e., each I can give birth to another I with rate b or it can die with rate k, at any time. Ignoring the stochasticity, the average behavior of the system is described by exponential birth and death. The population n¯(t) can be determined as follows:

dn¯(t)dt=(bk)n¯(t)n¯(t)=e(bk)t,
(A2)

where we have assumed that the initial size of the population is one. As this is a one-step process, the probability of finding n copies of I in the sample at time t satisfies the following master equation:

dPn(t)dt=k(n+1)Pn+1(t)+b(n1)Pn1(t)(k+b)nPn(t).
(A3)

The factor of n is needed because the birth or death could happen to anyone. Equation (A3) can be solved by an ansatz of the form Pnfn for n1, which together with the initial condition Pn(0)=δn,1 gives us the solution as

Pn(t)=n¯(1k/b)2(n¯1)(n¯k/b)(n¯1n¯k/b)n.
(A4)

The distribution can be used to calculate the first two moments,

n(t)=nnPn(t)=n¯(t)=e(bk)t,
(A5)
Δn2=[nn]2=(b+kbk)e(bk)t[e(bk)t1],
(A6)

which reveal more interesting features about the system. First, it is reassuring that the average population size behaves according to the mean-field description above that predicted exponential growth or decay. A quantity of interest is

Δn2n¯=(b+kbk)[e(bk)t1],
(A7)

which probes whether number fluctuations follow a characteristic Poisson behavior. In the long time limit, we have

Δn2n¯={,b>k,k+bkb,b<k,
(A8)

which shows that while a decaying population that corresponds to b<k has a Poisson behavior, a growing population corresponding to b>k has giant number fluctuations, which can be characterized via

Δnn¯=b+kbk1e(bk)t,
(A9)

which leads to

Δnn¯=b+kbk
(A10)

in the long time limit. In other words, the fluctuations scale with the average population size when b>k and with the square root of the average population size when b<k.

The above solution allows us to calculate the extinction probability of the population P0(t), which is an absorbing state. We find

P0(t)=1n=1Pn(t)=kbe(bk)t1e(bk)tk/b,
(A11)

which is a very interesting result. When k>b, n¯0 at long times, and we obtain P0=1. It is no surprise that extinction at long times is a certainty when the death rate is larger than the birth rate. However, when k<b, n¯ at long times, and we obtain P0=k/b, a result that is in contradiction with the prediction of the average behavior of the system, which is exponential growth. Therefore, number fluctuations could completely annihilate an exponentially growing population.

Here, we derive quantitative estimates that allow us to compare the timescale of the extinction process to that of the deterministic peak in the SIR model. This is conceptually interesting in its own right, but it also allows us to meaningfully differentiate between “real” maxima and random transient peaks in the number of infected individuals in sub-populations that experience extinction.

In the pure birth–death process, the fraction of extinction events, 0ϕx1, that have already happened by time t can easily be calculated from Eq. (A11) as

ϕx(t)=P0(t)limtP0(t)=11k/be(bk)tk/b.
(B1)

This equation can be inverted to yield the time tx by which a fraction ϕx of extinction events have happened,

tx(ϕx)=1bklog(1ϕxk/b1ϕx).
(B2)

On the other hand, we can also estimate the fraction of non-extinct populations, 0ϕc1, that will still be below a cutoff size nc at time t,

ϕc(t,nc)=n=1nc1Pn(t)=1e(bk)tk/be(bk)t1(11k/be(bk)tk/b)nc.
(B3)

Evaluating ϕc(tx(ϕx),nc), therefore, yields the fraction of populations still below nc when a fraction ϕx of extinction events have already happened. This expression can be inverted to yield the simple relationship

nc(ϕx,ϕc)=1+log(1ϕc)logϕx,
(B4)

giving the number of infected individuals below which a fraction ϕc of non-extinct populations will still be at the time when a fraction ϕx of populations destined for extinction have already reached the extinct state.

In order to estimate the effect of extinction in our numerical simulations (cf. Fig. 4), we detect the maximum number of infected individuals in each sub-population (independent of their timing) and compare the sum of these numbers to our estimate Imaxext from the main text. In the sub-populations that experience random extinction of the epidemic, the detected numerical maxima will in reality be transient fluctuations before extinctions. These contribute more and more as R0=b/k1 when the deterministic peak value Nγ=N(1(1+logR0)/R0)20 decreases and the extinction probability 1/R0 increases. Using the estimates above, we can exclude these false maxima based on their timing by only considering those maxima for which

tmax>tx(ϕx)
(B5)

and simultaneously ensuring that

nc(ϕx,ϕc)<Nγ
(B6)

is fulfilled. ϕx and ϕc play the role of accuracy parameters. The first condition ensures that false maxima are excluded with probability ϕx, while the second one ensures that a pure birth–death process would not have reached the deterministic SIR peak by the same time with probability ϕc. Note that the latter is a conservative estimate, as growth in the SIR model is significantly slowed before reaching its peak compared to a pure birth–death process. In Fig. 4, we use a value of ϕx=ϕc=0.99 to exclude 99% of false maxima and still detect more than 99% of deterministic SIR maxima except for the data points marked as unreliable, for which Eq. (B6) is not fulfilled, and, therefore, the extinction process and the deterministic SIR peak are not clearly separated in time. Conversely, this also means that for all other parameters (i.e., larger R0=b/k), extinction usually happens well before the deterministic SIR dynamics reaches its peak.

It is worth emphasizing that, in the limit bk and small populations, the distinction between an initial stochastic phase and a deterministic time course becomes meaningless since γ eventually becomes order 1/N and the mean extinction time diverges. At this point, the dynamics throughout will be dominated by random growth of the number of infected individuals and stochastic fluctuations will continue to contribute, even as the number of susceptibles decreases, eventually ending the epidemic (i.e., during and beyond the maximum). In addition, the assumption that there is no depletion of susceptibles in the early phase (and thus the equivalence to a pure birth–death process) breaks down. However, in this study, we are interested in the regime where even sub-populations are still large and while b is sufficiently close to k to yield a significant extinction probability k/b, it is large enough to lead to a significant deterministic outbreak peak. Therefore, we do not investigate this regime.

The fact that the early phase of the dynamics in the SIR model (when SN and I is small) corresponds to a simple birth–death process also allows us to obtain an analytical estimate for the peak time distributions of the sub-populations. This can be readily adapted from a similar calculation performed on an equivalent problem in evolution, where the dynamics of a small mutant sub-population with a given selective advantage can likewise be understood as a birth–death branching process,26 for which the transition from the initial stochastic regime where extinction is still possible to the deterministic regime of exponential growth corresponds to the establishment of the mutation in the population (which precedes fixation).

We obtain an approximation for the establishment time distribution of the disease in a sub-population as

PSIRest(τ)=k(1k/b)exp((bk)τkbe(bk)τ),
(C1)

where we have corrected for an additional minus sign missing from Ref. 26. The variation in the timing of the later deterministic dynamics is due entirely to fluctuations in this initial stochastic phase. To compare this analytical approximation with our simulation results for the peak time in the main text, we plot the non-normalized, unconditional distribution,

P(tpeak)=[1(k/b)n]PSIRest(tpeak+τ¯tpeak),
(C2)

which is diminished by a factor [1(k/b)n] [from Eq. (C1)] accounting for the probability of extinction in a population with initially n infected individuals and has its mean shifted to the measured mean peak time tpeak. Here,

τ¯τ=1bkln(γkb),
(C3)

where γ=1.7810724 is the exponential of Euler’s constant.

We note that simply shifting the mean of the distribution is justified because the dynamics is predominantly identical in different sub-populations once they are in the deterministic regime, while only lagging by a random time span τ. This simple argument depends on the assumption that stochastic fluctuations can be ignored before deviations from exponential behavior (i.e., saturation effects) have to be considered for the deterministic dynamics. This is true for the scenarios we consider in the SIR model since our sub-populations still consist of thousands of individuals and we are explicitly focusing on cases where b is not arbitrarily close to k.

For estimating the peak reduction effect due to desynchronization of the sub-populations, it is convenient to work with the normalized equations for s=S/N and i=I/N, which read

s˙=bsi,
(D1a)
i˙=bsiki.
(D1b)

When i reaches its peak γ=i(tpeak), new infections and recovery balance according to Eq. (D1b) and s(tpeak)=k/b. Based on this known value, we use the following ansatz for s:

s(t)=kb(1+ε(t)),
(D2)

with ε(tpeak)=0. Since we are interested in the regime where there is a substantial extinction probability k/b, s(tpeak) is also still of order 1. Together with the fact that i˙(tpeak)=0 by definition, we expect from Eq. (D1a) that the lowest (linear) order of ε will suffice to describe the dynamics around the peak; i.e., ε(t)ε1(ttpeak) (conversely, we expect this approximation to break down when bk). Substituting the ansatz into Eq. (D1a) yields ε1=bγ or

ε(t)bγ(ttpeak).
(D3)

With this, we can obtain an approximation for i around the peak. From Eq. (D1b), we know that

dlog(i)dt=bs(t)k=kε(t),
(D4)

which can easily be solved. Together with the condition i(tpeak)=γ, we obtain

i(t)γexp(12bkγ(ttpeak)2).
(D5)

Now that we have an approximation for i(t) near the peak, we can calculate how these time courses add up across individual sub-populations by assuming that they all have the shape (D5), with the peak time tpeak stochastically distributed according to Eq. (C2). Defining

i¯(t)=limNs1Nsj=1Nsi(tpeakj;t),
(D6)

where each i(tpeakj;t) represents a time course as in Eq. (D5) with tpeakj drawn from the distribution (C2) for each j, we obtain an average superposition of many sub-populations in the limit Ns,

i¯(t)=dtpeakPSIRest(tpeak+τ¯)i(tpeak;t).
(D7)

Note that [as compared to Eq. (C2)] we use here the normalized distribution, without the diminishing factor due to extinction, in order to extract the reduction strictly due to desynchronization. We have also set tpeak to 0 without loss of generality, as a different value would simply shift i¯(t) by the corresponding time.

The integral in Eq. (D7) cannot be integrated in a closed form. We, therefore, replace PSIRest by a normal distribution N(0,σ2) with the same variance σ2=π2/[6(bk)2]. It is useful to note that, as for the normal distribution, the variance completely determines the shape of the Gumbel distribution in Eq. (C1), which means that the systematic error introduced by this replacement is parameter independent. Finally, we can calculate

i¯(t)=dtpeakexp(tpeak2/(2σ2))2πσi(tpeak;t)=γαexp(12bkγt2α2),
(D8)

with

α=1+π2bkγ6(bk)2.

The maximum of the resulting time course occurs at t=0 (due to our arbitrary choice of the mean for tpeak) and is i¯(0)=γ/α. Since the expected peak value without desynchronization is γ, desynchronization reduces this peak value by a factor of α1. According to Eq. (the definition of α above), α itself depends on γ, which in turn is a function of R0=b/k. Using the well known approximation γ=1[1+log(R0)]/R0,20 which is valid as long as SN initially, we rewrite α as

α=1+π2[R01log(R0)]6(R01)2.
(D9)

While we expect the quantitative estimate to be less accurate toward higher R0 (see above), we note that the important limits

limR01α=1,
(D10)
limR011α=1212+π20.7407
(D11)

exist. The first one signifies that there is no peak reduction due to desynchronization for R0, consistent with the disappearance of the stochastic phase at the beginning of the dynamics. The second limit indicates a finite reduction by a factor 0.7407 toward R0=b/k=1. Since the timescales of both the stochastic fluctuations and the deterministic peak behavior diverge for R01 (and are ill-defined for R0=1), this means that they must exhibit identical scaling behavior in order for neither of them to dominate. In between the two extremes, 1/α increases monotonically with R0, which implies that the maximum reduction that can be achieved by desynchronization is about 26% and is reached close to R0=1. It is important to note, however, that several assumptions even about the deterministic time course (for example, the value of γ) break down when R0 is so close to 1 that γ becomes of order 1/N; therefore, a fully stochastic treatment would be needed to fully capture this regime. This does not limit the validity of the results in the regime we are interested in, i.e., where sub-populations still exhibit clear deterministic outbreaks (or extinction).

1.
W. O.
Kermack
and
A. G.
McKendrick
, “
A contribution to the mathematical theory of epidemics
,”
Proc. R. Soc. Lond. Ser. A
115
,
700
721
(
1927
).
2.
J.
Dehning
,
J.
Zierenberg
,
F. P.
Spitzner
,
M.
Wibral
,
J. P.
Neto
,
M.
Wilczek
, and
V.
Priesemann
, “
Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions
,”
Science
15
,
eabb9789
(
2020
).
3.
H.
Heesterbeek
,
R. M.
Anderson
,
V.
Andreasen
,
S.
Bansal
,
D.
De Angelis
,
C.
Dye
,
K. T. D.
Eames
,
W. J.
Edmunds
,
S. D. W.
Frost
,
S.
Funk
,
T. D.
Hollingsworth
,
T.
House
,
V.
Isham
,
P.
Klepac
,
J.
Lessler
,
J. O.
Lloyd-Smith
,
C. J. E.
Metcalf
,
D.
Mollison
,
L.
Pellis
,
J. R. C.
Pulliam
,
M. G.
Roberts
,
C.
Viboud
, and
Isaac Newton Institute IDD Collaboration
, “
Modeling infectious disease dynamics in the complex landscape of global health
,”
Science
347
,
990
993
(
2015
).
4.
L.
Hufnagel
,
D.
Brockmann
, and
T.
Geisel
, “
Forecast and control of epidemics in a globalized world
,”
Proc. Natl. Acad. Sci. U.S.A.
101
,
15124
15129
(
2004
).
5.
J. T.
Matamalas
,
A.
Arenas
, and
S.
Gómez
, “
Effective approach to epidemic containment using link equations in complex networks
,”
Sci. Adv.
4
,
eaau4212
(
2018
).
6.
R.
Pastor-Satorras
,
C.
Castellano
,
P.
Van Mieghem
, and
A.
Vespignani
, “
Epidemic processes in complex networks
,”
Rev. Mod. Phys.
87
,
925
979
(
2015
).
7.
Metapopulation Dynamics: Empirical and Theoretical Investigations, edited by M. Gilpin and I. Hanski (Academic Press, London, 1991).
8.
J.
Dushoff
and
S.
Levin
, “
The effects of population heterogeneity on disease invasion
,”
Math. Biosci.
128
,
25
40
(
1995
).
9.
A.
Lunelli
and
A.
Pugliese
, “
Final attack ratio in SIR epidemic models for multigroup populations
,”
Ric. Mat.
67
,
49
68
(
2018
).
10.
V.
Andreasen
and
F. B.
Christiansen
, “
Persistence of an infectious disease in a subdivided population
,”
Math. Biosci.
96
,
239
253
(
1989
).
11.
H.
Andersson
and
T.
Britton
, Stochastic Epidemic Models and Their Statistical Analysis, Springer Lecture Notes in Statistics Vol. 151 (Springer-Verlag, New York, 2000).
12.
M.
Hartfield
and
S.
Alizon
, “
Introducing the outbreak threshold in epidemiology
,”
PLoS Pathog.
9
,
e1003277
(
2013
).
13.
R.
Huerta
and
L. S.
Tsimring
, “
Contact tracing and epidemics control in social networks
,”
Phys. Rev. E
66
,
056115
(
2002
).
14.
L. J. S.
Allen
, “
A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis
,”
Infect. Dis. Modell.
2
,
128
142
(
2017
).
15.
I.
Tiwari
,
P.
Sarin
, and
P.
Parmananda
, “
Predictive modeling of disease propagation in a mobile, connected community using cellular automata
,”
Chaos
30
,
081103
(
2020
).
16.
M. J.
Keeling
, “
Metapopulation moments: Coupling, stochasticity and persistence
,”
J. Anim. Ecol.
69
,
725
736
(
2000
).
17.
B.
Grenfell
, “
(Meta)population dynamics of infectious diseases
,”
Trends Ecol. Evol.
12
,
395
399
(
1997
).
18.
The combination of traditional containment measures with relative regional isolation is the subject of a separate study.25
19.
D. T.
Gillespie
, “
Exact stochastic simulation of coupled chemical reactions
,”
J. Phys. Chem.
81
,
2340
2361
(
1977
).
20.
H.
Weiss
, “
The SIR model and the foundations of public health
,”
Mater. Mat.
3
,
1
17
(
2013
), ISSN: 1887-1097.
21.
P.
Whittle
, “
The outcome of a stochastic epidemic—A note on Bailey’s paper
,”
Biometrika
42
,
116
122
(
1955
).
22.
Note that, for simplicity of presentation, we are using the long time limit of the extinction probability to define p1ext rather than P0(t). This simplification implies that extinction happens fast enough on the timescale relevant for our problem. It can be justified by comparing the timescale of the extinction process to that of the infection peak in the SIR model (see further below in the main text and in  Appendixes A– D.
23.
L.
Ferretti
,
C.
Wymant
,
M.
Kendall
,
L.
Zhao
,
A.
Nurtay
,
L.
Abeler-Dörner
,
M.
Parker
,
D.
Bonsall
, and
C.
Fraser
, “
Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing
,”
Science
368
,
eabb6936
(
2020
).
24.
Q.
Li
,
X.
Guan
,
P.
Wu
,
X.
Wang
,
L.
Zhou
,
Y.
Tong
,
R.
Ren
,
K. S. M.
Leung
,
E. H. Y.
Lau
,
J. Y.
Wong
,
X.
Xing
,
N.
Xiang
,
Y.
Wu
,
C.
Li
,
Q.
Chen
,
D.
Li
,
T.
Liu
,
J.
Zhao
,
M.
Liu
,
W.
Tu
,
C.
Chen
,
L.
Jin
,
R.
Yang
,
Q.
Wang
,
S.
Zhou
,
R.
Wang
,
H.
Liu
,
Y.
Luo
,
Y.
Liu
,
G.
Shao
,
H.
Li
,
Z.
Tao
,
Y.
Yang
,
Z.
Deng
,
B.
Liu
,
Z.
Ma
,
Y.
Zhang
,
G.
Shi
,
T. T. Y.
Lam
,
J. T.
Wu
,
G. F.
Gao
,
B. J.
Cowling
,
B.
Yang
,
G. M.
Leung
, and
Z.
Feng
, “
Early transmission dynamics in Wuhan, China, of novel coronavirus—Infected pneumonia
,”
N. Engl. J. Med.
382
,
1199
1207
(
2020
).
25.
P.
Bittihn
,
L.
Hupe
,
J.
Isensee
, and
R.
Golestanian
, “Local measures enable COVID-19 containment with fewer restrictions due to cooperative effects,” medRxiv 2020.07.24.20161364 (2020).
26.
M. M.
Desai
and
D. S.
Fisher
, “
Beneficial mutation-selection balance and the effect of linkage on positive selection
,”
Genetics
176
,
1759
1798
(
2007
).