# Sandpiles on Multiplex Networks

###### Abstract

We introduce the sandpile model on multiplex networks with more than one type of edge and investigate its scaling and dynamical behaviors. We find that the introduction of multiplexity does not alter the scaling behavior of avalanche dynamics; the system is critical with an asymptotic power-law avalanche size distribution with an exponent on duplex random networks. The detailed cascade dynamics, however, is affected by the multiplex coupling. For example, higher-degree nodes such as hubs in scale-free networks fail more often in the multiplex dynamics than in the simplex network counterpart in which different types of edges are simply aggregated. Our results suggest that multiplex modeling would be necessary in order to gain a better understanding of cascading failure phenomena of real-world multiplex complex systems, such as the global economic crisis.

###### pacs:

89.75.Hc## I Introduction

Complex network theory has provided a useful tool for studying the structure and the dynamics of complex systems newman (); havlin (). Various dynamical processes occuring on networks have also been extensively studied dynamics (). An important class of dynamical processes is the so-called cascading failure dynamics, addressing how the failure of nodes induces other nodes’ failure over the network. Various real-world phenomena such as blackouts of electrical power grids motter02 (); sachtjen (), Internet congestion cascades motter02 (); ejlee (), or global economic crises kyumin () can be modeled under the framework of cascading failure on networks. Prototypical models of cascading failure include the sandpile model btw (), the threshold cascade model of behavioral adoptions watts02 (), and the capacity-based overload cascade model motter02 ().

Most studies thus far have focused on networks with a single type of link, the so-called simplex networks. In many real-world complex systems, nodes participate in more than one type of interaction with other nodes in the system; People in society have friendship, family ties, professional collaborative relations, etc. Countries in the global economy interact through various financial and political channels; they are multiplex networks. Different interaction channels do not always operate in the same way, nor do they operate in a completely autonomous way. Therefore, the multiplexity of networks introduces another nontrivial structural factor in complex networks, the dynamic consequence of which is to be understood. In this perspective, there have been recent studies related to the multiplexity of networks in the form of interdependency between network layers buldyrev10 (), interacting networks leicht (), and multi-relational structure of online social networks szell (). These studies showed that the multiplexity presents nontrivial structural patterns and connectivity leicht (); szell (); kyumin11 () and can drastically change the robustness properties of interdependent network systems buldyrev10 (). This suggests that the simplex network framework would generically be incomplete and thus one needs to consider the multiplex network framework for a better understanding of complex systems.

Each type of interaction in multiplex networks defines a network layer; for example, a social network may comprise a friendship layer, a family tie layer, etc. A change of a node’s state may be determined by intralayer dynamics; for example a person in a social network may adopt a behavior (using a smartphone app, say) by the peer pressure from his/her friends or by a demand due to the work environment alone. However, the potential influence of his/her state change is not confined within that layer; even though you adopted the app due to your work demand, your usage of the app may spread to your friends or family as you may recommend it through your friendship and family tie layers. In this way, the layers of a network get coupled and become interdependent. Such a multiplex coupling is generically not equal to a simple aggregation of link types, thereby leading to nontrivial dynamical consequences.

In this paper, we apply multiplex modeling to a simple model of cascading failure, the sandiple model, and investigate its dynamical impact. We begin by reviewing the sandpile dynamics on random, simplex networks in Sec. II. The multiplex sandpile model is introduced in Sec. III, and the principal results will be presented in Sec. IV. Other model variants are considered in Sec. V, and the paper concludes in Sec. VI.

## Ii Sandpiles on random networks

The sandpile model was originally introduced by Bak, Tang, and Wiesenfeld (BTW) in 1987 btw (). Since then, this model has been an archetypical theoretical model for investigating cascading failure dynamics and self-organized criticality soc (). The dynamics proceeds as follows. i) Each node is assigned a prescribed threshold height, usually taken to be its coordination number . ii) At each step, a grain is added at a randomly chosen node, say , raising the node’s height by one, . iii) If the height of the node becomes equal to or exceeds its threshold, that is , the node becomes unstable and topples grains to its neighbors, such that and for being the neighbor of . iv) If this toppling event induces instability in some neighbors, they also topple in the same way; a cascade of toppling events (an avalanche) occurs until there remain no more unstable nodes in the system. v) Repeat ii)–iv). BTW showed that without a fine-tuning of the external parameters, this model could exhibit a scale invariance in the form of a power-law distribution of avalanche sizes (the number of nodes participating in a given avalanche),

