key: cord-0631492-0codtgi3
authors: Radicchi, Filippo; Bianconi, Ginestra
title: Epidemic plateau in critical SIR dynamics with non-trivial initial conditions
date: 2020-07-29
journal: nan
DOI: nan
sha: b13528cdf3f6d3bb1b765208ca3cc2a27ec5d644
doc_id: 631492
cord_uid: 0codtgi3

Containment measures implemented by some countries to suppress the spread of COVID-19 have resulted in a slowdown of the epidemic characterized by time series of daily infections plateauing over extended periods of time. We prove that such a dynamical pattern is compatible with critical Susceptible-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 containment measures to an ongoing epidemic, however, has the effect to make the system enter in its critical regime with a number of infected individuals potentially large. We describe how such non-trivial starting conditions affect the critical behavior of the SIR model. We perform a theoretical and large-scale numerical investigation of the model. We show that the expected outbreak size is an increasing function of the initial number of infected individuals, while the expected duration of the outbreak is a non-monotonic function of the initial number of infected individuals. Also, we precisely characterize the magnitude of the fluctuations associated with the size and duration of the outbreak in critical SIR dynamics with non-trivial initial conditions. Far from heard immunity, fluctuations are much larger than average values, thus indicating that predictions of plateauing time series may be particularly challenging.

At the onset of the COVID-19 pandemic, world-wide time series of the number of infected individuals have displayed an exponential growth. Such a behavior is well predicted by standard epidemic frameworks [1] . In slightly later stages, however, time series have exhibited non-trivial dynamical patterns. Many papers have attempted to model observed behaviors and to determine the role of containment measures [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] . The common and reasonable assumption is that containment measures implemented in the attempt of mitigating the outbreak have strongly influenced the unfolding of the epidemic. Unfortunately, this a setting where modeling attempts are particularly challenging. The effective implementation of containment measures imposed by authorities rely on people personal judgements and adaptive behavior, and while epidemic spreading is a well-studied branch of mathematical biology [18] , statistical physics [19] and network science [20] [21] [22] [23] [24] [25] [26] , the modelling of adaptive behavior is only at its infancy [27] [28] [29] .

According to the data, in several countries, the slowdown of the epidemic spread is characterized by an almost flat time series of daily number of new infections. Moreover, the time series of the number of removed individuals display power-law growth instead of an exponential growth as a function of time [12, 14, 16, 17] . Here, we propose a theoretical interpretation of those features as the signature of the system being in (or near) its critical regime. Criticality is a fundamental property characterizing the dynamics of biological and socio-technical systems [30] [31] [32] . Our work consists of an in-depth investigation of a critical Susceptible-Infected-Removed (SIR) dynamics starting from a non-trivial initial configuration characterized by n 0 initially infected individuals. We interpret the emergence of the critical regime as the result of disease containment strategies, and the non-trivial initial condition as the configuration of the system when spreading becomes critical. In the typical setting considered in statistical mechanics [33, 34] , a single seed is generally used as the initial condition for critical SIR dynamics; the mapping of the critical SIR to the critical standard branching process allows for a full characterization of the spreading dynamics [35, 36] . The realistic assumption of having an initial number of infected individuals n 0 > 1 introduces an additional scale in the system affecting in a non-trivial manner the scaling properties of the SIR critical dynamics. While in other nonequilibrium systems a non-trivial initial condition may lead to a change of the critical exponent values [37] [38] [39] , in the critical SIR, the introduction of a non-trivial initial condition does not change the critical exponents that characterize the distribution of outbreak size and duration. However, it introduces lower exponential cutoffs in the distributions. As a result, the expected size and duration of the outbreak, as well as their standard deviations, have a non-trivial dependence on the initial condition n 0 . In this paper, we evaluate, by means of analytic arguments and large-scale simulations, the scaling of these quantities as functions of the population size N .

