Two golden times in two-step contagion models

TL;DR

This paper identifies two types of golden times in two-step contagion models, scaling differently with system size, and explores their universality and dependence on initial conditions and model specifics.

Contribution

It uncovers two distinct scaling behaviors of the golden time in two-step contagion models and analyzes their dependence on initial conditions and network size.

Findings

Two types of golden times scale as O(N^{1/3}) and O(N^{4})

The 1/3 exponent is universal across models with discontinuous transitions

The 4 exponents are model-dependent and close to 1/4

Abstract

The two-step contagion model is a simple toy model for understanding pandemic outbreaks that occur in the real world. The model takes into account that a susceptible person either gets immediately infected or weakened when getting into contact with an infectious one. As the number of weakened people increases, they eventually can become infected in a short time period and a pandemic outbreak occurs. The time required to reach such a pandemic outbreak allows for intervention and is often called golden time. Understanding the size-dependence of the golden time is useful for controlling pandemic outbreak. Here we find that there exist two types of golden times in the two-step contagion model, which scale as $O(N^{1/3})$ and $O(N^{\zeta})$ with the system size $N$ on Erd\H{o}s-R\'enyi networks, where the measured $\zeta$ is slightly larger than $1/4$. They are distinguished by the initial number of infected nodes, $o(N)$ and $O(N)$, respectively. While the exponent $1/3$ of the $N$-dependence of the golden time is universal even in other models showing discontinuous transitions induced by cascading dynamics, the measured $\zeta$ exponents are all close to $1/4$ but show model-dependence. It remains open whether or not $\zeta$ reduces to $1/4$ in the asymptotically large-$N$ limit.

Equations

$$\mathrm{S}+\mathrm{I}$$

$$\mathrm{W}+\mathrm{I}$$

$$\mathrm{I}$$

$$P_S(\ell )=(1-\kappa -\mu {)}^{\ell }.$$

$$P_W(\ell )=\sum _{n=0}^{\ell -1}(1-\kappa -\mu {)}^n\mu (1-\eta {)}^{\ell -n-1}.$$

$$m=i_{0}+\sum_{q=1}^{\infty}P_{d}(q)\Big{(}s_{0}f_{q}(u)+w_{0}g_{q}(u)\Big{)},$$

$$f_q(u)=\sum _{\ell =1}^q\left(\genfrac{}{}{0ex}{}{q}{\ell }\right)u^{\ell }(1-u{)}^{q-\ell }{P}_{IR}(\ell )$$

$$g_q(u)$$

$$u_{n+1}=i_{0}+\sum_{q=1}^{\infty}\dfrac{qP_{d}(q)}{z}\Big{(}s_{0}f_{q-1}(u_{n})+w_{0}g_{q-1}(u_{n})\Big{)},$$

$$u_{n+1}$$

$$P_{IR}^{(S\to I)}(\ell )=\sum _{m=0}^{\ell -1}(1-\kappa -\mu {)}^m\kappa =\frac{\kappa }{\kappa +\mu }[1-(1-\kappa -\mu {)}^{\ell }].$$

$$\displaystyle F^{(S\to I)}(u_{n})=\dfrac{\kappa}{\kappa+\mu}\Big{[}1-e^{-(\kappa+\mu)zu_{n}}\Big{]}.$$

$$u_{n+1}^{(S\to I)}=z\kappa u_n-a{u}_n^2+O(u_n^3),$$

$$F^{(S\to W\to I)}(u_n)=F(u_n)-F^{(S\to I)}(u_n)$$

$$F^{(S\to W\to I)}(u_n)$$

$$u_{n+1}^{(S\to W\to I)}=b{u}_n^2+O(u_n^3),$$

$$u_{n+1}=u_n+(b-a)u_n^2+O(u_n^3).$$

$$u_{n+1}=i_0+(1-i_0)F(u_n)$$

$$G(u)\equiv i_0+(1-i_0)F(u)-u=0$$

$$h_{n+1,\ell}=i_{\ell}+\sum_{q=0}^{\infty}\dfrac{z^{q}e^{-z}}{q!}\Big{(}s_{\ell}f_{q}(h_{n,\ell})+w_{\ell}g_{q}(h_{n,\ell})\Big{)},$$

$$s_{\ell }=s_0{e}^{-(\kappa +\mu )z{u}_{\ell -1}}$$

$$w_{\ell }=1-u_{\ell }-s_{\ell }.$$

$$h_{n+1,n^{\ast }}$$

$${ϵ}_{n+1}=d_0+(1+\delta d_1){ϵ}_n+(d_2+\delta d_2){ϵ}_n^2+O({ϵ}_n^3),$$

$$d_0$$

$$\delta d_1$$

$$d_2$$

$$\delta d_2$$

$$\dot{\epsilon}=d_{0}+(d_{2}+\delta d_{2})\bigg{(}\epsilon+\dfrac{\delta d_{1}}{2(d_{2}+\delta d_{2})}\bigg{)}^{2}-\dfrac{(\delta d_{1})^{2}}{4(d_{2}+\delta d_{2})}.$$

$$\dot{{ϵ}^{\prime }}=d_0+{d_2}^{\prime }{ϵ}^{\prime 2}$$