(1) |

hence, the name self-organized criticality btw ().

In Euclidean lattices in low dimensions, the BTW sandpile model exhibits nontrivial critical behaviors soc (). Its dynamics on random networks, however, can be understood more easily, following the mean-field behaviors, alstrom (); bonabeau (). Later, the model has been considered on the scale-free networks with a power-law degree distribution, , where the degree is the number of links of a node goh03 (); jkps04 (); physa04 (); physa05 (). For the BTW sandpile model with on scale-free networks, the unconventional mean-field behavior was obtained such that the avalanche size exponent is given by goh03 ()

(2) |

Here, we briefly present the analytical approach based on the branching process approximation leading to these results goh03 (); jkps04 (). Assuming that the avalanche proceeds without forming a loop, it can be mapped onto a branching process. In this mapping, the probability to make a -branch at each branching event is given by

(3) |

The first factor comes from the likelihood that a node with degree would gain a grain from its unstable neighbor, which is proportional to its degree. The second factor comes from the assumption that the timing at which the node gets a grain from its neighbor is random, such that the probability to become unstable (that is, at the moment) is . is given by the normalization as , and the mean branching number is obtained as

(4) |

Therefore, the branching process is critical. The asymptotic (large ) behavior of the avalanche size distribution can be obtained from the singular behavior near unity in the corresponding generating function jkps04 (). On scale-free networks, the branching probability decays algebraically, , which gives rise to a singular term in its generating function, leading to the asympototic power-law behavior of for as given by Eq. (2). Thus, the behavior of the branching probability provides an essential clue to the sandpile dynamics.

## Iii Sandpile model on multiplex networks

In this paper we generalize the BTW model on multiplex networks. In multiplex networks, each type of interaction functions as a network layer within which the dynamics takes place as if in the usual single network case. However, the state change of a node in one layer applies also to other layers, giving rise to couplings between layers’ dynamics. Such an interlayer coupling induces interdependency between layers and renders the multiplex dynamics rather nontrivial buldyrev10 (); bjkim (); brummitt (). To illustrate the effect of multiplexity in its simplest setting, we consider in the paper multiplex networks with two kinds of interactions or a duplex networks.

Specifically, the BTW sandpile model on duplex networks is run as follows: i) Each node is assigned the thresholds equal to its degree for each layer; that is, for and . ii) At each step, a grain is added to a randomly chosen node in a randomly chosen layer ; . iii) If the height of this node reaches or exceeds its threshold in any layer, that is, either or , the node becomes unstable, initiating an avalanche. The layer in which the threshold is exceeded is called the unstable layer. iv) Each unstable node topples on both layers; in the unstable layers, it topples grains, one for each neighbor, whereas it topples only grains in the other layer, giving one grain to each of a random subset of neighbors [Fig. 1(b)]. That is,

(5) |

and the height of all the nodes receiving a grain increases by one. v) If this toppling event causes other nodes to be unstable, these unstable nodes topple in the same way as in iv). vi) Repeat iv)–v) until there are no more unstable nodes in the system, and the avalanche ends. During each toppling, the grain is lost with a small probabilty, taken to be in this paper, which is necessary to prevent overloading the system. We call this multiplex dynamics model the 12 model.

In our 12 multiplex setting, a node’s failure (toppling) due to a specific layer dynamics induces its simultaneous failure in other layers regardless of its conditions in those layers. This setting is motivated by a number of real-world situations such as the global economy, in which one country’s recession due to an imbalance in the trade system can affect the country’s failure in a credit network, for example. The 12 model is a stylized model of such scenarios, aiming to address the effect of multiplexity. To highlight the effect of multiplexity in the 12 model, we compare it to a non-multiplex, or simplex, counterpart of the model, which we call the 12 model. In the 12 model, we ignore the individuality of different interactions and simply overlay two layers [Fig. 1(c)] to form a simplex aggregate network, on which the conventional BTW model is run. In the next section, we study and compare the avalanche dynamics of these two models.

## Iv Results

