Two golden times in two-step contagion models: A nonlinear map approach.

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. Using an approach based on a nonlinear mapping, 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^{ζ}) with the system size N on Erdős-Rényi networks, where the measured ζ 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 ζ exponents are all close to 1/4 but show model-dependence. It remains open whether or not ζ reduces to 1/4 in the asymptotically large-N limit. Our method can be applied to several models showing a hybrid percolation transition and gives insight into the origin of the two golden times.

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 [1]. 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 [2,3] and the susceptible–infected–susceptible (SIS) model [4]. Analytical and numerical studies of those models revealed that a continuous phase transition occurs on Erdős-Rényi (ER) random networks [5]. 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 [6,7]. 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 [8-15]. In another instance of modified SIR models, a network evolves by rewiring links at a certain rate during the spread of contagion [16]. 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 [17] 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 [18].

It was revealed [14,17] that when a disease starts spreading from a single node, the golden time scales as with 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 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 for the golden time is consistent with the characteristic size of the giant cluster generated in the SIR model [19,20], thus it is 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 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 to a value slightly larger than 1/4. We estimate this scaling behavior using the saddle-node bifurcation theory [21] and discuss the underlying mechanism.

Similar size dependencies of mean cascading time at a transition point were studied for other cascade dynamics models such as -core percolation [22-27] and cascading failure model on interdependent network (CFoIN) [28-32]. It was found [28,29] that in the CFoIN, the mean cascading time is proportional to or depending on the way of choosing the transition points. References [26,27] showed that the exponent 1/3 is also obtained in -core percolation. Thus, the scaling behavior of is robust. However, for the case, a different scaling behavior with  [33] 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 -core percolation and the threshold model [34,35]. We show that when the cascade starts from a fixed number of multiple seeds , the golden times for both models also become proportional to , where 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 in large- 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. III A. We show that the roles of the linear and nonlinear terms are well separated. In Sec. III B, 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 -core percolation and the threshold model and show that the numerical values of are slightly larger than . We also discuss the possibility of 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.

THE SWIR MODEL

The SWIR model is a generalization of the SIR model by including two states, a weakened state (denoted as ) and an infected state (), between the susceptible state () and removed state (), instead of a single infected state alone, as in the SIR model. Nodes in state are involved in the reactions and , which occur in addition to the reactions and in the SIR model. At each discrete time step , the following processes are performed. (i) All the nodes in state 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 , it changes its state in one of two ways: either to with probability or to with probability . If a neighbor is in the state , it changes to with probability , where , and are the contagion probabilities for the respective reactions. (iii) All nodes in the list change their states to . 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:

In an absorbing state, each node is in one of three states, the susceptible, weakened, or removed state. We define as the conditional probability that a node remains in state in the absorbing state, provided that it has neighbors in state and was originally in state S. This means that the node remains in state even though it has been in contact times with these neighbors in state before they change their states to . Thus, we obtain Next, is similarly defined as the conditional probability that a randomly selected susceptible node is in state after it contacts neighbors in state before they change their states to . The probability is given as Finally, is the conditional probability that a node has been infected in any state, either or , provided that it was originally in state and its neighbors are in state in the absorbing state. Using the relation , one can determine in terms of and .

On a network with a degree distribution , we consider the case in which the initial densities of susceptible, weakened, and infectious nodes are given as and , provided that . The order parameter , the density of nodes in state after the system falls into an absorbing state, is given using the local tree approximation as where and is the probability that an arbitrarily chosen edge leads to a node in state or in the absorbing state; however, the infection of that node was not made through the chosen edge. We define similarly to but at time step . The probability can be derived from as follows: where is the mean degree of the network and the factor is the probability that a node connected to a randomly chosen edge has degree . As converges to .

GOLDEN TIMES IN THE SWIR MODEL

The single-seed case

First, we consider the case in which the initial number of infectious nodes is ; that is, and in the thermodynamic limit. In this case, the SWIR model exhibits a mixed-order transition [14] at a transition point when the mean degree is larger than a critical value. The order parameter displays a discontinuous transition from to , whereas other physical quantities such as the outbreak size distribution exhibit a critical behavior. The behavior of the order parameter as a function of is schematically shown in Fig. 1(a).

FIG. 1.

Schematic plots of the order parameter versus for (a) (the single-seed case) and (b) (the multiple-seed case). For appropriate values of parameters , and , hybird percolation transitions occur. (a) Even though is increased, the order parameter remains zero up to a transition point . However, at , it remains at or jumps to . For (b), as is increased, gradually increases from to a finite value at . At remains at or jumps to .