$$\mathcal{T}={\int }_{-\infty }^{\infty }\frac{d{ϵ}^{\prime }}{d_0+{d_2}^{\prime }{ϵ}^{\prime 2}}\sim \frac{\pi }{\sqrt{d_0}},$$

$$⟨n_c⟩\sim N^{0.252\pm 0.001}$$

Full text

I Introduction

Epidemic spread of diseases and rumors and their control and containment have become a central issue in recent years as the real world becomes “smaller.” It is a general observation that there is a slow phase in the spreading process before the sudden pandemic outbreak review_epidemics . This slow period is called golden time as it allows for intervention, which is much more difficult after the disease becomes global. Modeling of epidemic spread with essential factors is necessary to control catastrophic outbreaks within this golden time. To this end, several epidemic models have been investigated on complex networks, for instance, the susceptible–infected–removed (SIR) model sir ; sir_newman and the susceptible–infected–susceptible (SIS) model sis . Analytical and numerical studies of those models revealed that a continuous phase transition occurs on Erdős-Rényi (ER) random networks ER . Thus, abrupt pandemic outbreaks on a macroscopic scale, which often occur in the real world, cannot be reproduced using those models.

Considerable effort has been devoted recently to construct mathematical models that exhibit a discontinuous epidemic transition at a finite transition point on complex networks. A natural way is to appropriately extend the conventional SIR and SIS models. For instance, an extended SIR model includes more than one infected state of different pathogens that are cooperatively activated in contagion: A person who is suffering from the flu can be more easily infected by pneumonia. This model is referred to as a cooperative contagion model grassberger_nphy . Similar instances include a two-step contagion process. A patient becomes weakened first and then becomes sick. This model is referred to as the susceptible–weakened–infected–removed (SWIR) model janssen ; janssen_spinoal ; grassberger_pre ; hasegawa1 ; hasegawa3 ; chung ; choi_2016 ; choi_multi . In another instance of modified SIR models, a network evolves by rewiring links at a certain rate during the spread of contagion sir_rewire . The rewiring takes into account the mobility of humans. Then, epidemic spread can be accelerated as the rewiring rate is increased, which can lead to a discontinuous transition representing the pandemic outbreak.

When diseases spread, we need to keep susceptible people separate from infected patients or vaccinate the susceptible people before the diseases spread on a macroscopic level. A recent study universal showed that for the SWIR model on ER networks, a system exhibits a long latent period (called a golden time) within which measures can be taken, beyond which the disease spreads explosively over the system at a macroscopic level. Estimating the golden time is important for the prevention of pandemic outbreaks. Moreover, it is necessary to get early-warning signals if a critical threshold is approached indicator .

It was revealed choi_2016 ; universal that when a disease starts spreading from a single node, the golden time $n_{c}$ scales as $n_{c}(N)\sim N^{\zeta}$ with $\zeta=1/3$ at the epidemic threshold. Here we reconsider this problem and represent the pattern of disease transmission using a nonlinear mapping. We show that the linear and nonlinear terms of the nonlinear mapping separately behave dynamically well. The linear term is responsible for one-step contagion without weakened states and the nonlinear term describes the two step contagion, which includes weakened state. Thus, the previous result of $N^{1/3}$ for the golden time is consistent with the characteristic size of the giant cluster generacolorted in the SIR model bennaim , thus it has got verified within this new framework. Next, we consider another case, which is the main concern of this paper, in which an epidemic starts to spread from endemic multiple seeds of $O(N)$ on ER networks also at the epidemic threshold. In this case, long latent period appears not immediately but after some characteristic time. Thus fluctuations induced by the stochastic process of disease transmission in the early time heavily affect the behavior during the latent period, which changes the measured exponent $\zeta$ to a value slightly larger than 1/4. We estimate this scaling behavior using the saddle-node bifurcation theory book and discuss the underlying mechanism.

Similar size dependences of mean cascading time at a transition point were studied for other cascade dynamics models such as $k$-core percolation kcore1 ; kcore2 ; kcore3 ; kcore4 ; kcore5 ; kcore_goltsev and cascading failure model on interdependent network (CFoIN) buldyrev ; zhou ; son ; baxter_mcc ; bianconi . It was found buldyrev ; zhou that in the CFoIN, the mean cascading time is proportional to $N^{1/3}$ or $N^{1/4}$ depending on the way of choosing the transition points. Refs. kcore4 ; kcore5 showed that the exponent 1/3 is also obtained in $k$-core percolation. Thus, the scaling behavior of $N^{1/3}$ is robust. However, for the $\zeta>1/3$ case, a different scaling behavior with $\zeta\approx 0.280$ grassberger_mcc was numerically obtained for a surface growth model effectively equivalent to the CFoIN.

Here we extend our formalism of nonlinear mapping used in the SWIR model to other models such as $k$-core percolation and the threshold model watts ; dodds . We show that when the cascade starts from a fixed number of multiple seeds $O(N)$, the golden times for both models also become proportional to $N^{\zeta}$, where $\zeta$ are estimated to be slightly larger than 1/4 within our simulation range and those values are different to each other, suggesting non-universal behavior. However, we cannot exclude the possibility $\zeta=1/4$ in large-$N$ limit. We shall discuss this point in Sec. IV.