The primary quantity of interest in sandpile dynamics is the avalanche size distribution , giving the likelihood that an avalanche of size is to occur in the steady state of the model btw (); soc (). We meaure the avalanche size in terms of the number of nodes becoming unstable during an avalanche. A related quantity is the number of toppling events during an avalanche, which scales similarly when the formation of loops during an avalanche is negligible, which is the case for sandpile dynamics on random graphs goh03 (). One can also measure the distribution of avalanche durations, which also follows a power law, whose exponent is related to the avalanche size exponent by a scaling relation soc (). In this work we focus exclusively on the avalanche size distribution.

### iv.1 On Duplex Erdős-Rényi Layers

We first consider the sandpile dynamics on multiplex networks with two Erdős-Rényi (ER) er () layers of nodes and the equal mean degree . Each node has degrees and in the first and the second layers, respectively, each of which is a Poisson random variable with the mean equal to . To obtain the avalanche size distribution and related statistics, we sampled successive avalanches after the initial transient period. For the simplex 12 model, the substrate network is nothing but an ER network with mean degree . Therefore, the 12 model is trivially critical, with the standard mean-field exponent bonabeau (); alstrom (), which is also confirmed by numerical simulations (Fig. 2).

For the multiplex 12 model, the sandpile dynamics would become more complicated due to layers’ coupling or interdependency. Even the criticality of the model is not obvious a priori. The numerical simulation of the 12 model, however, showed that its avalanche size distribution is essentially identical to that of the 12 model; follows a power law, that is, the model is still critical, with the same mean-field exponent (Fig. 2). Such a robust mean-field behavior was also reported recently for sandpile models on interdependent networks in somewhat different settings charlie ().

Although the two models show similar with the same exponent , the specific sandpile dynamics underlying this distribution are different, given the presence/absence of multiplexity. Compared to the simplex case that can be understood in terms of the framework presented in Sec. II, the multiplex dynamics follows a more complicated cascade process. To provide more detailed insight, we again employ the branching process approximation following Sec. II. For the simplex 12 case, the branching probability is given by, from Eq. (3),

(6) |

where is the degree distribution of the aggregate network of two layers. Equation. (6) is confirmed by numerical simulations [Fig. 3(a)]. For the multiplex 12 model, the associated branching process can be divided into two parts: one through the unstable layer (denoted layer 1) and the other through layer 2. The latter is the consequence of multiplexity and the layers’ interdependency. To obtain the branching probability into branches, one has to consider the composite branching of the -branch in layer 1 and the -branch in layer 2. Therefore, one can write

(7) |

The first summand gives the probability that a node with degree receives a grain and becomes unstable, making a -branch. The second summation gives the probability that the node has grains in layer 2 at the moment (thereby branching), assuming again the random timing of topplings. Summing the combined factor over , we have . Multiplexity is clearly manifested in the second factor, but it does play a role in the first factor, too. The unspecified factor addresses the probability that a node becomes unstable as it receives a grain, which was for simplex sandpile, Eq. (3). The assumption behind it was that all heights are equally likely when encountering a toppling, as the grains are successively piled. In the multiplex model, however, between two successive receipts of grains, a node can topple due to an instability in the other layer. Accordingly, the likelihood that a node becomes unstable as it receives a grain becomes smaller than , and we numerically find that it scales roughly as with in our simulation setting [Fig. 3(a), inset]. Therefore, the multiplex coupling induces truly nontrivial dynamics that cannot be viewed as a simple sum of two individual layers’ dynamics.

The overall is shown in Fig. 3(a) along with the simplex counterpart, showing that in multiplex sandpile dynamics, branching into smaller number of branches (toppling a smaller number of grains) occurs more frequently than in simplex sandpiles. Meanwhile, the mean branching number is in both cases; numerically, we measured and . This supports the robustness of self-organized criticality. Although takes a nontrivial form, its large- behavior is still dominated by the Poisson tail, which provides an explanation for the robustness of the avalanche size exponent.

To differentiate the cascade dynamics of the two models in more detail, we measured the frequency of topplings in nodes with and , and plotted the normalized difference in Fig. 4. One can see that most nodes fail more frequently in the 12 model than in the 12 model whereas the nodes with links in only one layer (either or ) fail more in the simplex dynamics. This result suggests higher vulnerability of most nodes in multiplex sandpile cascades than in simplex scenarios. This picture is in line with recent studies indicating higher vulnerability of interdependent networks buldyrev10 () and multiplex cascade models charlie ().

