Time-dependent heterogeneity leads to transient suppression of the COVID-19 epidemic, not herd immunity

Significance Epidemics generally spread through a succession of waves that reflect factors on multiple timescales. Here, we develop a general approach bridging across these timescales and demonstrate how to incorporate population heterogeneity into a wide class of epidemiological models. We demonstrate that a fragile state of transient collective immunity emerges during early, high-paced stages of the epidemic, leading to suppression of individual epidemic waves. However, this state is not an indication of lasting herd immunity: Subsequent waves may emerge due to stochastic changes in individual social activity. Parameters of transient collective immunity are estimated using empirical data from the COVID-19 epidemic in several US locations.

T he COVID-19 pandemic is nearly unprecedented in the level of disruption it has caused globally, but also, potentially, in the degree to which it will change our understanding of epidemic dynamics and the efficacy of various mitigation strategies. Ever since the pioneering works of Kermack and McKendrick (1), epidemiological models have been widely and successfully used to quantify and predict the progression of infectious diseases (2)(3)(4)(5)(6). More recently, the important role in epidemic spreading played by population heterogeneity and the complex structure of social networks has been appreciated and highlighted in multiple studies (7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25). However, integration of this conceptual progress into reliable, predictive epidemiological models remains a formidable task. Among the key effects of heterogeneity and social network structure are 1) the role played by superspreaders and superspreading events during initial outbreaks (12,13,16,(26)(27)(28)(29) and 2) a substantial reduction of the final size of epidemic (FSE) as well as the herd immunity threshold (HIT) compared to the homogeneous case (9, 10, 19-21, 23, 30). The COVID-19 pandemic has reignited interest in the effects of heterogeneity of individual susceptibility to the disease, in particular, to the possibility that it might lower both HIT and FSE (31)(32)(33)(34)(35). In studying epidemics in heterogeneous populations, it is important to emphasize the qualitative nature of two important timescales. First, overdispersion is dominated by short-term patterns of behavior and even accidental events rather than persistent population behavioral heterogeneity. Second, short-term overdispersion is generally assumed to have a limited impact on the long-term epidemic dynamics, being important primarily during early outbreaks dominated by superspreading events. In this paper, we attempt to provide a multiscale theory for epidemic progression and show that both overdispersion and persistent heterogeneity affect the overall progression of the COVID-19 epidemic. The significance of this multiscale perspective is that it provides a natural formalism to predict the occurrence and nature of successive epidemic waves, even when it might seem that a first wave has attained a state which could be mistaken for herd immunity.
There are several existing approaches to model the effects of heterogeneity on epidemic dynamics, each focusing on a different characteristic and parameterization. In the first approach, one stratifies the population into several demographic groups (e.g., by age), and accounts for variation in susceptibility of these groups and their mutual contact probabilities (2). While Significance Epidemics generally spread through a succession of waves that reflect factors on multiple timescales. Here, we develop a general approach bridging across these timescales and demonstrate how to incorporate population heterogeneity into a wide class of epidemiological models. We demonstrate that a fragile state of transient collective immunity emerges during early, high-paced stages of the epidemic, leading to suppression of individual epidemic waves. However, this state is not an indication of lasting herd immunity: Subsequent waves may emerge due to stochastic changes in individual social activity. Parameters of transient collective immunity are estimated using empirical data from the COVID-19 epidemic in several US locations. this approach represents many aspects of population dynamics beyond the homogeneous and well-mixed assumption, it clearly does not encompass the whole complexity of individual heterogeneity, interpersonal communications, and spatial and social structures. These details can be addressed in a second approach, where one analyzes the epidemic dynamics on real-world or artificial social networks (8,13,23,36,37). Through elegant mathematics, it is possible to obtain detailed results in idealized cases, including the mapping onto well-understood models of statistical physics such as percolation (9,38). As demonstrated in ref. 21, the FSE is a very robust property of the epidemic, insensitive to fine details of its dynamics (39) defined by 1) distribution of susceptibilities in the population (20,30,40) and 2) correlations between infectivity and susceptibility. Importantly, it does not depend on the part of heterogeneous infectivity that is not correlated with susceptibility. However, these approaches, so far, have been mostly limited to the analysis of the FSE and distribution of outbreak sizes on a static social network.
For practical purposes, it is desirable to predict the complete time-dependent dynamics of an epidemic, preferably by explicitly including heterogeneity into classical well-mixed mean-field compartmentalized models. This approach was developed some time ago in the context of epidemics on networks (10,23), and the resulting mean-field theory effectively reproduces the structure of heterogeneous well-mixed models extensively studied in the applied mathematics literature (19-22, 24, 30). The overall impact of heterogeneity occurs because, as the disease spreads, it preferentially immunizes the more susceptible individuals, so the remaining population is less susceptible, and spread is inhibited. This effect is further enhanced by a positive correlation between infectivity and susceptibility. In the context of static network models, this correlation is perfect, since both factors are proportional to the degree of individual nodes. Ref. 24 studied a hybrid model in which social heterogeneity represented by network degree was combined with a biological one. These approaches have been recently applied in the context of COVID-19 (31,32,35,41,42). The conclusion of these studies was that the HIT may be well below that expected in classical homogeneous models.
These approaches to heterogeneity delineate end-members of a continuum of theories: overdispersion describing shortterm, bursty dynamics (e.g., due to superspreader accidents), as opposed to persistent heterogeneity, which is a long-term characteristic of an individual and reflects behavioral propensity, for example, to socialize in large gatherings without prudent social distancing. Overdispersion is usually modeled in terms of a negative binomial branching process (12,13,16,(26)(27)(28). Strictly speaking, both persistent heterogeneity and short-term variations contribute to the overdispersion of individual reproduction number. However, we will see below that the former is likely to be a much weaker source of variation compared to the latter. It is also generally presumed that short-term overdispersion is uncorrelated in time and thus has no effect on epidemic dynamics. Indeed, large variations in an individual's infectivity would average out as long as they are not correlated with susceptibility. But, since the initial exposure and secondary infections are separated by a single generation interval (typically about 5 d for , the levels of individual social activity at those times are expected to be correlated, and (at least partially) reflect short-term overdispersion. How, then, can one understand the epidemic progression across various timescales, from early stages of a fast exponential growth to the final state of the herd immunity?
Below, we present a comprehensive yet simple theory that accounts for both social and biological aspects of heterogeneity, and predicts how, together, they modify early and intermediate epidemic dynamics, as well as global characteristics of the epidemic such as its HIT. Importantly, early epidemic dynamics may be sensitive to both persistent heterogeneity and short-term overdispersion. As a result, the apparent early-stage heterogeneity is typically enhanced compared to its long-term persistent level. This may lead to a suppression of the first wave of the epidemic due to reaching a novel state that we call transient collective immunity (TCI) determined by a combination of shortterm and long-term heterogeneity, whose threshold is lower than the eventual HIT. The implication is that the first wave of an epidemic ends due to a combination of both persistent heterogeneity and whatever mitigation constraints are imposed on the population. As the latter are relaxed by authorities or through behavioral changes associated with seasonal factors, subsequent waves can still occur. Thus, TCI is dramatically different from the idea of herd immunity due to heterogeneity.
Our starting point is a generalized version of the heterogeneous well-mixed theory in the spirit of ref. 10, but we use the age-of-infection approach (1) rather than compartmentalized susceptible, infectious, recovered (SIR)/susceptible, exposed, infectious, recovered (SEIR) models of epidemic dynamics (see, e.g., ref. 2). Similar to multiple previous studies, we first completely ignore any time dependence of individual susceptibilities and infectivities, focusing only on their long-term persistent components. This approach implicitly assumes that any short-term overdispersion (responsible for, e.g., the superspreading phenomenon) is uncorrelated in time and thus effectively averaged out. This is a valid assumption if the large short-term variations in individual infectivity are completely uncorrelated with an individual's susceptibility. However, this approximation is hard to justify in the case of COVID-19. Indeed, if the two are correlated on the timescale of a single generation interval (5 d), this will strongly affect the overall epidemic dynamics. Therefore, our initial analysis is eventually expanded to a more general theory accounting for the nonnegligible effects of short-term overdispersion. In the case of persistent heterogeneity, we demonstrate how the model can be recast into an effective homogeneous theory that can readily encompass a wide class of epidemiological models, including various versions of the popular SIR/SEIR approaches. Specific innovations that emerge from our analysis are the nonlinear dependence of the effective reproduction number Re on the overall population fraction S of susceptible individuals, and another nonlinear function Se that gives an effective susceptible fraction, taking into account preferential removal of highly susceptible individuals.
A convenient and practically useful aspect of this approach is that it does not require extensive additional calibration in order to be applied to real data. In the effort to make quantitative predictions from epidemic models, accurate calibration is arguably the most difficult step, but is necessary due to the extreme instability of epidemic dynamics in both growth and decay phases (43,44). We find that, with our approach, the entire effect of heterogeneity is, in many cases, well characterized by a single parameter which we call the immunity factor λ. It is related to the statistical properties of heterogeneous susceptibility across the population and to its correlation with individual infectivity. The immunity factor determines the rate at which Re drops during the early stages of the epidemic as the pool of susceptibles is being depleted: Re ≈ R0(1 − λ(1 − S )). Beyond this early linear regime, for an important case of gamma-distributed individual susceptibilities, we show that the classical proportionality, Re = R0S , transforms into a power law scaling relationship Re = R0S λ . This leads to a modified version of the result for the HIT, Heterogeneity in the susceptibility of individual members of the population has several different contributions: 1) biological, which takes into account differences in factors such as strength of immune response, genetics, age, and comorbidities; and 2) social, reflecting differences in the number and frequency of close contacts of different people. The immunity factor λ in our model combines these sources of heterogeneous susceptibility as well as its correlation with individual infectivity. As we demonstrate, under certain assumptions, the immunity factor is simply a product of social and biological contributions: λ = λs λ b . In our study, we leverage existing studies of real-life face-to-face contact networks (13,19,37,(45)(46)(47)(48) to estimate the social contribution to heterogeneous susceptibility, and the corresponding immunity factor λs . The biological contribution, λ b , is expected to depend on specific details of each infection.
To test this theory, we use empirical data for the COVID-19 epidemic to independently estimate the immunity factor λ. In particular, we apply our previously described epidemic model that features multichannel Bayesian calibration (43) to describe epidemic dynamics in New York City (NYC) and Chicago from the start of the epidemic in mid-March until the end of the first wave around June 15, 2020. This model uses high-quality data on hospitalizations, intensive care unit (ICU) occupancy, and daily deaths to extract the underlying Re (S ) dependence in each of two cities. In addition, we perform a similar analysis of data on individual states in the United States, using data generated by the model in ref. 49. Using both approaches, we find that the locations that were severely impacted by the COVID-19 epidemic show a more pronounced reduction of the effective reproduction number. This effect is much stronger than predicted by classical homogeneous models, suggesting a significant role of heterogeneity. The estimated immunity factor ranges between four and five. Importantly, this represents a transient value of the parameter λ observed on intermediate timescales and dependent on both persistent and short-term heterogeneity. Our estimates of the long-term value of the immunity factor defined by persistent heterogeneity only is considerably lower, about two. This difference explains why achieving the state of TCI after the first wave of the epidemic does not imply long-term herd immunity.
Finally, we integrate the persistent heterogeneity theory into our time-of-infection epidemiological model (43), and project possible outcomes of the second wave of the COVID-19 epidemic during the summer months in NYC and Chicago, using data up to the end of June 2020. By considering the worst-case scenario of a full relaxation of any currently imposed mitigation, we find that the results of the heterogeneity-modified model significantly modify the results from the homogeneous mode. In particular, based on our estimate of the immunity factor, our model predicts no second wave in NYC immediately after release of mitigations in June and up to September 2020, indicating that the TCI has likely been achieved there. Chicago, on the other hand, has not passed the TCI threshold that we infer, but the effects of heterogeneity would still result in a substantial reduction of the magnitude of the second wave there, even under the worst-case scenario. This, in turn, suggests that the second wave can be completely eliminated in such medium-hit locations, if appropriate and economically mild mitigation measures are adopted, including, for example, mask wearing, contact tracing, and targeted limitation of potential superspreading events, through limitations on indoor bars, dining, and other venues. We further investigate the issue of fragility of collective immunity in heterogeneous populations. By allowing rewiring of the social network with time, we demonstrate that the TCI may wane, much like an individual's acquired immunity may wane due to biological factors. This phenomenon would result in the emergence of secondary epidemic waves after a substantial period of low infection counts.

Theory of Epidemics in Populations with Persistent Heterogeneity
Following in the footsteps of refs. 10,19,[22][23][24]30, and 32, we consider the spread of an epidemic in a population of individuals who exhibit significant persistent heterogeneity in their suscep-tibilities to infection α. This heterogeneity may be biological or social in origin, and we assume these factors are independent: α = α b αs . Effects of possible correlations between α b and αs have been discussed in ref. 24. The biologically driven heterogeneous susceptibility α b is shaped by variations of several intrinsic factors such as the strength of individuals' immune responses, age, or genetics. In contrast, the socially driven heterogeneous susceptibility αs is shaped by extrinsic factors, such as differences in individuals' social interaction patterns (their degree in the network of social interactions). Furthermore, individuals' different risk perceptions and attitudes toward social distancing may further amplify variations in socially driven susceptibility heterogeneity. We only focus on susceptibility that is a persistent property of an individual. For example, people who have elevated occupational hazards, such as health care workers, typically have higher, steady values of αs . Similarly, people with low immune response, highly social individuals (hubs in social networks), or scofflaws would all be characterized by above-average overall susceptibility α.
In this work, we group individuals into subpopulations with similar values of α and describe the heterogeneity of the overall population by the probability density function (pdf) of this parameter, f (α). Since α is a relative measure of individual susceptibilities, without loss of generality, we set α ≡ ∞ 0 αf (α)d α = 1. Each person is also assigned an individual reproduction number Ri , which is an expected number of people that this person would infect in a fully susceptible population with α = 1. Accordingly, from each subpopulation with susceptibility α, there is a respective mean reproductive number Rα to which we refer as infectivity throughout this study. Any correlations between individual susceptibility and infectivity will significantly impact the epidemic dynamics. Such correlations are an integral part of most network-based epidemiological models, due to the assumed reciprocity in underlying social interactions, which leads to Rα ≈ α (9, 10, 23). In reality, not all transmissions involve face-to-face contacts, and biological susceptibility need not be strongly correlated with infectivity. Therefore, it is reasonable to expect only a partial correlation between α and Rα.
Let Sα(t) be the fraction of susceptible individuals in the subpopulation with susceptibility α, and let jα(t) = −Ṡα be the corresponding daily incidence rate per capita in that subpopulation. At the start of the epidemic, we assume everyone is susceptible to infection: Sα(0) = 1. The course of the epidemic is described by the following age-of-infection model: Here t is the physical time, and τ is the time since infection for an individual; . . . represents averaging over α with pdf f (α). J (t) is the force of infection, that is, per capita incidence rate in a fully susceptible subpopulation with α = 1. Rα is the previously introduced infectivity, that is, the mean reproductive number of the subpopulation with susceptibility α. K (τ ) is the pdf of the generation interval, which we assume to be independent of α, for the sake of simplicity. The homogeneous version of the age-ofinfection model was introduced in the early days of mathematical epidemiology in the classical paper by Kermack and McKendrick (1). It is based on the observation that the force of infection, J (t), is defined by the number of previously infected individuals. The contribution of each individual depends on his/her time since infection τ and is weighted by the infectivity profile K (τ ). As shown in ref. 50, the rate of the exponential growth of the epidemic can be inferred from the Laplace transform of K (τ ).

Tkachenko et al.
Time-dependent heterogeneity leads to transient suppression of the COVID-19 epidemic, not herd immunity Well-known compartmentalized models correspond to specific functional forms of K (τ ), such as the exponential distribution for the SIR model. According to Eq. 1, the fraction of susceptibles in the subpopulation with any given α can be expressed as Here The total susceptible fraction of the population is related to the moment generating function Mα of the distribution f (α) (i.e., the Laplace transform of f (α)) according to [

4]
Similarly, the effective reproductive number Re can be expressed in terms of the parameter Z , Note that, for Z = 0, this expression gives the basic reproduction number R0 = αRα / α . This result is reminiscent of the well-known result (23) R0 = kin kout / kin for epidemic spread on a directed network with in and out degrees kin (analogue of our susceptibility α) and kout (analogue of Rα). Note that, due to our choice of normalization α = 1, the prefactor in Eq. 5 can be omitted. Since both S and Re depend on time only through Z (t), Eqs. 4 and 5 establish a parametric relationship between these two important quantities during the time course of an epidemic. In contrast to the classical case when these two quantities are simply proportional to each other, that is, Re = SR0, the relationship in the present theory is nonlinear due to heterogeneity. Eq. 5 was derived by substituting Eq. 1 into Eq. 2. This leads to the renewal equation for J (t) of the same form as in a homogeneous problem, Furthermore, by averaging Eq. 1 over all values of α, one arrives at the following heterogeneity-induced modification to the relationship between the force of infection and incidence rate: Here is the effective susceptible fraction of the population, which is less than S , due to the disproportionate removal of highly susceptible individuals. Just as with Re , it is a nonlinear function of S , defined parametrically by Eqs. 4 and 8. Further generalization of this theory for the time-modulated age-of-infection model is presented in SI Appendix. There, we also discuss the adaptation of this approach for the important special case of a compartmentalized SIR/SEIR model. One of the striking consequences of the nonlinearity of Re (S ) is that the effective reproduction number could be decreasing at the early stages of an epidemic significantly faster than predicted by homogeneous models. Specifically, for 1 − S (t) Z (t) 1, one can linearize the effective reproduction number as We named the coefficient λ the immunity factor because it quantifies the effect that a gradual buildup of population immunity has on the spread of an epidemic. The classical value of λ is one, but it may be significantly larger in a heterogeneous case. By linearizing Eq. 5 in terms of 1 − S Z 1 and dividing the result by R0 = αRα , one gets λ = α 2 Rα αRα . [10] As one can see, the value of the immunity factor thus depends both on the statistics of susceptibility α and on its correlation with infectivity Rα. We previously defined the overall susceptibility as a combination of biological and social factors: α = αs α b . Here αs is a measure of the overall social connectivity or activity of an individual, such as the cumulative time of close contact with other individuals averaged over a sufficiently long time interval (known as node strength in network science). Since the contribution of interpersonal contacts to an epidemic spread is mostly reciprocal, we assume Rα ≈ αs . On the other hand, in our analysis, we neglect a correlation between biological susceptibility and infectivity, as well as between α b and αs . Under these approximations, the immunity factor itself is a product of biological and social contributions, λ = λ b λs . Each of them can be expressed in terms of leading moments of α b and αs , respectively, These equations follow from Eq. 10 in the limit Rα = const and Rα ≈ α, respectively. Although these equations resemble classical results for R0 in heterogeneous networks (8-10, 23), here they describe a completely different effect of suppression of Re in response to depletion of susceptible population S . That is why λs in Eq. 12 is proportional to the third moment of αs instead of the second moment in the case of R0 = αRα ≈ α 2 s . Note that the biological contribution to the immunity factor depends only on the coefficient of variation CV b of α b . On the other hand, the social factor λs depends on both the coefficient of variation CVs and the skewness γs of the distribution of αs . Due to our normalization, αs α b ≈ αs α b = α = 1.
The relative importance of biological and social contributions to the overall heterogeneity of α may be characterized by a single parameter χ. For a log-normal distribution of α b , χ appears as a scaling exponent between infectivity and susceptibility: Rα ≈ α χ (see SI Appendix for details). The corresponding expression for the overall immunity factor is λ = α 2+χ / α 1+χ . The limit χ = 0 corresponds to a predominantly biological source of heterogeneity, that is, where CVα is the coefficient of variation for the overall susceptibility. In the opposite limit χ = 1, population heterogeneity is primarily of social origin; hence λ ≈ λs is affected by both CVα and the skewness γα of the pdf f (α). The biological contribution λ b depends on specific biological details of the disease and thus is unlikely to be as universal and robust as the social one. For the COVID-19 epidemic, there is no strong evidence of a wide variation in attack rates unrelated to social activity, geographic location, or socioeconomic status. For instance, there is very little age variability in COVID-19 prevalence as reported by the NYC Department of Health (51) based on the serological survey that followed the first wave of the epidemic. Therefore, below, we will largely ignore possible biological heterogeneity, and focus on social heterogeneity. So far, our discussion has focused on the early stages of epidemics, when the Re (S ) dependence is given by a linearized expression Eq. 9. To describe the nonlinear regime, we consider a gamma-distributed susceptibility: f (α) ≈ α 1/η−1 exp(−α/η), where η = CV 2 α . In this case, according to Eqs. 4 and 5, Re , Se , and S are related by scaling relationships (see SI Appendix), and Re (S ) = R0S λ .

[14]
The exponent λ = 1 + (1 + χ)CV 2 α = 1 + (1 + χ)η coincides with the early-epidemics immunity factor defined in Eqs. 9 and 10 for a general case. Note that, without correlation (χ = 0), both scaling exponents would be the same; this result has been previously obtained for the SIR model in ref. 30 and more recently reproduced in ref. 41 in the context of COVID-19. The scaling behavior Re (S ) is shown in Fig. 1 for λ = 3 ± 1. While this range is arbitrary, it includes the empirical values of λ estimated below. This function is dramatically different from the classical linear dependence Re = SR0. To emphasize the importance of this difference, Eq. 14 immediately leads to a major revision of the classical result for the HIT 1 − S0 = 1 − 1/R0. S0 is the fraction of susceptible population at which the growth stops, while 1 − S0 is the relative size of the epidemic at that time. By setting Re = 1 in Eq. 14, we obtain Nonlinear modifications to homogeneous epidemiological models similar to Eqs. 13 and 14 have been proposed in the past as plausible descriptions of heterogeneous populations in various contexts. Specifically, they were used as empirical fits to simulations of the SIR model on small-world networks (19), as well as to the behavior of the Agent-Based Model on realistic urban contact networks (18). A conceptual explanation of the origin of a nonlinear relation between S and Re was proposed in refs. 19 and 30. However, the scaling law similar to Eqs. 13 and 14 has not been derived except in a special case of the SIR model without correlation between susceptibility and infectivity (30). As we were finalizing this paper for public release, a preprint by Aguas et al. (52) appeared online that independently obtained our Eqs. 14 and 15 for gamma-distributed susceptibilities. The same result has also been recently obtained in ref. 42. Our approach is more general: It provides the exact mapping of a wide class of heterogeneous well-mixed models onto homogeneous ones, and provides a specific relationship between the underlying statistics of α and Rα and the nonlinear functions Re (S ) and Se (S ). Of course, our methodology has the same limitations as the original heterogeneous wellmixed approximation (10). This approximation was shown to provide an adequate description for many classes of networks (19). Additional corrections may still arise, for example, due to clustering and other network structure not captured in its degree (α) distribution. Our focus on the gamma distribution is well justified by the observation that the social strength αs is approximately exponentially distributed, that is, it is a specific case of the gamma distribution with η = CV 2 α = 1 (see more discussion of this in the next section). A moderate biological heterogeneity would lead to an increase of the overall CVα, but the pdf f (α) will still be close to the gamma distribution family. From the conceptual point of view, it is nevertheless important to understand how the function Re (S ) would change if f (α) had a different functional form. In SI Appendix, we present analytic and numerical calculations for two other families of distributions: 1) an exponentially bounded power law f (α) ≈ e −α/α + /α q (q ≥ 1, with an additional cutoff at lower values of α), and 2) the log-normal distribution. In addition, we give an approximate analytic result that generalizes Eq. 14 for an arbitrary skewness of f (α). This generalization works remarkably well for all three families of distributions analyzed in this work. As suggested by Eqs. 11 and 12, as the distribution becomes more skewed, the range between the χ = 0 and χ = 1 curves broadens. For instance, for distributions dominated by a power law, f (α) ≈ 1/α q , with 3 < q < 4 and χ = 1, λ diverges even though CVα remains finite. This represents a crossover to the regime of so-called scale-free networks (2 ≤ q ≤ 3), which are characterized by zero epidemic threshold yet strongly self-limited dynamics: The epidemic effectively kills itself by immunizing the hubs on the network (10,23,53).