The paper is structured as follows: in Sec II we provide the theoretical interpretation of the plateau as a critical SIR dynamics starting from n 0 > 1 initial condition; in Sec III, by we perform a statistical mechanics investigations of the statistical properties of the critical SIR dynamics with non-trivial initial condition, supported by extensive numerical simulations of the process; finally in Sec IV we provide the concluding remarks. The Appendix describes the Gillespie algorithm used in this work to simulate the critical SIR dynamics.

We consider the Susceptible-Infected-Removed (SIR) model on a well-mixed population [19] [20] [21] [22] [23] [24] [25] [26] . At any point in time, individuals can be found in three possible states: susceptible, infected and removed. Susceptible individuals do not carry the disease but they can be infected; infected individuals carry the disease, and they can spread it to susceptible individuals; removed individuals are either removed or deceased, and they do no partecipate in the spreading dynamics. We indicate with λ the rate of infection, i.e., the expected number of spreading events occurring per unit of time. Without loss of generality, we set the recovery rate equal to one.

We start our discussion by focusing on the deterministic treatment of the SIR model on a well-mixed population with infinite size. If we indicate with s, i and r the fractions of susceptible, infected and removed individuals, respectively, we can write

Please note that s + i = r = 1. The critical dynamical regime is characterized by

If we start from an initial condition consisting of a fraction i(0) of infected individuals and a fraction r(0) = 0 of removed individuals, at the onset of the epidemic, i.e., t 1, we observe a different behavior depending on the value of λ. In the non-critical regime, i.e., λ = λ c , the deterministic equations for i and r read

Solutions of the above equations are

In essence, in the subcritical regime, i.e., λ < λ c , the number of infected individuals decays exponentially fast, and the number of removed individuals remains vanishing. In the supercritical regime, i.e., λ > λ c , the number of infected and removed individuals display an exponential increase. At criticality, i.e., λ = λ c , the deterministic equations for i and r are

leading to

Therefore, according to the deterministic approach, we should expect that the number of removed individuals at criticality increases linearly in time with a slope that is given by the initial condition i(0). The population size is N = 10 7 , and spreading is started from n0 = 1 seed. Different curves correspond to different values of λ = 1 + 2 −i with 2 ≤ i ≤ 9. b) Same as in panel a, but for λ = 1. Different curves correspond to different numbers of initially infected nodes n0 = 2 i , with 0 ≤ i ≤ 9. As λ approaches the critical value λc = 1 and n0 decreases toward one, we observe a plateauing of the time series.

From the deterministic Eqs. (1), it is evident that

The equation can be integrated to obtain the well-known solution [19] 

Using Eqs.(1), we can express the logarithmic derivative of the number of infected individuals as

where λ s(t) is the reproduction number. The former equations implies that the time series of infected individuals i(t) has a peak at t = t determined by

The fraction of susceptible individuals at the peak of the epidemic is given by s = s(t ) = 1/λ. By making the further assumption that the epidemic starts from a fraction i(0) = 1 − s(0) of infected individuals and zero removed individuals r(0) = 0 in Eq. (8), we obtain

Using Eqs. (10) and (11) in the first of Eqs. (1), we get

It follows that the second derivative of ln i is given by

where ρ is defined as

We note that ρ is zero, i.e., we reach a plateau, only for s(0) = 1 and λ = 1. This fact implies that, in the deterministic approach, a perfect plateau of the time series ln i is never achieved for λ > 1.

In the vicinity of the critical point, the time series of the infected individuals is still well described by a plateau. Developing the right hand side of Eq. (13) around λ = 1, s(0) = 1, we get s(0) 1 and

The above equation indicates that the conditions to have a near-plateau dynamics are having an infectivity rate λ close to one, and having the system as far as possible from heard immunity, i.e., 1 − s(0) 1. In summary, the near-critical state for λ 1 is a fragile state that can be characterized by a very slow dynamics if containment measures do not further decrease the infectivity λ below one (see Figure 1 ).