We are interested in how infected nodes spread as a function of the cascade step 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., , where is the mean degree, Eq. (10) is reduced as follows: We remark that on ER networks, in the limit becomes equivalent to obtained from Eq. (7).

We pick up the contribution of the reaction from Eq. (11) but neglect the contribution of the reaction . Then, the probability that a node becomes directly infected by infectious neighbors, which is denoted by , is given as Applying the formula for the Poisson degree distribution to Eq. (10), we obtain that Because the order parameter increases from , we assume that is small in the early time regime. Thus, where . Actually, the coefficient of the first-order term is the mean branching ratio in the early time regime. The critical branching (CB) process occurs at the transition point . However, the discrete mapping Eq. (14) at may be rewritten in the form of a saddle-node bifurcation, , where is a function of the continuous time variable and the overdot denotes differentiation with respect to it. Because is a stable fixed point for , and this point represents the fixed point of the SIR model, indicating a second-order transition.

Next, we consider the two successive reactions and , in which a susceptible node becomes infected in two steps and is eventually removed. Because a node can be infected either by the reaction or by the reactions and , the probability can be obtained using the relation as Again, using , we obtain that where . Here we note that the first-order term is absent. Combining Eqs. (14) and (17) at criticality, we obtain that Thus, . When , i.e., , the fixed point is stable, and thus a continuous transition occurs. Otherwise, the fixed point is unstable, and a discontinuous transition occurs. The condition for a discontinuous transition is consistent with previously obtained results [8,14].

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 [17], i.e., the mean branching ratio is almost unity, and the main contribution is that of the reaction . Thus, in the thermodynamic limit, always stays zero so that nonlinear terms in Eq. (18) do not appear. However, in finite systems, grows gradually and the nonlinear term becomes significant after a characteristic time . It was argued in Ref. [19] that for the SIR model at the epidemic threshold, the maximum size of outbreaks is proportional to in the mean field limit. When grows up to , 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 , because the fractal dimension of the CB tree is two. However, for the SWIR model, the coefficient of the nonlinear term Eq. (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 to that of the SIR model, which is what we observed in a previous work [14,17].

The multiple-seed case

Next, when the number of infectious nodes is , i.e., , and in the thermodynamic limit, it was shown [7,10-12,15] that there exists a critical value such that when , a hybrid phase transition occurs at a transition point , whereas when , a continuous transition occurs. Here we focus on the former case.

In the multiple-seed case, Eq. (10) becomes when the network is an ER network with mean degree . Fixed points of Eq. (19) satisfy the equation, and the smallest solution among them is the order parameter . We note that contains the parameters , and . As already shown in the single-seed case, for appropriately given values of and and . Thus when is sufficiently small, i.e., satisfies for values of near zero. Then increases continuously as is increased untill reaches a critical value such that and are satisfied. We note that depends on and . When is reduced to , the transition point of the single-seed case. exhibits a critical behavior as approaches and subsequently jumps from to another value 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 [15]. Accordingly, the reaction can occur even in the early time regime, as shown in Fig. 2(c) with red zigzag (lowest) curve. Moreover, the CB process appears not from the beginning but slightly after that indicated by an arrow at in Fig 2(c), at which the density of removed nodes is close to indicated in Fig. 1(b). From this step remains almost constant for a long time as shown in Fig. 2(d).

FIG. 2.

(a) For the single-seed case, plot of the branching ratios as a function of time step for each type of reaction at . (b) Evolution of the densities of removed nodes (red curve, top) and weakened nodes (blue curve, bottom) as a function of for the single-seed case. Data are obtained from a single realization of infinite outbreak in the SWIR model starting from a single seed () with reaction probabilities and for both (a) and (b). (c) Similar to (a) but for the multiple-seed case at . (d) Similar to (b) but for the multiple-seed case. For both (c) and (d), data are obtained using the parameters , and . The ER networks on which the simulations were performed have a size and mean degree . Legends “transferring to ” in (a) and (c) indicate the mean number of neighbors of an infected node that change their state from to at step . In (d), the characteristic time steps and , from and at which the CB process starts and ends, respectively, are marked.

However, in finite systems, due to the fluctuations arising in the stochastic process of epidemic spread, the densities of each species of nodes at can be different for each realizaton. Those fluctuations lead to different for each realization, where is the golden time at which begins to increase drastically as shown in Fig 2.

We denote the densities of each species of nodes at time as and , respectively. Then, on ER networks, for satisfies where and Moreover, using Eq. (10), the relation , and , we determine and .

We focus on the density fluctuations of each species at . We split the densities of each species of nodes into two parts: ( or ), where the first term represents densities of -species nodes obtained from Eqs. (19), (22), and (23) at . The second one is the deviation. Then for , Eq. (21) becomes We did not take into account the density fluctuations induced after , because they are negligible compared to those at . At , Eq. (24) has a nontrivial fixed point in the thermodynamic limit. Then Eq. (24) is rewritten with as where Here, . is similarly defined. Neglecting higher order terms of , Eq. (25) is rewritten in an alternative form, Because for large , the last term can be neglected compared to and Eq.(30) is rewritten simply as where and . We note that is a real number, while is a positive number.

The nonlinear mapping Eq. (31) includes several features: When reaches a fixed point ; when remains at zero; when , there arises the so-called bottleneck effect at  [21,25,28]. The time to pass through the bottleneck is calculated as which is approximately the time interval of the plateau region, i.e., . Because is much smaller than , which is the golden time for a single realization of the process. can have different values for different realizations, yielding a different . Thus, we need to take average of over different realizations to obtain .

We performed extensive numerical simulations at the transition point of the SWIR model starting from multiple seeds , and obtained that and as shown in Fig. 3. represents the ensemble average over only positive values of . For negative values of , Eq. (31) has a fixed point and does not diverge by repeating iterations. We remark that the exponent value is larger than . The numerical exponent values in Eqs. (33) and (34) are obtained with the data only within the range . The data beyond that range deviate from the tendency, which may be caused by too long passing time through too narrow bottle necks as the system size becomes large. The noise term was obtained at , at which a critical branching process starts as shown in Fig. 2(c).

FIG. 3.

For the multiple-seed case, (a) plot of the average cascade time step () and () versus system size at . The notation in indicates that only positive values of are considered in taking the average. is measured at . The guideline has a slope of 0.252. (b) Plot of (), (), and () versus . Data were obtained from ER networks of different sizes but with the same mean degree, . , and were used. Average is taken over more than realizations for each data point for .

Now we consider the distribution of obtained from different realizations but at the same for the system size , denoted as . We define the standard deviation of as where represents the average over all range of and . behaves as as shown in Fig. 4(d). If we assume that any moment of the distribution is determined by a single scale, so that then would scale as . However, this result is not consistent with the numerical result Eq. (34).

FIG. 4.

Plot of the probability distribution versus for different system sizes, (a), (b), and (c). The distribution is obtained from more than realizations for each . Plot of the standard deviation (d), (e), and the skewness (f) of as a function of the system size . Average is taken over more than realizations for each data point.

We check the -dependent behavior of . Figure 4(e) shows that behaves as asymptotically but the data points deviate in small region. This discrepancy mainly originates from the asymmetry of , which is caused by the multiplicative noise induced by the stochastic process. has a longer tail in its positive side than in the oppposite side as shown in Figs. 4(a)–4(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 defined as in Fig. 4(f). The above result suggests that the distribution remains asymmetric in any finite systems but becomes symmetric only in the limit . becomes a Gaussian distribution in that limit. The asymmetry of decreases because the ratio of the noise to the mean number of infected nodes becomes smaller for larger systems.

Due to those features, behaves differently from 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 contain heavy fluctuations, particularly in the large-system-size region. For much larger system sizes, are so close to the Gaussian distribution that one may think that behaves as , i.e., in the thermodynamic limit . However, it is a challenging task to verify that numerically.

When is naturally obtained as . Then, we do not need to take average over ensembles for sufficiently large because sample to sample fluctuations of become negligible compared to it. Then, Numerical result in Fig. 5 supports this prediction.

FIG. 5.

For the multiple-seed case, scaling plot of the average cascade time step versus for different system sizes . Data for different system sizes collapse well onto a single curve, indicating that for . Numerical simulations were performed on ER networks with mean degree and initial density of seeds . Average is taken over more than realizations for each data point.

-CORE PERCOLATION AND THE THRESHOLD MODEL

In our previous work [17], we showed that there exists universal mechanism of avalanche dynamics in the SWIR model, -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 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 -core percolation [22] and the threshold model [34] 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 . We note that the SWIR model and the two models above can be regarded as special cases of generalized epidemic process [8,35] with heterogeneous transmission probabilities. Thus, similar behaviors of golden time are expected. Let us begin with -core percolation.

-core percolation

In -core percolation, the dynamics starts with the removal of all nodes that have degree less than . These removals may decrease the degrees of the remaining nodes to less than . If such nodes exist, they are removed as well. This pruning process is repeated until no more such nodes remain. When the mean degree of the original network is larger than a threshold , a -core subgraph of size can exist.

Here we consider the avalanche dynamics of -core percolation [27] induced by failures of multiple nodes. We first construct a -core subgraph from an ER random graph with mean degree . After that, we perform an avalanche process by removing randomly selected nodes from the -core subgraph simultaneously. This removal may cause additional successive deletion of nodes. The avalanche size is defined as the total number of removed nodes in the avalanche process, and the duration time of the avalanche is the number of pruning steps. The avalanche size is classified as finite or infinite according to the criterion: If it is finite, the -core would still remains; if it is infinite, the -core would collapse to zero. The avalanche size depends on and . For sufficiently large , there exists a critical density such that an infinite avalanche can occur when in the thermodynamic limit.

In Fig. 6, we measure the mean cascade time step (golden time) of infinite avalanches at the transition point for different system sizes . We also measure , where is the noise measured at , at which a critical branching process starts. The definition of is presented in the Appendix. It was found that scales as and scales as . The distribution function of is shown in Fig. 7(a) for two different system sizes. Similar to the case of the SWIR model, is distributed asymmetrically and the distribution becomes more symmetric for larger system sizes. The skewness of is plotted as a function of the system size in Fig. 7(b). It shows a power-law decaying behavior. Those obtained facts lead to the measured exponent being larger than , supporting the numerical result.

FIG. 6.

For -core percolation with , (a) plot of the golden time () and (), and () versus for -core percolation with starting from multiples nodes of . is measured at and denotes the value of measured at . A solid (dashed) guideline has slope of 0.2655 (0.25). (b) Plot of (), (), (), and () versus . Data were obtained from ER networks of different sizes but with the same mean degree, , and . Average is taken over more than realizations for each data point.

FIG. 7.

(a) Plot of the probability density function of for two different system sizes and . (b) Plot of skewness of as a function of the system size . (c) Plot of the probability density function of for two different system sizes and . Average is taken over more than realizations for each data point.

For comparison, is proportional to where is the noise measured at , which represents structural fluctuation at [Fig. 6(a)]. Since the multiplicative fluctuations of cascade dynamics are absent at , the distribution of do not change in shape for different system sizes [Fig. 7(c)]. Thus, it satisfies , which makes scale as .

The threshold model

Next we consider the threshold model [34], which was introduced to study the spread of cultural fads on social networks. Each node is assigned its threshold value and has one of two states, active or inactive. An inactive node surrounded by active neighbors and inactive neighbors changes its state to active when the fraction of active neighbors . For a given set of , 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. Here, we initially introduce active nodes in a system. At each generation, every inactive node whose number of active neighbors is identified and changes its state to active. For convenience, we choose a single threshold value for all nodes on ER networks. Then the critical mean degree is determined as a function of and . We performed simulations with and . Then the critical point is determined as in the thermodynamic limit. The mean cascade time step of infinite outbreaks, , is obtained numerically as in Fig. 8. Thus, the measured exponent is also larger than .

FIG. 8.

(a) Plot of the average cascade time step for the threshold model starting from initial multiple active nodes. Guideline has a slope of 0.263. (b) Plot of () and () versus . Data were obtained from ER networks of different sizes but with the same . Average is taken over more than realizations for each data point.

PUZZLE IN CFOIN

A similar size dependency of the golden time was addressed in the CFoIN. Zhou et al. [29] 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 of each realization, the mean golden time scales as . However, when a single mean-field transition point is taken for all realizations, the golden time scales as . The authors presented the hand-waving argument that in finite systems of size , individual follows a standard Gaussian distribution having the mean value and the standard deviation proportional to  [28,29]. Using a formula similar to Eq. (31) with following a Gaussian distribution having the standard deviation , they obtained . However, the author of Ref. [33] investigated the scaling relation of the golden time numerically using a different algorithm, and obtained the exponent different from . 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, 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 at will follow a Gaussian distribution in a similar way to the -core percolation case. However, because the dynamics proceeds from stochastically, noises are accumulated during the avalanche dynamic process. In this case, for a given network at , noises of obtained at do not form a regular Gaussian distribution but do an asymmetric distribution, and the observed scaling of is not . We think that the result obtained in Ref. [33] 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. [29] and [33] originates from the times at which the distribution of the fluctuation is measured.

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 with . This scaling behavior is universal for cascade-induced dynamic models such as the threshold model, -core percolation and the CFoIN.

When the dynamics starts from multiple seeds of , we measured a change in the value of to 0.252 in the SWIR model. In this case, a long CB process does not appear from , but it does at some characteristic time . During the time until , clusters of ever infected nodes merge and form a cluster of size . 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 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 with slightly larger than . 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 approaches . This problem could not be ultimately solved by our study. A very precise analysis of corrections to scaling would be needed for it, which was not possible in spite of our massive numerical efforts.

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