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: α=α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:−dSαdt=jα(t)=αSα(t)J(t)[1]J(t)=∫0∞RαK(τ)jα(t−τ)dτ.[2]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-of-infection 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(τ). 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 asSα(t)=exp(−αZ(t)).[3]Here Z(t)≡∫0tJ(t′)dt′. 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 toS(t)=∫0∞f(α)e−αZ(t)dα=Mα(−Z(t)).[4]Similarly, the effective reproductive number Re can be expressed in terms of the parameter Z,Re(t)=1⟨α⟩∫0∞αRαf(α)e−αZ(t)dα.[5]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=⟨kinkout⟩/⟨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,J(t)=∫0∞dτK(τ)Re(t−τ)J(t−τ).[6]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:dSdt=−SeJ.[7]HereSe(t)=∫0∞αf(α)e−αZ(t)dα=−dMα(−Z(t))dZ[8]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 asRe≈R0(1−λ(1−S)).[9]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λ=α2Rαα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,λb=αb2αb2=1+CVb2[11]λs=⟨αs3⟩⟨αs⟩⟨αs2⟩=1+CVs2(2+γsCVs)1+CVs2.[12]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α⟩≈αs2. Note that the biological contribution to the immunity factor depends only on the coefficient of variation CVb 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, λ≈λb=1+CVα2, 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),Se(S)=S1+η[13]andRe(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 obtain1−S0=1−1R01/λ.[15]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 well-mixed 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Re/R vs. S dependence for gamma-distributed susceptibility with λ=3±1 (light blue area). The dashed line shows the classical homogeneous result, Re=R0S. The dark blue region corresponds to the range 2<R0<3.5 representative for COVID-19.

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⇓⇓–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,R0=R⟨α2⟩+δai2¯.[16]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λeff(t)=λ∞+11−S(t)∫0∞δλ(t,t′)J(t−t′)dt′[17]λ∞=⟨α3⟩+αiδai2¯⟨α2⟩+δai2¯[18]δλ(t,t′)=δai2(t)δai(t−t′)¯⟨α2⟩+δai2¯.[19]Constant λ∞ reflects suppression of Re due to the buildup of the long-term collective immunity. On the other hand, the time-dependent term δλ(t′) leads to an additional suppression of Re over intermediate timescales. This term has likely played a significant 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⇓⇓–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 Radio-frequency 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 λ=⟨αs3⟩/⟨αs2⟩=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 short-term overdispersion for a particular individual is proportional to the persistent value of that person’s social activity: δai2¯≈α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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

TCI threshold (dot-dashed), long-term heterogeneous HIT (solid), and homogeneous HIT (dotted) for various values of R0 (x axis). HIT (solid line) is determined by persistent heterogeneity. The corresponding immunity factor λ∞≈2 was estimated from Eq. 20 assuming strong short-term overdispersion (χ*≪1) and the exponential distribution of α (η=1). For transient behavior, λeff≈4 is assumed based on analysis of empirical data for COVID-19 epidemic in select locations.

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

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Correlation between the relative reduction in the effective reproduction number Re(t)/R0 (y axis) with the susceptible population S(t). (A) The progression of these two quantities for NYC and Chicago, as given by the epidemiological model described in ref. 43. (B) The scatter plot of Re(t0)/R0 and S(t0) in individual states of the United States, evaluated in ref. 49 (t0 is the latest date covered in that study).

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 hardest-hit 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).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Projections of daily deaths under the hypothetical scenario in which any mitigation is completely eliminated as of June 15, 2020, for (A) NYC and (B) Chicago. Different curves correspond to different values of the transient immunity factor λeff=1 (blue), 3 (red), 4 (green), and 5 (black lines). The model described in ref. 43 was fully calibrated on daily deaths (circles), ICU occupancy, and hospitalization data up to the end of May. See SI Appendix for additional details, including CIs.

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 self-isolation 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,Ṡa=−αSaJ−1τs(Sa−S).[21]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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Effect of social rewiring on the epidemic dynamics. The time course of an epidemic in a heterogeneous SIR model with R0=2.5 and λ=3. During the first 100 d, a mitigation factor μ=0.7 is applied. Social networks gradually rewire with a time constant τs=150 d. The figure shows multiple waves. (Inset) Re(t) plotted as a function of S(t). Solid black line shows the homogeneous limit reached after multiple waves.