This paper is organized as follows: We first introduce the SWIR model and set up the evolution equation of the epidemic dynamics in Sec. II. Next, we derive a nonlinear mapping for the epidemic spread from a single seed in Sec. IIIA. We show that the roles of the linear and nonlinear terms are well separated. In Sec. IIIB, we derive a similar nonlinear mapping for the multiple-seed case, and show how the multiplicative feature of the fluctuations of epidemic spreading affects scaling of the golden time. In Sec. IV, we obtain the golden times of the multiple-seed case for $k$-core percolation and the threshold model and show that the numerical values of $\zeta$ are slightly larger than $1/4$. We also discuss the possibility of $\zeta=1/4$ in the thermodynamic limit. In Sec. V, we discuss the origin of the puzzle in view of nonlinear dynamics theory. A summary is presented in Sec. VI.

II The SWIR model

The SWIR model is a generalization of the SIR model by including two sates, a weakened state (denoted as $W$) and an infected state ($I$), between the susceptible state ($S$) and recovered state ($R$), instead of a single infected state $I$ alone, as in the SIR model. Nodes in state $W$ are involved in the reactions $S+I\to W+I$ and $W+I\to 2I$, which occur in addition to the reactions $S+I\to 2I$ and $I\to R$ in the SIR model. At each discrete time step $n$, the following processes are performed. (i) All the nodes in state $I$ are listed in random order. (ii) The states of the neighbors of each node in the list are updated sequentially as follows: If a neighbor is in state $S$, it changes its state in one of two ways: either to $I$ with probability $\kappa$ or to $W$ with probability $\mu$. If a neighbor is in the state $W$, it changes to $I$ with probability $\eta$, where $\kappa$, $\mu$, and $\eta$ are the contagion probabilities for the respective reactions. (iii) All nodes in the list change their states to $R$. This completes a single time step, and we repeat the above processes until the system reaches an absorbing state in which no infectious node is left in the system. The reactions are summarized as follows:

[TABLE]

In an absorbing state, each node is in one of three states, the susceptible, weakened, or recovered state. We define $P_{S}(\ell)$ as the conditional probability that a node remains in state $S$ in the absorbing state, provided that it has $\ell$ neighbors in state $R$ and was originally in state S. This means that the node remains in state $S$ even though it has been in contact $\ell$ times with these $\ell$ neighbors in state $I$ before they change their states to $R$. Thus, we obtain

[TABLE]

Next, $P_{W}(\ell)$ is similarly defined as the conditional probability that a randomly selected susceptible node is in state $W$ after it contacts $\ell$ neighbors in state $I$ before they change their states to $R$. The probability $P_{W}(\ell)$ is given as

[TABLE]

Finally, $P_{IR}(\ell)$ is the conditional probability that a node has been infected in any state, either $I$ or $R$, provided that it was originally in state $S$ and its $\ell$ neighbors are in state $R$ in the absorbing state. Using the relation $P_{S}(\ell)+P_{W}(\ell)+P_{IR}(\ell)=1$, one can determine $P_{IR}(\ell)$ in terms of $P_{S}$ and $P_{W}$.

On a network with a degree distribution $P_{d}$, we consider the case in which the initial densities of susceptible, weakened, and infectious nodes are given as $s_{0}$, $w_{0},$ and $i_{0}$, provided that $s_{0}+w_{0}+i_{0}=1$. The order parameter $m$, the density of nodes in state $R$ after the system falls into an absorbing state, is given using the local tree approximation as

[TABLE]

where

[TABLE]

[TABLE]

and $u$ is the probability that an arbitrarily chosen edge leads to a node in state $R$ or $I$ but not infected through the chosen edge in the absorbing state. We define $u_{n}$ similarly to $u$ but at time step $n$. The probability $u_{n+1}$ can be derived from $u_{n}$ as follows:

[TABLE]

where $z\equiv\sum_{q}qP_{d}(q)$ is the mean degree of the network and the factor $qP_{d}(q)/z$ is the probability that a node connected to a randomly chosen edge has degree $q$. As $n\to\infty$, $u_{n}$ converges to $u$.

III Golden times in the SWIR model

III.1 The single-seed case

First, we consider the case in which the initial number of infectious nodes is $o(N)$; that is, $i_{0}=w_{0}=0$ and $s_{0}=1$ in the thermodynamic limit. In this case, the SWIR model exhibits a mixed-order transition choi_2016 at a transition point $\kappa_{c}$ when the mean degree is larger than a critical value. The order parameter displays a discontinuous transition from $m(\kappa_{c})=0$ to $m_{0}$, whereas other physical quantities such as the outbreak size distribution exhibit a critical behavior. The behavior of the order parameter $m(\kappa)$ as a function of $\kappa$ is schematically shown in Fig. 1(a).

We are interested in how infected nodes spread as a function of the cascade step $n$ when the order parameter jumps. As a particular case, when the network is an ER network having a degree distribution that follows the Poisson distribution, i.e., $P_{d}(q)=(q+1)P_{d}(q+1)/z=z^{q}e^{-z}/q!$, where $z$ is the mean degree, Eq. (10) is reduced as follows:

[TABLE]