From now on, we assume that the system is in the critical regime. We further assume that spreading dynamics is started from n 0 > 1 initial seeds. The two assumptions serve to rationalize two main features of real time series. First, time series are characterized by long temporal windows of almost flat behavior. This is a signature of criticality. Second, plateaus are observed only after initial growths in the number of infected individuals, meaning that the critical regime is reached only after that containment strategies have effectively changed the spreading dynamics of the disease. Whereas critical properties of the SIR model are well understood for spreading processes initiated by n 0 = 1 individual, we are not aware of existing studies dealing with non-trivial initial conditions consisting of n 0 > 1 seeds. How do the properties of the critical dynamics change with n 0 ? What is the be- Please note that all the above questions cannot be answered with a purely deterministic approach. SIR outbreak sizes and durations obey probability distributions that are well peaked around their expected value only if the system is off criticality. However, the very fact that the system is assumed to be in the critical regime implies that fluctuations have a dominant role in the determination of the properties of the dynamical system. In Figure  2 for example, we display time series representative for the critical regime of the dynamics. Ground-truth time series are obtained by simulating the SIR stochastic dynamics (see Appendix A for details). They are compared with the deterministic expectation value obtained by integrating Eqs. (1) . We note that some realizations of the process are more persistent and more pervasive in the population than what predicted by the expected value.

From here on, we abandon the deterministic SIR equations and we embrace a stochastic approach. Critical SIR dynamics starting from a single initial seed, i.e., n 0 = 1, is known to be characterized by extremely large fluctuations of the outbreak size and duration. These fluctuations can be quantified by leveraging the mapping between critical SIR in a well-mixed population and the mean-field branching process. In the following sections, we first review results valid for n 0 = 1. Then, we focus our attention on the non-trivial case n 0 > 1.

A. Critical dynamics with n0 = 1 initial seed If the initial condition is such that only one node is in the infected state while all other nodes are in the susceptible state, the critical SIR model gives rise to outbreaks that follow the statistics of a critical branching process [35, 36] corrected by some scaling functions F T (N/N T ) and F R (N/N R ) that implement the effective cutoff caused by finite-size effects [33, 34] . Here, N is the size of the system; N T and N R are instead parameters that determine when the cutoff takes place. Specifically, the distribution P (T ) of the duration T of an outbreak follows the law

while the size of the outbreak R follows the distribution

The cutoff sizes N T and N R have been derived in Refs. [33, 34] . They are given by

From the expressions for P (T ) and P (R) given by Eqs. (16) and (17), respectively, and further assuming a sharp cutoff, it is easy to deduce that the scaling with the system size N of the average outbreak size R , the average duration T , and the standard deviations σ R and σ T [33, 34] obey

We observe that all the above quantities are subextensive, as they all grow sub-linearly with the system size. The expected critical outbreak size R grows as the system size to the power of 1/3. However, the standard deviation associated to the outbreak size, i.e., σ R , grows with increasing system size much faster than R . This fact indicates that it is very challenging to make predictions if the dynamics is critical. Similarly, the outbreak duration is characterized by large fluctuations in the large population limit. We note that the exponents 2 and 3/2 of the distribution P (T ) and P (R) are the critical mean-field exponents. These exponents are universal and are observed for many critical spreading processes [40] . They characterize the critical SIR on network topologies too as long as the underlying network has a homogeneous degree distribution. In power-law networks, these exponents can deviate from their mean-field values as investigated in Refs. [40, 41] . We have seen in Sec. II that the deterministic approach predicts a linear increase of the number of removed individuals with time. However, such a prediction is not accurate for the ground-truth dynamics; accounting for stochastic effects correctly predicts a quadratic growth of the number of removed individuals in time when the epidemic starts with a single initial seed. To this end, the number of removed individuals grows in time as a power-law

where z is a dynamical critical exponent, and t is the expectation value of the time necessary to observe R removed individuals. The value of the dynamical critical exponent can be obtained in different ways [36] . Here, we present the derivation of the value of the dynamical exponent based on Langevin-like equations for the dynamics.

