Power-law population heterogeneity governs epidemic waves.

We generalize the Susceptible-Infected-Removed (SIR) model for epidemics to take into account generic effects of heterogeneity in the degree of susceptibility to infection in the population. We introduce a single new parameter corresponding to a power-law exponent of the susceptibility distribution at small susceptibilities. We find that for this class of distributions the gamma distribution is the attractor of the dynamics. This allows us to identify generic effects of population heterogeneity in a model as simple as the original SIR model which is contained as a limiting case. Because of this simplicity, numerical solutions can be generated easily and key properties of the epidemic wave can still be obtained exactly. In particular, we present exact expressions for the herd immunity level, the final size of the epidemic, as well as for the shape of the wave and for observables that can be quantified during an epidemic. In strongly heterogeneous populations, the herd immunity level can be much lower than in models with homogeneous populations as commonly used for example to discuss effects of mitigation. Using our model to analyze data for the SARS-CoV-2 epidemic in Germany shows that the reported time course is consistent with several scenarios characterized by different levels of immunity. These scenarios differ in population heterogeneity and in the time course of the infection rate, for example due to mitigation efforts or seasonality. Our analysis reveals that quantifying the effects of mitigation requires knowledge on the degree of heterogeneity in the population. Our work shows that key effects of population heterogeneity can be captured without increasing the complexity of the model. We show that information about population heterogeneity will be key to understand how far an epidemic has progressed and what can be expected for its future course.

1 Introduction

Diseases that spread by transmission between individuals can give rise to epidemic waves that pass through a population [1, 2]. One infected person can infect several others who are susceptible to the infection, characterized by the basic reproduction number R0, initially typically generating an exponential growth of the number of infections. The number of infections reaches a peak and later dies down when there is a sufficient number of individuals that have gained immunity after they recovered from the infection so that further growth is hampered. The fraction of immune individuals reached at the point when the epidemic starts to recede is called herd immunity [1, 3, 4].

There are big uncertainties as to when and why an epidemic reaches its peak and the levels of herd immunity required [5]. Simple models of infections dynamics predict that for an initially fast growing epidemic most of the population will become infected before the epidemic dies down [1, 3, 6]. It was noted early by William Farr when investigating smallpox and other epidemics that epidemics appear to follow a general time course in the form of a skewed bell shaped curve [7, 8]. They first grow fast, reach a peak and then die down quickly, typically much before the majority of a population has been affected. The fact that an epidemic dies down is usually attributed to the fact that there exists some degree of immunity in the population [9]. The uncertainty about when the peak of an epidemic is reached and why an epidemic dies out even if there remains a large number of still susceptible individuals reveals that the factors that limit an epidemic are not well understood. Furthermore, the effectiveness and impact of mitigation measures such as social distancing to counter a fast growing epidemic are not known.

Simplified models of infection dynamics, such as the classic Susceptible-Infected-Removed (SIR) model have been used for a long time to describe the dynamics of epidemics spreading through a population [1, 3, 6, 10]. Such models capture key features of the epidemic as a nonlinear wave with qualitative properties that match observed bell-shaped dynamics of epidemic waves. However, more quantitatively, such models exhibit the robust feature that a quickly growing epidemic does not stop unless the majority of a susceptible population has reached immunity after going through the infection [1]. This raises the question whether important factors are missing in these simple and elegant models. To understand at what conditions and at what levels epidemic waves become self-limiting and die down remains an important challenge. This aspect is also key to understand the role and effectiveness of social distancing measures to influence dynamics of an epidemic wave [10-12].

Simple epidemic models treat the population as effectively consisting of identical individuals. However, individuals in a population can differ widely. The importance of population heterogeneity was put forward to understand smallpox epidemic which could not be captured by simple models [13]. Such heterogeneity has been taken into account by adding details such as introducing several compartments to a model [14] or by introducing distributions of susceptibility [13, 15, 16] or infectiousness [15-17]. It was suggested that population heterogeneity reduces effective herd immunity levels [13, 16, 18, 19].

In this paper, we present a generalization of the SIR model that takes into account effects of population heterogeneity. We show here that effects of heterogeneity can be added without losing the simplicity of the SIR model and keeping its mathematical structure. We introduce a single new parameter, the susceptibility exponent α, which characterizes a generic power-law heterogeneity in the distribution of infection susceptibilities of the population. Power laws are often found in nonlinear and complex systems [20-23]. In the present context, power laws could be expected for example based on a variability of immune responses of different individuals which could imply a wide variability in the efficiency of the transmission of an infection [24, 25]. Furthermore, population heterogeneity could be relevant at very different scales, from the immune response of cells to the behaviors of individuals that affect infection rates. Such as broad range of relevant scales could give rise to approximately scale free properties or power laws.

In the heterogeneous SIR model proposed here, the qualitative behaviors of the epidemic wave are unchanged. However, as a function of the parameter α, the wave can become self-limited at much lower levels of infected individuals as compared to the classic SIR model. In the limit of large α we recover the classic SIR model of homogeneous populations. For smaller α we find that the number of infections at the peak and the cumulative number of infections after the epidemic has passed can be strongly reduced. Our work has implications for the concept of herd immunity and clarifies that herd immunity cannot be discussed independently of population heterogeneity.

We discuss the dynamics of the SARS-CoV-2 pandemics using the heterogeneous SIR model applied to data on reported infection numbers and COVID-19 associated deaths in Germany [26]. We estimate parameter values including the susceptibility exponent α and show that the time course observed in Germany is compatible with different scenarios ranging from a homogeneous population strongly affected by mitigation to a self-limited epidemic wave in a heterogeneous population where social distancing measures play a minor role.

II. The Susceptible-Infected-Removed model

The Susceptible-Infected-Removed (SIR) model captures key features of a spreading epidemic as a mean field theory based on pair-wise interactions between infected and susceptible individuals. This model captures generic and robust features without aiming to describe specific details. In the presence of I infected individuals in a population of N individuals, the infection can be transmitted to susceptible individuals. They stay infectious during an average time γ−1 after which they no longer contribute to infections. The number of susceptible individuals S and the number of infected individuals I obey where the dots denote time derivatives. The infection rate is denoted β and can in general depend on time t. This time dependence could correspond to seasonal changes or mitigation measures [10, 12, 27]. The infection rate can be modulated by the average infection susceptibility which we introduce to capture effects of population heterogeneity discussed below. In the classical SIR model . The number of removed (or recovered) individuals is given by N − S − I and the cumulative number of infections is C = N − S. A key parameter is the basic reproduction number which denotes the average number of new infections generated by an infected individual. The growth rate of infections is , where is a time dependent reproduction number.

The time course of an epidemic is often provided as the number of new cases per day. This corresponds to the rate of new infections per unit time with and R = J/(γI).

A. Infection dynamics in homogeneous populations

In the simple case of a homogeneous population, all individuals have the same degree of susceptibility, x = 1 and the population average of x is independent of time. This is the classic SIR model. An example for a solution to these equations for homogeneous population and constant β is given in Fig 1(a) and 1(b). The corresponding time dependent reproduction number is presented in Fig 1(c). The number of infections first grows exponentially with growth rate

Fig 1

Effects of population heterogeneity on the dynamics of SIR models.

Examples for the time course of fraction of susceptible S/N (green), fraction of infected I/N (orange) and fraction of the cumulative number of infected C/N (blue) in a SIR model with N total individuals. (a)-(c) Homogeneous SIR model with R0 = 2.5 and γ = 0.13 day−1. (d)-(f) Heterogeneous SIR model with same R0 and γ and with α = 0.1. (a) and (d) show time course as linear plot, (b) and (e) show semi logarithmic plots of the same variables. (c) and (f) show the normalized time dependent reproduction number R(t)/R0 (yellow) and the average susceptibility (purple) as a function of time. The dotted lines in (a),(b),(d) and (e) indicate the herd immunity level C. Other parameters: N = 8107 individuals and I0 = 10 initially infected.