We remark that on ER networks, $u_{n}$ in the limit $n\to\infty$ becomes equivalent to $m$ obtained from Eq. (7).

We pick up the contribution of the reaction $S+I\to 2I$ from Eq. (11) but neglect the contribution of the reaction $W+I\to 2I$. Then, the probability that a node becomes directly infected by $\ell$ infectious neighbors, which is denoted by $P_{IR}^{(S\to I)}(\ell)$, is given as

[TABLE]

Applying the formula for the Poisson degree distribution to Eq. (10), we obtain that

[TABLE]

Because the order parameter increases from $m=0$, we assume that $u_{n}$ is small in the early time regime. Thus,

[TABLE]

where $a\equiv\kappa(\kappa+\mu)z^{2}/2$. Actually, the coefficient $z\kappa$ of the first-order term is the mean branching ratio in the early time regime. When the critical branching (CB) process occurs, the mean branching ratio becomes unity, so the transition occurs at $\kappa_{c}=1/z$. On the other hand, the discrete mapping (14) at $\kappa_{c}$ may be rewritten in the form of a saddle-node bifurcation, ${\dot{u}}^{(S\to I)}=-au^{2}$, where $u$ is a function of the continuous time variable $n$ and the overdot denotes differentiation with respect to it. Because $a>0$, $u^{*}=0$ is a stable fixed point for $u\geq 0$, and this point represents the fixed point of the SIR model, indicating a second-order transition.

Next, we consider the two successive reactions $S+I\to W+I$ and $W+I\to 2I$, in which a susceptible node becomes infected in two steps and eventually recovers. Because a node can be infected either by the reaction $S+I\to 2I$ or by the reactions $S+I\to W+I$ and $W+I\to 2I$, the probability $f^{(S\to W\to I)}(u_{n})$ can be obtained using the relation

[TABLE]

as

[TABLE]

Again, using $u_{n}\ll 1$, we obtain that

[TABLE]

where $b\equiv{\mu\eta}z^{2}/2$. Here we note that the first-order term $O(u_{n})$ is absent. Combining Eqs. (14) and (17), we obtain that

[TABLE]

Thus, ${\dot{u}}=(b-a)u^{2}$. When $b-a<0$, i.e., $\mu\eta<\kappa_{c}^{2}+\kappa_{c}\mu$, the fixed point $u^{*}=0$ is stable, and thus a continuous transition occurs. Otherwise, the fixed point $u^{*}=0$ is unstable, and a discontinuous transition occurs. The condition $\mu\eta>\kappa_{c}^{2}+\kappa_{c}\mu$ for a discontinuous transition is consistent with previously obtained results grassberger_pre ; choi_2016 .

When contagion starts from a single infectious node, its spread in the early time regime is governed by the linear term of Eq. (18). It proceeds in the form of a CB tree universal , i.e. the mean branching ratio $(u_{n+1}-u_{n})/(u_{n}-u_{n-1})$ is almost unity, and the main contribution is that of the reaction $S+I\to 2I$. Thus in the thermodynamic limit, $u_{n}$ always stays zero so that nonlinear terms in Eq. (18) do not appear. On the other hand, in finite systems, $u_{n}$ grows gradually and the nonlinear term $(b-a)u_{n}^{2}$ becomes significant after a characteristic time $n_{c}(N)$. It was argued in bennaim that for the SIR model at the epidemic threshold, the maximum size of outbreaks is proportional to $N^{2/3}$ in the mean field limit. When $u_{n}$ grows up to $O(N^{2/3})$, the nonlinear terms in Eq. (14) suppresses further growth of the cluster, leading to a subcritical branching process. This means that the CB process driven by Eq. (14) persists up to $O(N^{1/3})$, because the fractal dimension of the CB tree is two. On the other hand, for the SWIR model, the coefficient of the nonlinear term (18) is positive, and the nonlinear term enhances further increase of removed nodes. The CB process turns into a supercritical process, leading to a pandemic outbreak. Accordingly, the golden time, the duration of the CB process, scales similarly as $\sim N^{1/3}$ to that of the SIR model, which is what we observed in a previous work choi_2016 ; universal .

III.2 The multiple-seed case

Next, when the number of infectious nodes is $O(N)$, i.e., $i_{0}>0$, $s_{0}=1-i_{0}$, and $w_{0}=0$ in the thermodynamic limit, it was shown janssen_spinoal ; hasegawa1 ; hasegawa3 ; choi_multi that there exists a critical value $i_{0}^{(c)}$ such that when $i_{0}<i_{0}^{(c)}$, a hybrid phase transition occurs at a transition point $\kappa_{c}$, whereas when $i_{0}=i_{0}^{(c)}$, a continuous transition occurs. Here we focus on the former case.

In the multiple-seed case, Eq. (10) becomes

[TABLE]

when the network is an ER network with mean degree $z$. Fixed points of Eq. (19) satisfy the equation,

[TABLE]