Role of Short-Term Variations in Social Activity
Short-term overdispersion in transmission is commonly presumed to have no effect on the overall epidemic dynamics, aside from the early outbreak often dominated by superspreaders. This would indeed be the case if overdispersed transmission were completely uncorrelated with individual susceptibility. But, since the timescale for an individual's infectivity (about 2 d) is comparable to a single generation interval (about 5 d) for the COVID-19 epidemics, ignoring such correlations appears unreasonable. We therefore developed a generalization of the theory presented in the previous section, that takes into account a time dependence of individual susceptibilities and infectivities, as well as temporal correlations between them. The theory is presented, using a path integral formulation, in SI Appendix. Here we present several important results directly related to the transient suppression of an epidemic and differentiate these effects from herd immunity.
Since fast variations are primarily caused by bursty dynamics of social interactions (54)(55)(56)(57), and since heterogeneous biological susceptibility appears subdominant in the context of COVID-19, we set α b = 1 for all individuals. So α has a purely social origin. Let ai (t) = αi + δai (t) be the time-dependent susceptibility of a person, which we associate with a variable level of social activity. Here δai (t) represents the time-dependent deviation of ai (t) from its persistent long-term average αi . Note that index i labels individuals rather than population groups. The level of social activity quantified by ai (t) also determines individual infectivity a(t)R around time t. Interestingly, even the classical result for the basic reproduction number in a heterogeneous system, R0 = R α 2 , needs to be modified due to correlated short-term variations in social activity, Here the bar. . . represents averaging over individual members of the population indexed by i, in contrast with . . . , averaging over all subgroups with various values of persistent heterogeneity α.
In the time-dependent generalization of our theory, Re and S no longer have a fixed functional relationship between them. Instead, this relationship becomes nonlocal in time. For instance, our result for the suppression of Re at the early stages of the epidemic is still formally valid, but the effective value of immunity factor λ becomes time dependent, and Eqs. 9 and 10 become