Starting from an initial fraction of infected individuals i(0) = n 0 /N and a fraction r(0) = 0 of removed individuals we write di dt

where η(t) is an uncorrelated white noise with E(η(t)) = 0 and E(η(t)η(t )) = δ(t − t ) and c is a constant. At criticality, i.e., λ = λ c = 1, thus, assuming t 1 and

We now perform a simple scaling analysis of this stochastic equations as usually done in non-equilibrium statistical mechanics, e.g., Refs. [37, 42, 43] . If we rescale time as

and define the scaling exponents z, α, for r, i

as the exponents that leave the SIR critical dynamics unchanged. The SIR stochastic Eqs.

from which we can derive the scaling exponents

In summary, in the critical dynamical regime, if the spreading is started from a single initial seed, we expect that the average number of removed individuals grows quadratically with time.

B. Critical dynamics with n0 > 1 initial seeds Critical SIR dynamics started from the non-trivial initial condition n 0 > 1 differs from the critical SIR dynamics started from n 0 = 1 seed. To include an explicit dependence on the parameter n 0 in the scaling of Eq. (20), we correct it by introducing the scaling function F (x), where x = n 0 /t(R). We impose that

with

According to the deterministic SIR equations for t 1 and i(0) 1, the number of removed individuals grows linearly in time with a slope i(0) = n 0 /N . Thus, we deduce that

This value is well supported by extensive numerical results (see Figure 3 ) which confirm that there is a crossover between linear and quadratic dependence of R on t(R).

In the simulation of the SIR model, it is natural to study the behavior of R as a function of t(R). However, in real epidemic time series, the number of infected individuals is measured over constant time intervals. The two ways of monitoring the evolution of the process, i.e., R versus t(R) rather than R(t) versus time t, may lead to the observation of different scaling exponents. The discrepancy is due to the stochastic nature of the spreading process. The phenomenon is apparent from the results of Figure 3 : depending on the type of measurement performed on the system, the power-law increase of the number of removed individuals as a function of time can be described by a continuous range of exponents ranging from ξ = 1 to ξ ∼ 2.5. We can therefore write

where ξ is a decreasing function of n 0 , and h(t, N ) is a modulating function expressing the deviation from the pure power-law behavior. The ansatz of the above equation is compatible with the power-law scaling of the empirical time series of removed individuals as a function of time observed in countries where containment measures have been implemented extensively [12, 14, 16, 17] .

C. Distribution of avalanche durations and sizes for the critical SIR model initiated by n0 > 1 seeds

In this section we investigate the statistical properties of the distribution of outbreak duration and size for the critical SIR dynamics in a well mixed-population when the initial condition is non-trivial, i.e., n 0 > 1. Scaling arguments suggest the following expression for the distribution P (T ) of the critical outbreak duration T P (T ) ∼ n 0 T −2 F T (N/N T , n 0 /T ). (33) The above scaling function is a natural modification of Eq. (16) by assuming that n 0 scales like time. In particular, the distribution P (T ) is characterized by a lower cutoff depending on n 0 . This fact is intuitive as an outbreak with a larger number of initially infected individuals is not expected to reach the absorbing state faster than an outbreak started by a single seed (see Figure 4a ). We note that, in the critical SIR dynamics, the dependence on n 0 does not lead to a change of the critical exponent values, as for example observed in other non-equilibrium phase transitions [37] [38] [39] . In Figure 5a , we display the function

and we demonstrate that the scaling function w T (n 0 /T ) for n 0 N T can be approximated as

The scaling behavior, valid for n 0 N 1/3 , can be justified by assuming that each of the n 0 seeds generates an independent outbreak obeying the statistics of the critical branching process. A critical avalanche started from a single infected individual has a duration T following the power-law distribution π(T ) ∼ T −2 [35, 36] . Thus, assuming independence among the n 0 avalanches, we can estimate the probability P (T ) as the probability that among all n 0 outbreaks the last outbreak to get extinguished is extinguished at time T . Therefore in the infinite population limit we obtain

