Epidemic plateau in critical SIR dynamics with non-trivial initial conditions
EEpidemic plateau in critical SIR dynamics with non-trivial initial conditions
Filippo Radicchi
Center for Complex Networks and Systems Research, Luddy School of Informatics,Computing, and Engineering, Indiana University, Bloomington, IN 47408
Ginestra Bianconi
The Alan Turing Institute, 96 Euston Rd, London NW1 2DB, United KingdomSchool of Mathematical Sciences, Queen Mary University of London, London, E1 4NS, United Kingdom
Containment measures implemented by some countries to suppress the spread of COVID-19 haveresulted in a slowdown of the epidemic characterized by time series of daily infections plateauingover extended periods of time. We prove that such a dynamical pattern is compatible with criticalSusceptible-Infected-Removed (SIR) dynamics. In traditional analyses of the critical SIR model,the critical dynamical regime is started from a single infected node. The application of containmentmeasures to an ongoing epidemic, however, has the effect to make the system enter in its criticalregime with a number of infected individuals potentially large. We describe how such non-trivialstarting conditions affect the critical behavior of the SIR model. We perform a theoretical andlarge-scale numerical investigation of the model. We show that the expected outbreak size is anincreasing function of the initial number of infected individuals, while the expected duration of theoutbreak is a non-monotonic function of the initial number of infected individuals. Also, we preciselycharacterize the magnitude of the fluctuations associated with the size and duration of the outbreakin critical SIR dynamics with non-trivial initial conditions. Far from heard immunity, fluctuationsare much larger than average values, thus indicating that predictions of plateauing time series maybe particularly challenging.
I. INTRODUCTION
At the onset of the COVID-19 pandemic, world-widetime series of the number of infected individuals havedisplayed an exponential growth. Such a behavior iswell predicted by standard epidemic frameworks [1]. Inslightly later stages, however, time series have exhibitednon-trivial dynamical patterns. Many papers have at-tempted to model observed behaviors and to determinethe role of containment measures [2–17]. The commonand reasonable assumption is that containment measuresimplemented in the attempt of mitigating the outbreakhave strongly influenced the unfolding of the epidemic.Unfortunately, this a setting where modeling attemptsare particularly challenging. The effective implementa-tion of containment measures imposed by authorities relyon people personal judgements and adaptive behavior,and while epidemic spreading is a well-studied branchof mathematical biology [18], statistical physics [19] andnetwork science [20–26], the modelling of adaptive be-havior is only at its infancy [27–29].According to the data, in several countries, the slow-down of the epidemic spread is characterized by an al-most flat time series of daily number of new infections.Moreover, the time series of the number of removed indi-viduals display power-law growth instead of an exponen-tial growth as a function of time [12, 14, 16, 17]. Here,we propose a theoretical interpretation of those featuresas the signature of the system being in (or near) its crit-ical regime. Criticality is a fundamental property char-acterizing the dynamics of biological and socio-technicalsystems [30–32]. Our work consists of an in-depth inves-tigation of a critical Susceptible-Infected-Removed (SIR) dynamics starting from a non-trivial initial configurationcharacterized by n initially infected individuals. We in-terpret the emergence of the critical regime as the resultof disease containment strategies, and the non-trivial ini-tial condition as the configuration of the system whenspreading becomes critical. In the typical setting con-sidered in statistical mechanics [33, 34], a single seed isgenerally used as the initial condition for critical SIR dy-namics; the mapping of the critical SIR to the criticalstandard branching process allows for a full characteri-zation of the spreading dynamics [35, 36]. The realisticassumption of having an initial number of infected in-dividuals n > n .In this paper, we evaluate, by means of analytic argu-ments and large-scale simulations, the scaling of thesequantities as functions of the population size N .The paper is structured as follows: in Sec II we providethe theoretical interpretation of the plateau as a criticalSIR dynamics starting from n > a r X i v : . [ phy s i c s . s o c - ph ] J u l by extensive numerical simulations of the process; finallyin Sec IV we provide the concluding remarks. The Ap-pendix describes the Gillespie algorithm used in this workto simulate the critical SIR dynamics. II. THE THEORETICAL INTERPRETATIONOF THE PLATEAU
We consider the Susceptible-Infected-Removed (SIR)model on a well-mixed population [19–26]. At any pointin time, individuals can be found in three possible states:susceptible, infected and removed. Susceptible individu-als do not carry the disease but they can be infected; in-fected individuals carry the disease, and they can spreadit to susceptible individuals; removed individuals are ei-ther removed or deceased, and they do no partecipate inthe spreading dynamics. We indicate with λ the rate ofinfection, i.e., the expected number of spreading eventsoccurring per unit of time. Without loss of generality, weset the recovery rate equal to one.We start our discussion by focusing on the determinis-tic treatment of the SIR model on a well-mixed popula-tion with infinite size. If we indicate with s , i and r thefractions of susceptible, infected and removed individu-als, respectively, we can write dsdt = − λsi,didt = λsi − i,drdt = i. (1)Please note that s + i = r = 1. The critical dynamicalregime is characterized by λ c = 1 . (2)If we start from an initial condition consisting of a frac-tion i (0) of infected individuals and a fraction r (0) = 0of removed individuals, at the onset of the epidemic, i.e., t (cid:28)
1, we observe a different behavior depending on thevalue of λ . In the non-critical regime, i.e., λ (cid:54) = λ c , thedeterministic equations for i and r read didt = λsi − i (cid:39) ( λ − i,drdt = i. (3)Solutions of the above equations are i ( t ) = i (0) e ( λ − t ,r ( t ) = i (0) λ − e ( λ − t . (4)In essence, in the subcritical regime, i.e., λ < λ c , thenumber of infected individuals decays exponentially fast, and the number of removed individuals remains vanish-ing. In the supercritical regime, i.e., λ > λ c , the numberof infected and removed individuals display an exponen-tial increase. At criticality, i.e., λ = λ c , the deterministicequations for i and r are didt = ( s − i (cid:28) ,drdt = i , (5)leading to i ( t ) (cid:39) i (0) ,r ( t ) (cid:39) i (0) t. (6)Therefore, according to the deterministic approach, weshould expect that the number of removed individuals atcriticality increases linearly in time with a slope that isgiven by the initial condition i (0). -7 -3 (a)(b) FIG. 1: a) The time series of the number of infected individ-uals I ( t ) are plotted close to the critical point as predictedby the deterministic SIR equations. The population size is N = 10 , and spreading is started from n = 1 seed. Differ-ent curves correspond to different values of λ = 1 + 2 − i with2 ≤ i ≤
9. b) Same as in panel a, but for λ = 1. Differentcurves correspond to different numbers of initially infectednodes n = 2 i , with 0 ≤ i ≤
9. As λ approaches the criti-cal value λ c = 1 and n decreases toward one, we observe aplateauing of the time series. From the deterministic Eqs. (1), it is evident that dids = − λs . (7)The equation can be integrated to obtain the well-knownsolution [19] s ( t ) + i ( t ) − λ ln s ( t ) = s (0) + i (0) − λ ln s (0) . (8)Using Eqs.(1), we can express the logarithmic deriva-tive of the number of infected individuals as d ln idt = λs − , (9)where λ s ( t ) is the reproduction number.The former equations implies that the time series ofinfected individuals i ( t ) has a peak at t = t (cid:63) determinedby λs ( t (cid:63) ) = 1 . (10)The fraction of susceptible individuals at the peak of theepidemic is given by s (cid:63) = s ( t (cid:63) ) = 1 /λ . By making thefurther assumption that the epidemic starts from a frac-tion i (0) = 1 − s (0) of infected individuals and zero re-moved individuals r (0) = 0 in Eq. (8), we obtain i (cid:63) = i ( t (cid:63) ) = 1 − λ − λ ln( λs (0)) . (11)Using Eqs. (10) and (11) in the first of Eqs. (1), we get dsdt (cid:12)(cid:12)(cid:12)(cid:12) t = t (cid:63) = − λs (cid:63) i (cid:63) = − i (cid:63) (12)It follows that the second derivative of ln i is given by d idt (cid:12)(cid:12)(cid:12)(cid:12) t = t (cid:63) = i ( t (cid:63) ) d ln idt (cid:12)(cid:12)(cid:12)(cid:12) t = t (cid:63) = − λ ( i (cid:63) ) = − ρ (cid:63) , (13)where ρ (cid:63) is defined as ρ (cid:63) = − λ ( λ − − ln( λs (0))) . (14)We note that ρ (cid:63) is zero, i.e., we reach a plateau, only for s (0) = 1 and λ = 1. This fact implies that, in the de-terministic approach, a perfect plateau of the time seriesln i is never achieved for λ > λ = 1, s (0) = 1, we get s (0) (cid:39) ρ (cid:63) (cid:39) − λ (cid:20) ln( s (0)) + 12 ( λ − (cid:21) (cid:39) (cid:20) − s (0) + 12 ( λ − (cid:21) . (15)The above equation indicates that the conditions to havea near-plateau dynamics are having an infectivity rate λ close to one, and having the system as far as possiblefrom heard immunity, i.e., 1 − s (0) (cid:28)
1. In summary,the near-critical state for λ (cid:39) λ belowone (see Figure 1). III. SIR CRITICAL DYNAMICS WITHNON-TRIVIAL INITIAL CONDITION
From now on, we assume that the system is in the crit-ical regime. We further assume that spreading dynamicsis started from n > n = 1 individual, we are not awareof existing studies dealing with non-trivial initial condi-tions consisting of n > n ? What is the be-havior of the expected duration of the outbreak? Whatabout the expected size of the outbreak? What abouttheir fluctuations? (a)(c)(e) (f)(d)(b) FIG. 2: We show three examples of time series for the num-ber of infected individuals I ( t ) (panels a, c and e) and thecorresponding number of removed individuals R ( t ) (panels b,d and f) for a critical SIR dynamics with non-trivial initialcondition. The time series are obtained by simulating thestochastic SIR dynamics at criticality (with λ = 1) on a well-mixed population with identical parameters: initial numberof infected individuals n = 128 and population size N = 10 .The dashed lines indicate the corresponding deterministic pre-dictions. Please note that all the above questions cannot be an-swered with a purely deterministic approach. SIR out-break sizes and durations obey probability distributionsthat are well peaked around their expected value only ifthe system is off criticality. However, the very fact thatthe system is assumed to be in the critical regime impliesthat fluctuations have a dominant role in the determina-tion of the properties of the dynamical system. In Figure2 for example, we display time series representative forthe critical regime of the dynamics. Ground-truth timeseries are obtained by simulating the SIR stochastic dy-namics (see Appendix A for details). They are comparedwith the deterministic expectation value obtained by in-tegrating Eqs. (1). We note that some realizations of theprocess are more persistent and more pervasive in thepopulation than what predicted by the expected value.From here on, we abandon the deterministic SIR equa-tions and we embrace a stochastic approach. Critical SIRdynamics starting from a single initial seed, i.e., n = 1,is known to be characterized by extremely large fluctu-ations of the outbreak size and duration. These fluctu-ations can be quantified by leveraging the mapping be-tween critical SIR in a well-mixed population and themean-field branching process. In the following sections,we first review results valid for n = 1. Then, we focusour attention on the non-trivial case n > A. Critical dynamics with n = 1 initial seed If the initial condition is such that only one node isin the infected state while all other nodes are in thesusceptible state, the critical SIR model gives rise tooutbreaks that follow the statistics of a critical branch-ing process [35, 36] corrected by some scaling functions F T ( N/N (cid:63)T ) and F R ( N/N (cid:63)R ) that implement the effectivecutoff caused by finite-size effects [33, 34]. Here, N is thesize of the system; N (cid:63)T and N (cid:63)R are instead parametersthat determine when the cutoff takes place. Specifically,the distribution P ( T ) of the duration T of an outbreakfollows the law P ( T ) ∼ T − F T ( N/N (cid:63)T ) , (16)while the size of the outbreak R follows the distribution P ( R ) ∼ R − / F R ( N/N (cid:63)R ) . (17)The cutoff sizes N (cid:63)T and N (cid:63)R have been derived inRefs. [33, 34]. They are given by N (cid:63)T = T , N (cid:63)R = R / . (18)From the expressions for P ( T ) and P ( R ) given byEqs. (16) and (17), respectively, and further assuminga sharp cutoff, it is easy to deduce that the scaling withthe system size N of the average outbreak size (cid:104) R (cid:105) , theaverage duration (cid:104) T (cid:105) , and the standard deviations σ R and σ T [33, 34] obey (cid:104) R (cid:105) ∼ N / , σ R ∼ N / , (cid:104) T (cid:105) ∼ ln N, σ T ∼ N / . (19) We observe that all the above quantities are sub-extensive, as they all grow sub-linearly with the systemsize. The expected critical outbreak size (cid:104) R (cid:105) grows asthe system size to the power of 1 /
3. However, the stan-dard deviation associated to the outbreak size, i.e., σ R ,grows with increasing system size much faster than (cid:104) R (cid:105) .This fact indicates that it is very challenging to makepredictions if the dynamics is critical. Similarly, the out-break duration is characterized by large fluctuations inthe large population limit. We note that the exponents2 and 3 / P ( T ) and P ( R ) are thecritical mean-field exponents. These exponents are uni-versal and are observed for many critical spreading pro-cesses [40]. They characterize the critical SIR on networktopologies too as long as the underlying network has a ho-mogeneous degree distribution. In power-law networks,these exponents can deviate from their mean-field valuesas investigated in Refs. [40, 41].We have seen in Sec. II that the deterministic approachpredicts a linear increase of the number of removed in-dividuals with time. However, such a prediction is notaccurate for the ground-truth dynamics; accounting forstochastic effects correctly predicts a quadratic growthof the number of removed individuals in time when theepidemic starts with a single initial seed. To this end,the number of removed individuals grows in time as apower-law R = t z , (20)where z is a dynamical critical exponent, and t is theexpectation value of the time necessary to observe R re-moved individuals. The value of the dynamical criticalexponent can be obtained in different ways [36]. Here, wepresent the derivation of the value of the dynamical expo-nent based on Langevin-like equations for the dynamics.Starting from an initial fraction of infected individuals i (0) = n /N and a fraction r (0) = 0 of removed individ-uals we write didt = ( λs − i + c √ iη ( t ) ,drdt = i, (21)where η ( t ) is an uncorrelated white noise with E ( η ( t )) =0 and E ( η ( t ) η ( t (cid:48) )) = δ ( t − t (cid:48) ) and c is a constant. Atcriticality, i.e., λ = λ c = 1, thus, assuming t (cid:28) i (0) (cid:28)
1, we have λs − (cid:39) − i (0). We can thereforewrite didt = c √ iη ( t ) ,drdt = i. (22)We now perform a simple scaling analysis of this stochas-tic equations as usually done in non-equilibrium statisti-cal mechanics, e.g., Refs. [37, 42, 43]. If we rescale timeas t → bt (23)and define the scaling exponents z, α, for r, ir → b z t,i → b α t, (24)(25)as the exponents that leave the SIR critical dynamicsunchanged. The SIR stochastic Eqs. (22) read b α − didt = cb α/ − / η ( t ) ,b z − drdt = b α i, (26)from which we can derive the scaling exponents α = 1 (27) z = 2 . (28)In summary, in the critical dynamical regime, if thespreading is started from a single initial seed, we expectthat the average number of removed individuals growsquadratically with time. B. Critical dynamics with n > initial seeds Critical SIR dynamics started from the non-trivial ini-tial condition n > n = 1 seed. To include an explicit de-pendence on the parameter n in the scaling of Eq. (20),we correct it by introducing the scaling function F ( x ),where x = n /t ( R ). We impose that R (cid:39) t ( R ) z F (cid:18) n t ( R ) (cid:19) (29)with F ( u ) ∼ (cid:26) u (cid:29) u β if u (cid:28) t (cid:28) i (0) (cid:28)
1, the number of removed individuals growslinearly in time with a slope i (0) = n /N . Thus, wededuce that β = 1 . (31)This value is well supported by extensive numerical re-sults (see Figure 3) which confirm that there is a cross-over between linear and quadratic dependence of R on t ( R ).In the simulation of the SIR model, it is natural tostudy the behavior of R as a function of t ( R ). However,in real epidemic time series, the number of infected in-dividuals is measured over constant time intervals. Thetwo ways of monitoring the evolution of the process, i.e., R versus t ( R ) rather than R ( t ) versus time t , may leadto the observation of different scaling exponents. The -3 -2 -1 -1 (b)(a) FIG. 3: Dynamical properties of the critical SIR process withnon-trivial initial condition. a) We plot the number of re-moved individuals R versus the expected time ¯ t ( R ) requiredto observe R removed individuals. The different curves indi-cate different initial conditions, from bottom to top, n = 2 i with 0 ≤ i ≤
9. The population size is N = 10 . Thedashed lines are guides to the eye and correspond to linearand quadratic growth of R versus ¯ t ( R ), respectively. b) Ex-pected number ¯ R ( t ) of removed individuals as a function oftime t . Data are the same as in panel a. The dashed lines areguides to the eye and correspond to power-law growth of ¯ R versus t with exponent ξ = 1 and ξ = 2 .
5, respectively. discrepancy is due to the stochastic nature of the spread-ing process. The phenomenon is apparent from the re-sults of Figure 3: depending on the type of measurementperformed on the system, the power-law increase of thenumber of removed individuals as a function of time canbe described by a continuous range of exponents rangingfrom ξ = 1 to ξ ∼ .
5. We can therefore write R ( t ) (cid:39) n t ξ h ( t, N ) (32)where ξ is a decreasing function of n , and h ( t, N ) is amodulating function expressing the deviation from thepure power-law behavior. The ansatz of the above equa-tion is compatible with the power-law scaling of the em-pirical time series of removed individuals as a function oftime observed in countries where containment measureshave been implemented extensively [12, 14, 16, 17]. C. Distribution of avalanche durations and sizesfor the critical SIR model initiated by n > seeds In this section we investigate the statistical propertiesof the distribution of outbreak duration and size for thecritical SIR dynamics in a well mixed-population whenthe initial condition is non-trivial, i.e., n >
1. Scalingarguments suggest the following expression for the distri-bution P ( T ) of the critical outbreak duration TP ( T ) ∼ n T − F T ( N/N (cid:63)T , n /T ) . (33)The above scaling function is a natural modification ofEq. (16) by assuming that n scales like time. In particu-lar, the distribution P ( T ) is characterized by a lower cut-off depending on n . This fact is intuitive as an outbreakwith a larger number of initially infected individuals isnot expected to reach the absorbing state faster than anoutbreak started by a single seed (see Figure 4a). Wenote that, in the critical SIR dynamics, the dependenceon n does not lead to a change of the critical exponentvalues, as for example observed in other non-equilibriumphase transitions [37–39]. In Figure 5 a , we display thefunction w T ( n , T ) = − ln F T ( N/N (cid:63)T , n /T ) F T ( N/N (cid:63)T , /T ) (34)and we demonstrate that the scaling function w T ( n /T )for n (cid:28) N (cid:63)T can be approximated as w T ( n , T ) (cid:39) − n − T . (35)The scaling behavior, valid for n (cid:28) N / , can be jus-tified by assuming that each of the n seeds generatesan independent outbreak obeying the statistics of thecritical branching process. A critical avalanche startedfrom a single infected individual has a duration T fol-lowing the power-law distribution π ( T ) ∼ T − [35, 36].Thus, assuming independence among the n avalanches,we can estimate the probability P ( T ) as the probabilitythat among all n outbreaks the last outbreak to get ex-tinguished is extinguished at time T . Therefore in theinfinite population limit we obtain P ( T ) (cid:39) n [1 − π c ( T )] n − π ( T ) , (36)where π c ( T ) is the probability that an outbreak gener-ated by a single infected individual is not extinguishedat time T , with π c ( T ) (cid:82) ∞ T π ( x ) dx (cid:39) T . By assuming1 (cid:28) n (cid:28) N / , we get P ( T ) (cid:39) n exp (cid:18) − n − T (cid:19) . (37)Finally we note that while the scaling behavior describedin Eq. (35) has strong numerical confirmation for T (cid:28) N (cid:63)T for values of T ∼ N (cid:63)T the scaling function w T ( n , T )signals a dependence of the cutoff on n (see Figure 5).Scaling arguments suggest that the distribution P ( R )of critical outbreak size R should obey P ( R ) ∼ n R − / F R ( N/N (cid:63)R , n , R ) , (38)where the function F R ( N/N (cid:63)R , n , R ) implements a lowercutoff dependent exponentially on n (see Figure 4b). -5 -5 (b)(a) FIG. 4: The rescaled distribution of outbreak duration P ( T )(panel a) and outbreak size P ( R ) (panel b) are plotted fora well-mixed population of N = 10 individuals and initialnumber of infected individuals n equal to 2 i , with 0 ≤ i ≤ n increases, and an upper cutoff whose value does not stronglydepend on n . (a)(b) FIG. 5: The function w T ( n , T ) and w R ( n , R ) defined inEq. (34) and Eq. (39) are plotted for a population of N = 10 individuals and different values of n = 2 i with 0 ≤ i ≤ In Figure 5 b , we show the function w R ( n , R ) = − ln F R ( N/N (cid:63)R , n , R ) F R ( N/N (cid:63)R , , R ) (39)which, for n (cid:28) N (cid:63)R ( R, n ), can be approximated as w R ( n , R ) (cid:39) − n R − n / R . (40)This scaling function indicates that, for large values of R , R scales like n . For small values of R , it is possibleto observe some corrections would be required to fullydescribe the scaling. We notice that, in the first order in n , the normalization constant of the distribution P ( R )is independent of n . A way to interpret the result is byconsidering the infinite population limit approximatingthe distribution P ( R ) as the convolution of the n sizesof independent outbreak events. In this limit we have P ( R ) = (cid:90) dωe iωR [ F ( ω )] n , (41)where F ( ω ) = (cid:80) r Π( r ) e − iωr is the generating functionof the distribution Π( r ) of avalanches sizes of SIR crit-ical dynamics starting from a single seed. Assumingin first approximation that Π( r ) is a pure power-lawΠ( r ) ∼ r − / it follows that the logarithm of the generat-ing function behaves, for small ω , as ln F ( ω ) ∼ √ ω . Theresult, together with Eq. (41), indicates that R shouldscale as n for n (cid:28) N / . For more more details on theinfinite population limit we refer the reader to Ref. [44]. D. Statistical properties of the critical outbreakstarted by n > seeds We performed large-scale simulations of the criticalSIR model to address fundamental questions regardingthe distributions of duration T and size R of outbreaksstarted by a non-trivial initial condition n >
1. In Fig-ure 6, we display the average values and the standarddeviations of both T and R for a large system composedof N = 10 individuals. We display the moments of thedistributions as a function of the number of initial seeds n . The main outcomes are as follows. The expectedsize (cid:104) R (cid:105) is a growing function of n . The expected du-ration (cid:104) T (cid:105) is a non-monotonic function of n , displayinga single peak. The standard deviations σ T and σ R alsodisplay a peak as a function of n . The coefficient of vari-ation σ T / (cid:104) T (cid:105) and σ R / (cid:104) R (cid:105) are monotonically decreasingwith n . Therefore we conclude that the role of criticalfluctuation is very relevant for the critical dynamics, al-beit in relative terms is most severe for the SIR dynamicsstarting from a single initial seed. E. Scaling analysis of (cid:104) T (cid:105) and σ T We make the ansatz that the average duration (cid:104) T (cid:105) canbe described by (cid:104) T (cid:105) (cid:39) G ( n | α, β, H, K ) , (42)where the function G ( n | α, β, H, K ) is given by G ( n | α, β, H, K ) = n α H n β K . (43)Here, α and β are, in the large population limit, inde-pendent of N . On the contrary, H and K are dependent
06 (see Figure 9). We note thatthe logarithmic scaling of H is expected from the knownscaling of (cid:104) T (cid:105) for the SIR critical model starting from asingle initial seed. The exponents α and β are given by α = 0 . ± . , β = 1 . ± . . (48)The standard deviation of the outbreak duration σ T can be described in the same exact way as (cid:104) T (cid:105) . Theansatz σ T (cid:39) G ( n | α, β, H, K ) , (49)leads to the data collapse shown in Figure 7b. The scalingparameters H and K obey the scaling relations K (cid:39) a K N δ K , H (cid:39) a H N δ H . (50)with δ H = 0 . ± . b H = 1 . ± . δ K = − . ± . a K = 2 . ± . δ H (cid:39) / n = 1 the scaling reduces to the well-known scaling for the critical SIR model starting from asingle initial seed. Moreover, the exponent α and β for σ T are given by α = 0 . ± . , β = 1 . ± . . (51) -3 -1 n K -3 -2 -1 < T > H - K / -3 -1 n K -2 -1 T H - K / (a) (b) FIG. 7: Finite-size scaling analysis for the duration T of thecritical SIR started from n initial seeds. a) Data describingthe average value of the outbreak duration (cid:104) T (cid:105) for differentsystem sizes N are collapsed on a unique universal curve usingthe scaling of Eq. (45). Data are shown for population sizes N ranging from N = 10 to 10 . Each data point is obtainedby simulating the SIR process 10 times. b) Same as in panela, but for the standard deviation σ T . F. Scaling analysis of (cid:104) R (cid:105) and σ R The ansatz for (cid:104) R (cid:105) is slightly different from the oneappearing in Eq. (43), as it includes an additional loga-rithmic correction (cid:104) R (cid:105) (cid:39) n α H n β (ln n ) γ K . (52)We take H = 32 N / , K = N − / (53)and α = 1 , β = 0 . . (54) N H N -3 -2 K N H N -2 K (a) (c)(b) (d) FIG. 8: a) Scaling parameter H as a function of the systemsize N for the average duration (cid:104) T (cid:105) of critical outbreaks. Datapoints are the same as those of Figure 7. The line displaysthe best fit of the data points with Eq. (47). b) Same as inpanel a, but for the scaling parameter K . c) Same as in panela, but for the standard deviation σ T . The orange line is thebest fit of the data points with Eq. (50). d) Same as in panelc, but for the scaling parameter K . and we perform a fit of the exponent γ . As Figure 9aand Figure 10a demonstrate, the function gives rise toexcellent data fits as long as the exponent γ is γ = a γ ln N + b γ (55)with a γ = 0 . ± .
003 and b γ = 0 . ± . (cid:104) R (cid:105) can be rescaled and the data obtainedfor different N collapsed (see Figure 10). This is done bynoting that (cid:104) R (cid:105) y ( n ) = g ( x ( n ) | α, β ) , (56)where x ( n ) = n (ln n ) γ/β K /β ,y ( n ) = H − K α/β (ln n ) αγ/β , (57)and the function g ( x | α, β ) is given by Eq. (44). The ex-pected size of the outbreak (cid:104) R (cid:105) does not display a maxi-mum as a function of n , i.e., it is a monotonous increas-ing function of n . The standard deviation σ R can beinstead fitted using the same ansatz as σ T , i.e. σ R (cid:39) G ( n | α, β, H, K ) . (58)The best estimates of the parameters are α = 0 . ± . , β = 0 . ± . , (59)and K (cid:39) a K N δ K , H (cid:39) a H N δ H . (60)with δ H = 0 . ± . b H = 1 . ± . δ K = − . ± . a K = 0 . ± . n n H n -3 -2 -1 K (a)(b)(c) FIG. 9: a) The fitting parameter γ , i.e., Eq.(52), as a functionof the system size N . Data points are the same as those ofFigure 7. The line displays the best fit of the data points withEq. (55). b) Scaling parameter H as a function of the systemsize N for the standard deviation σ R of critical outbreaks.The solid line corresponds to the scaling function of Eq.(60).c) Same as in panel b, but for the scaling parameter K . Thesolid line corresponds to the scaling function of Eq.(60). -3 n K -2 -1 R H - K / -8 -4 x(n ) -7 -5 -3 < R > y ( n ) FIG. 10: Finite-size scaling analysis for the size R of thecritical SIR started from n initial seeds. a) Data describingthe average value of the outbreak size (cid:104) R (cid:105) for different systemsizes N are collapsed on a unique universal curve using thescaling of Eq. (56). Data are the same as of Figure 7. b)Same as in panel a, but for the standard deviation σ R . IV. CONCLUSIONS
Motivated by the current COVID-19 pandemic,we have investigated the critical properties of theSusceptible-Infected-Removed (SIR) dynamics in wellmixed-populations starting from non-trivial initial condi-tions consisting of n > n > n > n of the initial seed set increases. Moreover, numer-ical results indicate that, as the initial number of infectedindividuals n increases, the expected size of the outbreakincreases while the expected duration first increases andthen decreases, displaying a maximum. Using scaling ar-guments, and extensive numerical simulations we havededuced the scaling of the maximum duration and thecorresponding number of initially infected individuals. Acknowledgements
We thank P.L. Krapivsky, Geza Odor and R. M. Ziff forinteresting discussions. F.R. acknowledges support fromthe National Science Foundation (CMMI-1552487). [1] R. M. Anderson, B. Anderson, and R. M. May,
Infec-tious diseases of humans: dynamics and control (Oxforduniversity press, 1992).[2] S. Maslov and N. Goldenfeld, arXiv:2003.09564 (2020).[3] G. N. Wong, Z. J. Weiner, A. V. Tkachenko, A. El- banna, S. Maslov, and N. Goldenfeld, arXiv preprintarXiv:2006.02036 (2020).[4] L. Ferretti, C. Wymant, M. Kendall, L. Zhao, A. Nur-tay, D. G. Bonsall, and C. Fraser, Science (2020),10.1126/science.abb6936. [5] G. Bianconi, H. Sun, G. Rapisardi, and A. Arenas, arXivpreprint arXiv:2007.05277 (2020).[6] D. Fanelli and F. Piazza, Chaos, Solitons & Fractals ,109761 (2020).[7] T. Carletti, D. Fanelli, and F. Piazza, arXiv preprintarXiv:2005.11085 (2020).[8] A. Bianconi, A. Marcelli, G. Campi, and A. Perali, arXivpreprint arXiv:2004.04604 (2020).[9] S. Bradde, B. Cerruti, and J.-P. Bouchaud, arXivpreprint arXiv:2006.09829 (2020).[10] P. Maheshwari and R. Albert, arXiv preprintarXiv:2006.09189 (2020).[11] A. Arenas, W. Cota, J. Gomez-Gardenes, S. G´omez,C. Granell, J. T. Matamalas, D. Soriano-Panos, andB. Steinegger, MedRxiv (2020).[12] A. L. Ziff and R. M. Ziff, MedRxiv preprint (2020),10.1101/2020.02.16.2002382.[13] G. Bianconi and P. L. Krapivsky, arXiv preprintarXiv:2004.03934 (2020).[14] M. Nekovee, medRxiv (2020).[15] O. Valba, V. Avetisov, A. Gorsky, and S. Nechaev,arXiv:2003.12290 (2020).[16] B. Blasius, arXiv:2004.00940 (2020).[17] A. Brandenburg, arXiv preprint arXiv:2002.03638(2020).[18] J. D. Murray, Mathematical biology: I. An introduction ,Vol. 17 (Springer Science & Business Media New York,2007).[19] P. L. Krapivsky, S. Redner, and E. Ben-Naim,
A kineticview of statistical physics (Cambridge University Press,Cambridge, 2010).[20] A.-L. Barab´asi et al. , Network science (Cambridge Uni-versity Press, Cambridge, 2016).[21] M. Newman,
Networks (Oxford University Press, Oxford,2010).[22] G. Bianconi,
Multilayer networks: structure and function (Oxford University Press, Oxford, 2018).[23] A. Barrat, M. Barthelemy, and A. Vespignani,
Dynami-cal processes on complex networks (Cambridge universitypress, 2008).[24] S. N. Dorogovtsev,
Lectures on complex networks , Vol. 24(Oxford University Press, Oxford, 2010).[25] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, andA. Vespignani, Rev. Mod. Phys. , 925 (2015).[26] M. A. Porter and J. P. Gleeson, Frontiers in Applied Dy-namical Systems: Reviews and Tutorials (2016).[27] M. Nanni, G. Andrienko, C. Boldrini, F. Bonchi, C. Cat-tuto, F. Chiaromonte, G. Comand´e, M. Conti, M. Cot´e,F. Dignum, et al. , arXiv preprint arXiv:2004.05222(2020).[28] T. Gross, C. J. D. DLima, and B. Blasius, Physical re-view letters , 208701 (2006).[29] T. Gross and H. Sayama, in Adaptive networks (Springer,2009) pp. 1–8.[30] T. Mora and W. Bialek, Journal of Statistical Physics , 268 (2011).[31] M. A. Munoz, Reviews of Modern Physics , 031001(2018).[32] J. P. Gleeson and R. Durrett, Nature communications ,1 (2017).[33] E. Ben-Naim and P. L. Krapivsky, Phys. Rev. E ,050901 (2004).[34] E. Ben-Naim and P. Krapivsky, Eur. Phys. J. B , 145(2012). [35] S. Zapperi, K. B. Lauritsen, and H. E. Stanley, Physicalreview letters , 4071 (1995).[36] K. B. Lauritsen, S. Zapperi, and H. E. Stanley, PhysicalReview E , 2483 (1996).[37] M. Henkel, H. Hinrichsen, and S. L¨ubeck, Non-equilibrium phase transitions: Absorbing Phase Transi-tions , Vol. 1 (Springer, 2008).[38] H. Janssen, B. Schaub, and B. Schmittmann, Zeitschriftf¨ur Physik B Condensed Matter , 539 (1989).[39] H. Hinrichsen and G. ´Odor, Physical Review E , 311(1998).[40] F. Radicchi, C. Castellano, A. Flammini, M. A. Mu˜noz,and D. Notarmuzi, arXiv preprint arXiv:2002.11831(2020).[41] K.-I. Goh, D.-S. Lee, B. Kahng, and D. Kim, Physicalreview letters , 148701 (2003).[42] A.-L. Barab´asi and H. E. Stanley, Fractal concepts insurface growth (Cambridge university press, 1995).[43] J. Marro and R. Dickman,
Nonequilibrium phase tran-sitions in lattice models (Cambridge University Press,2005).[44] P. Krapivsky, (private communication) (2020).[45] D. T. Gillespie, Journal of computational physics , 403(1976). Appendix A: Stochastic SIR dynamics on wellmixed populations