[19]
Constant λ∞ reflects suppression of Re due to the buildup of the long-term collective immunity. On the other hand, the timedependent term δλ(t ) leads to an additional suppression of Re over intermediate timescales. This term has likely played a signif-icant role in shaping the transient self-limiting dynamics during the first wave of COVID-19 epidemic in some hard-hit locations. Note that, according to Eq. 17, δλ(t ) is being averaged with the weight proportional to the force of infection J (t − t ), since 1 − S (t) ≈ ∞ 0 J (y − t )dt . Since δλ(t, t ) decreases with time difference t , its effect on λ eff should be the strongest during the initial period of fast exponential growth. The initial suppression of the epidemic is caused by the combined effect of mitigation measures and both terms in λ eff . Since λ eff > λ∞, the population may reach the state of TCI earlier than the actual long-term herd immunity determined by persistent heterogeneity. However, this state is fragile and may wane with time. Specifically, as J (t) drops after the first wave, the second term in Eq. 17 gradually decays, bringing λ eff (t) closer to λ∞. According to Eq. 19, it is the correlation time of bursty social activity δa(t) that sets the timescale over which this TCI state deteriorates, and the new epidemic wave may get ignited. The relationship between this relaxation time and the duration of a single epidemic wave also determines the typical value of λ eff during that wave.
Despite a large number of empirical studies of contact networks (54)(55)(56)(57), information about the temporal correlations in α(t) or its proxies remains limited. On the other hand, much more is known about parameters of persistent heterogeneity. Recently, real-world networks of face-to-face communications have been studied using a variety of tools, including Radiofrequency identification (RFID) devices (45), Bluetooth and Wi-Fi wearable tags, and smartphone apps (46,47), as well as census data and personal surveys (13,37,48). Despite coming from a wide variety of contexts, the major features of contact networks are remarkably robust. In particular, pdfs of both the degree (the number of contacts per person) and the node strength plotted in log-log coordinates appear nearly constant, followed by a sharp fall after a certain upper cutoff. This behavior is generally consistent with an exponential distribution in fs (αs ) (19,46,48), f (α) ≈ e −α/ α . That sets the value of η = CV 2 α ≈ 1. If not for short-term overdispersion, that would yield λ = α 3 s / α 2 s = 3!/2! = 3 according to Eq. 12. However, with temporal effects taken into account, the buildup of long-term collective immunity is determined by λ eff (∞) = λ∞. In order to estimate it, we make a simple model assumption that the shortterm overdispersion for a particular individual is proportional to the persistent value of that person's social activity: δa 2 i ≈ αi . This leads to λ∞ = 1 + η(1 + χ * ). [20] Here χ * = α 2 /ai (t) 2 is a parameter that measures the relative strength of persistent heterogeneity and the overdispersion on the timescale of a single generation interval. Note that, formally, we recover our original result for λ in the purely persistent case, with χ * replacing the parameter χ that originally quantified the correlation between infectivity and susceptibility. By assuming the limit of strong short-term overdispersion (χ * 1), we obtain λ∞ ≈ 2. This estimate is consistent with numerical simulations of the agent-based epidemiological model on urban contact networks carried out in ref. 18. As shown in SI Appendix, the very same value of λ∞ should be used as a scaling exponent for long-term behavior of Re (S ). Therefore, HIT is set by Eq. 15 with λ = λ∞ ≈ 2. Its value is plotted vs. R0 in Fig. 2, along with the homogeneous result, and the estimated threshold of TCI. To estimate the corresponding transient immunity factor λ eff , we analyze the empirical data for the first wave of the COVID-19 epidemic, below.