where π c (T ) is the probability that an outbreak generated by a single infected individual is not extinguished at time T , with π c (T )

Finally we note that while the scaling behavior described in Eq. (35) has strong numerical confirmation for T N T for values of T ∼ N T the scaling function w T (n 0 , T ) signals a dependence of the cutoff on n 0 (see Figure 5 ). Scaling arguments suggest that the distribution P (R) of critical outbreak size R should obey

where the function F R (N/N R , n 0 , R) implements a lower cutoff dependent exponentially on n 0 (see Figure 4b ). The distributions display a lower cutoff that increases as n0 increases, and an upper cutoff whose value does not strongly depend on n0. In Figure 5b , we show the function

which, for n 0 N R (R, n 0 ), can be approximated as

This scaling function indicates that, for large values of R, R scales like n 2 0 . For small values of R, it is possible to observe some corrections would be required to fully describe the scaling. We notice that, in the first order in n 0 , the normalization constant of the distribution P (R) is independent of n 0 . A way to interpret the result is by considering the infinite population limit approximating the distribution P (R) as the convolution of the n 0 sizes of independent outbreak events. In this limit we have

where F(ω) = r Π(r)e −iωr is the generating function of the distribution Π(r) of avalanches sizes of SIR critical dynamics starting from a single seed. Assuming in first approximation that Π(r) is a pure power-law Π(r) ∼ r −3/2 it follows that the logarithm of the generating function behaves, for small ω, as ln F(ω) ∼ √ ω. The result, together with Eq. (41) , indicates that R should scale as n 2 0 for n 0 N 2/3 . For more more details on the infinite population limit we refer the reader to Ref. [44] . We performed large-scale simulations of the critical SIR model to address fundamental questions regarding the distributions of duration T and size R of outbreaks started by a non-trivial initial condition n 0 > 1. In Figure 6 , we display the average values and the standard deviations of both T and R for a large system composed of N = 10 8 individuals. We display the moments of the distributions as a function of the number of initial seeds n 0 . The main outcomes are as follows. The expected size R is a growing function of n 0 . The expected duration T is a non-monotonic function of n 0 , displaying a single peak. The standard deviations σ T and σ R also display a peak as a function of n 0 . The coefficient of variation σ T / T and σ R / R are monotonically decreasing with n 0 . Therefore we conclude that the role of critical fluctuation is very relevant for the critical dynamics, albeit in relative terms is most severe for the SIR dynamics starting from a single initial seed.

We make the ansatz that the average duration T can be described by

where the function G(n 0 |α, β, H, K) is given by

Here, α and β are, in the large population limit, independent of N . On the contrary, H and K are dependent on the population size. By introducing the function

we observe that it is possible to rescale the curves obtained for different values of N by performing the transformation H −1 K α/β G(n 0 |α, β, H, K) = g(n 0 K 1/β |α, β). (45) This expression allows us to perform a data collapse of the data obtained for T at different values of n 0 and different population size N (see Figure 7a ). We observe that, if we start from a non-trivial initial condition n 0 , the expected duration of the outbreak T a maximum at

The scaling parameters H and K obey the scaling relation

with a H = 0.58 ± 0.08, b H = −2.0 ± 0.9 δ K = −0.37 ± 0.04 a K = 0.75 ± 0.06 (see Figure 9 ). We note that the logarithmic scaling of H is expected from the known scaling of T for the SIR critical model starting from a single initial seed. The exponents α and β are given by α = 0.78 ± 0.03, β = 1.10 ± 0.1.

The standard deviation of the outbreak duration σ T can be described in the same exact way as T . The ansatz

leads to the data collapse shown in Figure 7b . The scaling parameters H and K obey the scaling relations