and the smallest solution among them is the order parameter $m(\kappa)$. We note that $G(u)$ contains the parameters $\kappa$, $\mu$, $\nu$, $z$ and $i_{0}$. As already shown in the single-seed case, for appropriately given values of $\mu$ and $\eta$, $G(0)=i_{0}$, $G^{\prime}(0)=(1-i_{0})(z\kappa-1)$ and $G^{\prime\prime}(0)=(1-i_{0})(b-a)>0$. Thus when $i_{0}$ is sufficiently small, i.e., $i_{0}<i_{0}^{(c)}$, $m(\kappa)$ satisfies G^{\prime}\big{(}m(\kappa)\big{)}<0 for values of $\kappa$ near zero. Then $m(\kappa)$ increases continuously as $\kappa$ is increased untill $\kappa$ reaches a critical value $\kappa_{c}$ such that G\big{(}m(\kappa_{c})\big{)}=G^{\prime}\big{(}m(\kappa_{c})\big{)}=0 and G^{\prime\prime}\big{(}m(\kappa_{c})\big{)}>0 are satisfied. We note that $\kappa_{c}$ depends on $i_{0}$ and $z$. When $i_{0}=0$, $\kappa_{c}$ is reduced to $1/z$, the transition point of the single-seed case. $m(\kappa)$ exhibits a critical behavior as $\kappa$ approaches $\kappa_{c}$ and subsequently jumps from $m(\kappa_{c})=m_{d}$ to another value $m_{u}$ as represented in Fig. 1(b). Thus the transition is hybrid.

We notice that at a transition point for the multiple-seed case, an infected node can be in contact with a node that was weakened by a different infectious root choi_multi . Accordingly, the reaction $W+I\to 2I$ can occur even in the early time regime, as shown in Fig. 2(c) with red zig-zag (lowest) curve. Moreover, the CB process appears not from the beginning but slightly after that indicated by an arrow at $n^{*}$ in Fig 2(c), at which the density of recovered nodes $r_{n^{*}}$ is close to $m_{d}$ indicated in Fig. 1(b). From this step $n^{*}$, $r_{n}$ remains almost constant for a long time as shown in Fig. 2(d).

However, in finite systems, due to the fluctuations arising in the stochastic process of epidemic spread, the densities of each species of nodes at $n^{*}$ can be different for each realizaton. Those fluctuations affect $n_{c}$, which can be also different for different realizations, where $n_{c}$ is the golden time, from which $r_{n}$ increases drastically.

We denote the densities of each species of nodes at a certain time step $\ell$ as $s_{\ell},w_{\ell},i_{\ell},$ and $r_{\ell}$, respectively. Then on ER networks, for $n>\ell$, $h_{n,\ell}\equiv u_{n}-r_{\ell}$ satisfies

[TABLE]

where

[TABLE]

and

[TABLE]

Moreover, using Eq. (10), the relation $i_{n}=u_{n+1}-u_{n}$, and $u_{n}=i_{n}+r_{n}$, we determine $i_{\ell}$ and $r_{\ell}$.

We focus on the density fluctuations of each species at $n^{*}$. We split the densities of each species of nodes into two parts: $x_{n^{*}}+\delta x_{n^{*}}$ ($x=s,w,i$ or $r$), where the first term represents densities of $x$-species nodes in thermodynamic limit at $n^{*}$, and the second one is the deviation. Then for $n>n^{*}$, Eq. (21) becomes

[TABLE]

We did not take into account the density fluctuations induced after $n^{*}$, because they are negligible compared to those at $n^{*}$. At $\kappa=\kappa_{c}$, Eq. (24) has a nontrivial fixed point $h_{d,n^{*}}=m_{d}-r_{n^{*}}$ in the thermodynamic limit. Then Eq. (24) is rewritten with $\epsilon_{n}=u_{n}-m_{d}$ as

[TABLE]

where

[TABLE]

Neglecting higher order terms of $\epsilon$, Eq. (25) is rewritten in an alternative form,

[TABLE]

Because $\delta i_{n^{*}}\sim\delta s_{n^{*}}\sim\delta w_{n^{*}}\sim\delta r_{n^{*}}\ll 1$ for large $N$, the last term can be neglected compared to $d_{0}$ and Eq.(30) is rewritten simply as

[TABLE]

where ${d_{2}}^{\prime}=d_{2}+\delta d_{2}$ and $\epsilon^{\prime}=\epsilon+\delta d_{1}/2{d_{2}}^{\prime}$. We note that $d_{0}$ is a real number, while $d_{2}^{\prime}$ is a positive number.

The nonlinear mapping Eq.(31) includes several features: When $d_{0}<0$, $\epsilon^{\prime}$ reaches a fixed point $\epsilon_{*}^{\prime}=-\sqrt{|d_{0}|/d_{2}^{\prime}}$; when $d_{0}=0$, $\epsilon^{\prime}$ remains at zero; when $d_{0}>0$, there arises the so-called bottleneck effect at $\epsilon^{\prime}=0$ book ; kcore3 ; buldyrev . The time step to pass through the bottleneck is calculated as

[TABLE]

which is approximately the time interval of the plateau region, i.e, $n_{c}-n^{*}$. Because $n^{*}$ is much smaller than $n_{c}$, $n_{c}\approx\cal{T}$, which is the golden time for a single realization of the process. $d_{0}$ can have different values for different realizations, yielding a different $n_{c}$. Thus, we need to take average of $n_{c}$ over different realizations to obtain $\langle n_{c}\rangle$.

