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

^{a}Center for Functional Nanomaterials, Brookhaven National Laboratory, Upton, NY 11973;^{b}Department of Physics, University of Illinois at Urbana–Champaign, Urbana, IL 61801;^{c}Department of Bioengineering, University of Illinois at Urbana–Champaign, Urbana, IL 61801;^{d}Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana–Champaign, Urbana, IL 61801;^{e}Department of Civil Engineering, University of Illinois at Urbana–Champaign, Urbana, IL 61801;^{f}Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, IL 61801

See allHide authors and affiliations

Edited by Andrea Rinaldo, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, and approved March 2, 2021 (received for review July 28, 2020)

## 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.

## Abstract

Epidemics generally spread through a succession of waves that reflect factors on multiple timescales. On short timescales, superspreading events lead to burstiness and overdispersion, whereas long-term persistent heterogeneity in susceptibility is expected to lead to a reduction in both the infection peak and the herd immunity threshold (HIT). Here, we develop a general approach to encompass both timescales, including time variations in individual social activity, and demonstrate how to incorporate them phenomenologically into a wide class of epidemiological models through reparameterization. We derive a nonlinear dependence of the effective reproduction number

The 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⇓⇓⇓–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⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–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⇓⇓–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⇓⇓⇓–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 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 short-term, 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⇓–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 COVID-19), 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 short-term 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

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

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:

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

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⇓–24, 30, and 32, we consider the spread of an epidemic in a population of individuals who exhibit significant persistent heterogeneity in their susceptibilities to infection α. This heterogeneity may be biological or social in origin, and we assume these factors are independent:

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,

Let

According to Eq. **1**, the fraction of susceptibles in the subpopulation with any given α can be expressed as**5** can be omitted. Since both S and **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, **5** was derived by substituting Eq. **1** into Eq. **2**. This leads to the renewal equation for **1** over all values of α, one arrives at the following heterogeneity-induced modification to the relationship between the force of infection and incidence rate:**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 *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

We previously defined the overall susceptibility as a combination of biological and social factors: **10** in the limit **12** is proportional to the third moment of

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 *SI Appendix* for details). The corresponding expression for the overall immunity factor is

So far, our discussion has focused on the early stages of epidemics, when the **9**. To describe the nonlinear regime, we consider a gamma-distributed susceptibility: **4** and **5**, *SI Appendix*),**9** and **10** for a general case. Note that, without correlation (**14** immediately leads to a major revision of the classical result for the HIT **14**, we obtain**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 **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

Our focus on the gamma distribution is well justified by the observation that the social strength *SI Appendix*, we present analytic and numerical calculations for two other families of distributions: 1) an exponentially bounded power law **14** for an arbitrary skewness of **11** and **12**, as the distribution becomes more skewed, the range between the

## 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⇓⇓–57), and since heterogeneous biological susceptibility appears subdominant in the context of COVID-19, we set

In the time-dependent generalization of our theory, **9** and **10** become

Note that, according to Eq. **17**, **17** gradually decays, bringing **19**, it is the correlation time of bursty social activity

Despite a large number of empirical studies of contact networks (54⇓⇓–57), information about the temporal correlations in **12**. However, with temporal effects taken into account, the buildup of long-term collective immunity is determined by

As shown in *SI Appendix*, the very same value of **15** with

## 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⇓–60). We used these data to perform multichannel calibration of our model (43), which allows us to infer the underlying time progression of both *A*. In both cases, a sharp drop of **14**, with the slope corresponding to transient immunity factor *SI Appendix*, we repeated our analysis in which *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:

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 *B*, with *SI Appendix*, Fig. S3, we present an extended version of this analysis for the 10 hardest-hit states and the District of Columbia, which takes into account the overall time progression of

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

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

## 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