As the number of susceptible decreases, the epidemic reaches a peak number of infected Imax = I(t) at time t = t with and R(t) = 1. At this peak, a fraction S/N = 1/R0 of individuals remain susceptible. The cumulative number of infections C at the maximum of I thus obeys

Eq (6) is the classic herd immunity level which is the fraction of immune individuals in the population beyond which the epidemic can no longer grow. Finally the epidemic dies down exponentially with rate where is the fraction of susceptible individuals that remain after long times, see Appendix A. Here W(z) denotes Lambert W-function [28]. The total fraction of infections over the course of the epidemic is C∞/N = 1 − S∞/N.

For a classic SIR model with homogeneous population we have for R0 = 2.5, a herd immunity level C/N of 60% of the population, see Fig 1(a) and 1(b). After the infection has passed C∞/N ≃ 89% of the population have been infected, see Fig 2(a) and 2(b) (green lines). The fraction of the population that become infected increase for larger R0. The SIR model thus suggests that for R0 > 2 the epidemic wave exceeds a majority of the population before the epidemic begins to die out.

Fig 2

Fraction of susceptible individuals at long times.

(a) Fraction S∞/N of susceptible individuals that remain at long times as a function of the basic reproduction number R0 for different degrees of population heterogeneity characterized by the values of α. The limit α → ∞ corresponds to the classic case with homogeneous populations (green). In the limit α → 0 populations are most heterogeneous (blue). (b) Fraction S/N of susceptible individuals as a function of R0 at the peak where the number of infected is maximal for different α. (c) Ratio of infections after the peak S − S∞ and infections before the peak N − S as a function of R0.

B. Infection dynamics with population heterogeneity

Not all individuals are the same and for some susceptible individuals the probability of infection per time is lower than for others. This can be captured by a distribution of susceptibilities x [13, 15, 16]. We denote s(x)dx the number of individuals with susceptibility between x and x+ dx. The total number of susceptible individuals is then . For each sub-population s(x) with susceptibility x, the number of susceptible individuals decreases as which for the whole population implies Eq (1) with average susceptibility which is in general time dependent.

This time dependence can be discussed by introducing the variable τ that increases monotonically during the epidemic wave. It starts with τ = 0 and is defined via the equation , which implies that it reaches a final value when the epidemic has decayed. Therefore τ can be interpreted as a measure of how far the epidemic has progressed. Eq (9) can then be written as ∂ s = −xs, and the number of susceptible individuals is where s0(x) is the initial susceptibility distribution at time t = t0 with average , see Appendix B.

C. Infection dynamics with generic power law heterogeneity

The dynamics of epidemic waves depends on the shape of the initial distribution s0(x). Here, we consider distributions that have the special property of shape invariance under the dynamics of epidemics. This property is satisfied by a gamma distribution which is governed by a power-law at small x, characterized by the exponent α > 0, and a cut off at large x, see Eq (C1) in Appendix C. The distribution s0(x) has average and variance 1/α. Indeed, we have , where the time dependence enters via , see Appendix C. This shape invariance implies that the gamma distribution is maintained at all times and is not merely an initial condition. Furthermore, starting with any initial distribution that exhibits a power law s0(x)∼x−1+ at small x, the distribution will converge for large τ to the shape invariant gamma distribution, which therefore is an attractor of the dynamics, see Appendix B. Note that in the limit of large α, we recover the classic SIR model for a homogeneous population. For small α, the population is strongly heterogeneous.

By inserting the distribution given in Eq (12) into Eq (11) we obtain as detailed in Appendix C. The average susceptibility is which starts from for τ = 0 and decreases for increasing τ, thus dampening the epidemic. We can now express the dynamics given in Eqs (1) and (2) as two dynamic equations for I(t) and τ(t) which read with initial values I(0) = I0, τ(0) = 0 and S(0) = N − I0. The number of susceptible individuals at time t is then given by S(t) = S(τ(t)). An example of a time course of this model for α = 0.1 is shown in Fig 1(d), 1(e) and 1(f).

We can discuss how the shape of the epidemic wave depends on the parameter α. The epidemic starts out with exponential growth of infected individuals at rate λ0 = γ(R0 − 1), with R0 = β/γ. The time dependent reproduction number is

When the reproduction number drops to R = 1, the number of infected reaches a maximum

Beyond the herd immunity level given by the cumulative number of infections at the maximum of I the reproduction number R drops below 1 and the epidemic dies down. In Eqs (17)–(19) we have considered the limit of small I0/N for simplicity. The general derivation is given in Appendix C. In the limit of large α, these expressions converge to those obtained for the homogeneous SIR model, see Appendix C. The remaining fraction of susceptible individuals at the peak and after the epidemic has passed is shown as a function of α in Fig 2(a) and 2(b). This reveals that as α is reduced, the fraction of the population reached by the epidemic decreases and can become very small for small α. At the same time the infections are more spread out over time and a larger fraction occurs after the peak when α is reduced, see Fig 2(c).

An important case is a strongly heterogeneous population. For small α ≪ 1, we obtain simple analytical expressions for the behavior of the system, see Appendix E. In this limit we have and C/N ≃ α ln R0. An important quantity is the rate J of new cases per time. For small α it takes the maximal value

The final number of susceptible individuals is given by where for small α the average susceptibility after the infection has passed is . Here W−1(z) denotes the −1 branch of Lambert W function. We finally have for small α

A key result is that for small α the herd immunity level can be much below the classical value suggested by the SIR model. For example for R0 = 2.5 and α = 0.1, we have Imax/N ≃ 2.8%, and a fraction C/N ≃ 8% of infected individuals required for herd immunity, much lower than is usually suggested. The total number of infected at long times is C∞/N ≃ 14%, see Appendix C. The reason for the small amplitude of the epidemic wave in a heterogeneous population (see Fig 1(d) and 1(e)) as compared to the amplitude for a homogeneous population (see Fig 1(a) and 1(b)) is the drop of the average infection susceptibility (see Fig 1(f)). This drop occurs because in a heterogeneous population the most susceptible individuals are first removed, thereby lowering the average susceptibility.

III. Application to the SARS-CoV-2 epidemic in Germany

We analyze the dynamics of the SARS-CoV-2 epidemic in Germany using public data provided by the Robert Koch institute [26]. These daily reports provide the numbers of reported positive tests for each day, but also the dates of reporting of those infections which later turn out as fatal. The total number of new reported infections per day (red symbols) are shown in Fig 3(a) together with the number of reported infections per day that were later fatal (blue symbols), which we denote . Both sets of data can be interpreted as proxies for the rate J of new cases per day up to an unknown factor. They show qualitatively similar behavior, a rapid growth and a decline after passing a maximum. However there are quantitative differences, in particular the growth rates at early and late times, given by the slopes of the data in a single logarithmic plot are different, see Fig 3(b). The number of new cases per day that are later fatal is related to the number of new infections per day as , where f denotes the infection fatality rate, the fraction of infections that are fatal, which we consider to be constant for simplicity.

Fig 3

(a) Daily new SARS-CoV-2 infections reported in the early months of 2020 in Germany. The number of new reported infections per day (red symbols) is shown together with the number of reported infections per day for those cases with later fatal outcome (blue symbols). (b) Semi logarithmic representation of the same data. The dashed and solid lines represent linear and cubic fits to the data in specific time intervals. They are used to estimate the initial and final growth rates λ0 and λ∞ as well as and at the maximum of the rate of new cases Jmax. We find λ0 ≃ 0.269 day−1 (0.336 day−1), λ∞ ≃ −0.068 day−1 (−0.038 day−1), A2 ≃ −10−2 day−2 (−0.9110−2 day−2) and A3 ≃ 6.810−4 day−3 (7.510−4 day−3) for the fatal cases (for all reported cases).

A. Comparison to heterogeneous SIR model