We performed extensive numerical simulations at the transition point $\kappa_{c}\approx 0.11494875096512$ of the SWIR model starting from multiple seeds $i_{0}=0.002$, and obtained that

[TABLE]

and

[TABLE]

as shown in Fig. 3. $\langle\cdots\rangle_{+}$ represents the ensemble average over only positive values of $d_{0}$. Otherwise, $\epsilon^{\prime}$ does not diverge by repeating iterations. We remark that the exponent value is larger than $1/4$. The numerical exponent values in Eqs. (33) and (34) are obtained with the data only within the range $N<10^{8}$. The data beyond that range get out of the trend abruptly, which may be caused by too long passing time through too narrow bottle necks as the system size becomes large. The noise term $d_{0}$ was obtained at $n^{*}\approx 80$ in Fig. 2, at which a critical branching process starts.

Now we consider the distribution of $d_{0}$ obtained from different realizations but at the same $n^{*}$ for the system size $N$, denoted as $Q_{N}(d_{0})$. We define the standard deviation $\sigma_{N}$ of $Q_{N}(d_{0})$ as

[TABLE]

where $\langle\cdots\rangle$ represents the average over all range of $d_{0}$ and $\langle d_{0}\rangle>0$. $\sigma_{N}$ behaves as $\sim N^{-1/2}$ as shown in Fig. 4(d). If we assume that any moment of the distribution $Q_{N}(d_{0})$ is determined by the single scale, so that

[TABLE]

then it would behave as $N^{1/4}$. However, this result is not consistent with the numerical result (34).

We check the $N$-dependent behavior of $\sqrt{\langle d_{0}^{2}\rangle_{+}}$. Fig. 4(e) shows that $\sqrt{\langle d_{0}^{2}\rangle_{+}}$ behaves as $N^{-1/2}$ asymptotically but the data points deviate in small $N$ region. This discrepancy mainly originates from the asymmetry of $Q_{N}(d_{0})$, which is caused by the multiplicative noise induced by the stochastic process. $Q_{N}(d_{0})$ has a longer tail in its positive side than in the oppposite side as shown in Fig. 4(a)$-$(c). As the system size becomes larger, it becomes not only narrower but also more symmetric. To quantify this asymmetric feature of the distribution, we measure the skewness of $Q_{N}(d_{0})$ defined as

[TABLE]

in Fig. 4(f). The above result suggests that the distribution remains asymmetric in any finite systems but becomes symmetry only in the limit $N\to\infty$. $Q_{N}(d_{0})$ becomes a Gaussian distribution in that limit. The asymetry of $Q_{N}(d_{0})$ decreases because the ratio of the noise to the mean number of infected nodes becomes smaller for larger systems.

Due to those features, \big{\langle}d_{0}^{-1/2}\big{\rangle}_{+} behaves differently from $\sigma_{N}^{-1/2}$ within our numerical range; however, it is not certain yet how it would be in the thermodynamic limit because our simulation data (Fig. 3) of \big{\langle}d_{0}^{-1/2}\big{\rangle}_{+} contain heavy fluctuations, particularly in the large-system-size region. For much larger system sizes, $Q_{N}(d_{0})$ are so close to the Gaussian distribution that one may think that \big{\langle}d_{0}^{-1/2}\big{\rangle}_{+} behaves as $\sigma_{N}^{-1/2}$, i.e., $\sim N^{1/4}$ in the thermodynamic limit $N\to\infty$. However, it is a challenging task to verify that numerically.

When $\kappa>\kappa_{c}$, $d_{0}$ is naturally obtained as d_{0}=(1-i_{0})({{\partial f(u_{n},\kappa)}/{\partial\kappa}})\big{|}_{m_{d},\kappa_{c}}(\kappa-\kappa_{c}). Then, we do not need to take average over ensembles for sufficiently large $\kappa-\kappa_{c}$ because sample to sample fluctuations of $d_{0}$ become negligible compared to it. Then,

[TABLE]

Numerical result in Fig. 5 supports this prediction.

IV $K$-core percolation and the threshold model

In our previous work universal , we showed that there exists universal mechanism of avalanche dynamics in the SWIR model, $k$-core percolation, the threshold model and the CFoIN, when an avalanche starts from a single seed. Due to the universal mechanism, the golden time scales as $N^{1/3}$ in a universal way for those models. During that study, we found that numerical simulations for the CFoIN with large system sizes require long computational times and memory space, so that numerical results with limited ensemble average were not neat. Based on such experience, here we limit our interest on to the golden time problem with multiple seeds to $k$-core percolation kcore1 and the threshold model watts besides the SWIR to check the universal behavior. We find that for both models, the exponent values of the golden time are also measured to be slightly larger than 1/4. We note that the SWIR model and the two models above can be regarded as special cases of generalized epidemic process grassberger_pre ; dodds with heterogeneous transmission probabilities. Thus, similar behaviors of golden time are expected. Let us begin with $k$-core percolation.

IV.1 $k$-core percolation