with δ H = 0.17 ± 0.01, b H = 1.7 ± 0.2 δ K = −0.36 ± 0.03 a K = 2.3 ± 0.7 (see Figure 9 ). We note that δ H 1/6. Therefore, for n 0 = 1 the scaling reduces to the wellknown scaling for the critical SIR model starting from a single initial seed. Moreover, the exponent α and β for σ T are given by α = 0.50 ± 0.05, β = 1.0 ± 0.1.

(51) 

The ansatz for R is slightly different from the one appearing in Eq. (43), as it includes an additional logarithmic correction

We take

and α = 1, β = 0.5. and we perform a fit of the exponent γ. As Figure 9a and Figure 10a demonstrate, the function gives rise to excellent data fits as long as the exponent γ is

with a γ = 0.053 ± 0.003 and b γ = 0.30 ± 0.06. The function R can be rescaled and the data obtained for different N collapsed (see Figure 10 ). This is done by noting that R y(n 0 ) = g(x(n 0 )|α, β),

where x(n 0 ) = n 0 (ln n 0 ) γ/β K 1/β ,

and the function g(x|α, β) is given by Eq. (44) . The expected size of the outbreak R does not display a maximum as a function of n 0 , i.e., it is a monotonous increasing function of n 0 . The standard deviation σ R can be instead fitted using the same ansatz as σ T , i.e.

σ R G(n 0 |α, β, H, K).

The best estimates of the parameters are α = 0.57 ± 0.03, β = 0.95 ± 0.05,

and

with δ H = 0.48 ± 0.02, b H = 1.8 ± 0.5 δ K = −0.25 ± 0.03 a K = 0.49 ± 0.3 (see Figure 9 ). The corresponding data collapse is shown in Figure 10b . Figure 7 . b) Same as in panel a, but for the standard deviation σR.

Motivated by the current COVID-19 pandemic, we have investigated the critical properties of the Susceptible-Infected-Removed (SIR) dynamics in well mixed-populations starting from non-trivial initial conditions consisting of n 0 > 1 infected individuals. Although the modeling framework oversimplifies the real-world scenario, the setting is realistic in two main respects. First, the plateauing time series observed in empirical data are compatible with the critical dynamical regime. Second, the initial condition n 0 > 1 is representative for a critical regime reached, thanks to effective containment measures, after a significant community transmission took already place. We have shown that a non-trivial initial condition n 0 > 1 introduces another typical scale on the dynamics inducing a lower cutoff in the distribution of the duration and size of critical outbreaks. The critical dynamics is characterized by very strong fluctuations, but the presence of a non-trivial initial condition mitigates the role of the fluctuations. In particular, while for a single initial seed the standard deviation on the outbreak size and duration is much larger than the corresponding expectation values, the relative error diminishes as the size n 0 of the initial seed set increases. Moreover, numerical results indicate that, as the initial number of infected individuals n 0 increases, the expected size of the outbreak increases while the expected duration first increases and then decreases, displaying a maximum. Using scaling arguments, and extensive numerical simulations we have deduced the scaling of the maximum duration and the corresponding number of initially infected individuals.

Infectious diseases of humans: dynamics and control

Mathematical biology: I. An introduction

A kinetic view of statistical physics

Multilayer networks: structure and function

Dynamical processes on complex networks

Adaptive networks

Nonequilibrium phase transitions: Absorbing Phase Transitions

Fractal concepts in surface growth

Nonequilibrium phase transitions in lattice models

The critical dynamics is obtained for λ = 1. The dynamics is iterated until the number of infected individuals is zero

We thank P.L. Krapivsky, Geza Odor and R. M. Ziff for interesting discussions. F.R. acknowledges support from the National Science Foundation (CMMI-1552487).

The critical SIR dynamics in a well-mixed population of N individuals is simulated with the following implementation of the Gillespie algorithm [45] . We indicate with S(t), I(t) and R(t) respectively the number of susceptible, infected and removed individuals as a function of time t. We start from the initial condition of I(0) = n 0 , S(0) = N − n 0 and R(0) = 0. At each elementary step, the algorithm proceeds as follows:(i) Time increases by the amount ∆twhere ∆t is given by