The calculated number of new cases per day J obtained as solution to Eqs (15) and (16) for a heterogeneous population and scaled by the factor f to match the data of fatal cases are shown in Fig 4(a)–4(d) as solid blue lines. These lines are shown together with the number of new fatal cases per day as blue symbols. The factor f was determined such that the cumulated cases per day matches the cumulative number of cases C on June 15. The time axis is chosen such that the model matches the data. From a fit of the model to the data we obtain the parameter estimates R0 ≃ 2.67 and γ ≃ 0.146 day−1. Good fits to the data are found for a range of α sufficiently small, about α < 0.2. The resulting infection fatality rates f vary as α is changed. Using α = 0.05 corresponds to an infection fatality rate f ≃ 0.13%. It could be larger or smaller if a different value of α was used. This would not significantly affect the quality of the fit as long as α < 0.2. The calculated time courses I(t), J(t) and C(t) corresponding to these fits to are shown in (e). The dependence of the average susceptibility on time and the function τ(t) are shown in (f). The increase of τ with time represents the advance of the epidemic. It reaches a final value at long times as discussed in Appendix C. The inset in (f) shows the shape of the distribution of susceptibility in the population at different stages characterized by different values of τ.

Fig 4

Time course of the SARS-CoV-2 epidemic in Germany (symbols) compared to solutions of the heterogeneous SIR model (lines).

(a) and (b) Data on daily reported infections (red) and on reported infections with later fatal outcome (blue) as logarithmic and linear plots. (c) and (d) same data and model solutions as in (a) and (b) but for cumulative numbers of cases. The horizontal dashed lines indicate scaled herd immunity levels. Parameter values for the model solution are R0 = 2.67 (R0 = 3.91), γ = 0.146 (γ = 0.069), α = 0.05 and N = 8107 for the fatal cases (for all cases). The case fatality rate that corresponds to this solution is f = 0.13% (f = 0.11%). (e) Time courses of the fraction of infected I/N (blue), the new cases per day J/N (red) and the fraction of cumulative cases C/N (yellow) for R0 = 2.67 and γ = 0.146. (f) Time course of the average susceptibility (blue), where R is the time dependent reproduction number and of (red) for the solution shown in (e). Inset: distributions of susceptibility in the population for different values of τ.

It is surprising that the model fits the data of fatal cases with just two fit parameters while yielding a reasonable infection fatality rate. This is further clarified when using the fit values of R0 and γ to calculate λ0 ≃ 0.24 day−1, slightly smaller than the estimate given in Fig 3(b). Using Eq (22), we also find λ∞ ≃ −0.069 day−1, very close to the estimate from the data. The data of all reported cases can also be captured by the model for small α, see Fig 4(a) and 4(c) red symbols and red lines. This fit is not as close and the parameter values are different, see Fig 4. Our comparison of the model to the data shows that the model captures the time course of fatal cases surprisingly well for the case of strong heterogeneity for infection fatality rates that fall in the range of estimates from immunological studies [29-33].

B. Quantification of the shape of the epidemic wave

In order to understand how the shape of the wave of infections constrains the possible parameter values of R0, γ and α, we consider in addition to the initial growth rate λ0 and the final decay rate λ∞ two coefficients describing the epidemic dynamics near its peak, using the expansion where the linear term disappears by definition at the maximum Jmax = J(t) at time t. The coefficients and can be obtained for the homogeneous and heterogeneous SIR model, see Appendix C, D and E. Fig 5 shows dimensionless combinations of these values as a function of R0 for different α ranging from the homogeneous case α → ∞ to strongly heterogeneous with α → 0 as solid lines of different color. The values obtained from the fits shown in Fig 3 are indicated as dashed lines together with shaded regions corresponding to estimated uncertainty ranges of these values.

Fig 5

Role of population heterogeneity for the behavior of the generalized SIR model.

Plots of various dimensionless ratios of parameters characterizing the shape of the infection wave for different values of α. Here the limit α → ∞ corresponds to the homogeneous SIR model, the limit α → 0 to the strongly heterogeneous case. Here λ0 and λ∞ denote the initial and final growth rate, and describe the shape of the wave at the maximum of new cases per day J. The horizontal dashed lines correspond to estimates from fits shown in Fig 3, the shaded regions indicate uncertainty ranges, see Appendix D.

We find that the ratio , which is independent of γ depends only weakly on α. We can therefore use it to estimate R0, see Fig 5. Using A2 ≃ −0.01 day−2 and λ∞ ≃ −0.07 day−1 determined from the data of fatal cases, we have leading to the estimate R0 ≃ 2.5, see Fig 5. This estimate can now be used to infer bounds on α. The ratio λ∞/λ0 ≃ −0.3 is only consistent with R0 ≥ 2.6 and the lower value corresponds to the limit of small α, see Fig 5. This reveals that α ≪ 1 must be small and that the classic SIR model with homogeneous population is not consistent with this data. We can now estimate γ using the small α limit. For R0 ≃ 2.6, we have γ ≃ 2λ∞ ≃ 0.14 day−1, see Fig 5(d). The data does not provide information about the true total number of infections. Therefore the precise value of α remains unknown. We can use estimates from immunological studies estimating the number of infections [29, 31, 32] to determine α. This suggests a range of about 0.01 < α < 0.15, corresponding to 0.65%>f > 0.04%. Fig 5 also shows the estimated ranges for data on all reported cases in red. For this case the inferred values of R0 is larger and the consistency with the data is less strong.

C. Effects of mitigation and social distancing measures

During an epidemic conditions can change over time. For example, mitigation by social distancing measures, quarantining or seasonal changes could affect how quickly an infection spreads on average from person to person. Given that such changes are global, they may be captured by a time-dependence of the rate β(t) [10, 12, 27]. In the following, we discuss scenarios of mitigated epidemics, starting from a reference point with an initial infection rate β0 prior to mitigation. We use this reference to define the herd immunity C of the population via Eq (19). The herd immunity level depends on the basic reproduction number R0 = β0/γ and on the population heterogeneity α. For immunity levels above herd immunity, C ≥ C, the population is stable after mitigation measures are completely relaxed and β is restored to its original value β0.

We examine three different scenarios with a comparable total number of infections. These scenarios are shown in Fig 6. They are characterized by different levels of immunity relative to herd immunity at July 1 and thus differ in the future behaviors beyond this time. Starting in all scenarios with R0 = 2 and using γ = 0.24, the model follows the initial growth at rate λ0 of the reported cases. If β is kept constant, β = β0, the model deviates from the data at later times, see dotted lines in Fig 6. If β is permitted to change in time, almost any reported time course could be described by the model. We use the data to infer a time course of β(t) such that the model follows the data, see Appendix G. The inferred values of β are shown as circles in Fig 6(c), 6(f) and 6(i). In order to fit the model to the data in different mitigation scenarios, we use a piecewise linear modulation of β. The time dependence of β(t) that resulted from these fits are shown in Fig 6(c), 6(f) and 6(i) as solid lines. The value of β decreases sharply at the onset of mitigation. After this decrease it stays roughly constant or increases at constant rate, thus relaxing mitigation. The magnitude of maximal mitigation and the two slopes of β(t) were used as fit parameters.

Fig 6

Scenarios of mitigation.

(a)-(c) Early mitigation by strong reduction of β for a homogeneous population (large α limit). The new cases per day are shown in (a) as symbols. A fit of the mitigated model is shown as solid lines. The solution for same parameter values R0 = 2 and γ = 0.24 but without mitigation is shown as dotted lines. The corresponding cumulative numbers of cases are shown in (b). Herd immunity levels corresponding to these solutions are indicated as horizontal dashes lines. The time dependence of β(t) are shown in (c) as solid lines. The time courses of β inferred from the data is shown as symbols. Mobility data indicating social activities in Germany relative to baseline values are shown in orange for comparison. (d)-(f) same plots as in (a)-(c) but for a moderate mitigation and heterogeneous population with R0 = 2, γ = 0.24 and α = 0.1. (g)-(i) Heterogeneous population with mild mitigation and release leading to almost herd immunity. Red symbols and lines correspond to the case of all reported infections, blue data and lines correspond to reported infections of fatal cases.