Here we first consider the avalanche dynamics of $k$-core percolation. First we construct a $k$-core subgraph from an ER random graph with mean degree $z$. When $z$ is larger than a threshold $z_{c}$, a $k$-core subgraph of size $O(N)$ can exist. After this step, $\rho_{0}N$ nodes are removed simultaneously. There may exist some nodes that have degree less than $k$. In this case, those nodes are removed repeatedly until no more such nodes remain. The avalanche size can be either finite or infinite depending on $z$ and $\rho_{0}$. If it is finite, the $k$-core would still exist; If it is infinite, the $k$-core would collapse to zero. For sufficiently large $z$, there exists a critical density $\rho_{c}$ such that an infinite avalanche can occur when $\rho_{0}>\rho_{c}$ in the thermodynamic limit. In Fig. 6, we measure the mean cascade time step (golden time) $\langle n_{c}\rangle$ of infinite avalanches at the transition point for different system sizes $N$. We also measure \big{\langle}{d_{0}(n^{*})}^{-1/2}\big{\rangle}_{+} and \big{\langle}{d_{0}(0)}^{-1/2}\big{\rangle}_{+}, where $d_{0}(n^{*})$ is the noise measured at $n^{*}=60$, at which a critical branching process occurs. The definition of $d_{0}$ is presented in Appendix A. $d_{0}(0)$ is measured at $n=0$, which represents structural fluctuation at $n=0$. It was found that $\langle n_{c}\rangle$ scales as $N^{0.2655\pm 0.003}$ and \big{\langle}{d_{0}(n^{*})}^{-1/2}\big{\rangle}_{+} scales as $N^{0.262\pm 0.007}$. On the other hand, \big{\langle}{d_{0}(0)}^{-1/2}\big{\rangle}_{+} is proportional to $N^{0.25}$.

The distribution functions of $\sqrt{N}d_{0}(0)$ and $\sqrt{N}d_{0}(n^{*})$ are shown in Fig. 7 for two different system sizes. Similar to the case of the SWIR model, $d_{0}(n^{*})$ is distributed asymmetrically and the distribution becomes more symmetric for larger system sizes (Fig. 7(a)). Such factors make the exponent larger than 1/4. On the other hand, since the multiplicative fluctuations of cascade dynamics are absent at $n=0$, the distribution of $d_{0}(0)$ do not change in shape for different system sizes (Fig. 7(b)). Thus it satisfies $Q_{N}(d_{0}(0))=\sqrt{N}Q(\sqrt{N}d_{0}(0))$, which makes \big{\langle}{d_{0}(0)}^{-1/2}\big{\rangle}_{+} scale as $\sim N^{1/4}$. However, it is also uncertain that the value $\zeta$ remains unchanged for larger systems.

IV.2 The threshold model

Next we consider the threshold model, which was introduced to study the spread of cultural fads on social networks. Each node $i$ is assigned its threshold value $\phi_{i}$ and has one of two states, active or inactive. An inactive node $i$ surrounded by $m_{i}$ active neighbors and $k_{i}-m_{i}$ inactive neighbors changes its state to active when the fraction of active neighbors $m_{i}/k_{i}>\phi_{i}$. For a given set of $\{\phi_{i}\}$, the order parameter, the density of active nodes in an absorbing state, jumps and exhibits a hybrid phase transition at a critical value of the mean degree $z_{c}$. Here, we initially introduce $i_{0}N$ active nodes in a system. At each generation, every inactive node $i$ whose number of active neighbors $m_{i}>k_{i}\phi_{i}$ is identified and changes its state to active. For convenience, we choose a single threshold value $\phi$ for all nodes on ER networks. Then the critical mean degree $z_{c}$ is determined as a function of $\phi$ and $i_{0}$. We performed simulations with $\phi=0.18$ and $i_{0}=0.01$. Then the critical point is determined as $z_{c}=9.191\dots$ in the thermodynamic limit. The mean cascade time step of infinite outbreaks, $\langle n_{c}\rangle$, is obtained numerically as $\sim N^{0.263}$ in Fig. 8. Thus the measured exponent is also larger than $1/4$.

V puzzle in CFoIN

A similar size dependency of the golden time was addressed in the CFoIN. Zhou et al. zhou revealed that the choices of different types of transition points lead to different scaling behaviors of golden time in the CFoIN. They showed that when the golden time is measured at the transition point $p_{c}$ of each realization, the mean golden time scales as $\langle n_{c}\rangle\sim N^{1/3}$. On the other hand, when a single mean-field transition point $p_{c}^{\textrm{MF}}$ is taken for all realizations, the golden time scales as $\langle n_{c}^{\textrm{MF}}\rangle\sim N^{1/4}$. The authors presented the hand-waving argument that in finite systems of size $N$, individual $p_{c}$ follows a standard Gaussian distribution having the mean value $p_{c}^{\textrm{MF}}$ and the standard deviation proportional to $N^{-1/2}$ buldyrev ; zhou . Using a formula similar to Eq. (31) with $d_{0}$ following a Gaussian distribution having the standard deviation $\sim N^{-1/2}$, they obtained $\langle n_{c}^{\textrm{MF}}\rangle\sim N^{1/4}$. On the other hand, the author of Ref. grassberger_mcc investigated the scaling relation of the golden time numerically using a different algorithm, and obtained the exponent $\zeta\approx 0.28$ different from $1/4$. Thus, the two results are not consistent with each other and this discrepancy has remained as a puzzle in the cascade-induced discontinuous percolation.