### iv.2 On Scale-free Layers

We have also performed numerical simulations on multiplex networks with i) two scale-free (SF) layers and ii) one ER and one SF layer. For the SF layer, we used the static model goh01 () with the same mean degree and the degree exponent . For the simplex dynamics, the degree distribution of aggregate networks in both cases has a power-law tail with exponent , which leads to the predicted avalanche size exponent according to Eq. (2). The multiplex dynamics is still critical ( and ), having the same power-law avalanche size distribution in both cases [Figs. 5(a,b)]. The detailed cascade dynamics are different; the branching probability is larger for very small and very large , but smaller for intermediate in the multiplex dynamics than in simplex ones [Figs. 5(c,d)]. This result suggests an interesting picture that hubs becomes more vulnerable in multiplex sandpiles. In an ordinary (simplex) sandpile, hubs play the role of a shock-absorbers, by tolerating a large amount of stress (grains) without failing; in multiplex dynamics, hubs fail more often as there is an additional chance of failure due to its interdependency with a low degree node. Therefore, although the multiplex and the simplex dynamics do not differ in terms of scaling behaviors, the microscopic picture of cascade dynamics at the individual node level does show differences that may have implications for risk assessment and aversion strategy for cascading failure in real-world phenomena.

## V Modifications of the model

In this section, we consider other variants of the BTW sandpile model on multiplex networks. The first variant is the one with a different rule iii) of the 12 model in Sec. III. In the new model, a node fails (topples) if it is unstable in both layers, that is, for both . When a node fails, it topples grains in each layer. We call this version the AND-I model. We also consider another rule such that a node topples all its grains. We call this variant the AND-II model. Therefore, for these variants we have instead of Eq. (5),

(AND-I model), | |||||

(8) |

For the AND-II model, when a node topples more than grains in the layer , it first distributes grains to its neighbors and the remaining grains to randomly chosen neighbors. The scaling behavior is found to be still robust with respect to both variants of multiplex sandpile dynamics. They are critical ( and ) and have the same avalanche size exponent on duplex ER networks [Fig. 6(a)]. Interestingly, the branching probability for these two variants is found to be the same as that of the simplex 12 model [Fig. 6(b)].

The final variant we consider is to change rule iv) of the 12 model such that a node fails if it is unstable in any layer, but the unstable node topples grains in both layers even when its current height is less than the threshold in layer 2. That is, we have,

(9) |

for both , and the height of all the neighbors of increases by one. We symbolized this as the 12 model, to distinguish it from the original multiplex 12 model. In this version, grains are not conserved, but produced, during the avalanche. This setting is less physical, but can model situations in which the failure due to interdependency induces the full stress () rather than the current stress level. We found that in this model, the system quickly enters into an infinite avalanche, in which the number of grains increases abruptly and eventually the avalanche covers the whole system, after which a finite fraction of nodes topple ad infinitum [Figs. 7(a-c)]. The branching analysis indeed shows that the mean branching number becomes larger than one, increasing up to [Fig. 7(d)]. Therefore, the branching process becomes highly supercritical and can continue forever. The excess stress (grains) generated by interdependency regenerates further excess stress, and the system falls into the trap of an infinite sequence of cascading failure. The system is no longer critical. A similar absence of scaling had been observed for the conventional BTW model on Euclidean lattices without dissipation at boundaries grassberger ().

## Vi Summary and conclusion