In the case of early mitigation, Fig 6(a)–6(c), fast reduction of β suppresses the epidemic before any appreciable progress towards herd immunity was made. Mitigation needs to be strong and sustained to be compatible with the data. By July, the population reaches only about 6−7% of herd immunity in this case. Note that this is the only scenario of the classic SIR model with a homogeneous population (α → ∞) that could be compatible with the data.

For heterogeneous populations with α ≲ 0.2, scenarios with milder mitigation and with infection levels closer to herd immunity are compatible with the data, see Fig 6d and 6g. A case of moderate mitigation with α = 0.1 is shown in Fig 6d–6f. The population in this case reaches by July 1st about ∼45% of herd immunity. A sustained mitigation is needed to account for the data, albeit with smaller magnitude compared to the first case. If the epidemic starts slightly earlier (3 days for the case shown in Fig 6g–6i as compared to (d-f)), the population reaches ∼95% of herd immunity by July 1st. Here, mitigation has the effect to reduce the cumulative number of infections as compared to a non-mitigated case (C/N = 5.8% compared to 11% by July 1st). This reduction of cumulative infections C results from a reduction of the number of infectious individuals I at the point when herd immunity is reached. In the absence of mitigation, I reaches its maximum when C = C, whereas mitigation can reduce I to small numbers as herd immunity is reached, preventing further infections. The minimal number of infections that can be achieved by temporary mitigation is C, which is up to 50% smaller than the long lime limit C∞ in an unmitigated epidemic (see Fig 2c).

The scenarios of temporary reduction of β could capture the mitigation effects of social distancing measures. To relate the inferred time dependence of β to measures of social activity, we show in Fig 6(c), 6(f) and 6(i), β(t) together with mobility data from Ref. [34] for comparison, see Appendix H. This mobility data shows a sharp decline and a slow but steady return to the initial state roughly in line with inferred changes of β(t).

The three scenarios differ in the fraction of herd immunity they reach by July 1 and therefore in their future trajectories. However, β(t) was adjusted by a fitting procedure such that all scenarios are consistent with the data on reported infections. This reveals that it can be difficult to distinguish effects of heterogeneity leading to a time dependent average susceptibility from mitigation effects corresponding to time-dependent β. Indeed our analysis shows that changes in mitigation strength can be compensated to some degree by changes of heterogeneity described by α.

IV. Conclusions and perspectives

We have presented a generalization of the classic Susceptible-Infected-Removed model for epidemic waves, which adds one new parameter to the model that captures population heterogeneity by a power-law exponent α. This exponent describes the power law that characterizes the distribution of susceptibility in the population s(x) ∼ x−1+ for small x. A special case for such distributions is the gamma distribution. Gamma distributions have been used before to describe heterogeneous populations [13, 15–17] and an approach similar to the one presented here has been proposed in [15] where also the shape invariance of gamma distributions is mentioned. Here, we have shown that gamma distributions have the special properties that they are both shape invariant under the dynamics and attractors of the dynamics for power-law distributions. This implies that for each α there exists a class of distributions with the same power law at small x which share the same limiting dynamics and distribution. The generalization of the SIR model introduced here captures the effects of these power laws by the parameter α in a generic way. Note that this generalization does not change the simplistic nature of the SIR model and does not change its numerical or analytical complexity.

For α > 1, population heterogeneity is weak and in the limit of large α, one recovers the classic SIR model of homogeneous population, see Fig 1(a)–1(c). For α < 1 population heterogeneity plays a key role in limiting the peak of the epidemic wave. We show that as a result of strong population heterogeneity (small α), the wave peaks when only a small minority of individuals have been infected, see Fig 1(d)–1(f). The herd immunity level, the point where the epidemic dies down spontaneously becomes very small for small α, see Eq (19). Thus our model shows that for small α, an epidemic wave can die out after reaching only a small fraction of the population even though a majority of the population is still susceptible. In this case the population is stable with respect to introducing new infected individuals because the average susceptibility has dropped significantly, see Fig 1(f).

Many properties of the nonlinear wave in this generalized model can be obtained exactly as a function of α and in the limit of small α. Numerical solutions can be generated quickly and efficiently. In a heterogeneous population the average susceptibility stays almost constant at early stages of the epidemic where the number of new cases grows exponentially with rate λ0 = γ(R0 − 1). At this stage the dynamics is the same as in the classic model and independent of α. However, then drops rather quickly and the epidemic waves thus reaches its peak and dies down, see Fig 1(b) and 1(f). This sudden drop in average susceptibility results from a shift of the distribution of susceptibility. The most susceptible individuals are removed from the dynamics at higher rates than those with low susceptibility. This leads to a rapid reduction of the average susceptibility until it has dropped to a low value where the time dependent reproduction number R falls below 1, see Eq (17). The wave then dies down at rate λ∞ and the average susceptibility approaches a final value . Thus the qualitative behavior of the classic SIR model is unchanged and the key parameters, the recovery rate γ and the basic reproduction number R0 have the same values and properties. However, the power-law distribution of susceptibility can dramatically change the peak of the epidemic and alters the precise shape of the wave. The dynamics effectively shifts the edge of the susceptibility distribution, see inset in Fig 4(f), which changes the stability of the population from prone to an exponentially growing wave to a stably decaying wave without requiring a large number of infections.

The simple SIR model does not aim to capture details such as the population structure, the geography or human travel. In the spirit of statistical physics it is based on the idea that the collective behaviors of many individuals give rise to an emergent epidemic wave with robust and generic features that can be captured by a simplified model that focuses on key aspects. Here we show that power-law heterogeneity is a key factor that should not be left out. Note that our approach to take into account population heterogeneity can also be applied to extensions and generalizations of the SIR model such as the SEIR and the SIRS models [16, 35–37].

We apply our model to data on the SARS-CoV-2 epidemic in Germany in 2020. The data from Germany provides time stamps on reporting dates of infections and reporting dates of infections that later are fatal. Surprisingly, the data for fatal cases is well described by the heterogeneous SIR model with constant parameters and small α but not by the classic SIR model with constant parameters. In the case of SARS-CoV-2, immunological data suggests that only a minority of the population exhibits antibodies [29, 31, 32]. This is consistent with a fit of our model to the data using a small value of α. The data on all reported cases can also be captured by the model, but the fit is less convincing. Comparing the data on all reported cases to the data on the time course of cases that are fatal reveals some differences. Clearly the fatal cases represent a different sampling as these cases correspond to predominantly old individuals and therefore measure a different quantity. However, starting from all reported cases and then using the fatal outcome as a second criterion could reduce biases due to changes in testing rates, testing strategies as well as testing errors.

An epidemic wave does not progress under constant conditions but is subject to changes such as mitigation measures and seasonal effects. We use our model in comparison with the data from Germany to investigate different scenarios of mitigation that correspond to different level of immunity in the population. In the case of a homogeneous population the data can only be accounted for if mitigation is strong and suppresses the epidemic far below herd immunity Fig 6(a)–6(c). This scenario further requires mitigation to be sustained and it leads to a fragile and unstable state when mitigation measures are relaxed.

In the case of heterogeneous populations, intermediate scenarios are possible which stay below herd immunity or just reach herd immunity, see Fig 6(d)–6(i). In the latter example shown in Fig 6(g)–6(i), mitigation effectively reduces the total number of infections by keeping immunity just at herd immunity level, leading to a stable state where mitigation can be safely relaxed. This is a desirable outcome because the number of infections could be reduced by mitigation by up to 40% without the need of sustaining mitigation, see Fig 2(c).

When discussing rapidly evolving epidemics such as SARS-CoV-2, herd immunity is often not considered to be reachable as it is predicted to require an unacceptably high fraction of cumulative infections [38]. Interestingly, the picture changes dramatically if a strongly heterogeneous population is considered. In this case herd immunity can be reached rather quickly while a large majority of the population is still susceptible. This raises the question of what are the features that are variable and that give rise to heterogeneity and how widely they are expected to vary in the population. One possibility is that differences of susceptibility stem from differences in the abilities of immune systems of susceptible individuals to react to a new pathogen. In addition to adaptive immunity related to the presence of specific antibodies, many individuals show a T-cell response to SARS-CoV-2 [25]. This response could for example be due to less specific cross reactions related to earlier encounters with related viruses [24].