We recall that for the SWIR model, $i_{0}N$ seeds are selected at random. Thus, the dynamics started from those nodes can be different for each sample. Because these choices are random, the distribution of $d_{0}$ at $n=0$ will follow a Gaussian distribution in a similar way to the $k$-core percolation case. However, because the dynamics proceeds from $n=0$ stochastically, noises are accumulated during the avalanche dynamic process. In this case, for a given network at $n=0$, noises of $d_{0}$ obtained at $n^{*}$ do not form a regular Gaussian distribution but do an asymmetric distribution, and the observed scaling of $\langle n_{c}\rangle$ is not $N^{1/4}$. We think that the result obtained in Ref. grassberger_mcc shares the common origin with the one we have in the SWIR model. Therefore, we think that the puzzle arising between the results of Refs. zhou and grassberger_mcc originates from the times at which the distribution of the fluctuation is measured.

VI Summary

The SWIR model is a simple two-step contagion model, enabling us to understand the machanism underlying a pandemic outbreak. Using this model, we obtained the scaling behavior of the golden time with respect to the system size. Using the local tree approximation, we set up a nonlinear dynamic equation in the form of saddle-node bifurcation that represents the cascade dynamics of two-step contagion. When the epidemic dynamics starts from a single infected node, we showed that the linear and the nonlinear terms of the nonlinear mapping play distinct their roles. In the early time regime, the linear term governs a critical branching (CB) process. The CB tree can be regarded as a critical cluster in percolation. However, in the late time regime, the nonlinear term causes an explosive spread of epidemic disease. The golden time is determined by the finite-size effect on the linear term, which scales as $\sim N^{\zeta}$ with $\zeta=1/3$. This scaling behavior is universal for cascade-induced dynamic models such as the threshold model, $k$-core percolation and the CFoIN.

When the dynamics starts from multiple seeds of $O(N)$, we measured a change in the value of $\zeta$ to 0.252 in the SWIR model. In this case, a long CB process does not appear from $n=0$, but it does at some characteristic time $n^{*}$. During the time until $n^{*}$, clusters of ever infected nodes merge and form a cluster of size $O(N)$. The size fluctuates for different realizations. The fluctuations are induced by the stochastic process of disease transmission. We found that these fluctuations change the value of $\zeta$ from 1/3 to about 0.252 for the multiple-seed case. Due to the multiplicative noise of disease spread, the size distribution of those clusters over different realizations becomes asymmetric with a long tail in its positive region. It seems that due to such non-Gaussianity the golden time scales as $\sim N^{\zeta}$ with $\zeta$ slightly larger than $1/4$. However, this asymmetry decreases gradually in a power-law manner of the skewness function as the system size is increased. This leaves the possibility open that asymptotically the value of $\zeta$ approaches $1/4$. This problem could not be ultimately solved by our study. A very precize analysis of corrections to scaling would be needed for it, what was not possible in spite of our massive numerical efforts.

We also obtained the similar behavior, $\zeta>1/4$ for the two other cascade dynamics models, $k$-core percolation and the threshold model. On the basis of the numerical results of the SWIR model, the threshold model and $k$-core percolation, the exponent $\zeta$ seems to be non-universal for the multiple-seed case. However, as the $\zeta$ exponents for those models deviate only slightly from 1/4 and the simulation sizes are limited, the asymptotic universal behavior cannot be entirely excluded.

Acknowledgements.

This work was supported by the National Research Foundation of Korea by grant no. NRF-2014R1A3A2069005 and by H2020 FETPROACT-GSS CIMPLEX Grant No. 641191 (JK).

Appendix A Derivation of $d_{0}$ in $k$-core percolation

We consider $k$-core subgraph of a given ER network of size $N$ and mean degree $z$. If $z$ is larger than a critical value $z_{c}$, a $k$-core subgraph of size $M(z)N$ can exist. $M(z)$ was obtained analytically in kcore_goltsev . We define $P_{d}(q)$ as the probability that a node in the $k$-core subgraph has degree $q$. We consider the evolution of the avalanche process after removing $\rho_{0}N$ nodes from the $k$-core subgraph. We define $u_{n}$ as the probability that a node attached to the end of an randomly chosen edge of a network will have degree less than $k$ at time step $n$. Then the evolution of $u_{n}$ satisfies the following equation,

[TABLE]

where $i_{0}=\rho_{0}/M(z)$ and

[TABLE]

Taking similar steps as for the SWIR model, we now consider the avalanche process after a certain time step $m$, which is

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Here $Q_{m}$ denotes the probability that a node attached to the end of randomly chosen edge in the remaining graph at time step $m$ has degree $\ell$.

We now consider sample to sample fluctuations at the characteristic time $n^{*}$. At a transition point $\rho_{0}=\rho_{c}$, Eq. (39) has a nontrivial fixed point $u_{d}<1$. Defining $h_{d,n^{*}}\equiv(u_{d}-u_{n^{*}-1})/(1-u_{n^{*}-1})$ and $\epsilon_{n}=h_{n,n^{*}}-h_{d,n^{*}}$, Eq. (39) with fluctuations becomes

[TABLE]

which can take the form of Eq. (25) with

[TABLE]