key: cord-0431412-x7j60uyl
authors: Browne, Christopher A.; Amchin, Daniel B.; Schneider, Joanna; Datta, Sujit S.
title: Infection percolation: A dynamic network model of disease spreading
date: 2021-03-17
journal: nan
DOI: 10.3389/fphy.2021.645954
sha: 14de74a8247fbafc49a94f31aff22c1df34534b3
doc_id: 431412
cord_uid: x7j60uyl

Models of disease spreading are critical for predicting infection growth in a population and evaluating public health policies. However, standard models typically represent the dynamics of disease transmission between individuals using macroscopic parameters that do not accurately represent person-to-person variability. To address this issue, we present a dynamic network model that provides a straightforward way to incorporate both disease transmission dynamics at the individual scale as well as the full spatiotemporal history of infection at the population scale. We find that disease spreads through a social network as a traveling wave of infection, followed by a traveling wave of recovery, with the onset and dynamics of spreading determined by the interplay between disease transmission and recovery. We use these insights to develop a scaling theory that predicts the dynamics of infection for diverse diseases and populations. Furthermore, we show how spatial heterogeneities in susceptibility to infection can either exacerbate or quell the spread of disease, depending on its infectivity. Ultimately, our dynamic network approach provides a simple way to model disease spreading that unifies previous findings and can be generalized to diverse diseases, containment strategies, seasonal conditions, and community structures.

Epidemic spreads-such as the 1918 flu, HIV/AIDS, and COVID-19 pandemics-highlight the critical importance of infectious disease modeling in our everyday lives. Mathematical models can provide key insights into the process by which disease is transmitted between individuals, help to forecast how a disease will continue to spread through a population, and assess the efficacy of different interventions. Developing accurate, computationally-tractable, and generally-applicable models of disease spreading is therefore of critical importance to public health.

The dynamics of infectious disease spreading are often modeled using compartmental models [1, 2] that employ the framework of reaction kinetics. For example, in one representation, members of a population are divided into three intermixed groups: those who are susceptible to infection (S), currently infected (I), or recovered from infection (R), with transitions from S → I and I → R occurring at prescribed rates. This standard model, known as the SIR model or the SI model in the absence of recovery, can successfully predict the initial exponential and eventual logistic dynamics of the total amount of infection in a population for many diseases [3] [4] [5] [6] . Thus, the SIR model provides a useful approach for disease modeling that is well established.

However, a notable omission of this model is that it does not consider the discrete, spatially-separated interactions between members of a population. Instead, it assumes a well-mixed population for which the crucial dynamics of disease transmission between individu- † These authors contributed equally to this work. * To whom correspondence should be addressed: ssdatta@princeton.edu als are lumped into the basic reproduction number R 0 , a macroscopic parameter describing the mean number of secondary transmissions from each infection. In practice, this quantity is used as a fitting parameter-limiting projection capabilities when the interactions between individuals change even slightly [7] [8] [9] [10] [11] [12] . A recent example highlighting this deficiency is the spread of COVID-19 in China: containment policies imposing spatial barriers and suppressing individual interactions are thought to have hindered exponential growth of infection, yet this pivotal effect cannot be captured by the classic SIR model without invoking additional fitting parameters [13] .