Here we have focused on data from Germany until July 2020 because it provides detailed information that is not available in most countries. Furthermore, Germany has relatively few reported infections and deaths per capita. Our work shows that even this rather mild manifestation of the epidemic can be captured by a heterogeneous SIR model with mild mitigation. The data of newly infected cases with fatal outcome can even be explained in a strongly heterogeneous population without considering any mitigation effects at all. This implies that in order to quantify the effects of mitigation, population heterogeneity has to be taken into account. In order to disentangle effects from heterogeneity and from mitigation the combination of different types of information is important. For example, analyzing in different countries the circumstances under which different sero-prevalence levels or multiple epidemic waves are observed could be key to understand the roles of mitigation and heterogeneity.

Appendix A: Properties of the homogeneous SIR model

For a homogeneous population with , the SIR model given in Eqs (1 and 2) can be written as with S = (N − I0)exp(−τ), I(0) = I0 and τ(0) = 0. We therefore have , where R0 = β/γ. For constant β this implies

We can eliminate τ and find

The maximum Imax = I(τ) with I′(τ) = 0, where the prime denotes a τ-derivative, occurs for

Therefore we the maximal number of infected individuals reads

At the maximum Imax, we have dI/dS = 0, which implies

At long times, the infection dies out when I(τ∞) = 0 with (1 − I0/N)exp(−τ∞) = 1 − τ∞/R0 and S∞ = (N − I0)exp(−τ∞). We therefore have S∞/N = 1 + ln(S∞/(N − I0))/R0 and where W(z) denotes the 0-branch of Lambert W-function. The time dependent solution τ(t) can be obtained from Ndτ/I(τ) = βdt via

To discuss empirical data, we consider the time-course of the rate of new cases J = βIS/N. We have with

The maximum Jmax = J(τ) is reached for K(τ) = 0, which implies where W−1(z) denotes the −1 branch of the Lambert W-function with W(z)e = z [28]. At the maximum Jmax of J we have and therefore and finally

Near the maximum of the rate jmax with we have and . For the homogeneous SIR model, we have

Appendix B: Distributions of infection susceptibility in the population

For an initial distribution s0(x) of susceptible individuals with susceptibility x, we define . We can then write the dynamics of the epidemic spreading given in Eqs (1) and (2) as two equations for I(t) and τ(t) with initial values I(0) = I0, τ(0) = 0 and S(0) = N − I0. The number of susceptible at time t is then given by S(t) = S(τ(t)). Defining the cumulant-generating function Γ(τ) = −ln S(τ), we have and the nth cumulant of x is for n > 1 given by

The classic case with homogeneous population then corresponds to s0(x) = (N − I0)δ(x), see Appendix A.

Here we consider distributions which exhibit a power-law behavior for small x with s0 ∼ x−1+. For α > 0 the power law must be cut off at large x > x0 for the distribution to be normalizable. From ∂ s = −xs, we have s(x, τ) = s0(x)e− which for τ ≫ 1/x0 approaches s(x, τ)∼x−1+ e−. The moments of this distribution can be obtained from the cumulant generating function The cumulants of this distribution are 〈x〉 = ατ−(n − 1)!. Using , the susceptibility thus approaches for large τ the limiting distribution which is the gamma distribution.

Appendix C: The generalized SIR model with population heterogeneity

We have shown in Appendix B that for susceptibility distribution with a power law at small x the gamma distribution is an attractor of the dynamics. We therefore choose at time t = 0 a gamma distribution with average as initial condition. It is given by

