as the original 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. We find that in strongly heterogeneous populations the epidemic reaches only a small fraction of the population. This implies that the herd immunity level can be much lower than in commonly used models with homogeneous populations. 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.
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 R 0 , 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] [11] [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] or infectiousness [15, 16] . It was suggested that population heterogeneity reduces effective herd immunity levels [13, 15, 17, 18] .
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 powerlaw heterogeneity in the distribution of infection susceptibilities of the population. Power laws are often found in nonlinear and complex systems [19] [20] [21] [22] . 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 [23, 24] . Furthermore, population heterogeneity could be relevant at very different scales, from the 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 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 [25] . We estimate parameter values including the susceptibility exponent α and 
where the dots denote time derivatives,x is a dimensionless average susceptibility and the rate β(t) describes a probability per unit time and per person to become infected, which can in general depend on time t. This time dependence could correspond to seasonal changes or mitigation measures [10, 12, 26] . 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İ
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 J =Ċ = −Ṡ and R = J/(γI).
In the simple case of a homogeneous population, all individuals have the same degree of susceptibility, x = 1 and the population average of x isx = 1 independent of time. This is the classic SIR model. An example for a solution to these equations for homogeneous populationx = 1 and constant β is given in Fig. 1 (a) ,(b). The corresponding time dependent reproduction number is presented in Fig. 1 (c) . The number of infections first grows exponentially with growth rate
As the number of susceptible decreases, the epidemic reaches a peak number of infected 
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. Here W (z) denotes Lambert W-function, see Appendix A. The total fraction of infections over the course of the
For a classic SIR model with homogeneous population we have for R 0 = 2.5, a herd immunity level C I /N of 60% of the population, see Fig. 1 The fraction of the population that become infected increase for larger R 0 . The SIR model thus suggests that for R 0 > 2 the epidemic wave exceeds a majority of the population before the epidemic begins to die out. 
which for the whole population implies Eq. (1) with average susceptibilitȳ
which is in general time dependent.
This time dependence can be discussed by introducing the variable τ that is a measure for how far the epidemic has advanced. It increases monotonically with time asτ = βI/N . Eq. (9) can be then be written as ∂ τ s = −xs, and the number of susceptible individuals is
where s 0 (x) is the initial susceptibility distribution at time t = t 0 with averagex = 1, see
Appendix B.
C. Infection dynamics with generic power law heterogeneity
The dynamics of epidemic waves depends on the shape of the initial distribution s 0 (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 α, and a cut off at large x. The distribution s 0 (x) has averagex = 1 and variance 1/α. Indeed, we have
where the time dependence enters viax(t), 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 s 0 (x) ∼ x −1+α at small x, it 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.
For the choice (12) we have
The average susceptibility isx
which starts fromx = 1 for τ = 0 and decreases for increasing τ , thus dampening the epidemic. We can now express the dynamics given in Eqns (1) and (2) as two dynamic equations for I(t) and τ (t) which reaḋ 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 = γ(R 0 − 1),
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 maxi-
the reproduction number R drops below 1 and the epidemic dies down. In Eqns. (17)- (19) we have considered the limit of small I 0 /N for simplicity. In the limit of large α, these expressions converge to those obtained for the homogeneous SIR model, see Appendix A.
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 (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 I max /N α(ln R 0 + 1/R 0 − 1) and C I /N α ln R 0 . 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 isx ∞
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 R 0 = 2.5 and α = 0.1, we have I max /N 2.8%, and a fraction C I /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.
We analyze the dynamics of the SARS-CoV-2 epidemic in Germany using public data provided by the Robert Koch institute [25] . 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 J t rep (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 J f rep (t). 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 J f rep (t) is related to the number of new infections per day as J f rep (t) = Jf , where f denotes the infection fatality rate, the fraction of infections that are fatal, which we consider to be constant for simplicity. 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 R 0 and γ to calculate λ 0 0.24 day −1 , slightly smaller than the estimate given in In order to understand how the shape of the wave of infections constrains the possible parameter values of R 0 , γ 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 J max = J(t J ) at time t J .
The coefficients A 2 =J/J| t=t J and A 3 = ... shaded regions corresponding to estimated uncertainty ranges of these values.
We find that the ratio A 2 /λ 2 ∞ , which is independent of γ depends only weakly on α. We can therefore use it to estimate R 0 , see . 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 R 0 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 [27, 29, 30] 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 R 0 is larger and the consistency with the data is less strong.
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, 26] . 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 I of the population via Eq. (19) . The herd immunity level depends on the basic reproduction number R 0 = β 0 /γ and on the population heterogeneity α. For immunity levels above herd immunity, C ≥ C I , 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 In the case of early mitigation, Fig. 6 (a-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. 6 (d) and
(g). A case of moderate mitigation with α = 0.1 is shown in Fig. 6 (d-f ). 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. 6 (g-i) 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 I , 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 I , which is up to 50% smaller than the long lime limit C ∞ in an unmitigated epidemic (see Fig. 2c ). 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 susceptibilityx 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 α.
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, 16] . that only a minority of the population exhibits antibodies [27, 29, 30] . 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)-(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)-(i). In the latter example shown in Fig. 6 (g)-(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 [32] . 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 [24] . This response could for example be due to less specific cross reactions related to earlier encounters with related viruses [23] .
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.
For a homogeneous population withx = 1, the SIR model given in Eqns. (1-2) can be written asİ 
We can eliminate τ and find
The maximum I max = I(τ I ) with I (τ I ) = 0, where the prime denotes a τ -derivative, occurs for
Therefore we the maximal number of infected individuals reads
At the maximum I max , we have dI/dS = 0, which implies
At long times, the infection dies out when I(τ ∞ ) = 0 with
where W (z) denotes the 0-branch of Lambert W-function. The time dependent solution τ (t) can be obtained from N dτ /I(τ ) = βdt via
To discuss empirical data, we consider the time-course of the rate of new cases J = βIS/N . We haveJ/J = βK with
The maximum J max = J(τ J ) is reached for K(τ J ) = 0 which implies
where W −1 (z) denotes the −1 branch of the Lambert W -function with W (z)e W (z) = z. At the maximum J max of J we haveṠI +İS = 0 and therefore
and therefore
Near the maximum of the rate J max withJ = 0 we have A 2 =J/J| t=t J and A 3 = ...
For the homogeneous SIR model, we have
For an initial distribution s 0 (x) of susceptible individuals with susceptibility x, we define
We can then write the dynamics of the epidemic spreading given in Eqns (1) and (2) as two equations for I(t) and τ (t) 
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 averagex = 1 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 s 0 (x) δ(x − 1). We then have S(τ ) = (N − I 0 )(1 + τ /α) −α and S = −(N − I 0 )(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, τ ) = s 0 (x)e −τ x and thus
withx(τ ) = (1 + τ /α) −1 . The time-invariant distribution is then given by
The dynamic equation of the heterogeneous SIR model reaḋ
Defining I =İ/τ , we have
For constant β, we then have by integrating over τ
The maximum of I is reached for τ = τ I with I = 0 and thus
We thus obtain
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 
We then writeJ
We then have
The maximum of J occurs at τ = τ J with K(τ J ) = 0. We thus have A 2 = (J/J)| t=t j = βK
Plugging in S(τ ) = N (1 + τ /α) α for the heterogeneous SIR model with I 0 << N , yields
At the maximum in J, we have K = 0, yielding 2α + 1
Using the function F (z, ν) defined by F ν (1 − F ) = z, we can solve this for τ J : This allows us to compute A 2 and A 3
The parameters λ 0 , λ ∞ , A 2 and A 3 can be obtained from linear and cubic fits to the logarithm the number of daily reported cases J rep . For these fits, time intervals corresponding to initial exponential growth (T i ), peak T p and final decay T f need to be defined. We use the time point t rep m , where J rep 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 J rep = 0 prior to t rep m and after the first day J rep = 0 after t rep m are excluded. The fits in Fig. 3 and the dashed horizontal lines in Fig. 5 and 7 correspond to fits with ∆t = 19days. The shaded areas in Fig. 5 and 7 depict the range of parameter values one obtains for fits with 10days ≤ ∆t ≤ 20days. 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 I 0 /N
The maximum of I occurs at τ = τ I when I = 0 or τ I /α R 0 − 1. We thus have
Similarly, using C I /N = 1 − (1 + τ I /α) −α , we find for small α
At long times, we have
In the limit of small α, u = τ /α is finite. The limiting function u(t) for small α can be expressed as 
where i 0 = I/(αN ) in the limit α = 0. The number of susceptible then becomes S N 1 − α ln(1 + u) .
Finally we discuss the maximum of the rate of new cases J = J max . We haveJ/J = βK,
At the maximum of J, τ = τ J with
Definingx J = (1 + τ J /α) −1 we have in the limit of small ᾱ
The value of J at the maximum is
We determine A 2 =J/J = β 2K and A 3 = ... J /J = β 3K , withK/β = K I/N andK/β 2 = K I 2 /N 2 . We then find
We also have A 2 2 A 3 = γR 0xJ 2 (E14) 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/R 0 , we have λ =İ/I = −βS /N − γ and the reproduction number
where β 0 = β(t = 0) and R 0 = β 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τ = βI/N . Growth of infection number is no longer possible for τ > τ I with − S (τ I ) = 1 R 0 (F2)
Thus the condition τ > τ I defines herd immunity conditions where the epidemic can no longer grow. If mitigation sets in early, before τ = τ I , 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 τ > τ I , mitigation facilitates the decay of infections by reducing 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.48 days −1 at March 15. We use an iterative scheme to calculate the rate for the next day as ln β(i + 1) = ln β(i) + ln J obs (i + 1) − ln J obs (i) − β(i) + γ,
where ln J obs (i) = (1/7) 3 ∆t=−3 J rep (i + ∆t) 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 
ln I(i + 1) = ln I(i)+ β(i)(1 −
ln β(i + 1) = ln β(i)+ ln J obs (i + 1) J obs (i) − ln I(i + 1) I(i)
The starting value of τ (0), can be derived by inverting Eq. (C7) for I(τ (0)) = I(0).
Data concerning the changes in mobility of the population has been provided by Google [31] . 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 3Feb 6,