Stochastic social behavior coupled to COVID-19 dynamics leads to waves, plateaus and an endemic state
Alexei V. Tkachenko, Sergei Maslov, Tong Wang, Ahmed Elbanna, George N. Wong, Nigel Goldenfeld
SStochastic social behavior coupled to COVID-19 dynamics leads to waves, plateausand an endemic state
Alexei V. Tkachenko † , Sergei Maslov , , † , Tong Wang , ,Ahmed Elbanna , George N. Wong , and Nigel Goldenfeld , Center for Functional Nanomaterials, Brookhaven National Laboratory, Upton, NY 11973, USA Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Department of Civil Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Department of Bioengineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA (Dated: February 22, 2021)It is well recognized that population heterogeneity plays an important role in the spread of epi-demics. While individual variations in social activity are often assumed to be persistent, i.e. constantin time, here we discuss the consequences of dynamic heterogeneity. By integrating the stochasticdynamics of social activity into traditional epidemiological models we demonstrate the emergenceof a new long timescale governing the epidemic in broad agreement with empirical data. Our modelcaptures multiple features of real-life epidemics such as COVID-19, including prolonged plateaus andmultiple waves, which are transiently suppressed due to the dynamic nature of social activity. Theexistence of the long timescale due to the interplay between epidemic and social dynamics providesa unifying picture of how a fast-paced epidemic typically will transition to the endemic state.
The COVID-19 pandemic has underscored the promi-nent role played by population heterogeneity in epi-demics. On one hand, the observed transmission of in-fection is characterized by the phenomenon of super-spreading, in which a small fraction of individuals areresponsible for a disproportionally large number of sec-ondary infections [1–4]. On the other hand, accordingto multiple models, population heterogeneity is expectedto suppress the herd immunity threshold (HIT) and re-duce the final size of an epidemic [5–10]. In the con-text of COVID-19, this observation led to a controversialsuggestion that a strategy relying exclusively on quicklyreaching herd immunity might be a viable alternative togovernment-imposed mitigation [11]. However, the expe-rience of locations that had embraced that strategy hasexposed its flaws. While the first wave of infections inthose locations never reached the scale of an unmitigatedepidemic predicted by classical homogeneous models, italso failed to provide long-lasting protection against newwaves [12].Another puzzling aspect of the COVID-19 pandemicis plateau-like dynamics, where the incidence rate staysat an approximately constant level for a prolonged time[13–15]. These dramatic departures from predictions ofboth classical epidemiological models and their hetero-geneous extensions have led to a greater appreciation ofthe role played by human behavior in epidemic dynam-ics. In particular, one plausible mechanism that might beresponsible for both suppression of the early waves andplateau-like dynamics is that individuals modify their be-havior based on information about the current epidemio-logical situation [15–17]. Another possibility is that longplateaus might arise because of the underlying structureof social networks [14].Here we study epidemic dynamics accounting for ran-dom changes in levels of individual social activity. Wedemonstrate that this type of dynamic heterogeneity, even without knowledge-based adaptation of human be-havior (e.g. in response to epidemic-related news) [15–17], leads to a substantial revision of the epidemicprogression, consistent with the empirical data for theCOVID-19 pandemic. In a recent study [8] we havepointed out that population heterogeneity is a dynamicproperty that roams across multiple timescales. A strongshort-term overdispersion of the individual infectivitymanifests itself in the statistics of super-spreading events.At the other end of the spectrum is a much weaker per-sistent heterogeneity operating on very long timescales.In particular, it is this long-term heterogeneity thatleads to a reduction of the HIT compared to that pre-dicted by classical homogeneous models [7–10, 18]. How-ever, the epidemic dynamics is also sensitive to transienttimescales over which the bursty short-term social activ-ity of each individual crosses over to its long-term aver-age. By including the effects of dynamic heterogeneity,we demonstrate that a suppression of the early waves ofthe COVID-19 epidemic, even without active mitigation,does not signal achievement of long-term herd immunity.Instead, it is associated with Transient Collective Immu-nity, a fragile state which degrades over time as indi-viduals change their social activity patterns [8]. As wedemonstrate below, the first wave is generally followedeither by secondary waves or by long plateaus character-ized by a nearly constant incidence rate. In the context ofCOVID-19, both long plateaus and multi-wave epidemicdynamics have been commonly observed [13]. Accord-ing to our analysis, the number of daily infections duringthe plateau regime, as well as the individual wave trajec-tories, are robust properties of the epidemic and dependon the current level of mitigation, degree of heterogeneityand temporal correlations of individual social activity.Our work implies that, once the plateau-like dynamicsis established, the epidemic gradually evolves towards thelong-term HIT determined by persistent population het- a r X i v : . [ q - b i o . P E ] F e b erogeneity. However, reaching that state may stretch overa surprisingly long time, from months to years. On theselong timescales, both waning of individual biological im-munity and mutations of the pathogen become valid con-cerns, and would ultimately result in a permanent en-demic state of the infection. Such endemic behavior isa well-known property of most classical epidemiologicalmodels [19]. However, the emergence of the endemic statefor a newly introduced pathogen is far from being com-pletely understood [20–22]. Indeed, most epidemiologicalmodels would typically predict complete extinction of apathogen following the first wave of the epidemic, well be-fore the pool of susceptible population would be replen-ished. A commonly accepted, though mostly qualitative,explanation for the onset of endemic behavior of such dis-eases as measles, seasonal cold, etc., involves geographicheterogeneity: the pathogen may survive in other geo-graphic locations until returning to a hard hit area with adepleted susceptible pool [20, 21]. In contrast, our theoryprovides a simple and general mechanism that preventsan overshoot of the epidemic dynamics and thus natu-rally and generically leads to the endemic fixed point.The importance of temporal effects has been long rec-ognized in the context of network-based epidemiologicalmodels [23–26]. On one hand, available high-resolutiondata on real-world temporal contact networks allows di-rect modeling of epidemic spread on those networks.On the other hand, building upon successes of epidemicmodels on static unweighted networks [5, 27–29], a va-riety of their temporal generalizations have been pro-posed. Those typically involve particular rules for dis-crete or continuous network rewiring [24–26] such as e.g.in activity-based network models [17, 30, 31]. While im-portant theoretical results have been obtained for some ofthese problems, especially regarding the epidemic thresh-old, many open questions and challenges remain in thefield. In this paper, we start with a more traditional het-erogeneous well-mixed model, which is essentially equiv-alent to the mean-field description of an epidemic on anetwork [6, 29, 32], and include effects of time-variablesocial activity that modulates levels of individual suscep-tibilities and infectivities.The basic idea behind our model is represented in Fig.1. Each individual i is characterized by time-dependentsocial activity a i ( t ) proportional to his/her current fre-quency and intensity of close social contacts. This quan-tity determines both individual susceptibility to infectionas well as ability to infect others. The time evolution ofcontact frequency, and hence a i ( t ), is in principle mea-surable by means of proximity devices, such as RFID,Bluetooth, Wi-Fi, etc. [5, 23, 33, 34]. In our model wecombine a simple mathematical description of social dy-namics with the standard Susceptible-Infected-Removed(SIR) epidemiological model. Qualitatively it leads tolong-term epidemic dynamics fuelled by replenishmentof susceptible population due to changes in the level ofindividual social activity from low to high. Fig. 1(a)illustrates this process by showing people with low so- ab FIG. 1. Schematic illustration of our model in which each in-dividual is characterized by a time-dependent social activity.a) People with low social activity (depicted as socially isolatedfigures at home) occasionally increase their level of activity(depicted as a party). The average activity in the populationremains the same, but individuals constantly change their ac-tivity levels from low to high (arrows pointing up) and back(arrows pointing down). Individuals are colored according totheir state in the SIR epidemiological model: susceptible -green, infected - red, and removed - blue. The epidemic is fu-elled by constant replenishment of susceptible population withhigh activity due to transitions from the low activity state.b) examples of individual time-dependent activity a i ( t ) (solidlines), with different persistent levels α i (dot-dashed lines).S,I,R states of an individual have the same color code as in(a). Note that pathogen transmission occurs predominantlybetween individuals with high current activity levels. cial activity (depicted as socially isolated at home) oc-casionally increasing their level of activity (depicted as aparty). Fig. 1(b) represents the same dynamics in termsof individual functions a i ( t ). Note that each person ischaracterized by his/her own long-term average activitylevel α i (dot-dashed lines), but the transmission occurspredominantly between individuals with high levels of current social activity. This is because a i ( t ) determinesboth current susceptibility and individual infectivity of aperson. However, the secondary transmission is delayedwith respect to the moment of infection, by a time of theorder of a single generation interval τ g (around 5 daysfor COVID-19). Studies of real-world contact and inter-personal communication networks have shown that indi-vidual social activity is bursty and varies across multipletimescales—from seconds to years [35–38].For any individual i the value of a i ( t ) has a tendencyto gradually drift towards its persistent average level α i ,which itself varies within the population. In our model,we assign a single timescale τ s to this mean reversionprocess. This is of course a simplification of the multi-scale relaxation observed in real social dynamics. While τ s can be treated as a fitting parameter of our model,here we simply set it to be τ s = 30 days, several foldlonger than the mean generation interval of COVID-19, τ g = 5 days. Note that from the point of view of theepidemic dynamics, variations in activity on timescalesshorter than the mean generation interval may be safelyignored. For example, attending a single party wouldincrease an individual’s risk of infection but would notchange his/her likelihood of transmission to others 5 dayslater.The individual social activity a i ( t ) is governed by thefollowing stochastic equation:˙ a i ( t ) = α i − a i ( t ) τ s + η i ( t ) (1)Here η ( t ) is a short-time noise that gives rise to time-dependent variations in a i ( t ). We set (cid:104) η i ( t ) η i ( t (cid:48) ) (cid:105) = a i ( t ) τ s k δ ( t − t (cid:48) ), which results in a diffusion in the spaceof individual social activity with diffusion coefficient pro-portional to a i . This stochastic process is well known inmathematical finance as the Cox–Ingersoll–Ross (CIR)model [39] and has been studied in probability theorysince 1950s [40]. The major properties of this model are(i) reversion to the mean and (ii) non-negativity of a i atall times, both of which are natural for social activity.Furthermore, the steady state solution of this model ischaracterized by the Gamma-distributed a i . This is con-sistent with the empirical statistics of short-term overdis-persion of disease transmission manifesting itself in su-perspreading events [1, 3, 4]. More specifically, for agiven level of persistent activity α , this model generatesthe steady-state distribution of “instantaneous” values ofsocial activity a following gamma distribution with mean α and variance α/k : f α ( a ) ∼ a k α − e − k a .The statistics of super-spreader events is usually repre-sented as a negative binomial distribution, derived froma gamma-distributed individual reproduction number [1].The observed overdispersion parameter k ≈ . − . k . In addition, we assume persistentlevels of social activity α i to also follow a Gamma dis-tribution with another dispersion parameter, κ . In sev-eral recent studies of epidemic dynamics in populationswith persistent heterogeneity [8, 9, 41] it has been demon-strated that κ determines the herd immunity threshold.Multiple studies of real-world contact networks (summa-rized, e.g. in [6]) report an approximately exponentialdistribution of α , which corresponds to κ (cid:39)
1. Through-out this paper, we assume a more conservative value, κ = 2, i.e. coefficient of variation 1 /κ = 0 .
5, half waybetween the fully homogeneous case and that with ex-ponentially distributed α . For consistency with the re-ported value of the short-term overdispersion parameter[4], 1 /k ≈ /κ + 1 /k ≈
3, we set k = 0 . α , effec-tively diffuse in the space of their current social activity a . This leads to major modifications of the epidemicdynamics. For instance, the equation for the suscepti-ble fraction in classical epidemic models [19] acquires thefollowing form:˙ S α ( a, t ) = (cid:20) − aJ ( t ) + ak τ s ∂ ∂a + α − aτ s ∂∂a (cid:21) S α ( a, t )(2)Here S α ( a, t ) is the fraction of susceptible individualswithin a subpopulation with a given value of persistentsocial activity α and with current social activity a , atthe moment of infection, and J ( t ) is the current strengthof infection. Its time evolution can be described by anytraditional epidemiological model, such as e.g. age-of-infection, SIR/SEIR, etc [19]. 𝝉 𝒔 𝝉 𝒈 FIG. 2. Schematic representation of the mechanisms that leadto self-limited epidemic dynamics. In traditional epidemicmodels, the major factor is the depletion of susceptible pop-ulation. Government-imposed mitigation and/or behavioraladaptation to the perceived risk create another feedback loop(purple). Yet another mechanism, is due to the dynamic het-erogeneity of the attack rate parameterized by h ( t ) (black).On the one hand, the attack rate heterogeneity is being gener-ated by the current infection. On the other hand, it suppressesitself on the timescale of τ s due the reversion of individual so-cial activity towards the mean (black feedback loop). Eq. (2), is dramatically simplified by writing S α ( a, t ) ≡ e − Z ( t ) α − k h ( t ) a , implicitly defining Z ( t ) as ameasure of persistent heterogeneity of the attack rateand the transient heterogeneity parameter h ( t ). As Z in-creases, so does the difference in depletion of susceptiblesamong subpopulations with different α ’s, i.e. various lev-els of persistent social activity. On the other hand, h ( t )parametrizes the transient heterogeneity within each ofthese subpopulations. Both Z ( t ) and h ( t ) indicate thecurrent level of heterogeneity of the attack rate i.e. sus-ceptible population structure. In the long run, transientheterogeneity disappears due to the diffusion in the a -space, thus h ( t ) asymptotically approaches 0 as t → ∞ .We combine this ansatz with a general methodology [8]that provides a quasi-homogeneous description for a widevariety of heterogeneous epidemiological models. For aspecific case of SIR dynamics, we assign each person astate variable I i set to 1 when the individual is infectiousand 0 otherwise. Now, the activity-weighted fraction ofthe infected population is defined as I ( t ) = (cid:104) I i a i ( t ) (cid:105) / (cid:104) a i (cid:105) ,and the current infection strength is proportional to it: J ( t ) = R M ( t ) I ( t ) /τ g . Here M ( t ) is a time-dependentmitigation factor, which combines the effects of govern-ment interventions, societal response to the epidemic, aswell other sources of time modulation, such as e.g. sea-sonal forcing.Using this ansatz, the epidemic in a population withboth persistent and dynamic heterogeneity of individualsocial activity can be compactly described as a dynamicalsystem with only three variables: the susceptible popu-lation fraction S ( t ), the infected population fraction I ( t )(activity-weighted) which is proportional to strength ofinfection J ( t ), and the transient heterogeneity parameter h ( t ). In the ( S, I, h )-space, the dynamics is given by thefollowing set of differential equations: dIdt = JS λ (1 + h ) − Iτ g (3) dSdt = − JS /κ (1 + h ) (4) dhdt = Jk − h (1 + h ) τ s (5)As in the case of persistent heterogeneity without tempo-ral variations [8], the long-term herd immunity threshold,1 − ( R M ) − /λ , is determined by the immunity factor λ .The latter depends both on the short-term and persistentdispersion parameters: λ = (cid:0) κ − (cid:1) (cid:0) k − g + 2 κ − (cid:1) k − g + κ − (6)Here k g = k (1 + τ g /τ s ) is the dispersion parameter foractivity a ( t ) averaged over a timescale of a single genera-tion interval. For parameters k = 0 . κ = 2, τ s = 30days used throughout our study, k g = 0 . λ = 1 .
7, con-sistent with our earlier estimate of λ ∞ ≈ Time / days M ( t ) R a D a il y i n c i d e n c e p e r b HomogeneousHeterogeneous (persistent only)Heterogeneous (dynamic)ABM, homogeneousABM, heterogeneous (persistent only)ABM, heterogeneous (dynamic)0 50 100 150 200 250 300 350
Time / days A tt a c k r a t e HIT, heterogeneousHIT, homogeneous c FIG. 3. (a)-(c) Comparison of the epidemic dynamics withhomogeneous population (black curves), persistent popula-tion heterogeneity (brown curves), and with dynamic het-erogeneity (green curves): mitigation profile (a), daily inci-dence (b), and cumulative attack rate (c). While parametersin cases (b) and (c) correspond to the same herd immunitythreshold (HIT), the behavior is drastically different. In thepersistent model, the epidemic quickly overshoots above HITlevel. In the case of dynamic heterogeneity, the initial waveis followed by a plateau-like behavior with slow relaxationtowards the HIT. Note an excellent agreement between thequasi-homogeneous theory described by Eqs. (3- 5) (solidlines) and the Agent-Based Model with 1 million agents whosestochastic activity is given by Eq. (1) (shaded area = therange of 3 independent simulations). rate, quantified by parameter h ( t ). Due to the long-termrelaxation of h ( t ) this feedback loop limits the scale ofa single epidemic wave, but does not provide long-termprotection against new ones.As demonstrated below, the theory described by Eqs.(3)-(5) is in excellent agreement with simulations of theAgent-Based Model (ABM) in which social activities of 1million agents undergo stochastic evolution described byEq. (1) (compare solid lines with shaded areas in Figs.3, 4).Figure 3 illustrates a dramatic effect time-dependentheterogeneity has on the epidemic dynamics. It com-pares three cases: the classical homogeneous SIR model(black), the same model with persistent heterogeneity(brown), and the dynamic heterogeneity case considered Time / days M ( t ) R a D a il y i n c i d e n c e p e r b Early mitigationDelayed mitigationABM, early mitigationABM, delayed mitigation0 50 100 150 200 250 300 350
Time / days A tt a c k r a t e c FIG. 4. The time course of an epidemic with enhanced mitiga-tion during the first wave. (a) shows the M ( t ) R progressionfor two different strategies. In both cases, the enhanced mit-igation leads to a 50% reduction of M ( t ) R from 2 to 1. Inthe first scenario (early mitigation, blue curves), the reduc-tion lasted for only 15 days starting from day 27. In the sec-ond scenario (delayed mitigation, red curves), the mitigationwas applied on day 37 and lasted for 45 days. (b)-(c) showdaily incidence and cumulative attack rates for both strate-gies. As predicted, differences in the initial mitigation hadno significant effect on the epidemic in the long run: the twotrajectories eventually converge towards the universal attrac-tor. However, the early mitigation allows to suppress the peakof the infection, potentially reducing the stress on healthcaresystem. A delayed mitigation gives rose to a sizable secondwave. in this study (green). The latter two models share thesame HIT (green dashed line) which is reduced comparedto the homogeneous case (black dashed line). In the ab-sence of dynamic heterogeneity (black and brown) theinitial exponential growth halts once the respective HITis reached, but the overall attack rate “overshoots” be-yond that point, eventually reaching a significantly largerlevel, known as the final size of the epidemic (FSE). Im-portantly, in both these cases the epidemic has only asingle wave of duration set by the mean generation in-terval τ g multiplied by a certain R -dependent factor. Inthe case of dynamic heterogeneity (green), described byEqs. (3)-(5), the epidemic is transiently suppressed atthe level which is below even the heterogeneous HIT. As we argued in Ref. [8] this temporary suppression is due tothe population reaching the Transient Collective Immu-nity (TCI). That state originates due to the short-termpopulation heterogeneity being enhanced compared to itspersistent level. Stochastic contributions to social activ-ity responsible for this enhancement eventually averageout, leading to a slow degradation of the TCI state. Fig.3b illustrates that as the TCI state degrades, the dailyincidence rate develops an extended plateau on the greencurve. The cumulative attack rate shown in Fig. 3c re-laxes towards the HIT. This relaxation is characterizedby an emergent long time constant ˜ τ (cid:39) τ s /k > τ s .According to Eqs. (3-5) for a fixed mitigation level M ( t ), any epidemic trajectory would eventually convergeto the same curve, i.e. the universal attractor. The ex-istence of the universal attractor is apparent in Fig 4,where we compare two scenarios with different mitiga-tion strategies applied at early stages of the epidemic. Inboth cases, an enhanced mitigation was imposed leadingto a reduction of M ( t ) R by 50% from 2 to 1. In thefirst scenario (blue curves), the enhanced mitigation wasimposed on day 27 and lasted for 15 days. In the secondscenario (red curves), the mitigation was applied on day37 and lasted for 45 days. As predicted, this differencein mitigation has not had any significant effect on theepidemic in the long run: these two trajectories even-tually converged towards the universal attractor. How-ever, short- and medium- term effects were substantial.The early mitigation scenario (blue curve) resulted in asubstantial suppression of the maximum incidence dur-ing the first wave. Immediately following the release ofthe mitigation the second wave started and reached ap-proximately the same peak value as the first one. If theobjective of the intervention is to avoid the overflow ofthe healthcare system, this strategy would indeed help toachieve it. In contrast, the delayed mitigation scenario(red curve) turned out to be largely counterproductive.It did not suppress the peak of the first wave, but broughtthe infection to a very low level after it. Eventually, thatsuppression backfired as the TCI state deteriorated andthe epidemic resumed as a second wave, which is not asstrong as the first one.Since the late-stage evolution in our model is charac-terized by a long relaxation time ˜ τ , the possibility ofwaning of individual biological immunity or escape mu-tations of the pathogen accumulated over certain (pre-sumably, also long) time τ b , becomes a relevant effect.It can be incorporated as an additional relaxation term(1 − S ) /τ b in Eq. (4). The analysis of our equations,modified in this way, shows that the universal attractorleads to a fixed point corresponding to the endemic state.That point is located somewhat below the HIT and char-acterized by the finite residual incidence rate (1 − S ∞ ) /τ b and, respectively, by finite values of I and h . Here S ∞ is a susceptible population fraction in the endemic state,which is close to, but somewhat higher than that at theonset of the herd immunity. A similar endemic steadystate exists in most classical epidemic models (See [19]and references therein). However, in those cases, the epi-demic dynamics would not normally lead to that pointdue to the overshoot phenomenon. Instead, those mod-els typically predict a complete extinction of the diseasewhen the prevalence drops below one infected individual.This may happen before herd immunity is lost due towaning biological immunity and/or replenishment of thesusceptible population (e.g. due to births of immunolog-ically naive individuals). That is not the case when thetime-dependent heterogeneity is included.Note that for most pathogens the endemic point is notfixed but instead is subjected to periodic seasonal forc-ing in M ( t ). This leads to annual peaks and troughs inthe incidence rate. Our model is able to describe thisseasonal dynamics as well as transition towards it for anew pathogen (see Fig. 5). It captures the importantqualitative features of seasonal waves of real pathogens,e.g. three endemic coronavirus families studied in Ref.[42]. They are (i) sharp peaks followed by a prolongedrelaxation towards the annual minimum; (ii) a possibilityof multi-annual cycles due to parametric resonance. Time / years D a il y i n c i den c e pe r persistent heterogeneitydynamic heterogeneity h(t) J ( t ) s / k FIG. 5. Multi-year dynamics of a new pathogen. Effects ofwaning biological immunity with characteristic time τ b = 5yrs, and seasonal forcing are included (see SI for details). Inthe case of persistent heterogeneity without temporal varia-tions of social activity (brown solid line) the infection getsextinct following the initial wave of the epidemic. In con-trast, dynamic heterogeneity leads to an endemic state withstrong seasonal oscillations (green line). Insert: The epidemicdynamics in the ( J, h ) phase space. The black dotted line cor-responds to the universal attractor trajectory manifested e.g.as a plateau in green line in Fig. 3b. The attractor leads tothe endemic state (red point).
To understand the nature of the overall epidemic dy-namics, we focus on the behavior of variables J ( t ) and h ( t ). Their evolution is described by Eqs.(3) and (5),with R ∗ = R M ( t ) S ( t ) λ playing the role of a drivingforce. As a result of depletion of susceptible population,the driving force is gradually reduced, and the dynamicconverges towards a slow evolution along the universal attractor shown as a black dotted trajectory in ( h, J ) co-ordinates at the insert to Fig. 5. For initial conditionsaway from that trajectory (say, J ≈ h = 0), the linearstability analysis indicates that the epidemic dynamicshas a damped oscillatory behavior manifesting itself asa spiral-like relaxation towards the universal attractor.A combination of this spiral dynamics with a slow drifttowards the endemic state gives rise to the overall tra-jectory shown as the solid green line at the insert to Fig.5. The periodic seasonal forcing generates a limit cycleabout the endemic point (small green ellipse around thered point).More generally, any abrupt increase of the effective re-production number e.g. due to a relaxed mitigation, sea-sonal changes, etc. would shift the endemic fixed pointup along the universal attractor. According to Eqs. (3-5)this will once again trigger a spiral-like relaxation. It willmanifest itself as a new wave of the epidemic, such as thesecondary waves in Fig. 4b).A systematic calibration of our model and applica-tion to real-world epidemiological data is beyond thescope of this study. However, below we present a proof-of-principle demonstration that the progression in theCOVID-19 epidemic in the summer and fall of 2020 infour regions of the continental US: South, Northeast,Midwest and West can be well described by our theory.The time dependence of daily deaths per capita (a re-liable, albeit delayed measure proportional to the trueattack rate) is shown in Fig. 6bc for each of the regionsand fitted by our model with k = 0 . τ s = 30 days, κ = 2, together with IFR assumed to be 0 . M ( t ) follows regu-lar seasonal dynamics during that time period. This isreflected in a simple mitigation profile R M ( t ) shown inFig. 6a featuring a sharp relaxation of mitigation in earlysummer 2020 and gradual seasonal increase of the repro-duction number during the fall. Release of mitigationtriggered an immediate second wave in the South andWest (Fig. 6c) that had relatively low exposure duringthe first wave of the epidemic. In our model the waveis transiently suppressed when these regions reach theTCI state. Differences between the data and our modelobserved in the early fall may be tentatively attributedto government-imposed mitigation measures and/or toknowledge-based adaptation of the human behavior [15–17]. The late fall wave in all regions was triggered by sea-sonal changes in transmission. According to our modelthis wave was stabilized in mid-winter due to the popu-lation once again reaching the TCI state. Note a goodagreement between peak levels of this wave with our pre-dictions.Midwest and Northeast (Fig. 6b) exhibited similar be-havior except having a plateau instead of the summer R M ( t ) Northeast USMidwest USWest USSouth US
Jun 1 Jul 1 Aug 1 Sep 1 Oct 1 Nov 1 Dec 1 Dec 31
Date in 2020 D a il y dea t h s pe r , seasonalchange abc end oflockdowns FIG. 6. Fitting of the empirical data on COVID-19 epi-demic in Northeast (green), Midwest (blue), West (purple)and South (orange) of the USA. (a) The best-fit profiles R M ( t ), under the assumption that they are shaped by theend of lockdowns in the early summer of 2020, followed bygradual seasonal changes in the fall. Time dependence ofdaily deaths per capita for the Northeast and Midwest aswell as for South and and West US are shown in panels (b)and (c), respectively. Data points represent reported dailydeaths per 100,000 of population for each of the regions [13].Solid lines are the theoretical fits with our model. The follow-ing parameters of M ( t ) R curves were subject to variation:the mitigation level M ( t ) before and after the termination ofthe lockdown, the amplitude of the seasonal forcing, and thesummer-winter crossover time. Other parameters have beenfixed: k = 0 . κ = 2, τ s = 30 days, IF R = 0 . τ b = 5years. wave due to their higher levels of exposure during thefirst wave of the epidemic in spring 2020. Importantly,the transmission monotonically increases throughout the entire time window, yet our model captures the observedsecondary waves and plateaus. That behavior would beimpossible to explain using traditional epidemiologicalmodels well below the herd immunity.In conclusion, we have proposed a new theory integrat-ing the stochastic dynamics of individual social activityinto traditional epidemiological models. Our model de-scribes the so-called “zero intelligence” limit in whichthere is no feedback from the epidemic dynamics to so-cial activity e.g. mediated by the news. Hence our ap-proach is complementary to knowledge-based models ofRefs. [15–17]. The stochastic social activity in our ap-proach is described by the CIR model [39] which cap-tures the following important properties: (i) the activ-ity cannot be negative; (ii) for any given individual itreverses towards its long-term average value; (iii) it ex-hibits gamma-distributed short-term overdispersion (akasuperspreading) [1, 3, 4]. We mapped the overall epi-demic dynamics featuring heterogeneous time-varying so-cial activity onto a system of three differential equations,two of which generalize the traditional SIR model. Thethird equation describes the dynamics of the heterogene-ity parameter h ( t ), driven up by the current strength ofinfection J ( t ) and relaxing back to zero due to variablesocial activity.The emergent property of our theory is the new longtimescale τ s /k governing the relaxation towards eitherthe herd immunity or the endemic state of the pathogen.This behavior is in striking contrast to traditional epi-demiological models generally characterized by a largeovershoot beyond the herd immunity threshold leadingto a likely extinction of new pathogens. Our theory is ina good agreement with the empirical observation of longplateaus observed for many real-life epidemics includingCOVID-19 [13]. Dynamic heterogeneity also leads to atransient suppression of individual waves of the epidemicwithout reaching the long-term herd immunity [8]. Fi-nally, we demonstrated that our theory is in quantitativeagreement with the data describing secondary waves ofCOVID-19 epidemic in different regions of the USA. ACKNOWLEDGMENTS
This work was supported by the University of IllinoisSystem Office, the Office of the Vice-Chancellor for Re-search and Innovation, the Grainger College of Engineer-ing, and the Department of Physics at the University ofIllinois at Urbana-Champaign. This research was par-tially done at, and used resources of the Center for Func-tional Nanomaterials, which is a U.S. DOE Office of Sci-ence Facility, at Brookhaven National Laboratory underContract No. DE-SC0012704. [1] J. O. Lloyd-Smith, S. J. Schreiber, P. E. Kopp, W. M.Getz,
Nature , 355 (2005).[2] A. P. Galvani, R. M. May,
Nature , 293 (2005).[3] A. Endo, S. Abbott, A. Kucharski, S. Funk,
WellcomeOpen Research , 67 (2020).[4] K. Sun, et al. , Science , eabe2424 (2020).[5] R. Pastor-Satorras, C. Castellano, P. Van Mieghem,A. Vespignani,
Reviews of Modern Physics , 925(2015).[6] S. Bansal, B. T. Grenfell, L. A. Meyers, Journal of theRoyal Society Interface , 879 (2007).[7] M. G. M. Gomes, et al. , medRxiv (2020).[8] A. V. Tkachenko, et al. , MedRxiv (2020).[9] J. Neipel, J. Bauermann, S. Bo, T. Harmon, F. J¨ulicher,
PLoS one , e0239678 (2020).[10] T. Britton, F. Ball, P. Trapman, Science eabc6810 (2020).[11] C. Aschwanden,
Nature , 26 (2020).[12] L. F. Buss, et al. , Science , 288 (2021).[13] Data and visualization at http://91-divoc.com/pages/covid-visualization/ (2020).[14] S. Thurner, P. Klimek, R. Hanel,
Proceedings of the Na-tional Academy of Sciences , 22684 (2020).[15] J. S. Weitz, S. W. Park, C. Eksin, J. Dushoff,
Proceedingsof the National Academy of Sciences p. 202009911 (2020).[16] S. Funk, E. Gilad, C. Watkins, V. A. A. Jansen,
Pro-ceedings of the National Academy of Sciences , 6872(2009).[17] A. Rizzo, M. Frasca, M. Porfiri,
Phys. Rev. E , 042801(2014).[18] C. Rose, et al. , arXiv: (2020).[19] M. J. Keeling, P. Rohani, Modeling infectious diseases inhumans and animals (Princeton University Press, 2011).[20] N. D. Wolfe, C. P. Dunavan, J. Diamond,
Nature ,279 (2007).[21] A. Engering, L. Hogerwerf, J. Slingenbergh,
EmergingMicrobes & Infections , 1 (2013).[22] R. Pastor-Satorras, A. Vespignani, Phys. Rev. E ,066117 (2001).[23] M. Starnini, et al. , Social Informatics , G. L. Ciampaglia,A. Mashhadi, T. Yasseri, eds. (Springer InternationalPublishing, Cham, 2017), pp. 536–551.[24] E. Volz, L. A. Meyers,
Proceedings of the Royal SocietyB: Biological Sciences , 2925 (2007).[25] S. Bansal, J. Read, B. Pourbohloul, L. A. Meyers,
Jour-nal of biological dynamics , 478 (2010).[26] J. M. Read, K. T. D. Eames, W. J. Edmunds, Journal ofthe Royal Society Interface , 1001 (2008).[27] A. L. Lloyd, R. M. May, Science , 1316 (2001).[28] R. M. May, A. L. Lloyd,
Physical Review E , 066112(2001).[29] Y. Moreno, R. Pastor-Satorras, A. Vespignani, EuropeanPhysical Journal B , 521 (2002).[30] N. Perra, B. Gon¸calves, R. Pastor-Satorras, A. Vespig-nani, Scientific Reports (2012).[31] A. Vazquez, B. R´acz, A. Luk´acs, A.-L. Barab´asi, Phys.Rev. Lett. , 158702 (2007).[32] R. Pastor-Satorras, A. Vespignani, Physical Review Let-ters , 3200 (2001).[33] M. Salathe, et al. , Proceedings of the National Academy of Sciences , 22020 (2010).[34] L. Isella, et al. , Journal of Theoretical Biology , 166(2011).[35] D. Rybski, S. V. Buldyrev, S. Havlin, F. Liljeros, H. A.Makse,
Proceedings of the National Academy of Sciences , 12640 (2009).[36] A.-L. Barabasi,
Nature , 207 (2005).[37] J. Saram¨aki, E. Moro,
The European Physical Journal B , 164 (2015).[38] B. F. Nielsen, K. Sneppen, L. Simonsen, J. Mathiesen, medRxiv (2020).[39] J. C. Cox, J. E. Ingersoll, S. A. Ross, Econometrica ,385 (1985).[40] W. Feller, The Annals of Mathematics , 173 (1951).[41] R. Aguas, et al. , medRxiv (2020).[42] R. A. Neher, R. Dyrdak, V. Druelle, E. B. Hodcroft,J. Albert, Swiss Medical Weekly (2020).[43] S. Anand, et al. , The Lancet , 1335 (2020).[44] F. J. Angulo, L. Finelli, D. L. Swerdlow,
JAMA NetworkOpen , e2033706 (2021). SUPPLEMENTARY MATERIALSA. Epidemic dynamics with dynamic heterogeneity
Let a i ( t ) be a measure of individual’s social activity proportional to frequency and intensity of close contacts withother people around time t . We refer to it as (social) susceptibility to infection, but it also determines one’s potentialto infect others. In particular, infectivity of a person i infected at time t ∗ i is given by β i ( t ∗ i + τ ) = C i ( τ ) a i ( t ∗ + τ ) (S1)Here C i ( τ ) is person’s contagiousness level, at time τ after infection.Let j ( t ) be the fraction of infected individual, weighted proportionally to their current infectivity level, and M ( t )be mitigation factor that reflects government and social response to epidemics, seasonal effects etc. Their product, J ( t ) = M ( t ) j ( t ) is force of infection, i.e. a hypothetical incidence rate in fully susceptible homogeneous populationwith α = 1. Within the heterogeneous (but well-mixed) age-of-infection model, current value of j ( t ) is given by j ( t ) = (cid:104) β i ( t − t ∗ i ) (cid:105) i = (cid:90) ∞ (cid:104) C i ( τ ) a i ( t ) a i ( t − τ ) S i ( t − τ ) (cid:105) i J ( t − τ ) dτ (S2)Here S i ( t − τ ) is the state of an individual i (1 if susceptible, 0 otherwise). Since J ( t ) is by definition proportional to j ( t ), we obtain the quasi-homogeneous renewal equation: j ( t ) = (cid:90) ∞ K ( t, τ ) R e ( t − τ ) j ( t − τ ) dτ (S3)Here effective reproduction number R e and the probability density of the generation interval K ( τ ), are given by R e ( t ) = M ( t ) (cid:28) S i ( t ) (cid:90) ∞ a i ( t ) a i ( t + τ ) C i ( τ ) dτ (cid:29) i (S4) K ( t, τ ) = (cid:104) S i ( t ) a i ( t ) a i ( t + τ ) C i ( τ ) (cid:105) i (cid:10) S i ( t ) (cid:82) ∞ a i ( t ) a i ( t + τ ) C i ( τ ) dτ (cid:11) i (S5) B. Stochastic model for social activity
It is well known that social interactions are “bursty”. That is to say, individual social activity has both (nearly)permanent and significant time-dependent contributions: a i ( t ) = α i + δa i ( t ) (S6)Without loss of generality we set the population-averaged permanent and instantaneous susceptibility to 1: (cid:104) a i ( t ) (cid:105) i = (cid:104) α i (cid:105) i = 1. Beyond its average value, the overall statistics of instantaneous α ( t ) is properly defined only if that quantityis average over specified time window δt . Naturally, its variation will gradually decrease as the time widow increases.Individual reproductive number, R i , for COVID-19 epidemic is (in)famously over-dispersed. This is a result of super-spreading, when a majority of secondary infections are caused by a small fraction of index cases. The overdispersionreflects (i) variation of peak contagiousness level among individuals and (ii) dispersion of a i ( t ) which is effectivelyaveraged over a timescale of the peak infection period (approximately 2 days).Importantly, according to Eq.(S4), the reproductive number depends on correlations of a i across a time scale of asingle generation interval (on average, 4 to 5 days for COVID 19). Thus, any variations in a i ( t ) that do not persistover that timescale would be averaged out. Here we introduce a simple model to account for temporal variation ofsocial activity. In this model, a i may vary on short time scale, relax to the persistent value for a given individual overcertain relaxation time, τ s : ˙ a i = α i − a i τ s + η i ( t ) (S7)Here η ( t ) is short-time noise that gives rise to time dependent fluctuations. We set (cid:104) η i ( t ) η i ( t (cid:48) ) (cid:105) = a i ( t ) τ s k δ ( t − t (cid:48) ), whichgives rise to individual diffusion in a i space with diffusion coefficient proportional to a i . The evolution of populationwith a given value of persistent activity α in that space is given by the following Fokker-Plank Equation:˙Ψ α ( a, t ) = 1 k τ s ∂ ( a Ψ( a, t )) ∂a + 1 τ s ∂ (( a − α )Ψ α ( a, t )) ∂a (S8)The steady state solution to this equation gives a probability density function (pdf) for a , which turns out to be acommonly used gamma distribution: Ψ α ( a, t ) = f α ( a ) = a αk − e − k a α αk Γ( αk ) (S9)Note that the statistics of superspreading events is commonly modeled assuming the very same distribution forindividual reproduction number, R i . This gives a strong empirical support to the chosen model, in particular to thechoice of diffusion coefficient to be proportional to α . It also allows us to partially calibrate the model. Reporteddispersion parameter associated with superspreading events for COVID 19 is in the range of 0 . . k is expected to be larger than k , i.e. has a smaller dispersion. This is because variationsof a ( t ) over the timescale shorter than a single generation interval would be averaged out according to Eq. (S10),while the superspreading statistics effectively probes it over a shorter time interval of the infectivity peak in a singleindividual. The latter could be further enhanced by a variation of individual contagiousness e.g. due to biologicalfactors.It is well known that the mean reproduction number R in a heterogeneous population depends on the secondmoment of distribution of α (in network epidemic models it is related to individual degree). However, there is animportant modification to that result for time-dependent a ( t ): R = (cid:90) ∞ (cid:104) a i ( t ) β i ( t + τ ) (cid:105) i dτ = R (cid:104) α i (cid:105) + (cid:90) ∞ (cid:104) C i ( τ ) δa i ( t ) δa i ( t + τ ) (cid:105) i dτ (S10)Here R = (cid:104) (cid:82) C i ( τ ) dτ (cid:105) i is the net infection transmission probability of an average person. This gives: R = (cid:104) α i (cid:105) + µk − (S11)Here factor µ is related to the Laplace transform of average contagiousness profile, K ( τ ) = (cid:104) C i ( τ ) (cid:105) i /R . µ = (cid:90) ∞ K ( τ ) e − τ/τ s dτ (S12)Note that, according to Eq.(S5), generation interval pdf K ( τ ) is close, but not identical to K ( τ ): K ( τ ) = (cid:18) e − τ/τ s − µk (cid:104) α i (cid:105) + µ (cid:19) K ( τ ) (S13)Specifically, for the case of SIR model, i.e. K ( τ ) ∼ e − τ/τ , one obtains: µ = 1 / (1 + τ /τ s ) ≈ / (1 + τ g /τ s ). Heremean generation interval τ g is given by τ g = τ ( (cid:104) α i (cid:105) + µ k − ) (cid:104) α i (cid:105) + µk − ≈ τ (cid:18) − τ τ s (1 + k (cid:104) α i (cid:105) ) (cid:19) (S14)In this SIR case, one can assign each person a state variable I i set to 1 when the individual is infectious and 0otherwise. This allows to describe the epidemic dynamics in terms of activity-weighted fraction of the infectedpopulation, I ( t ) = (cid:104) I i a i ( t ) (cid:105) / (cid:104) a i (cid:105) . Note that variable j ( t ) and hence the strength of infection are proportional to it: J ( t ) = M ( t ) j ( t ) = R M ( t ) I ( t ) τ g (S15) C. Mapping on quasi-homogeneous dynamic system
Let S α ( a, t ) be the fraction of susceptibles among the sub-population with persistent activity level α and givenvalue of a , at time t . Change of function S α ( a, t ) is driven by two effects: (i) depletion of susceptible population dueto infection and (ii) diffusion of individual in α space. By substituting Φ α ( a, t ) = f α ( a ) S α ( t ) into Fokker-Plank Eq.(S8), and adding the infection term with rate − a ( t ), we obtain an evolution equation for S α ;˙ S α ( a, t ) = − aS α ( a, t ) J ( t ) + ak τ s ∂ S α ( a, t ) ∂a + (cid:18) α − aτ s (cid:19) ∂S α ( a, t ) ∂a (S16)This equation can be solved by using the following ansatz: S α ( a, t ) = exp [ − Z ( t ) α − k h ( t ) a ] (S17)Here Z ( t ) is a measure of persistent heterogeneity: the larger it is, the more is the difference in depletion of susceptibleamong subpopulations with different α ’s, i.e. various average levels of social activity. On the other hand, h ( t ) param-eterizes the transient heterogeneity within each of these subpopulations. In the long run, this type of heterogeneitydisappears due to evolution in a -space, thus h ( t ) asymptotically approaches 0 as t → ∞ . Substituting Eq. (S17) intoEq. (S16) results in simple equations for both Z ( t ) and h ( t ):˙ h = J ( t ) k − h ( t )(1 + h ( t )) τ s (S18)˙ Z = k h ( t ) τ s (S19)The renewal equation Eq. (S3) for j ( t ) completes our quasi-homogeneous description of the epidemic dynamics.However, to fully close this system of equations, one needs to express the effective reproduction number, R e , in termsof the functions M ( t ), Z ( t ) and h ( t ). This is done by substituting the ansatz, Eq. (S17), into Eq. (S4). We performthis calculation in two steps, by first finding the effective number R α for a sub-population with average level of activity α , followed by averaging over persistent heterogeneity. This gives R α = (cid:90) ∞ a ( α + µ ( a − α )) f α ( a ) e − Z ( t ) α − k h ( t ) a da = αR (cid:0) α + µk − + h (1 − µ ) (cid:1) e − ˜ Zα (1 + h ( t )) (S20)Here ˜ Z = Z + k ln(1 + h ) (S21)Note that ˙˜ Z = J ( t )1 + h ( t ) (S22)The averaging over persistent heterogeneity, under the assumption that α obeys the gamma distribution, p ( α ) ∼ α κ − e − κα , yields R e ( t ) = M ( t ) (cid:90) ∞ R α p ( α ) dα = χ + (cid:0) − χ )(1 + k h ( µ − − (cid:1) (cid:16) κ − ˜ Z ( t ) (cid:17) R M ( t ) (cid:16) κ − ˜ Z ( t ) (cid:17) κ (1 + h ( t )) (S23)Here χ = 1 + κ − κ − + µk − (S24)Similarly, we calculate S , which ends up having the same form as in the model with persistent heterogeneity [8]: S ( t ) = (cid:90) ∞ (cid:90) ∞ p ( α ) f α ( a ) e − Z ( t ) α − k h ( t ) a dadα = 1 (cid:16) κ − ˜ Z ( t ) (cid:17) κ (S25)By comparing Eqs.(S23) and (S25) we obtain R e in terms of S and h : R e ( t ) = R M ( t ) S λ q χ ( S, h )(1 + h ( t )) (S26)Here q χ ( S, h ) = (1 − χ ) (cid:0) k h ( µ − − (cid:1) S − χ/κ + χS (1 − χ ) /κ ≈ λ = 1 + 1 + χκ = (cid:0) κ − (cid:1) (cid:0) µk − + 2 κ − (cid:1) µk − + κ − (S28)Note that for most practical purposes, one can set q χ ( S, h ) = 1. The parameter λ is the “immunity factor” thatemerged in the context of our earlier study of persistent heterogeneity [8]. In that case, λ = 1 + 2 /κ appears asa scaling exponent in relationship between effective reproduction number R e ( t ) and the fraction of the susceptiblepopulation S ( t ): R e = R M S λ . Eq. (S26) generalizes that result.Eqs.(S3), (S18), S26, (S22) give the full description of the epidemic dynamics in heterogeneous system. For aparticular case of SIR model ( K ( τ ) ∼ e − τ/τ g ), we obtain a 3D dynamical system, in terms of variables I ( t ), S ( t ) and h ( t ): dIdt = JS λ (1 + h ) − Iτ g (S29) dSdt = − JS /κ (1 + h ) (S30) dhdt = Jk − h (1 + h ) τ s (S31)Here, infection strength J ( t ) is proportional to the activity-weighted fraction of susceptible population I ( t ), and themitigation profile M ( t ), as given by Eq. (S15). Eq. (S30) was derived by combining Eq. (S22) and Eq. (S25).Alternatively, after substituting result of integration of Eq. (S22) into (S25), one gets the explicit formula for S ( t ): S ( t ) = (cid:18) κ − (cid:90) t −∞ J ( t (cid:48) ) dt (cid:48) h ( t (cid:48) ) (cid:19) − κ (S32) D. Waves and plateaus
According to Eq. (S29), the combined driving force of the epidemic is R ∗ = R M ( t ) S λ ( t ). It includes both theeffects of mitigation M ( t ) and suppression associated with the build up of the long-term herd immunity. First, weassume R ∗ to be fixed or change very slowly (adiabatically), i.e. on the timescales longer than τ s . In that case, J ( t )and h ( t ) trail the driving force R ∗ ( t ), staying close to the corresponding adiabatic fixed point ( J ∗ , h ∗ ) in their 2Dphase space: h ∗ = √ R ∗ − J ∗ = k h ∗ (1 + h ∗ ) τ s (S34)The stability of this adiabatic fixed point, and the more rapid epidemic dynamics can be described by linearizing Eqs.(S29, S31) and (S18) around ( J ∗ , h ∗ ), i.e. by assuming h ( t ) = h ∗ + δh ( t ) and J ( t ) = J ∗ + δJ ( t ): ddt (cid:18) δhδJ (cid:19) = 1 τ s (cid:18) − (1 + 2 h ∗ ) τ s /k − k γh ∗ (cid:19) (cid:18) δhδJ (cid:19) (S35)The eigenmodes of this linearized system are both stable, but the rates have substantial imaginary components: r ± = − h ∗ τ s ± i (cid:115) γh ∗ τ s − (1 + 2 h ∗ ) τ s (S36)This indicates that relaxation towards point ( J ∗ , h ∗ ) has a pronounced oscillatory character. The period of theoscillations is T ≈ π (cid:114) τ s γh ∗ ≈ π (cid:115) τ s γ ( √ R ∗ −
1) (S37)The amplitude of the oscillations decays with the time constant 2 τ s / (1+2 h ∗ ). This oscillatory behavior would manifestitself as multiple epidemic waves. In reality, the dynamics is more complicated since rapid changes of M ( t ), e.g. dueto seasonal effects, government and societal response to the epidemic, would additionally modulate it.The assumption of R ∗ = R M ( t ) S λ ( t ) being fixed is not, of course, realistic. In particular, the mitigation factor M ( t ) may have both slow and fast variations. On top of that, the dependence of R ∗ on S ( t ) creates a negativefeedback suppressing the forcing on the long run. For a constant mitigation M , there is a line of fixed points( J, S, h ) = (0 , S, S ≤ S = ( R M ) − /λ . Here 1 − S represents the long-term herd immunity threshold(HIT) for a given mitigation level M . There is one particular solution ( ˜ J ( t ) , ˜ S ( t ) , ˜ h ( t )) corresponding to all threevariables slowly evolving in such a way that R e stays close to 1 at all times, eventually reaching the HIT point,(0 , S ). As follows from the above stability analysis, this solution acts as an attractor, with any trajectory in ( J, S, h )space converging towards it, unless perturbed by variations in mitigation M ( t ). To construct that solution, we setthe growth rate for I ( t ) in Eq.(S29) to 0, and use Eqs (S30)-(S31) to calculate the corresponding evolution of h ( t ): (cid:18) τ s + ˜ τ − τ s (1 + ˜ h ( t )) ν (cid:19) ˙˜ h = ˜ h ( t )(1 + ˜ h ( t )) (S38)˜ h ( t ) ≈ ν )( e ( t − t ) / ˜ τ −
1) (S39)Here ν = λκ , and ˜ τ = τ s (cid:18) ν ( R M ) ν k (cid:19) (cid:39) τ s k (S40)Remarkably, under assumption of strong overdispersion, k (cid:28)
1, the emergent timescale ˜ τ is significantly longerthan the social rewiring time, τ s . This long timescale corresponds to a slow process of individuals trapped in the lowactivity state, a ( t ) ≤ k , transitioning to the high activity level a ≥
1. Respective evolution of ˜ S ( t ) and ˜ J ( t ) are givenby: ˜ S ( t ) = S (cid:16) h ( t ) (cid:17) /λ (S41)˜ J ( t ) ≈ (cid:18) τ s − τ (cid:19) k ˜ h ( t )(1 + ˜ h ( t )) (S42) E. Waning of biological immunity
Our equations could be easily modified to account for the waning of biological immunity. This adds a new term inEq. (S30) which becomes: dSdt = − JS /κ (1 + h ) + 1 − Sτ b (S43)Here τ b is the lifetime of biological immunity, which we set to 5yrs throughout this work. The last term τ b (1 − S )describes the rate at which the recovered population (fraction 1 − S ) reverts back to the susceptible state. The endemicsteady state can be found by setting time derivatives Eqs. (S29),(S31) and (S43), to 0. Under the assumption that τ b (cid:29) ˜ tau , the endemic point in ( S, J, h ) is given by J en ≈ − S HIT τ b S /κHIT (S44) h en ≈ ˜ τ J en = ˜ ττ b − S HIT S /κHIT (S45) S en = S HIT (1 + h en ) /λ ≈ S HIT (S46)Here S HIT = R M /λ corresponds to the HIT. F. Seasonal forcing
Seasonal effects are commonly described as simple sin-shaped modulation of reproductive number [42]. In this work,we used a combination of sigmoidal functions to model transition between “winter” and “summer” values of M ( t ): M s ( t ) = 1 + σ ∞ (cid:88) n =0 (cid:20) − tanh (cid:18) t − t spring + nT ∆ (cid:19) + tanh (cid:18) t − t fall + nT ∆ (cid:19)(cid:21) (S47)Here T = 1yr, time parameters t spring < t fall and ∆ determine the timing of and sharpness of winter-summer-wintertransitions. σ determines the amplitude of seasonal changes. In particular, σ = 0 .
25 in Fig.2, and ranges between0 .
25 and 0 .
35 in our fits of epidemic dynamics for different US regions, Fig. 6.
G. Implementation of the agent-based model
All simulations for the agent-based model use 1 million agents and 3 simulation replicates. For each agent in thesimulation, at each time step, the social activity follows the stochastic dynamics described in Eq. 1. After that, theoverall force of infection is computed using J ( t ) = R M ( t ) τ g (cid:104) a ( t ) (cid:105) i N (cid:88) i a i I i , (S48)where I i is binary and used to denote whether or not the agent is infectious, N is the number of agents in the simulation.For a susceptible agent i , the chance of being infected in one simulation step is a i ( t ) J ( t ) dt which is proportional to theforce of infection, his/her activity a i ( t ), and dt - the length of the time step used in our simulations. For an infectiousagent, the probability of recovering from the infectious state in one simulation step is γdt . When the waning ofbiological immunity is ignored, recovered agents will always stay in the recovered state and cannot be infected again. TABLE I. Set of fixed model parameters k κ τ g τ s τ b ..