Here (α − 1)! denotes Euler’s gamma function. Note that in the limit of large α, this approaches the homogeneous case with s0(x)≃δ(x − 1). We then have S(τ) = (N − I0)(1 + τ/α)− and S′ = −(N − I0)(1 + τ/α)−(1+, where the prime denotes a derivative with respect to τ. As time evolves, the shape of the distribution s(x, t) is time independent. Indeed, with ∂ s = −xs we have s(x, τ) = s0(x)e− and thus with . The time-invariant distribution is then given by

The dynamic equation of the heterogeneous SIR model read

Defining , we have

For constant β, we then have by integrating over τ

The maximum of I is reached for τ = τ with I′ = 0 and thus

We thus obtain

The herd immunity level C/N is obtained by using C = N − S, where S is derived using Eq (13) together with Eq (C8). It follows that

The epidemic ends at long times for I(τ∞) = 0, for which with S∞ individuals that remain susceptible. This quantity obeys

We then find Where the function F(z, ν) is defined as the inverse of the function x(1 − x) via the condition F(1 − F) = z. Finally, using Ndτ/I(τ) = βdt, the time dependent solution τ(t) can be written as

Appendix D: Dynamics of the rate of daily new cases

Data on the dynamics of the epidemic typically provides information about new cases reported per day. We therefore consider the rate of nex cases J = −βIS′/N, where the prime denotes a τ derivative. Using I′ = −S′ − N/R0, we have

We then write with

We then have

The maximum of J occurs at τ = τ with K(τ) = 0. We thus have and .

Plugging in S(τ) = N(1 + τ/α) for the heterogeneous SIR model with I0 ≪N, yields

At the maximum in J, we have K = 0, yielding

Using the function F(z, ν) defined by F(1 − F) = z, we can solve this for τ:

This allows us to compute A2 and A3

The parameters λ0, λ∞, A2 and A3 can be obtained from linear and cubic fits to the logarithm the number of daily reported cases Jrep. For these fits, time intervals corresponding to initial exponential growth (T), peak T and final decay T need to be defined. We use the time point , where Jrep reaches its maximum as a reference point relative to which the intervals are given by:

These time intervals are further reduced depending on the used data and Δt such that all time points before the last day with Jrep = 0 prior to and after the first day Jrep = 0 after are excluded. The fits in Fig 3 and the dashed horizontal lines in Figs 5 and 7 correspond to fits with Δt = 19days. The shaded areas in Figs 5 and 7 depict the range of parameter values one obtains for fits with 10 days ≤ Δt ≤ 20days.

Fig 7

Normalized coefficient at the peak of new cases per day.

(a) The ratio as a function of R0 is shown for different values of α. (b) The ratio as a function of R0 for the same values of α. The dashed lines represent the values inferred from the data shown in Fig 3 for all cases (red) and fatal cases (blue). The shaded colored regions correspond to the uncertainties of the fits to the data.

Appendix E: Small α limit for heterogeneous populations

For small α the system reaches a well defined limiting dynamics that can be expressed analytically. We start from I(τ) for small I0/N which for small α becomes

The maximum of I occurs at τ = τ when I′ = 0 or τ/α ≃ R0 − 1. We thus have

Similarly, using C/N = 1 − (1 + τ/α)−, we find for small α

At long times, we have I(τ∞) = 0, where (1 + τ∞/α)− = 1 − τ∞/R0. For small α this implies α ln(1 + τ∞/α)≃τ∞/R0 and thus and , where

In the limit of small α, u = τ/α is finite. The limiting function u(t) for small α can be expressed as where i0 = I/(αN) in the limit α = 0. The number of susceptible then becomes

Finally we discuss the maximum of the rate of new cases J = Jmax. We have , where

At the maximum of J, τ = τ with

Defining we have in the limit of small α

The value of J at the maximum is

We determine and , with and . We then find

We also have and

Appendix F: Mitigation in the heterogeneous SIR model

We now consider the case where the rate of infections β(t) becomes time dependent because of overall changes of conditions such as seasonal effects or measures of social distancing. Using I′ = −S′+ N/R0, we have and the reproduction number where β0 = β(t = 0) and R0 = β0/γ. The epidemic can be mitigated by a reduction of β over time. However if the mitigation is relaxed the epidemic can grow again. As the epidemic advances, τ increases as . Growth of infection number is no longer possible for τ > τ with

Thus the condition τ > τ defines herd immunity conditions where the epidemic can no longer grow. If mitigation sets in early, before τ = τ, the epidemic is slowed and it takes more time to reach herd immunity. in this case a new wave starts after mitigation is relaxed. If mitigation occurs for τ > τ, mitigation facilitates the decay of infections by reducing λ∞ = −(β∞/β0)R0 S′(τ∞) − γ as compared to the value λ∞ = −R0 S′(τ∞) − γ without mitigation.

Appendix G: Inferring β(t) from reported cases

For a given time course of infections, there always exists a function β(t) such that the SIR model follows this time course. We first consider the classic SIR model. A change in the rate of new infections J = βIS/N can be decomposed in three different contributions,

In the case of an early mitigation, S ≈ N and thus . Together with Eq (2), we find

This provides a differential equation for ln β if ln J(t) is given, which does not require knowledge of the amplitude of J. We infer β(t) for each day, using the initial value β(0) = 0.48days−1 at March 15. We use an iterative scheme to calculate the rate for the next day as where is a running average over seven days of the number of reported cases.

For the two scenarios of a later mitigation, the heterogeneous SIR model was considered with J = βI(1 − I0/N)(1 + τ/α)−(. We then have

Again, this equation can be used to construct an iterative scheme to infer β(t). For given initial number of infected individuals on March 15 I(0), we can iteratively obtain the subsequent values as

The starting value of τ(0), can be derived by inverting Eq (C7) for I(τ(0)) = I(0).

Appendix H: Mobility data

Data concerning the changes in mobility of the population has been provided by Google [34]. The data reports the changes compared to a baseline of visits and length of stay at different places. The baseline depends on the specific day of the week and refers to the median value, for the corresponding day of the week, during the 5-week period Jan 3–Feb 6, 2020. Fig 8 shows these changes for Germany for a representative number of categories. These categories are defined in [34] as follows: “Grocery and pharmacy: Mobility trends for places like grocery markets, food warehouses, farmers markets, specialty food shops, drug stores, and pharmacies. Transit stations: Mobility trends for places like public transport hubs such as subway, bus, and train stations. Retail and recreation: Mobility trends for places like restaurants, cafes, shopping centers, theme parks, museums, libraries, and movie theaters. Residential: Mobility trends for places of residence. Workplaces: Mobility trends for places of work. The residential category shows a change in duration while the other categories measure a change in total visitors.”

Fig 8

Mobility changes for a representative set of commonly visited places in Germany up to July 1 2020 from [34].

3 Sep 2020

PONE-D-20-24977

Power-Law Population Heterogeneity Governs Epidemic Waves

PLOS ONE

Dear Dr. Jülicher,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Comments of reviewers are valuable and give you an opportunity to improve the paper.

Please submit your revised manuscript by Oct 18 2020 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

We look forward to receiving your revised manuscript.

Kind regards,

Vygintas Gontis, Ph.D.

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Thank you for inlcuding your funding statement;

"The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript."

At this time, please address the following queries:

Please clarify the sources of funding (financial or material support) for your study. List the grants or organizations that supported your study, including funding received from your institution.

State what role the funders took in the study. If the funders had no role in your study, please state: “The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.”

If any authors received a salary from any of your funders, please state which authors and which funders.

If you did not receive any funding for this study, please state: “The authors received no specific funding for this work.”

Please include your amended statements within your cover letter; we will change the online submission form on your behalf.

3. In your Data Availability statement, you have not specified where the minimal data set underlying the results described in your manuscript can be found. PLOS defines a study's minimal data set as the underlying data used to reach the conclusions drawn in the manuscript and any additional data required to replicate the reported study findings in their entirety. All PLOS journals require that the minimal data set be made fully available. For more information about our data policy, please see http://journals.plos.org/plosone/s/data-availability.

Upon re-submitting your revised manuscript, please upload your study’s minimal underlying data set as either Supporting Information files or to a stable, public repository and include the relevant URLs, DOIs, or accession numbers within your revised cover letter. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories. Any potentially identifying patient information must be fully anonymized.

Important: If there are ethical or legal restrictions to sharing your data publicly, please explain these restrictions in detail. Please see our guidelines for more information on what we consider unacceptable restrictions to publicly sharing data: http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions. Note that it is not acceptable for the authors to be the sole named individuals responsible for ensuring data access.

We will update your Data Availability statement to reflect the information you provide in your cover letter.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: In the paper under review, the authors give a detailed study of the classical SIR model for epidemics in the case when the heterogeneity of the susceptible population is allowed. The main focus is made on the case where the initial susceptibility distribution has gamma density with parameter $\\alpha$. This, in particular, leads to the time-dependent reproduction number $R(t)$, decreasing as hyperbolic function, i.e.\\ $R(t)=(\\bar x(t))^{1+\\alpha} R_0$. The authors provide quantitative properties of the model, such as the herd immunity level, the final size of epidemics, etc. An important conclusion is that, in the heterogeneous population, the herd immunity level can be much lower than in homogeneous case (typically 60\\%).

In my opinion, the subject and results of the paper are very interesting, the paper is well written and I recommend it for publishing in the journal after minor revision. Attached, please find a list of comments.

Reviewer #2: In this manuscript, the authors propose a generalized SIR model taking into account the heterogeneity of the population, i.e., a distribution of the susceptibility of being infected, in terms of a parameter alpha, which is the power-law exponent of the susceptibility distribution when small values of alpha are considered. In other words, when alpha -> infinity, the classical SIR model is recovered and for the case of alpha -> 0, the heterogeneity of the population is incorporated into the SIR model. The key-result of their work is that the population herd immunity is earlier achieved in the case of a heterogeneous population, implying in a lower number of infected people and fatalities. The authors employ their model to analyze the case of the Covid-19 spread in Germany and discuss the importance of taking the population's heterogeneity in the effectiveness of mitigation actions, such as social distancing. Below, I raise some points to be addressed:

- In panel (a) and (d) of Fig. 1, it is not quite clear to me why the number of infected people is lower for the heterogeneous SIR model. Also, what is the justification for using the specifically values of R0 = 2.5, gamma = 0.13 day^-1 and alpha = 0.1? I suggest that the authors make these points clear;

- For future works, I believe it would be interesting to explore the very same analysis here employed in terms of the heterogeneity of the population in the light of the SIRS model, since it would be interesting to analyse how the population's susceptibility distribution is affected in the case where Recovered people can become susceptible of being infected again. Perhaps it would interesting to mention this in the main text;

- On page 6, the authors mention about the Lambert W function and points to a more detailed discussion about it in Appendix A. At this point, I believe it is worth adding a couple of references for clarity both in the main text and in Appendix A about the Lambert W function, just for the sake of completeness;

- The authors should correct on page 3, third paragraph, first sentence, the typo "the" twice in the sentence "In the heterogeneous SIR model proposed here, the qualitative behaviors of the the epidemic wave are unchanged.";

- I suggest that the first sentence of page 8 "Eq. (9) can be then be written as" to be corrected to "Eq. (9) can then be written as...";

- On page 8, the authors write "The dynamics of epidemic waves depends on the shape of the initial distribution s0(x). Here, we consider distributions that have the special property of shape invariance under the dynamics of epidemics. This property is satisfied by a gamma distribution". I suggest that the authors state if only the gamma distribution satisfies such a condition and, if it is not the case, then it would be interesting to add a sentence about the consideration of other distributions as well;

- In Fig 4 panel (a), I suggest that the authors state why their solutions of the heterogeneous SIR model do not incorporate the initial increase in the number of infected people and the small increase between Jun and Jul for both the number of infections and fatalities. The same holds for panel (c) in the case of the initial growth of both the cumulative number of infected people and fatalities. In panel (f), it would be interesting to discuss what is the meaning of tau saturating at the specific value of ~0.2, as well as about the meaning of the average susceptibility x saturating close to the \\tau curve? What does this mean? I suggest that the authors briefly discuss about these points;

- Based on their discussions, the achievement of the herd immunity is key regarding the fade of the disease spread. Based on their discussion about and heterogeneous population, do the authors have any suggestions for public policies in order to minimize the number of infections and, consequently, the number of fatalities?;

- On page 18, the authors write "We show that as a result of strong population heterogeneity (small alpha), the wave peaks when only a small minority of individuals have been infected, see Fig. 1 (d)-(f)." This is indeed true. However, upon analysing panels (a), (d), and (g) of Fig. 6, one notices that the model solution fit of the data (red solid line) present lower reported cases decrease rate than in the scenario without mitigation (red dotted line). This is particularly true after June. How can this be explained? It seems that, although the maximum is earlier achieved, the decrease rate is lower than in the case without mitigation. I suggest that the authors include a few sentences to discuss about this;

- Section "Discussion" seems more like "Conclusions and Perspectives";

- I believe that the number of references could be improved in this work since there are a lot of discussions throughout the manuscript that deserves more important references.

In summary, the work is relevant since in reality not everyone is equally susceptible to being infected and thus the authors' consideration of a susceptibility distribution among the population is solid and can indeed improve the understanding of the epidemics dissemination. I do recommend publication after minor revisions.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: nbbhj.pdf

Click here for additional data file.

9 Sep 2020

Response to Referee 1

>In the paper under review, the authors give a detailed study of the classical SIR

>model for epidemics in the case when the heterogeneity of the susceptible

>population is allowed. The main focus is made on the case where the initial

>susceptibility distribution has gamma density with parameter $\\alpha$.

>This, in particular, leads to the time- dependent reproduction number $R(t)$,

>decreasing as hyperbolic function, i.e.\\ $R(t)=(\\bar x(t))^{1+\\alpha} R_0$.

>The authors provide quantitative properties of the model, such as the herd

>immunity level, the final size of epidemics, etc. An important conclusion is that,

>in the heterogeneous population, the herd immunity level can be much lower

>than in homogeneous case (typically 60\\%).

>

>In my opinion, the subject and results of the paper are very interesting, the

>paper is well written and I recommend it for publishing in the journal after

>minor revision. Attached, please find a list of comments.

We thank the referee for a careful reading of our manuscript

>1. In the models for SARS-CoV-2 epidemics, the SEIR model is in widespread

>use, where E stands for the ”exposed” compartment. I wonder, if it is possible

>to relate it with the generalized SIR model in (1)–(2) with time- varying

>parameters?

In our work we focus on the simpler SIR model where the distinction between

exposed and infected is not made. The additional “exposed” state E in the SEIR

model introduces a time delay between exposure to a pathogen and the

onset of infectiousness. This extra delay due to the exposed state could be

captured by an exponential memory kernel in Eq (2) with an extra relaxation time

describing the delay. However, the main properties of the epidemic wave such as

herd immunity levels are only weakly affected by the extra short delay in the SEIR

model.

In our work we focus on the SIR model because it captures essential features of

an epidemics in a minimal model. In our revised manuscript we now relate to the

SEIR model in the discussion.

>2. The model in eqs. (1)–(2) is not very clearly explained. The classical SIR

>model assumes three components – susceptible (S), infected (I) and recov-

>ered (R). Whereas eqs. (1)–(2) of the paper does not include part R. Also,

>it is slightly unclear the meaning of x ¯. It would be useful to explain what does

>it mean ”dimensionless average susceptibility”. Usually the coefficient at IS/N

>is called ”infection rate” and can admit values larger (or smaller) than 1.

We use the symbol R for the time dependent reproduction number rather than

the number of recovered individuals. Note that the number of recovered

individuals is given by #recovered = N-S-I. We now add this information. In the

revised manuscript we have rewritten the explanation of the model in order to

be more clear about beta and \\bar x and we now call beta the infection rate.

>3. p. 5. It is mentioned that time-varying β could correspond to seasonal

>changes or mitigation measure. For illustrative purposes, it would be nice to

>see concrete forms of β(t) in such cases.

While in the first parts of the manuscript we consider constant beta, we

discuss in section III.C different mitigation scenarios where we choose specific

forms of beta(t) in order to discuss effects of mitigation during this year’s

SARS-Cov2 epidemic, see Fig. 6.

>4. It is well-known that SIR-type compartment models are rather sensitive to

>initial conditions. It would would be good to see some discussion on this

>sensitivity. At least, how the graphs will change if you take other values than

>I0 = 10.

The initial condition I0 is not important for our work and a change in I0 would

essentially only shift the absolute time of the initial point if the time of the

wave peak is known. Once the wave is progressing and reaches its peak, its

behaviour depends only very weakly on the initial conditions if the initial number

of infected is small compared to the population size. Note that in the appendices

we do provide exact expressions for many quantities as a function of I0/N,

revealing the role of initial conditions. However for all plots and results shown

in the manuscript the initial number I0 is completely unimportant. Because

the plots for other choices of I0 would not be distinguishable from the ones

shown we prefer not to show plots with other values of I0.

Another reason not to discuss different I0 is that at the early time of the

epidemics when I(t) went from say 10 to 20 we have absolutely no information

about the epidemic wave and therefore it does not to help discussing wether

it started with I0=2 very early or with I0=15 several days later.

>5. p. 7, l. -9. Since s depends on both x and t, I would recommend to write

>S(t) = ∫ \\int_0^∞ dxs(x, t). Similarly, in eq. (10), as the right-hand side

>depends on t, better to write x ¯(t) = (1/S(t))\\int_0^∞ dx xs(x, t).

We have changed these expression as the referee suggests.

>6. p. 7, l. -2. The meaning of the variable τ (”a measure for how far the

>epidemic has advanced”) should be explain better.

The variable tau emerges from a trick to solve the equations. It does not have

a physical meaning but it has some clear and important properties.

In the revised manuscript we now explain better the variable tau.

>7. p, 8, eq. (12). Strictly speaking, eq. (12) gives a density function, not

>distibution. Also, you may wish to add that α > 0 and to write a precise form

>s_0(x)=…(The symbol ’∼’ usually means asymptotic equivalence.) Also, I

>wonder if it would be possible to introduce additional flexibility in the initial

>susceptibility by introducing two-parameter gamma distribution with

>density s_0(x)=…..

We agree with the referee and have added that alpha >0. Note that we

thought about all these points carefully and we have good reasons to present

the work in the way we did. We have used the term distribution in the context

of our aim to make the paper more easily accessible for a broader readership

including non-experts for whom the term probability distribution is usually

better known. For this reason we have also kept the equations in the main

part of the manuscript rather light, but we provide all the precise and full

expressions in the Appendices for expert readers. The appendices carry a

lot of substance and should not be seen as secondary.

Note that the precise form of the initial density function s_0 is provided in

Appendix C in Eq. (C1). In our revised manuscript we now refer to the

appendix to clarify what the normalization factor is.

We have checked that the full 2-parameter gamma distribution does not

provide any further flexibility. The reason is that we choose without of loss

of generality \\bar x=1 at the initial time point. With the second parameter

of the gamma distribution we can change this initial value but this change

can be absorbed by changing the infection rate beta. So in fact one can

see beta as the second parameter. But since beta already exists in the

classical SIR model we keep it and only add one new parameter alpha

describing the ratio of variance and mean.

>Minor comments:

>p. 3, l. -11. Delete ’the’.

Done

>p. 8, eq. (13). Please add more details how this equality is obtained.

We added some information and now refer to Appendix C for details.

p. 9. Please add more details how eqs. (18) and (19) are obtained.

We added a reference to the Appendix C where the details are provided.

>p. 22. It seems that the Lambert W -function should be defined by

>W (z)eW (z) = z.

Done

Response to Referee 2

>In this manuscript, the authors propose a generalized SIR model taking

>into account the heterogeneity of the population, i.e., a distribution of the

>susceptibility of being infected, in terms of a parameter alpha, which is

>the power-law exponent of the susceptibility distribution when small values

>of alpha are considered. In other words, when alpha -> infinity, the

>classical SIR model is recovered and for the case of alpha -> 0, the

>heterogeneity of the population is incorporated into the SIR model. The

>key-result of their work is that the population herd immunity is earlier

>achieved in the case of a heterogeneous population, implying in a lower

>number of infected people and fatalities. The authors employ their model

>to analyze the case of the Covid-19 spread in Germany and discuss the

>importance of taking the population's heterogeneity in the effectiveness

>of mitigation actions, such as social distancing.

>Below, I raise some points to be addressed:

>- In panel (a) and (d) of Fig. 1, it is not quite clear to me why the number

>of infected people is lower for the heterogeneous SIR model. Also, what

>is the justification for using the specifically values of R0 = 2.5,

>gamma = 0.13 day^-1 and alpha = 0.1? I suggest that the authors

>make these points clear;

The figure just shows the fact that the number of infected people is lower

in the heterogeneous SIR model to clearly show this point. The reasons are

the lowered herd immunity levels given by Eq (19) resulting from the drop

in \\bar x as shown in Fig. 1 f. In appendix C we provide an exact analysis

of the nonlinear dynamics that reveals these surprising properties.

In order to explain this better, we now clarify after Eq (22) that the drop

of \\bar x is the reason for a reduced herd immunity level and resulting

lower infection numbers.

The parameter values are here just for illustrative purposes. The qualitative

behaviors do not depend on the parameter choice within broad ranges.

However the values chosen are rounded versions of values we found are

relevant to the current SARS-Cos2 epidemics are typical values used in the

current epidemics. gamma = 0.13 corresponds to individuals being

infectious during one week, and this parameter does not affect the shape

but only the duration of the wave. alpha=0.1 was chosen so that the

difference between panels (a) and (d) is clearly visible but such that the

maximum of I(t) can still be seen in (d).

>- For future works, I believe it would be interesting to explore the very

>same analysis here employed in terms of the heterogeneity of the population

>in the light of the SIRS model, since it would be interesting to analyse how

>the population's susceptibility distribution is affected in the case where

>Recovered people can become susceptible of being infected again.

>Perhaps it would interesting to mention this in the main text;

There are many open and interesting questions that one can address

with our approach. We agree that it will be very interesting to look into

effects where recovered individuals becomes infected again. In the

revised manuscript we now mention the SIRS model in the discussion.

>- On page 6, the authors mention about the Lambert W function and

>points to a more detailed discussion about it in Appendix A. At this point,

>I believe it is worth adding a couple of references for clarity both in the

>main text and in Appendix A about the Lambert W function, just for the

>sake of completeness;

We have added a reference on the Lambert W function.

>- The authors should correct on page 3, third paragraph, first sentence,

>the typo "the" twice in the sentence "In the heterogeneous SIR model

>proposed here, the qualitative behaviours of the the epidemic wave are

>unchanged.";

Thanks - corrected

>- I suggest that the first sentence of page 8 "Eq. (9) can be then be

>written as" to be corrected to "Eq. (9) can then be written as…";

Thanks - done

>- On page 8, the authors write "The dynamics of epidemic waves depends

>on the shape of the initial distribution s0(x). Here, we consider distributions

>that have the special property of shape invariance under the dynamics of

>epidemics. This property is satisfied by a gamma distribution". I suggest that

>the authors state if only the gamma distribution satisfies such a condition and,

>if it is not the case, then it would be interesting to add a sentence about the

>consideration of other distributions as well;

To our knowledge the gamma distribution is the only distribution that is shape

invariant under the dynamics. However to be cautious we avoid claiming this

fact as we do not have a proof. Our work and our results do not depend on

the gamma distribution being the only shape invariant distribution.

>- In Fig 4 panel (a), I suggest that the authors state why their solutions of

>the heterogeneous SIR model do not incorporate the initial increase in the

>number of infected people and the small increase between Jun and Jul for

>both the number of infections and fatalities. The same holds for panel (c) in

>the case of the initial growth of both the cumulative number of infected

>people and fatalities.

In Fig. 4 we simply compare fits of our model to data. The fits show that

the data based on reported cases that lead to deaths (blue), which probably

measures more reliably the progression of serious illnesses, is better

captured by the model than simply using the reported number of reported

cases (red). A possibility for the difference to the data early and in June is

that the blue data gives a better representation of the epidemics than the

red data. However we have refrained form making such statements because

all available data have problems and there are many reasons why there

could be serious biases and errors in the data.

In our manuscript we note on p. 13 that it is surprising that the model fits

both types of data rather well even though we only have three time

independent parameters (note the classical SIR model cannot account for

this data).

>In panel (f), it would be interesting to discuss what is the meaning of tau

>saturating at the specific value of ~0.2, as well as about the meaning of

>the average susceptibility x saturating close to the \\tau curve? What does

>this mean? I suggest that the authors briefly discuss about these points;

The final value of tau is explained in Appendix C where we show that tau

reaches a fixed final value when the epidemics dies out that is described

exactly by the implicit equation (C11) in the revised manuscript. In our

revised manuscript we now mention in the main text on p. 7 that tau

reaches a final value and we add a reference to appendix C when discussing

panel (f) of Fig. 4. Note that the relationship between tau and \\bar x is

given in Eq. (14). The similarity of \\bar x and \\tau at the end of the wave

is a coincidence and consistent with Eq. (14).

>- Based on their discussions, the achievement of the herd immunity is key

>regarding the fade of the disease spread. Based on their discussion about

>and heterogeneous population, do the authors have any suggestions for

>public policies in order to minimize the number of infections and, consequently,

>the number of fatalities?;

We have on purpose refrained from a discussion of public policy. Here we

want to focus on the concepts and the science. The science presented here

is clear and rigorous. Implications for policy are much less rigorous and depend

on many other factors and can be coloured by opinions. We think that our

work has many implications for public policy and that it will stimulate and be

useful for future discussions but we do not want to weaken our work by

adding elements that are uncertain.

>- On page 18, the authors write "We show that as a result of strong

>population heterogeneity (small alpha), the wave peaks when only a small

>minority of individuals have been infected, see Fig. 1 (d)-(f)." This is indeed

>true. However, upon analysing panels (a), (d), and (g) of Fig. 6, one notices

>that the model solution fit of the data (red solid line) present lower reported

>cases decrease rate than in the scenario without mitigation (red dotted line).

>This is particularly true after June. How can this be explained? It seems that,

>although the maximum is earlier achieved, the decrease rate is lower than

>in the case without mitigation. I suggest that the authors include a few

>sentences to discuss about this;

It is correct that in the case of the epidemics that is dying out because of

mitigation (Fig. 6), we find that the rate of decay of cases is slower than without

mitigation. This is similar to the idea to “flatten the curve”, i.e. mitigation

reduces the maximal number of infections but broadens the wave and thus

makes it slower. However this feature is not completely general and we prefer

not to enter this discussion.

>- Section "Discussion" seems more like "Conclusions and Perspectives";

We have changed the discussion to “Conclusions and Perspectives”

>- I believe that the number of references could be improved in this work

>since there are a lot of discussions throughout the manuscript that deserves

>more important references.

We have now added reference [28] on Lamberts W function and Ref. [35-37]

for the SEIR and the SIRS model. We think that all statements that need backing

by references have been referenced.

>In summary, the work is relevant since in reality not everyone is equally

>susceptible to being infected and thus the authors' consideration of a susceptibility

>distribution among the population is solid and can indeed improve the

>understanding of the epidemics dissemination. I do recommend publication after

>minor revisions.

We thank the referee for a careful reading of our manuscript

Submitted filename: responsetoreferees.pdf

Click here for additional data file.

11 Sep 2020

Power-Law Population Heterogeneity Governs Epidemic Waves

PONE-D-20-24977R1

Dear Dr. Jülicher,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Vygintas Gontis, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Congratulations, After careful analysis of the article revised version and Response to Reviewers, we have decided to accept the paper for the publication in Plos One.

Reviewers' comments:

17 Sep 2020

PONE-D-20-24977R1

Power-Law Population Heterogeneity Governs Epidemic Waves

Dear Dr. Jülicher:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Vygintas Gontis

Academic Editor

PLOS ONE