Application to the COVID-19 Epidemic
The COVID-19 epidemic reached the United States in early 2020, and, by March, it was rapidly spreading across multiple states. The early dynamics was characterized by a rapid rise in the number of cases, with doubling times as low as 2 d. In response to this, the majority of states imposed a broad range of mitigation measures including school closures, limits on public gatherings, and stay-at-home orders. In many regions, especially the hardest-hit ones like NYC, people started to practice some degree of social distancing even before government-mandated mitigation. In order to quantify the effects of heterogeneity on the spread of the COVID-19 epidemic, we apply the Bayesian age-of-infection model described in ref. 43 to NYC and Chicago (see SI Appendix for details). For both cities, we have access to reliable time series data on hospitalization, ICU room occupancy, and confirmed daily deaths due to COVID-19 (51,(58)(59)(60). We used these data to perform multichannel calibration of our model (43), which allows us to infer the underlying time progression of both S (t) and Re (t). The fits for Re (S ) for both cities are shown in Fig. 3A. In both cases, a sharp drop of Re that occurred during the early stage of the epidemic is followed by a more gradual decline. For NYC, there is an extended range over which Re (S ) has a constant slope in logarithmic coordinate. This is consistent with the power law behavior predicted by Eq. 14, with the slope corresponding to transient immunity factor λ eff = 4.5 ± 0.05. Chicago exhibits a similar behavior but over a substantially narrower range of S . This reflects the fact that NYC was much harder hit by the COVID-19 epidemic. Importantly, the range of dates we used to estimate the immunity factor corresponds to the time interval after state-mandated stay-at-home orders were imposed, and before the mitigation measures began to be gradually relaxed. The signatures of the onset of the mitigation and of its partial relaxation are clearly visible on both ends of the constant-slope regime. To examine the possible effects of variable levels of mitigation on our estimates of λ eff , in SI Appendix, we repeated our analysis in which Re (t) was corrected by Google's community mobility report in these two cities (61) (see SI Appendix). Although the range of data consistent with the constant slope shrank somewhat, our main conclusion remains unchanged. This provided us with a lower-bound estimate for the transient immunity factor: λ eff = 4.1 ± 0.1.
To test the sensitivity of our results to details of the epidemiological model and choice of the region, we performed an alternative analysis based on the data reported in ref. 49. In that study, the COVID-19 epidemic was modeled in each of the 50 US states and the District of Columbia. Because of the differences in population density, level of urbanization, use of public transport, etc., different states were characterized by substantially different initial growth rates of the epidemic, as quantified by the basic reproduction number R0. Furthermore, the time of arrival of the epidemic also varied a great deal between individual states, with states hosting major airline transportation hubs being among the earliest ones hit by the virus. As a result of these differences, at any given time, the infected fraction of the population differed significantly across the United States (49). We use state level estimates of Re (t), R0, and S (t) as reported in ref. 49 to construct the scatter plot Re (t0)/R0 vs. S (t0) shown in Fig. 3B, with t0 chosen to be the last reported date in that study, May 17, 2020. By performing the linear regression on these data in logarithmic coordinates, we obtain the fit for the slope λ eff = 5.3 ± 0.6 and for S = 1 intercept around 0.54. In SI Appendix, Fig. S3, we present an extended version of this analysis for the 10 hardesthit states and the District of Columbia, which takes into account the overall time progression of Re (t) and S (t), and gives similar estimate λ eff = 4.7 ± 1.5. Both estimates of the immunity factor based on the state data are consistent with our earlier analysis of NYC and Chicago. In light of our theoretical picture, this value of this transient immunity factor, λ eff 4, is set by the pace of the first epidemic wave in the United States. As expected, it exceeds our estimate of λ∞ ≈ 2 associated with persistent heterogeneity and responsible for the long-term herd immunity.
We can now incorporate this transient level of heterogeneity into our epidemiological model, and examine how future projections change as a result of this modification. This is done by plugging scaling relationships given by Eqs. 13 and 14 into the force of infection and incidence rate equations of the original model. These equations are similar to Eqs. 6 and 7, but also include time modulation due to the mitigation and a possible seasonal forcing (see SI Appendix for more details). After calibrating the model by using the data streams on ICU occupancy, hospitalization, and daily deaths up to the end of May, we explore a hypothetical worst-case scenario in which any mitigation is completely relaxed as of June 15, in both Chicago and NYC. In other words, the basic reproduction number R0 is set back to its value at the initial stage of the epidemic, and the only factor limiting the second wave is the partial or full TCI, Re = R0S λ . The projected daily deaths for each of the two cities under this (unrealistically harsh) scenario are presented in Fig. 4 for various values of λ. For both cities, the homogeneous model (λ = 1, blue lines) predicts a second wave which is larger than the first one, with an additional death toll of around 35, 000 in NYC and 12, 800 in Chicago. The magnitude of the second wave is greatly reduced by heterogeneity, resulting in no second wave in either of the two cities for λ = 5 (black lines). Even for a modest value λ = 3 (red lines), which is less than our estimate, the second wave is dramatically reduced in both NYC and Chicago (by about 90% and 70%, respectively).
Note that our predictions about the second wave in NYC and Chicago have been made based on the data up to June 10, 2020 and extended up to early September 2020. The real epidemic dynamic in both cities during this time interval was consistent with the "no second wave" scenario shown in Fig.  4. However, one must be warned against using it to put form bounds on λ eff , since we considered the worst-case scenario of full release of mitigation to prepandemic levels. In reality, some mitigation measures, for example, mask wearing and social distancing, stayed in place. Ultimately, second waves broke out in both cities in the late fall. The mechanisms leading to gradual degradation of the TCI state are described in the next section.

Fragility of TCI
One of the consequences of the bursty nature of social interactions is that the state of TCI gradually wanes due to changes of individual social interaction patterns on timescales longer than a single generation interval. This may be viewed as a slow rewiring of social networks. In the context of the COVID-19 epidemic, individual responses to mitigation factors such as stay-at-home orders may differ across the population. When mitigation measures are relaxed, individual social susceptibility αs inevitably changes. The impact of these changes on collective immunity depends on whether each person's αs during and after the mitigation are sufficiently correlated. For example, the TCI state would be compromised if people who practiced strict selfisolation compensated for it by an above-average social activity after the first wave of the epidemic has passed.
To illustrate the effects of postmitigation rewiring of social networks, we consider a simple modification of the heterogeneous model with no persistent heterogeneity (α = 1 for everyone) and exponentially distributed instantaneous levels of social activity ai (t). This corresponds to λ eff (0) = 3 and λ eff (∞) = λ∞ = 1. In this model, each individual completely changes the set of his/her social connections at some timescale τs . These changes destroy heterogeneity, giving rise to gradual relaxation of susceptible fraction Sa toward its overall mean value S . To model this, we modify Eq. 1 to include a simple relaxation term, Epidemiological models with rewiring of underlying social networks have been studied before (62) (see ref. 63 for a review), but under a constraint that the individual level of social activity quantified by network degree is preserved. In contrast, the dynamics described above stems from the individual level of social activity αs changing in time.
We simulate the full heterogeneous SIR model (see SI Appendix for details) in which Eq. 1 was replaced with Eq. 21. Fig. 5 shows the results of this simulation, where the first wave of the epidemic is mitigated, thereby reducing the effective reproduction number R0 = 2.5. During the course of the mitigation, R0 is multiplied by µ = 0.7. After the mitigation measures have been lifted at the end of the first wave, the population is positioned slightly below the TCI threshold, preventing the immediate start of the second wave. However, gradual rewiring of the social network with time constant τs = 150 d ultimately results in the second and even the third wave of the epidemic (Fig. 5). Fig. 5, Inset shows Re (t) plotted as a function of S (t) in this epidemic. Note that each of the waves follows the power law relationship between Re (t) ≈ S (t) λ predicted by Eq. 14. Since constant rewiring eliminates correlations in individual social activity on scales longer than τs , the epidemic stops after multiple waves bring the total fraction of infected individuals close to the unmodified (homogeneous) HIT 1/R0. Note, however, that, in this case, there is almost no overshoot, and thus the final size of the epidemic is reduced compared to the case of a purely homogeneous and unmitigated epidemic.

Discussion
In this work, we have demonstrated how the interplay between short-term overdispersion and persistent heterogeneity in a population leads to dramatic changes in epidemic dynamics on multiple timescales, transient suppression of the epidemic during its early waves all the way up to the state of long-term herd immunity. First, we developed a general approach that allows  for the persistent heterogeneity to be easily integrated into a wide class of traditional epidemiological models in the form of two nonlinear functions Re (S ) and Se (S ), both of which are fully determined by the statistics of individual susceptibilities and infectivities. Furthermore, Re (S ) is largely defined by a single parameter, the immunity factor λ, introduced in our study. Like susceptibility itself, λ has two contributions: biological and social (Eqs. 11 and 12).
We then expanded our approach to include effects of time dependence of individual social activity, and, in particular, of likely correlations over the timescale of a single generation interval. While our results for purely persistent heterogeneity confirmed and corroborated that HIT would be suppressed compared to the homogeneous case, addition of temporal variations led to a dramatic revision of that simple narrative. Both persistent heterogeneity and short-term overdispersion contributions lead first to a slowdown of a fast-paced epidemic, and to its medium-term stabilization. However, this state of TCI is fragile and does not constitute long-term herd immunity. HIT is indeed suppressed, but only due to the persistent heterogeneity. This suppression is significantly weaker than the initial stabilization responsible for the TCI state reached after the first wave of a fast-paced epidemic.
Among other implications of the TCI phenomenon is the suppression of the so-called overshoot. Namely, it is well known that most models predict that an epidemic will not stop once HIT is passed, ultimately reaching a significantly larger cumulative attack rate, FSE. Multiple prior studies (9,10,20,21,30,34) have shown that FSE would be suppressed by persistent heterogeneity, similarly to HIT. In SI Appendix, we present a simple result that unifies several previously studied limiting cases, and gives an explicit equation for the FSE for the gammadistributed susceptibility and variable level of its correlation with infectivity. However, because of the transient suppression of the early waves of the epidemic discussed in this work, the overshoot effect would be much weaker or essentially eliminated. For instance, our simple rewiring model demonstrates how the epi-demic, after several waves, ultimately reaches HIT level, but does not progress much beyond it (Fig. 5). The FSE result may still be used, but primarily as an estimate for the size of the first wave of an (unmitigated) epidemic. In that case, the transient value of immunity factor λ eff should be assumed.
By applying our theory to the COVID-19 epidemic, we found evidence that the hardest-hit areas, such as NYC, have likely passed TCI threshold by the end of the first wave, but are less likely to have achieved real long-term herd immunity. Other places that had intermediate exposure, such as, for example, Chicago, while still below the TCI threshold, have their effective reproduction number reduced by a significantly larger factor than predicted by traditional epidemiological models. This gives a better chance of suppressing the future waves of the epidemic in these locations by less disruptive measures than those used during the first wave, for example, by using masks, social distancing, contact tracing, control of potential superspreading events, etc. However, similar to the case of NYC, transient stabilization of the COVID-19 epidemic in Chicago will eventually wane. As for the permanent HIT, although suppressed compared to classical value, it definitely has not yet been passed in those two locations.
In a recent study (35), the reduction of HIT due to heterogeneity has been illustrated using a toy model. In that model, 25% of the population was assumed to have their social activity reduced by 50% compared to a baseline, while another 25% had their social activity elevated twofold. The rest of the population was assigned the baseline level of activity. According to Eq. 12, the immunity factor in that model is λ = 1.54. For this immunity factor, Eq. 15 predicts HIT at S0 = 64%, 55%, and 49%, for R0 = 2, 2.5, and 3, respectively. Despite the fact that the model distribution is not gamma shaped, these values are in a very good agreement with the numerical results reported in ref. 35: S0 = 62.5%, 53.5%, and 47.5%, respectively.
Thus there is a crucial distinction between the persistent heterogeneity, short-term variations correlated over the timescale  The transient and long-term suppression coefficients Re/R 0 are calculated using λ eff = 4 and λ∞ = 2, respectively. Fraction of susceptible population S as of early June 2020 is estimated from the cumulative reported death count per capita, assuming the infection fatality rate of 0.7% (67). This estimate is supplemented by seroprevalence data for late spring-early summer 2020. of a single generation interval, and overdispersion in transmission statistics associated with short-term superspreading events (12,13,16,(26)(27)(28). In our theory, a personal decision to attend a large party or a meeting would only contribute to persistent heterogeneity if it represents a recurring behavioral pattern. On the other hand, superspreading events are shaped by shorttime variations in individual infectivity (e.g., a person during the highly infectious phase of the disease attending a large gathering). Hence, the level of heterogeneity inferred from the analysis of such events (12,27) would be significantly exaggerated and should not be used to estimate the TCI threshold and HIT. Specifically, the statistics of superspreading events is commonly described by the negative binomial distribution with dispersion parameter k estimated to be in the range 0.1 to 0.3 for COVID-19 (28,64,65). This is much stronger overdispersion than the value k = 1 estimated from persistent heterogeneity, based on the exponential distribution of α. Thus, persistent heterogeneity is a weaker source of variation compared to shortterm variations. According to ref. 12, this is consistent with the expected value of the individual-level reproduction number Ri drawn from a gamma distribution with the shape parameter k 0.1 . . . 0.3. This distribution has a very high coefficient of variation, CV 2 = 1/k 3 . . . 10. In the case of a perfect correlation between individual infectivity and susceptibility α, this would result in an unrealistically high estimate of the immunity factor: λ = 1 + 2CV 2 = 1 + 2/k 7 . . . 20. For this reason, according to our perspective and calculation, the final size of the COVID-19 epidemic may have been substantially underestimated in ref. 31. Similarly, the degree of heterogeneity assumed in other recent studies (32,52) is considerably larger than our estimates. Based on our analysis, the value of the immunity factor λ depends on the pace of the epidemic and on the timescale under consideration. We estimated its long-term value (responsible for the permanent HIT) as λ∞ ≈ 2. However, the transient values are expected to be higher, especially during the first several waves of COVID-19 in select locations, characterized by large growth rates. Our analysis of the empirical data in NYC and Chicago indicates that the slowdown of the epidemic dynamics in those locations was consistent with λ ≈ 4. In Table 1, we present our estimates of the factor by which Re is transiently suppressed as a result of depletion of susceptible population in selected locations in the world, as of early June 2020, as well as the predicted long-term suppression related to acquisition of a partial herd immunity.
Finally, we summarize the assumptions and limitations of our study. First, we assume a long-lasting biological immunity of recovered individuals. Second, our approach is based on the well-mixed approximation in which geographic heterogeneity as well as nontrivial properties of the contact network (clustering, degree-degree correlations, etc.) are ignored. In addition, our description of transient epidemic dynamics is based on the approximation of a constant λ eff , while gradual degradation of TCI with time is illustrated using a simplified model, Eq. 21. A generalization of the present model including explicit description of stochastic social activity is needed (see ref. 66). Furthermore, additional calibration based on long-term empirical data is required before our approach can be used for reliably guiding policy decisions during an epidemic.
Population heterogeneity manifests itself at multiple scales. At the most coarse-grained level, individual cities or even countries can be assigned some level of susceptibility and infectivity, which inevitably vary from one location to another, reflecting differences in population density and its connectivity to other regions. Such spatial heterogeneity will result in self-limiting epidemic dynamics at the global scale. For instance, hard-hit hubs of the global transportation network, such as NYC during the COVID-19 epidemic, would gain full or partial TCI, thereby limiting the spread of infection to other regions during future waves of the epidemic. This might be a general mechanism that ultimately limits the scale of many pandemics, from the Black Death to the 1918 influenza.
Data Availability. All study data are included in the article and SI Appendix.