To address this issue, sophisticated extensions of this model have been developed to explicitly incorporate spatial variations in spreading [2, . For example, in one approach, different subpopulations-each modeled using different SIR dynamics-are coupled together [2, [36] [37] [38] . While powerful, this approach still employs lumped parameters for each subpopulation that do not explicitly consider differences between discrete individuals. Percolation theory provides another powerful way to explore disease spreading throughout a social network composed of discrete, spatially-separated individuals, yielding insights into the onset of disease spreading and the size of disease outbreaks [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] 39] . Nevertheless, models based on this approach also suffer from limitations: they typically either do not consider the dynamics of disease spreading and only treat the final stage of infection, or they also describe the dynamics of disease transmission using an ad hoc macroscopic parameter that aggregates the influence of random and uncorrelated individual interactions. However, transmission is known to depend sensitively on the full history of infection, on specific individual behaviors [40] including social distancing [13, 23, 24, 32, 41, 42] , and on interactions between different social networks [39] . Thus, the ability to accurately predict the temporal evolution of active infections in an overall population, as well as the full spatiotemporal features of disease spreading, remains limited.

Here, we build on this previous work to develop an infection percolation framework to model infectious disease spreading by applying SIR dynamics to discrete interactions between members of a population. Our framework explicitly considers the full spatiotemporal history of infection as well as individual interactions in describing the spatial variation and individual dynamics of disease transmission. Simulations employing this framework show that disease spreads through a social network as a traveling wave of infection, followed by a traveling wave of recovery-consistent with previous predictions obtained using analytical SIR modeling [19, 43, 44] . Analysis of these waves reveals general features of the total number of infections, maximal number of active infections, and the temporal evolution of active infections in a population, and clarifies how these features are determined by the interplay between disease transmission and recovery-consistent with the results of previous percolation simulations [45] . Our framework therefore provides a simple way to unify different findings previously obtained using disparate modeling approaches, and helps to clarify how disease spreading manifests for different diseases and containment strategies. Finally, as an example of how our framework can help to go beyond typical models of well-mixed populations that do not incorporate possible correlations in individual behavior, we demonstrate how disease spreading is strongly altered in a spatially heterogeneous social network.

Our approach is inspired by dynamic network modeling of fluid-driven transport in heterogeneous media, which similarly seeks to predict spatiotemporal features of spreading in complex settings [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] . We represent members of a population, all of whom are initially susceptible to infection, by sites of a network that are connected by bonds representing pairwise interactions between them (Methods). We describe the intrinsic infectivity of a given disease throughout this network by the parameter P * , where a high value of P * characterizes a highly infectious disease. Crucially, each interaction between individuals ij is characterized by its own barrier to disease transmission P th,ij [58] ; this threshold explicitly describes the propensity of individual i, if infected, to transmit the disease to the susceptible individual j, and depends on individual behaviors such as social distancing between i and j [13, 23, 32, 42] . Furthermore, if individuals i and j are infected and susceptible, respectively, the duration of disease transmission from i to j is given by ∆τ ij ; at the population scale, the mean of all the ∆τ ij can be thought of as the inverse of a macroscopic disease transmission rate. However, in what follows, we will describe our results in terms of the time durations, since our model utilizes the time duration, not rate, in calculations of different transmission events between individuals.

Previous studies indicate that disease is transmitted between individuals more rapidly, with shorter ∆τ ij , for stronger infections P * and reduced barriers to transmission P th,ij [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] . We therefore propose the ansatz that disease is only transmitted if P th,ij < P * , with a transmission time ∆τ ij = τ 0 /f (P * − P th,ij ), where f monotonically increases from 0 to 1 and τ 0 is a characteristic minimal transmission time. This feature contrasts with previous percolation-based models, which often assume that the individual disease transmission times and barriers to disease transmission are either constant or drawn from two independent distributions. Instead, by explicit linking the individual ∆τ ij and P th,ij , our work provides a way to analyze the competing roles of individual disease infectivity and susceptibility in influencing populationscale disease spreading, as we show below.

Finally, as in the classic SIR model, we introduce the recovery duration τ r,i after which individual i, if infected, transitions to a recovered state and can no longer infect other individuals. The mean of this quantity can be thought of as the inverse of a macroscopic infection recovery rate. However, in what follows, we will describe our results in terms of the recovery duration-again, because our model utilizes the time duration, not rate, in calculations of different recovery events.

In nondimensional form, these control parametersthe infectivity, individual barriers to infection, disease transmission times, and recovery duration-are then represented by the parametersP * ≡ P * /P th,max ,P th,ij ≡ P th,ij /P th,max , ∆τ ij ≡ ∆τ ij /τ 0 = 1/f (1 −P th,ij /P * ), andτ r,i ≡ τ r,i /τ 0 , respectively, where P th,max is the maximal value of P th,ij in the population. Thus, this framework parses the role of a macroscopic transmission rate or reproduction number R 0 into appropriate parameters that are globally constant for a particular disease and parameters that vary by individual, providing a systematic way to incorporate measurable spreading parameters to model populations with spatial heterogeneities.

As a first step toward exploring the spatiotemporal features of disease spreading in this framework, we construct two-dimensional (2D) simulations implementing these deterministic rules. We represent the social network as a static square lattice with N t = 10 4 sites, though exploring more complex networks with small world and scale-free features [20, [70] [71] [72] will be an important direction for future work. The disease is introduced at the central site at timeτ ≡ τ /τ 0 = 0. For simplicity, we take the barriers to disease transmissioñ P th,ij to be undirected, withP th,ij =P th,ji ; however, directed transmission is known to dramatically change disease spreading in some networks, and will be useful to explore in future implementations [40] . TheP th,ij have randomly assigned values that are chosen from a uni-form distribution-and kept fixed for all computations performed for a given distribution-though we later test other distributions as well. We describe the individual dynamics of disease transmission using the simplest linear function, f = 1 −P th,ij /P * , and take recovery to be an intrinsic property of a given disease, withτ r,i set to a constant valueτ r throughout the lattice network.

Infection percolation in a recovery-free population. We first investigate the recovery-free case with τ r → ∞. For a disease with low infectivity,P * 1, only a small subset of individual interactions lead to transmission; thus, a large portion of the population becomes inaccessible to infection, and disease spreading is localized (Movie S1)-reminiscent of subcritical bond percolation. For sufficiently high infectivity, however, a sufficient number of individual interactions permit transmission for the disease to spread throughout the populationreminiscent of supercritical bond percolation. The example ofP * = 0.6 (Movie S2) is shown at an intermediate timeτ = 100 in Figure 1a . For this infectivity, the disease spreads in a spatially heterogeneous, ramified pattern, reminiscent of supercritical bond percolation near the percolation threshold. This heterogeneous spreading leads to the formation of discrete clusters of bypassed individuals who remain uninfected (white in Fig.  1a ). For diseases of higher infectivity, disease spreading is more compact: the leading front of the infected population becomes more smooth, resulting in fewer and smaller uninfected clusters, as shown for the example of P * = 0.7 in Fig. 1b and Movie S3. Clearly, disease infectivity strongly influences the spatiotemporal features of spreading.

To gain further insight into disease spreading, we repeat these simulations for a broad range ofP * , and characterize the infection growth dynamics by measuring the time-dependent infected fraction of the population, φ. We again observe a bifurcation in infection behavior. For low infectivity, the disease spreads slowly (light curves in Fig. 1c) , ultimately only infecting a total fraction φ t 1 (points withP * < 0.5 in Fig. 1d ). By contrast, for sufficiently high infectivity, the disease spreads rapidly with φ ∼τ 2 (dark curves in Fig. 1c ), ultimately infecting nearly the entire population (points forP * > 0.5 in Fig.  1d ). The abrupt onset of rapid spreading throughout the population at a critical infectivityP c,0 ≈ 0.5 again suggests that disease spreading can be described as a dynamic process of percolation through infected bonds, consistent with previous calculations [17, 18, [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] 39] . This suggestion is further confirmed by the value ofP c,0 , which coincides with the critical probability of bond percolation on the 2D square lattice [34, 73, 74] . We find no discernible difference in these results when the simulations are performed with lattice networks that are up to two orders of magnitude larger (Fig. S1 ), suggesting that our results are not strongly sensitive to finite-size effects.

Why does infection growth show quadratic scaling in time forP * >P c,0 ? At the leading front of the infected population, new individuals are infected over a range of transmission times ∆τ ij = 1/(1 −P th,ij /P * ) ∈ [1, ∞), leading to heterogeneous disease spreading. However, as P * increases, a greater proportion of disease transmission between individuals occurs in the shortest possible time τ 0 , corresponding to ∆τ ij = 1, resulting in more compact spreading (Fig. 1b) . Hence, the leading front of the infected population spreads radially outward on the 2D social network lattice at a maximal rate of 1 new individual per τ 0 , therefore spanning an overall region with a maximal radius ofτ individuals. The maximal infected fraction of the population is then given directly by the area of this infected region: φ ≈ πτ 2 /N t , consistent with a previous result in percolation theory [75] and yielding the quadratic scaling shown in Fig. 1c .

As this infected region spreads, it eventually reaches the boundary and spans the entire population at a shortest possible time ofτ f,0 ≈ N t /2 = 71 (SI text)-in good agreement with the onset of the plateau in φ atτ ≈ 100 for the highestP * shown in Fig. 1c . AsP * decreases, we expect that the infected region spreads at a slower rate, as a greater proportion of disease transmission occurs at ∆τ ij > 1-also in good agreement with the variation of the curves in Fig. 1c . The variability in these curves reflects the increasing importance of the variability in the individual barriers; different simulations employing different choices of the randomly-chosenP th,ij would yield slightly different dynamics. However, in all cases, we expect that decreasingP * leads to slower transmission as a greater proportion of disease transmission occurs at ∆τ ij > 1.

Infection percolation in a population with recovery. How do these results change when infected individuals can recover? To address this question, we perform the same simulations as in Fig. 1a -d, but withτ r = 4. This modification markedly alters disease spreading. For P * = 0.6, the spread of disease is quickly quenched by recovery [31] , and only a few individuals are ever infected (Fig. 1e , Movie S4). Thus, even above the critical infectivity for percolation in the recovery-free case,P c,0 ≈ 0.5, recovery gives rise to subcritical spreading behavior. For a higher infectivity ofP * = 0.7, the disease does continue to spread, eventually reaching the boundaries of the population as in supercritical bond percolation. However, unlike the recovery-free case, it does so in a spatially heterogeneous, ramified pattern ( Fig. 1f , Movie S5). Close inspection of the spatiotemporal pattern of infection and recovery reveals the underlying cause: recovery of infected individuals before they can transmit the disease shields clusters of individuals who would have otherwise been infected, as exemplified by the large uninfected region in the top right of Fig. 1f , which was heavily infected in the recovery-free case shown in Fig. 1b . Together, these results hint that recovery suppresses infection percolation. The competition between infection and recovery also drastically alters the time evolution of the (c) Growth of infected fraction φ over timeτ is hindered for diseases with low infectivity, but exhibits a generic quadratic scaling over time as indicated by the triangle, which eventually plateaus near unity for diseases with infectivity above a critical value. Infection growth is slower asP * decreases. (d) Total infected fraction φt exhibits an abrupt increase above a critical infectivity,Pc,0 ≈ 0.5, indicating a percolation transition to an epidemic. For a population with recovery afterτr = 4: (e) A disease with intermediate infectivityP * = 0.6 remains localized and does not spread, while (f ) a disease with higher infectivityP * = 0.7 spreads in a spatially heterogeneous, ramified pattern.

(g) Growth of infected fraction φ over timeτ is hindered for diseases with low infectivity, but exhibits a generic quadratic scaling at intermediate times and linear scaling at longer times as indicated by the triangles, before reaching a peak value φp before dropping rapidly to zero for diseases with infectivity above a critical value. Infection growth is slower asP * decreases.

(h) Total infected fraction φt and peak infected fraction φp, which represent the total infected fraction at the end of a given simulation and the maximal infected fraction during the simulation, respectively, both exhibit an abrupt increase above a critical infectivity,Pc ≈ 0.65, that is larger than the recovery-free case-indicating that recovery suppresses the percolation transition to an epidemic. Images in a-b, e-f are for the same timeτ = 100. All data are for a uniform distribution of individual interaction barriersP th,ij ∈ [0, 1]. Each curve in c,g and each point in d,h represents a separate simulation, all of which preserve the same lattice network structure of theP th,ij and only varyP * andτr, enabling us to systematically test the influence of these two key parameters on disease spreading.

infected fraction φ. For low infectivity, recovery is sufficiently fast to quench the spread of disease (light curves in Fig. 1g ). As a result, the total fraction of the population ever infected φ t 1 (points withP * < 0.65 in Fig. 1h ). By contrast, for sufficiently high infectivity, the disease initially spreads rapidly, first with φ ∼τ 2 as in the recovery-free case and then with φ ∼τ 1 (dark curves in Fig. 1g) . As time progresses φ eventually reaches a peak value φ p before dropping rapidly to zero as the entire population recovers. Both φ t and φ p increase precipitously above the critical infectivityP c ≈ 0.65 >P c,0 (Fig. 1h) , again indicating that recovery suppresses infection percolation. Again, the variability in the curves in Fig. 1g reflects the variability in the individual barriers; different simulations employing different choices of the randomly-chosenP th,ij would yield slightly different dynamics. However, in all cases, we expect that decreasing P * leads to slower transmission as a greater proportion of disease transmission occurs at ∆τ ij > 1.

To further explore the competition between infection and recovery, we inspect the spatiotemporal patterns of both using simulations performed at several different values ofP * andτ r , with snapshots all taken atτ = 75 shown in Fig. 2 . For a large value ofτ r = 10, the critical infectivity for percolation isP c ≈ 0.6; forP * slightly above this value, disease spreading is heterogeneous, with recovered individuals again shielding clusters of individuals who would have otherwise been infected FIG. 2. The competition between infection and recovery determines the onset and pattern of disease spreading. For a population with a longer recovery durationτr = 10.0, (a) a disease with lower infectiv-ityP * = 0.7 spreads in a ramified pattern, while diseases with higher infectivitiesP * = 0.9 and 1.1 (bc) spread in a more compact manner. The leading front of the infected population is trailed by an inner region of recovery, leading to the formation of a circular traveling pulse of infection. For a shorter recovery duration ofτr = 3.0, (de) the threshold infectivity for appreciable disease spreading is larger, and (e-f ) the pulse of infection is thinner. These effects are even more pronounced for the shortest recovery duration ofτr = 2.1, as shown in (g-i). The thick green line indicates the transition to infection percolation. All images are shown for the same timeτ = 75 and all data are for a uniform distribution of individual interaction bar-riersP th,ij ∈ [0, 1] preserving the same lattice network structure of theP th,ij throughout, thus demonstrating how varyingP * andτr influences disease spreading. (Fig. 2a) . For even higherP * , disease spreading becomes more compact (Figs. 2b-c) . Moreover, the leading front of the infected population is trailed by an inner compact region of recovery, leading to the formation of a wide circular pulse of infection that travels outward through the population (dark blue region in Figs. 2b-c). This feature is notably similar to the traveling pulses of infection predicted by extensions of the classic SIR model [19, 43, 44] -highlighting the ability of our framework to reproduce previously-reported continuum-scale phenomena. We observe similar behavior at smaller values of τ r , but shifted to increased values ofP * (Figs. 2d-i): recovery can increasingly quench the spread of disease as it becomes faster relative to infection transmission. As a result, the critical infectivity for percolation,P c , shifts to higher values (thick green line in Fig. 2) , further confirming that recovery suppresses infection percolation. Recovery also strongly impacts the number of infected individuals: asτ r decreases, the thickness of the pulse of infected individuals decreases (compare Figs. 2c, f, i).

These results demonstrate the key influence of recovery on disease spreading; we therefore examine the underlying model further to develop and test analytical relations that confirm the internal consistency of our simulations. We first focus on the observation that faster recovery can quench the spread of disease. At the leading front of the infected population, new individuals are infected only when the infection transmission time ∆τ ij = 1/(1 −P th,ij /P * ) is shorter than the recovery timeτ r ; otherwise, an infected individual recovers from the disease before they are able to transmit it to a neighbor. Thus, only individual interactions with P th,ij <P * (1 −τ −1 r ) can transmit disease. For a disease to continually spread throughout the population, a sufficient number of these interactions must permit disease transmission, as given by the critical probability of bond percolationP c,0 ; that is, infection percolates only wheñ P c,0 <P * (1 −τ −1 r ). The critical infectivity in a population with recovery is then directly given by the relatioñ

which increases withτ r : faster recovery suppresses infection percolation. This dependence is in good agreement with the shift in the critical infectivity shown in Fig. 1d , h and by the thick green line in Fig. 2 .

To further confirm the internal consistency of our simulations, we explore values ofτ r spanning nearly two orders of magnitude. For each simulation, we vary the (d-f ) Data for the rescaled peak infected fractionφp ≡ φp/(φp)max show reasonable alignment when plotted as a function of the shifted disease infectivityP * −Pc; insets show that the theoretical φ(τp) calculated using Eq. 2 (dashed lines) provides a reasonable approximation to the maximal peak fraction (φp)max determined forP * = 2. Intriguingly, all the datasets appear to approximately converge asτr increases. All three distributions spanP th,ij ∈ [0, 1]; the normal distribution is centered at 0.5 and has standard deviation = 0.25, while the bimodal is constructed from two normal distributions centered at 0 and 1, both with standard deviation = 0.25. Each point represents a separate simulation, all of which preserve the same lattice network structure of theP th,ij for a given distribution and only varyP * andτr, enabling us to systematically test the influence of these two key parameters on disease spreading. disease infectivityP * and identify the critical infectivitỹ P c at which the total fraction of the population ever infected φ t abruptly increases, similar to the curves shown in Fig. 1d , h. Consistent with our expectation, all the data for differentτ r show similar growth when plotted as a function of the shiftedP * −P c (Fig. 3a) . Moreover, the variation of the critical infectivityP c withτ r shows excellent agreement with Eq. 1 (dashed line, Fig. 3a  inset) .

Finally, we test the generality of this relation by exploring two other distributions ofP th,ij : a normal distribution, representing a population with random and uncorrelated interactions whose variation is distributed about a single mean value, and a bimodal distribution, representing a population with distinct high-risk and low-risk subpopulations arising from, for example, differences in compliance with public health interventions. Our central findings that the onset and dynamics of disease spreading are regulated by the interplay between disease transmission and recovery are consistent across the different distributions tested. Specifically, for all the distributions tested, we observe a similar abrupt increase in φ t withP * above a critical infectivityP c , with excellent alignment of all the data for differentτ r (Figs. 3a-c) . In all cases, the value of P c determined shows close agreement with Eq. 1 (Figs. 3a-c insets). A bimodal population in which half of the contacts between members maintain high barriers to transmission while the other interspersed connections do not-such as through strong or weak social distancing, respectively-hasP c slightly less than this relation, indicating that it is slightly more susceptible to infection. with rescaled timeτ /τr shows general quadratic to linear scaling as indicated by the triangles, followed by a drop to zero at τ f /τr ≈ Nt/2/τr + 1. Data shown are for the largest disease infectivity tested,P * = 2. Distributions of individual interaction barriers and colors indicating different recovery durationsτr are the same as in Fig. 3 . Each point represents a separate simulation, all of which preserve the same lattice network structure of theP th,ij for a given distribution and only varyτr.

General scaling of infection growth dynamics. We next focus on the observation that, for highP * , the leading front of the infected population is trailed by an inner region of recovery. At the leading front of the infected population, the shortest possible disease transmission time is again ∆τ ij = 1; therefore, the leading front spans an overall region with a maximal radius ofτ individuals, as in the recovery-free case. For short timesτ <τ r , these individuals have not yet recovered; hence, we expect that φ ≈ πτ 2 /N t again in this regime, in agreement with the short-time quadratic scaling shown in Fig. 1g as well as the recovery-free quadratic scaling shown in Fig. 1c . At longer timesτ ≥τ r , infected individuals begin to recover, forming an inner region of recovery. The leading front of this region spreads at the same rate as the leading front of the infected population; however, it only spans a maximal radius ofτ −τ r individuals. The maximal infected fraction of the population is then given by the area between the leading fronts of the infected and recovered popula- (Fig. 1g) . Furthermore, as the infected region spreads, it eventually reaches the boundary, followed by the growing region of recoveryτ r later. We therefore expect that the infected fraction φ will peak at a timeτ p ≈ √ N t /2 before dropping rapidly and reaching zero atτ f ≈ N t /2 +τ r , corresponding toτ p ≈ 50 andτ f ≈ 75, respectively, for τ r = 4 (SI text)-in good agreement with the results shown in Fig. 1g . Together, these calculations yield a general expression for the full time evolution of the infected fraction of a population in the limit of highP * well above the threshold P c , in excellent agreement with the data shown in Fig. 1g :

AsP * decreases, we again expect that the infected region spreads at a slower rate, prolongingτ p -also in good agreement with the different curves shown in Fig. 1g .

Eq. 2 also provides an approximation of the peak infected fraction, φ p = φ(τ p ), in the limit of highP * ; for the case ofτ r = 4, we estimate φ p ≈ 0.1, consistent with the data shown in Fig. 1g . To further examine this relation, we analyze the results of all the simulations described in Fig. 3 with varyingτ r . For each simulation, we determine the maximal value of φ p , (φ p ) max , atP * = 2, the highest value tested. Consistent with our expectation, φ p = φ(τ p ) calculated using Eq. 2 (dashed line, Fig. 3d inset), provides a reasonable approximation to the measured (φ p ) max (data points, Fig. 3d inset) . Furthermore, all the data forφ p ≡ φ p /(φ p ) max show reasonable alignment when plotted as a function of the shiftedP * −P c (Fig. 3d) , with some deviation for the smallestτ r likely arising from geometric effects not taken into account in our simple estimate of φ p = φ(τ p ). We again test the generality of these results by exploring their applicability to a normal distribution and a bimodal distribution ofP th,ij ; the results are closely similar in all three cases (Figs. 3d-f) . Taken together, all of the results shown in Figs. 1-3 support the validity of our scaling relations.

As a final exploration of these relations, we use all of our simulations to directly test the general expression for the infected fraction φ given by Eq. 2. All of the data collapse onto this scaling curve for all values ofτ r and distributions ofP th,ij tested (Figs. 4a-c) . Furthermore, the peak in φ and its eventual drop to zero occur at values ofτ /τ r close to the predicted valuesτ p /τ r ≈ √ N t /2τ r andτ f /τ r ≈ N t /2/τ r + 1, respectively. Thus, the close agreement between all the data and the theoretical prediction confirm the internal consistency of our simulations with the scaling relations. We note that the analytical expression represents the solution for the high infectivity limit with highP * ; however, as indicated by the data in Figs. 1g and 4a-c, it provides a good approximation for the spreading dynamics of a broad range of diseases withP * P c . Intriguingly, similar quadratic power-law scalings have been reported for the initial regional epidemic spreads of COVID-19 [13] , though explaining the full dynamics of this pandemic involves additional complexities that we do not consider here.

Extension to spatially heterogeneous networks. While our simple implementation thus far considered distributions ofP th,ij that were homogeneously mixed over the social network lattice, our framework allows for spatially heterogeneous networks to be constructed to model disease spreading in specific community structures [8, 9, 76] . As an example of this idea, we test disease spreading from the boundary between two stratified subpopulations, each on a separate lattice with a different meanP th,ij , representing differences in susceptibility to infection that can arise from differences in containment strategies or in socioeconomic factors [3, [8] [9] [10] . Depending on the viral infectivityP * , we find that stratification can either markedly exacerbate or hinder the overall amount of infection, and the total amount of infection in each subpopulation, compared to the homogeneouslymixed case (Fig. 5) . Specifically, for diseases of intermediate infectivity, disease spreading is stronger in the stratified population-the total infected fraction φ t is two orders of magnitude larger than in the homogeneouslymixed case-due to the earlier onset of infection percolation in the more susceptible subpopulation. However, for diseases of higher infectivity, disease spreading is slightly stronger in the homogeneously-mixed population; in the stratified case, the less susceptible subpopulation buffers against continued spread of disease. These results exemplify how our framework provides a way to directly assess the critical role played by community structure in the spread of disease. More extreme forms of heterogeneity-for example with small-world connections between subpopulations, or power law distributions in susceptibility-can easily be incorporated within this framework, and will be an interesting direction for future research. For example, as suggested in previous work [44] , we anticipate that long-range connections between regions with disparate susceptibilities will increasingly fragment traveling waves of infection.

The framework presented here provides a direct way to merge the temporal dynamics underlying the classic SIR model with a network representation of all the discrete interactions between members of a population. We demonstrate this principle using 2D lattice simulations implementing a simplified form of this framework. Our simulations reveal that, for diverse diseases and populations, disease spreading can be understood as a process of infection percolation through a social network. By testing different recovery durations and distributions of individual interaction parameters, we find that the onset and dynamics of spreading are determined by the interplay between disease transmission and recovery at the scale of individual interactions. This finding thus builds on the rich body of previous work exploring disease spreading through the lens of percolation theory [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] .

Guided by these insights, we develop a scaling theory that yields general predictions for the total number of infections, maximal number of active infections, and the temporal evolution of active infections in a population. Importantly, our scaling theory clarifies how these predictions can be applied to different diseases, with varying infectivitiesP * and recovery durationsτ r , and different populations, with varying distributions of the individual barriers to interactionP th,ij e.g. due to different implementations of containment strategies like social distancing.

Our simulations represent a first step toward implementing this framework. Thus, we have necessarily made a number of simplifying assumptions that can be relaxed in future extensions of this work. For example, our simulations consider a social network represented by a static, 2D, square lattice-while in reality, social networks are dynamic and have more complex structures. While previous work suggests that the assumption of a static network may not greatly alter disease spreading [34] , our framework can be implemented on dynamically-changing networks using a time-dependent matrix of individual interactions (P th,ij ). Furthermore, while previous work suggests that disease spreading is well-approximated by spreading on a 2D network covering the Earth's surface [22, 77, 78] , networks of different dimensionality d can also be implemented-which would result in a more general form of Eq. 2 with φ ∼τ d at short times and φ ∼τ d−1 at longer times. Previous work also suggests that social networks can have broadly-distributed degree distributions [20, 72] , small world connections [70, 71] , or multiple layers [39] unlike the fixed 2D connectivity of our lattice; moreover, disease transmission can occur among groups, not just pairs of individuals, which would require consideration of a network with higher-order interactions [79] . Incorporating these features into our framework will be a valuable direction for future work. Indeed, the random long-range connections that arise in small-world networks could disrupt the traveling waves of infection and recovery we observed in our simple network, or even seed new traveling waves of infection and recovery, potentially leading to richer dynamics. Our treatment of individual interactions can also be extended in future work. For example, the barriers to disease transmissionP th,ij need not be undirected, and the recovery durationsτ r,i need not be constant throughout the network. Furthermore, our representation of the disease transmission function f = 1 −P th,ij /P * represents the simplest linear function, and can instead be replaced by a more complex form that incorporates the sophisticated dynamics of transmission specific to different diseases [80] .

Our framework could also enable straightforward assessment of the efficacy of different public health policies. For example, the implementation of strong social distancing results in an increase in P th,max , leading to a reduction inP * and hence a reduction in the peak infected fraction of the population, φ p (Figs. 3d-f)-consistent with previous work [10, 13, 81] . Alternatively, the development of better treatments shortening the recovery durationτ r also leads to a reduction in the maximal peak infected fraction of the population, (φ p ) max (Figs. 3d-f, insets); it also hastens the transition to slow linear scaling and eventual decline of infection growth (Figs. 4a-c)again consistent with previous work [41, 82] . The influence of other factors that are documented to impact disease spreading-e.g. seasonality of infectivity [11, 12, [83] [84] [85] [86] [87] , heterogeneity in community susceptibility [8, 9] , and targeted vaccination [10] -can also be evaluated through appropriate modifications toP * andP th,ij . For example, targeted vaccination yielding perfect immunity is typically implemented by removing nodes from the network [10, 88] ; in our framework, this could equivalently be accomplished by setting P th,ij > P * , therefore preventing further infection of a node, while imperfect immunization [89, 90] could alternatively be implemented through smaller increases in P th,ij .

Implementation of the network model. We implement the dynamic network model in MATLAB. To define each 2D square lattice of 100 by 100 individuals ("nodes"), we specify node locations and an adjacency matrix characterizing the connectivity of the network. To ensure only one horizontally-oriented border between the two subpopulations in the stratified population shown in Fig. 5 , we employ a periodic boundary condition in the horizontal direction by connecting the first and last nodes in each row for all simulations; the nodes at the top and bottom boundaries do not have such conditions, ensuring that the vertical direction is non-periodic. For the bonds between nodes ("edges"), we randomly assign the values of the interaction barrier P th,ij ∈ [0, 1] from a given distribution, as specified in the main text. Each simulation has a specified disease infec-tivityP * and recovery durationτ r . From the values of P th,ij andP * , we compute the discrete infection transmission times ∆τ ij = 1/(1 −P th,ij /P * ) for each edge. The simulations shown in Figs. 1-2, 3a ,d, and 4a are for the exact same lattice with the exact same configuration ofP th,ij taken from a uniform distribution. Similarly, the simulations shown in Figs. 3b,e, and 4b are for the exact same lattice with the exact same configuration ofP th,ij taken from a normal distribution, and the simulations shown in Figs. 3c,f, and 4c are for the exact same lattice with the exact same configuration ofP th,ij taken from a bimodal distribution.

To perform each simulation, we use a modified invasion percolation algorithm based on the method described by Masson et al. [91] . We start atτ = 0 by introducing the disease at the central node of the lattice, with all the other nodes specified as being susceptible (S). Then, for the next and each successive time step of the simulation, we use a binary tree structure to sort all edges in contact with infected nodes and find the most susceptible edge ij-the edge with the minimal infection transmission time ∆τ ij . The next node to become infected, j, is then the node that is connected to the infected region through this most susceptible edge. This target node is then specified as being infected (I), and its time of infection is specified by adding the time increment ∆τ ij to the overall elapsed timeτ . We also decrease the remaining transmission times for all edges in contact with infected nodes by this time increment. New edges made available for infection by node j are added to the binary tree; because the tree was mostly sorted in the last step, subsequent sorts are time-efficient. We incorporate recovery by identifying all infected nodes for which at leastτ r has elapsed since infection, setting its infection transmission time to all nodes connected to it as being equal to ∞, and specifying the node as being recovered (R).

Size dependence. Our simulations occur on a fixed finite size network. All simulations presented are conducted on networks with N = 10 4 nodes. We verify that our results are not considerably influenced by finite-size effects by repeating simulations for five network sizes. Figure 6 shows that the critical growth of the total infected fraction φ t aboveP * > 0.5, corresponding to Fig. 1d in the main text, is insensitive the system size, even when the number of nodes is increased by two orders of magnitude (N = 10 4 to 10 6 ). We hence anticipate our system size N = 10 4 to be sufficiently large to capture the general dynamics of this spread.

FIG. 6. Critical growth of the total infected fraction φt above a critical infectivityP * > 0.5 for a recovery-free population. Network size is varied logarithmically from N = 10 4 to 10 6 with no discernible difference (some curves lie beneath N = 10 4 ).

Estimate of characteristic disease spreading timescales. For the case of a recovery-free population, we calculate the shortest possible timeτ f,0 at which the infected fraction plateaus in the limit of high disease in-fectivityP * . As time progresses, the disease spreads radially outward. Because we consider a square network comprising N t nodes in total, √ N t on a side, the leading edge of the circular infected region first reaches the boundary of the population whenτ ≈ √ N t /2. However, the total infected fraction can continue to grow: it only plateaus when it spans the entire 2D network, including its corners. This occurs when the radius of the infected region is equal to half the diagonal of the square network, τ f,0 ≈ ( √ N t /2) × √ 2 = N t /2. For the case of a population with recovery durationτ r , we extend this calculation to estimate the time at which the infected fraction of the population will peak,τ p , as well as the time at which the infected fraction of the population reaches zero after all individuals recover,τ f , again in the limit of high disease infectivity. As time progresses, the disease spreads radially outward in a circular infected region, followed by an inner circular region of recovery that spreads at the same rate but is delayed bỹ τ r . We consider two separate regimes: the "thin pulse" regime withτ r < √ N t (1/ √ 2−1/2), and the "thick pulse" regime withτ r > √ N t (1/ √ 2 − 1/2). For a thin pulse, as in the recovery-free case, the leading edge of the infected region first reaches the boundary of the population whenτ ≈ √ N t /2 (Fig. 7a) . At this time, the total infected fraction is nearly maximal, and we therefore approximateτ p ≈ √ N t /2. As time progresses, the leading edge of the region of recovery then first reaches the boundary of the population at a timẽ τ ≈ √ N t /2 +τ r (Fig. 7b) . Both regions continue to spread into the corners of the square boundary, and the leading edge of the infected region eventually reaches the corners at a timeτ ≈ N t /2 as in the recovery-free case (Fig. 7c) . Subsequently, the region of recovery con- tinues to grow; the total infected fraction continues to decrease, eventually reaching zero when the region of recovery has reached the corners of the square boundary, τ f ≈ N t /2 +τ r .

For a thick pulse, the leading edge of the infected region again first reaches the boundary of the population whenτ ≈ √ N t /2 (Fig. 7d) . The leading edge of the infected region then reaches the corners of the square boundary at a timeτ ≈ N t /2 as in the recoveryfree case (Fig. 7e) . Thus, the time at which the infected fraction is maximal is between these two times: √ N t /2 τ p N t /2. As time progresses, the region of recovery then continues to grow, eventually first reaching the boundary of the population atτ ≈ √ N t /2 +τ r (Fig. 7f) . Subsequently, the region of recovery continues to grow; the total infected fraction continues to decrease, eventually reaching zero when the region of recovery has reached the corners of the square boundary, τ f ≈ N t /2 +τ r . For our simulations with N t = 10 4 , the transition between the thin and thick pulse regimes occurs atτ r ≈ 21; therefore, our analysis of the example system withτ r = 4 presented in the main text is in the thin pulse regime, withτ p ≈ √ N t /2 ≈ 50 andτ f ≈ N t /2 +τ r ≈ 75 as reported in the main text. Together with Eq. 2, these estimates provide a universal scaling for the peak infection time, φ p = φ(τ p ) (Figs. 3d-f insets).

Movie S1. Sequence of infection for a disease with low infectivityP * = 0.3, showing that disease spreading is quickly localized. This simulation is without recovery.

Movie S2. Sequence of infection for a disease with intermediate infectivityP * = 0.6, showing that the disease spreads in a spatially heterogeneous, ramified pattern, leading to the formation of discrete clusters of bypassed individuals who remain uninfected. Infected individuals are shown in dark blue, uninfected susceptible individuals are shown in white. This simulation is without recovery.

Movie S3. Sequence of infection for a disease with higher infectivityP * = 0.7, showing that the disease spreads in a more compact pattern, with a smoother leading edge, leading to the formation of fewer and smaller clusters of bypassed individuals. Infected individuals are shown in dark blue, uninfected susceptible individuals are shown in white. This simulation is without recovery.

Movie S4. Sequence of infection for a disease with intermediate infectivityP * = 0.6, showing that recovery causes disease spreading to be quickly localized. Infected individuals are shown in dark blue, recovered individuals are shown in light green, uninfected susceptible individuals are shown in white. This simulation is withτ r = 4.

Movie S5. Sequence of infection for a disease with higher infectivityP * = 0.7, showing that the disease spreads continually, but recovery causes the disease to spread in a spatially heterogeneous, ramified pattern, leading to the formation of discrete clusters of bypassed individuals who remain uninfected. Infected individuals are shown in dark blue, recovered individuals are shown in light green, uninfected susceptible individuals are shown in white. This simulation is withτ r = 4.

It is a pleasure to acknowledge Navid C.P.D. Bizmark for stimulating discussions. This work was supported by startup funds from Princeton University. This material is also based upon work supported by the National Science Foundation Graduate Research Fellowship Program (to C.A.B.) under Grant No. DGE1656466. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. 

Simple compartmental models for disease transmission

Mathematical models in population biology and epidemiology

Mathematical models of infectious disease transmission

An application of the theory of probabilities to the study of a priori pathometry.-Part I

A contribution to the mathematical theory of epidemics

Fractional order SIR model with constant population

Dynamics of an SIR model with vaccination dependent on past prevalence with high-order distributed delay

The effect of heterogeneity in measles vaccination on population immunity

Measles outbreak in a community with very low vaccine coverage, the Netherlands

Optimal vaccination schemes for epidemics among a population of households, with application to variola minor in Brazil

Multinational influenza seasonal mortality study group: Multinational impact of the 1968 Hong Kong influenza pandemic: evidence for a smoldering pandemic

Influenza pandemics of the 20th century

Effective containment explains subexponential growth in recent confirmed COVID-19 cases in China

The mathematical theory of infectious diseases and its applications (Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE

Spatial contact models for ecological and epidemic spread

On the critical behavior of the general epidemic process and dynamical percolation

Universality classes for spreading phenomena: A new model with fixed static but continuously tunable kinetic exponents

Geometrical cluster growth models and kinetic gelation

Travelling waves and spatial hierarchies in measles epidemics

Infection dynamics on scalefree networks

Epidemic outbreaks in complex heterogeneous networks

Percolation and epidemics in a two-dimensional small world

Spread of epidemic disease on networks

Superspreading and the effect of individual variation on disease emergence

Epidemic size and probability in populations with heterogeneous infectivity and susceptibility

Second look at the spread of epidemics on networks

On analytical approaches to epidemics on networks

The abundance threshold for plague as a critical percolation phenomenon

Effects of epidemic threshold definition on disease spread statistics

Epidemic threshold for the susceptible-infectious-susceptible model on random networks

Message passing approach for general epidemic models

The effect of heterogeneity on invasion in spatial epidemics: from theory to experimental evidence in a model system

Heterogeneity in susceptible-infectedremoved (SIR) epidemics on lattices

Shift of percolation thresholds for epidemic spread between static and dynamic smallworld networks

Unification of theoretical approaches for epidemic spreading on complex networks

Spatio-temporal synchronization of recurrent epidemics

Modeling the spatiotemporal epidemic spreading of covid-19 and the impact of mobility and social distancing interventions

Spatial structure in disease transmission models

Epidemic model with isolation in multilayer networks

Unexpected epidemic thresholds in heterogeneous networks: The role of disease transmission

Comparing nonpharmaceutical interventions for containing emerging epidemics

Moving infections: individual movement decisions drive disease persistence in spatially structured landscapes

Spatial-temporal dynamics in nonlocal epidemiological models

Analytic models for sir disease spread on random spatial networks

Epidemics spread in heterogeneous populations

Introduction to percolation theory

A fast two-step algorithm for invasion percolation with trapping

Controlling capillary fingering using pore size gradients in disordered media

Dynamics of invasion percolation

Temporal development of invasion percolation

Percolation and invasion percolation as fractal growth problems

Epidemics and percolation in small-world networks

disorder, and percolation in epidemics

Percolation on heterogeneous networks as a model for epidemics

Epidemics, disorder, and percolation

Percolation and epidemics in random clustered networks

Epidemic spreading and risk perception in multiplex networks: a self-organized percolation method

Epidemic models and percolation

Epidemic variability in complex networks

Appropriate models for the management of infectious diseases

Nonmarkovian infection spread dramatically alters the susceptible-infected-susceptible epidemic threshold in networks

Estimation of the infection time and latency distribution of AIDS with doubly censored data

Modeling the progression of HIV infection

A stochastic model for the assessment of the transmission pathways of heart and skeleton muscle inflammation, pancreas disease and infectious salmon anaemia in marine fish farms in Norway

Risk factors for pancreas disease (pd) outbreaks in farmed atlantic salmon and rainbow trout in Norway during

Transmission dynamics of pancreas disease (pd) in a Norwegian fjord: aspects of water transport, contact networks and infection pressure among salmon farms

Highly sensitive real-time in vivo imaging of an influenza reporter virus reveals dynamics of replication and spread

Efficient reverse genetics reveals genetic determinants of budding and fusogenic differences between nipah and hendra viruses and enables real-time monitoring of viral spread in small animal models of henipavirus infection

Visualizing real-time influenza virus infection, transmission and protection in ferrets

Collective dynamics of 'small-world'networks

Navigation in a small world

Scale-free networks: a decade and beyond

The critical probability of bond percolation on the square lattice equals 1/2

Invasion percolation on a random lattice

ZHURNAL EKSPERIMENTALNOI I TEORETICH-ESKOI FIZIKI

Self-similar community structure in a network of human interactions

Saprotrophic invasion by the soil-borne fungal plant pathogen rhizoctonia solani and percolation thresholds

Empirical evidence of spatial thresholds to control invasion of fungal parasites and saprotrophs

Abrupt phase transition of epidemic spreading in simplicial complexes

Epidemiological implications of the contact network structure for cattle farms and the 20-80 rule

Social distancing strategies against disease spreading

Convalescent plasma as a potential therapy for COVID-19

Pathogen adaptation to seasonal forcing and climate change

Absolute humidity and the seasonal onset of influenza in the continental United States

Improvement of disease prediction and modeling through the use of meteorological ensembles: human plague in Uganda

Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period

Seasonal dynamics of recurrent epidemics

Statistical physics of vaccination

Imperfect vaccines and herd immunity to HIV

Modeling malaria vaccines II: Population effects of stagespecific malaria vaccines dependent on natural boosting

A fast algorithm for invasion percolation

Competing Interests: The authors declare no competing interests.