To illustrate the effects of postmitigation rewiring of social networks, we consider a simple modification of the heterogeneous model with no persistent heterogeneity (**1** to include a simple relaxation term,

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 *Inset* shows **14**. Since constant rewiring eliminates correlations in individual social activity on scales longer than

## 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 **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 gamma-distributed 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 epidemic, 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

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, **12**, the immunity factor in that model is **15** predicts HIT at

Thus there is a crucial distinction between the persistent heterogeneity, short-term variations correlated over the timescale of a single generation interval, and overdispersion in transmission statistics associated with short-term superspreading events (12, 13, 16, 26⇓–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 short-time 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

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 **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*.

## Acknowledgments

We gratefully acknowledge discussions with Mark Johnson at Carle Hospital. The calculations we have performed would have been impossible without the data kindly provided by the Illinois Department of Public Health through a Data Use Agreement with Civis Analytics. This work was supported by the University of Illinois System Office, the Office of the Vice-Chancellor for Research and Innovation, the Grainger College of Engineering, and the Department of Physics at the University of Illinois at Urbana–Champaign. Z.J.W. is supported in part by the US Department of Energy (DOE) Computational Science Graduate Fellowship, provided under Grant DE-FG02-97ER25308. A.E. acknowledges partial support by NSF CAREER Award 1753249. This work made use of the Illinois Campus Cluster, a computing resource that is operated by the Illinois Campus Cluster Program in conjunction with the National Center for Supercomputing Applications and which is supported by funds from the University of Illinois at Urbana–Champaign. This research was partially done at, and used resources of, the Center for Functional Nanomaterials, which is a US DOE Office of Science Facility, at Brookhaven National Laboratory under Contract DE-SC0012704.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: oleksiyt{at}bnl.gov or maslov{at}illinois.edu.

Author contributions: A.V.T., S.M., A.E., and N.G. designed research; A.V.T., S.M., A.E., G.N.W., Z.J.W., and N.G. performed research; A.V.T., S.M., A.E., G.N.W., Z.J.W., and N.G. analyzed data; and A.V.T., S.M., A.E., G.N.W., Z.J.W., and N.G. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2015972118/-/DCSupplemental.

- Copyright © 2021 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY).

## References

- ↵
- ↵
- M. J. Keeling,
- P. Rohani

- ↵
- K. Rock,
- S. Brand,
- J. Moir,
- M. J. Keeling

- ↵
- J. Ma

- ↵
- ↵
- G. Chowell

- ↵
- A. L. Lloyd,
- R. M. May

- ↵
- R. M. May,
- A. L. Lloyd

- ↵
- M. E. J. Newman

- ↵
- Y. Moreno,
- R. Pastor-Satorras,
- A. Vespignani

- ↵
- Z. Dezső,
- A. L. Barabási

- ↵
- ↵
- ↵
- ↵
- ↵
- M. Small,
- C. Tse,
- D. M. Walker

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- W. Gou,
- Z. Jin

- ↵
- Y. Kim,
- H. Ryu,
- S. Lee

- ↵
- ↵
- ↵
- A. Endo,
- S. Abbott,
- A. Kucharski,
- S. Funk

- ↵
- J. C. Miller

- ↵
- ↵
- L. Hébert-Dufresne,
- B. M. Althouse,
- S. V. Scarpino,
- A. Allard

- ↵
- M. G. M. Gomes et al.

- ↵
- P. V. Brennan,
- L. P. Brennan

- ↵
- F. Ball

- ↵
- T. Britton,
- F. Ball,
- P. Trapman

- ↵
- J. C. Stack,
- S. Bansal,
- V. S. A. Kumar,
- B. Grenfell

- ↵
- ↵
- ↵
- ↵
- ↵
- C. Rose et al.

- ↵
- J. Neipel,
- J. Bauermann,
- S. Bo,
- T. Harmon,
- F. Jülicher

- ↵
- G. N. Wong et al.

- ↵
- M. Castro,
- S. Ares,
- J. A. Cuesta,
- S. Manrubia

- ↵
- ↵
- B. F. Nielsen,
- K. Sneppen,
- L. Simonsen,
- J. Mathiesen

- ↵
- G. L. Ciampaglia,
- A. Mashhadi,
- T. Yasseri

- M. Starnini et al.

- ↵
- ↵
- H. J. T. Unwin et al.

- ↵
- ↵NYC Department of Health and Mental Hygiene, Data from “Covid-19 County Cases, Tests, and Death by Day.” GitHub. https://github.com/nychealth/coronavirus-data. Accessed 20 June 2020.
- ↵
- R. Aguas et al.

- ↵
- ↵
- ↵
- G. Kossinets,
- D. J. Watts

- ↵
- D. Rybski,
- S. V. Buldyrev,
- S. Havlin,
- F. Liljeros,
- H. A. Makse

- ↵
- J. Saramäki,
- E. Moro

- ↵Illinois Department of Public Health, Data from “COVID-19 statistics.” Illinois Department of Public Health. https://www.dph.illinois.gov/content/covid-19-county-cases-tests-and-deaths-day. Accessed 20 June 2020.
- ↵The City, Data from “Covid-19 New York City Data.” GitHub. https://github.com/thecityny/covid-19-nyc-data. Accessed 20 June 2020.
- ↵The City, Data from “Coronavirus in New York City.” https://projects.thecity.nyc/2020_03_covid-19-tracker/. Accessed 20 June 2020.
- ↵
- J. Fitzpatrick,
- K. DeSalvo

*The Keyword*(2020). https://www.blog.google/technology/health/covid-19-community-mobility-reports?hl=en. Accessed 14 February 2021. - ↵
- ↵
- ↵
- Z. Susswein,
- S. Bansal

- ↵
- K. Sun et al.

- ↵
- A. V. Tkachenko et al.

- G. Meyerowitz-Katz,
- L. Merone

- Centers for Disease Control and Prevention

- Presidenza del Consiglio dei Ministri - Dipartimento della Protezione Civile, Data from “COVID-19.” GitHub. https://github.com/pcm-dpc/COVID-19. Accessed 20 June 2020.
- London Assembly, Data from “Coronavirus (COVID-19) Death.” https://data.london.gov.uk/dataset/coronavirus–covid-19–deaths. Accessed 15 December 2020.
- H. Ward et al.

- A. R. Demonbreun et al.

- M. Runesson, Data from “Covid-19/Corona data from Stockholm county project.” GitHub. https://github.com/mrunesson/covid-19. Accessed 15 December 2020.
- X. C. Dopico et al.

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Population Biology