We have studied a number of variants of the Bak-Tang-Wiesenfeld sandpile model on multiplex networks. To account for the multiplexity of real-world complex systems, we introduced the multiplex sandpile model in which the network layers are interdependent so that a node fails if it becomes unstable in any layer (called the 12 model). We found that this multiplex dynamics did not alter the scaling behavior of the mean-field-type sandpile dynamics, yet the detailed pattern of the microscopic cascade dynamics at the individual node level was affected by the multiplexity. Compared to the simplex counterpart, in multiplex dynamics, higher-degree nodes such as hubs in scale-free layers are more likely to fail. This finding resonates with the recent finding of higher vulnerability of scale-free interdependent networks than Poisson-random ones under mutual percolation buldyrev10 (). The robustness of scaling behaviors is further tested by considering other variants of multiplex dynamics. We could alter the scaling behavior only by breaking the conservation law, introducing excess stress due to interdependency, which drives the system supercritical. In this sense, the scaling properties of the BTW sandpile model is strongly robust on networks. In order to alter its scaling properties, one needs to alter the asymptotic form of , which is strongly constrained by the underlying substrate network structure. Therefore, random multiplex couplings between finite number of layers would not likely induce a nontrivial scaling behavior of BTW-type sandpile dynamics despite differences in the microscopic dynamics. We note, however, that for other classes of cascade dynamics, the multiplexity can play a critical role brummitt ().

###### Acknowledgements.

This work was supported in part by Mid-career Researcher Program (No. 2009-0080801) and Basic Science Research Program (No. 2011-0014191) through NRF grants funded by the MEST. K-ML is supported in part by Global Ph.D. Fellowship Program (No. 2011-0007174) through NRF, MEST.## References

- (1) M. E. J. Newman, A.-L. Barabási, and D. J. Watts, The structure and dynamics of networks (Princeton University Press, Princeton, 2006).
- (2) R. Cohen and S. Havlin, Complex networks: Structure, robustness and function (Cambridge University Press, Cambridge, 2010).
- (3) A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical processes on complex networks (Cambridge University Press, Cambridge, 2008).
- (4) M. L. Sachtjen, B. A. Carreras, and V. E. Lynch, Phys. Rev. E 61, 4877 (2000).
- (5) A. E. Motter and Y.-C. Lai, Phys. Rev. E 66, 065102 (2002).
- (6) E. J. Lee, K.-I. Goh, B. Kahng, and D. Kim, Phys. Rev. E 71, 056108 (2005).
- (7) K.-M. Lee, J.-S. Yang, G. Kim, J. Lee, K.-I. Goh, and I.-M. Kim, PLoS ONE 6, e18443 (2011).
- (8) P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987); Phys. Rev. A 38, 364 (1988).
- (9) D. J. Watts, Proc. Natl. Acad. Sci. U.S.A. 99, 5766 (2002).
- (10) S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin, Nature 464, 1025 (2010).
- (11) E. A. Leicht and R. M. D’Souza, arXiv:0907.0894.
- (12) M. Szell, R. Lambiotte and S. Thurner, Proc. Natl. Acad. Sci. U.S.A. 107, 13636 (2010).
- (13) K.-M. Lee, J. Y. Kim, W.-k. Cho, K.-I. Goh, and I.-M. Kim, arXiv:1111.0107.
- (14) H. J. Jensen, Self-organized criticality (Cambridge University Press, Cambridge, 1998).
- (15) P. Alstrøm, Phys. Rev. A 38, 4905 (1988).
- (16) E. Bonabeau, J. Phys. Soc. Jpn. 64, 327 (1995)
- (17) K.-I. Goh, D.-S. Lee, B. Kahng, and D. Kim, Phys. Rev. Lett. 91, 148701 (2003).
- (18) D.-S. Lee, K.-I. Goh, B. Kahng, and D. Kim, J. Korean Phys. Soc. 44, 633 (2004).
- (19) D.-S. Lee, K.-I. Goh, B. Kahng, and D. Kim, Physica A, 338, 84 (2004).
- (20) K.-I. Goh, D.-S. Lee, B. Kahng, and D. Kim, Physica A 346, 93 (2005).
- (21) J. Um, P. Minnhagen, and B. J. Kim, Chaos 21, 025106 (2011).
- (22) C. D. Brummitt, K.-M. Lee, and K.-I. Goh, arXiv:1112.0093.
- (23) P. Erdős and A. Rényi, Publ. Math. Inst. Hung. Acad. Sci. 5, 17 (1960).
- (24) C. D. Brummitt, R. M. D’Souza, and E. A. Leicht, arXiv:1106.4499v1.
- (25) K.-I. Goh, B. Kahng, and D. Kim, Phys. Rev. Lett. 87, 278701 (2001).
- (26) P. Grassberger and S. S. Manna, J. Phys. (Paris) 51, 1077 (1990).