Classes of critical avalanche dynamics in complex networks

Dynamical processes exhibiting absorbing states are essential in the modeling of a large variety of situations from material science to epidemiology and social sciences. Such processes exhibit the possibility of avalanching behavior upon slow driving. Here, we study the distribution of sizes and durations of avalanches for well-known dynamical processes on complex networks. We find that all analyzed models display a similar critical behavior, characterized by the presence of two distinct regimes. At small scales, sizes and durations of avalanches exhibit distributions that are dependent on the network topology and the model dynamics. At asymptotically large scales instead --irrespective of the type of dynamics and of the topology of the underlying network-- sizes and durations of avalanches are characterized by power-law distributions with the exponents of the mean-field critical branching process.


I. INTRODUCTION
In this paper, we study seven stochastic models that are prototypical to describe the diffusion of some sort of "activity" in networks [2,3]. Specifically, our analysis includes: the competition-induced-criticality model (CIC) [4][5][6], the voter model (VOT) [7,8], the invasion process (IP) [9], link dynamics (LD) [10], the contact process (CP) [11], the susceptibleinfected-susceptible model (SIS) and the susceptible-infectedrecovered model (SIR) [3] (see Figure 1). In all these models, active nodes can pass the active status to their inactive neighbors in the network, and can return to the inactive status either spontaneously or by interacting with inactive neighbors. The rules that govern these transitions are model specific. All of the models are characterized by the existence of one or more absorbing states, at which all dynamics ceases and the system state remains frozen. Typically, when parameters are set such that the system is at the interface between the active and inactive/absorbing phases, critical behavior -characterized by power-law distributions of the sizes and durations of activity avalanches-emerges.
Many studies of spreading phenomena in social, technological and biological networks [12][13][14] reveal such a critical behavior. These observations have triggered interest in understanding the origin of criticality in specific dynamical models and its relationship with the underlying network architecture. The existing literature reveals that the statistical properties of avalanches in some of the above models may be dependent on the topology of the network on top of which the dynamics proceeds [5,[15][16][17], while some other authors suggest that critical features are independent of network topology [18]. Here, we aim at reconsidering the problem of critical scaling behavior on networks, for all the above dynamical models within a unified and consistent perspective.
In many cases, the dependence of avalanche statistics on the * filiradi@indiana.edu topology of the underlying network is theoretically explained by regarding the avalanche as the result of a simple branching process (BP) [19][20][21][22], where the out-degree distribution P(k out ) is identified with the distribution of offspring number within the network. Such a mapping is exact as long as the evolution of an avalanche does not substantially change the probability for an active node to find inactive neighbors to infect, e.g., when the substrate is a directed tree. The mapping usually fails for arbitrary networks: after some transient time an active site may find neighbors that are already active, so that the process is not merely branching out, but interfering with itself. Reasonably, this failure is more dramatic in undirected networks as the front of an avalanche can immediately move backwards and, as a consequence, break the equivalence with a simple branching process. According to the standard BP theory, sizes S and durations T of avalanches at criticality -i.e., when on average there is one offspring per active node-are distributed according to power laws where τ and α are avalanche critical exponents, and G S (S /S C ) and G T (T/T C ) are cut-off (scaling) functions, with the cutoff scales, S C and T C , depending only on system/network size right at the critical point [23,24]. Moreover, the average avalanche size scales with its duration as S ∼ T θ , where the exponent θ obeys the general scaling relationship θ = (α − 1)/(τ − 1) [25,26].
The values of τ and α may depend on the offspring distribution, i.e., the probability for an active node to activate a given number of new nodes. If the second moment of such a distribution is finite, then τ = 3/2, α = 2 and θ = 2. (2) Eq. (2) defines the so-called "standard" mean-field (MF) or "branching process" exponents. Thus, this type of scaling is expected to emerge for critical avalanches in the case in which the second moment of P(k out ) in the network is finite. In fact, these values of τ and α are extremely universal and robust; they emerge in many different types of propagation processes such as directed percolation, CP, VOT, SIS, SIR and many others, as long as the underlying pattern of connections is either a high-dimensional lattice or a sufficiently homogeneous network [24,[27][28][29][30] 1 . This super-universality can be rationalized using a Langevin equation for the density ρ of active sites where F is a noise-amplitude constant and ξ(t) a zero-mean Gaussian white noise, which is shared, as an effective meanfield description, by all the above mentioned models [31]. Observe that the square-root term in Eq. (3) accounts for "demographic" fluctuations and is a direct consequence of the central limit theorem [30,35].
On the other hand, if the second moment of the offspring distribution -or, equivalently P(k out )-diverges, then the critical exponents of the associated branching process differ from the standard MF ones of Eq. (2). In particular, for P(k out ) ∼ [k out ] −γ with 2 < γ < 3, one obtains γ-dependent exponents [5,36,37], Observe that these "anomalous" branching process (ABP) exponents converge to those of Eq. (2) in the limit γ → 3, i.e., when the second moment of P(k out ) becomes finite (with the 1 A particularly simple proof of the emergence of the standard exponents when the underlying tree is homogeneous with k out = 2 can be found in Ref. [31]. A more systematic derivation -for different types of underlying regular or random tree topologies-can be obtained within the generatingfunction formalism [29,32,33]; for instance, already back in 1949, Otter computed the solution when P(k out ) is a Poisson distribution [34].
Real-world networks often exhibit power-law degree distributions P(k) ∼ k −γ with 2 < γ < 3 [39,40], and studies concerning spreading processes often assume underlying scalefree network topologies [2,3]. A naive extension of the standard BP mapping to networks with diverging second moment of the degree distribution suggests that one should generically observe anomalous exponents. Does anomalous scaling hold for avalanches on real networks?
The current literature on the existence of anomalous avalanche scaling in scale-free networks reports conclusions that are often contradictory or difficult to reconcile with each other. For example, according to Larremore et al. critical avalanches in networks are always characterized by the standard MF exponents of Eq. (2) irrespective of the underlying network topology [18]. However, numerical evidence in support of such a claim is presented only for power-law networks with degree exponent γ > 3. Furthermore, this is in apparent contradiction with what reported for the CIC model on directed scale-free networks. In particular, Gleeson et al. employ a map onto an anomalous branching process to argue that one should expect anomalous scaling for 2 < γ < 3 and standard MF critical exponents for γ > 3 [5], but offered limited computational evidence in support of such a claim.
Also for the broadly studied SIR model the current state of understanding is not entirely clear. First, the model is never studied directly; rather, claims follow from SIR equivalence with bond percolation [41], according to which the distribution of SIR avalanche sizes can be deduced from the percolation cluster-size distribution. Theoretical claims on bond percolation in scale-free networks mostly regard undirected networks [15][16][17]. This is a very difficult setting to consider given that the percolation threshold vanishes [42]. A largescale numerical study of the percolation cluster size distribution in scale-free graphs is the one of Ref. [43], where critical exponents seem compatible with the standard ones of Eq. (2) for any γ > 2.
The goal of this paper is to provide a coherent picture for avalanche statistics in critical processes taking place on networks. The study consists in extensive numerical simulations, combined with analytical arguments, of the various avalanche models on a variety of networks, both synthetic and real.

II. MODELS: NETWORKS AND DYNAMICS
The models studied here are described in detail in SM [1]. Figure 1 illustrates schematically the mechanisms at the basis of the various models under consideration. As a substrate for the dynamics of activity, we assume in all cases a network composed of N nodes. The topology of the network is fully specified by its adjacency matrix A, whose generic element A i j = 1 if an edge from node i to node j exists and A i j = 0 otherwise. We assume that no selfloops are present in the network, i.e., A ii = 0 for all i. In the most general case, we consider directed networks, where A i j A ji . For simplicity, we further assume that the network is composed of only one strongly connected component, so that at least one directed path between any pair of nodes exists.
The state of the system at time t is denoted by the vector is a discretevalued variable representing the state of node i at time t. In all models except for the SIR, σ i (t) can assume two values: σ i (t) = 1, 0 indicating that the node is active, inactive respectively. In the SIR one can also have that σ i (t) 0, 1 meaning that node i is recovered and does not participate any longer in the dynamics.
All models are stochastic Markov processes where the elementary reactions that lead to changes in system configurations are triggered by the random selection of network elements, either nodes, edges or both. Propensities of the various reactions may depend on exogenous parameters whose values can be tuned to bring the system into different dynamical regimes. The state σ = 0 = (0, 0, . . . , 0) T is an absorbing configuration for all models. Additional absorbing configurations are present in some models. For example, in the SIR model, all configurations with no infected nodes, but arbitrary number of recovered nodes are absorbing; in some other models, such as the CIC, the configuration σ = 1 = (1, 1, . . . , 1) T is also absorbing at the critical point. The detailed definitions of all the considered models are reported in the SM [1].
We are interested in the critical regimes of the considered dynamical models. The criterion to achieve criticality is model specific. VOT, IP and LD have no parameters and are intrinsically critical. CIC critical point is achieved by set-ting model parameters to network-independent values. For CP a known network-independent value is a good approximation of it. The critical point of SIS is approximated by considering the inverse of the largest eigenvalue of the adjacency matrix of the graph [18,44]. For SIR, the critical regime is approximated relying on the value of the largest eigenvalue of the non-backtracking matrix of the graph [45,46] (see SM [1] for further details).
We consider avalanches initiated by a single randomlychosen node, j, so that the initial configuration is σ i (0) = 0 for all i j, and σ j (0) = 1. We follow the dynamics of each avalanche until the system reaches an absorbing configuration. We define the duration T as the number of time steps needed to enter into an absorbing configuration. We also define the size S as the number of elementary spreading events occurred during T . An elementary spreading event is the occurrence of an active node passing activation to other, not necessarily inactive, nodes.
All nodes may participate multiple times in an avalanche (i.e., they can be "re-activated") so that the network size N is not an upper bound for S ; the only exception to this rule is the SIR model where nodes can be activated only once. We focus on finite avalanches only, i.e., those avalanches that end up in absorbing configuration σ = 0. In the CIC model, for instance, we exclude avalanches that end in σ = 1, as they can be viewed as infinite avalanches. For the SIR model, we consider instead all avalanches.
We study, by means of extensive computational simulations, avalanche statistics for all the above-mentioned prototypical dynamical models on top of scale-free networks generated by one of two possible standard generative models. Both models produce uncorrelated random graphs with power-law degree distributions.
First, we consider undirected scale-free networks obtained via the uncorrelated configuration model [47] with degree distribution P(k) ∼ k −γ with support [4, √ N]. In our numerical analyses, we set γ = 2.1 and vary the network size N. The choice γ = 2.1 is expedient because it corresponds to a large gap between standard MF and ABP exponents, easing computational discrimination of scaling regimes. Results for γ = 2.5 are reported in the SM [1]. We generated a single graph instance of the model for every N, and used these graphs in all our analyses. We tested that choosing a particular instance of the graph model does not affect the statistics of avalanches. For every network and model, we simulated 10 6 avalanches seeded in a randomly chosen single node and measured the corresponding avalanche size and duration distributions, i.e., P(S ) and P(T ), respectively.
The first major result of our analyses (see Fig. 2) is thatwhen deployed on undirected scale-free networks-all activation models considered are characterized asymptotically (for large S and large T ) by standard MF exponents. This happens regardless of the fact that we have set 2 < γ < 3, i.e., for networks for which a strict analogy with branching process would suggest anomalous exponents. An exhaustive report of the results of our analysis is shown in the SM [1].
Second, we analyze directed scale-free networks constructed according to the model of Ref. [5]. This is a sim-  A more complicated scenario emerges when k out max = N − 1. In Fig. 3, we show the results only for few selected models. Results for all other models are reported in SM [1]. For IP and 3 In the generation of a graph instance, each node i is connected to k out i other nodes, chosen at random in the network, so that in-degrees obey a Poissonian distribution with average value equal to the average out-degree. CP, we still observe clear MF scaling, in all networks. In all other models, the distribution of the avalanche sizes displays a crossover from anomalous (for small sizes/durations) to MF exponents (for large sizes/durations). The crossover point increases with the network size, suggesting that only anomalous exponents should be present in the limit of asymptotically large networks 4 . 4 Let us stress that the observed scaling of P(T ) and power-law relation between S and T provide much less clear evidence of anomalous BP behavior even for the case k out max = N − 1 (see SM [1]). This issue is due to the finite size of the networks, and it is visible also in numerical results concerning pure BP with finite-size constraints (SM [1]). Clearer observations of anomalous critical exponents for P(T ) and the power-law relation between S and T can be obtained for γ = 2.5; such a choice of the γ value leads however to much less noticeable differences between anomalous and standard exponents for the distribution P(S ) (SM [1]).
In summary, our results provide strong support for MF exponents in all situations where a priori we expect standard BP behavior (i.e., finite second moment of the out-degree distribution). Settings for which one could predict a priori anomalous BP scaling (i.e. for scale-free networks with 2 < γ < 3) generate results that are much less cleancut. Anomalous exponents can at most be observed only in the regime of small avalanches for the distribution of the avalanche size. Strong deviations from the predicted anomalous power-law scalings are observed otherwise.

III. ANALYTICAL APPROACH
Two alternative types of approaches are frequently used to study avalanches in networks: the first one is, as discussed above, a branching process approximation for cascades of a small size at the beginning of the process, and the second one is an approximation by a dynamical system once the avalanche spreads to a significant fraction of the underlying network [48]. The mathematical approach we develop in what follows belongs to this second group.
Depending on their specific features, all dynamical models under consideration can be grouped into three classes described by different types of mathematical equations. These groups are: (i) IP and CP; (ii) CIC, VOT, and LD; (iii) SIS and SIR.
The first group is trivially described by dynamical processes that are insensitive to the out-degree sequence of the underlying network. In other words, anomalous propagation events in which a single active node propagates activity to an arbitrarily large number of nearest neighbors are simply not allowed by the dynamics. Henceforth, anomalous type of scaling is not expected to appear, even at the level of a naive mapping onto a branching process. Thus, IP and CP avalanches are expected to be always characterized by standard MF exponents (see SM [1]), in perfect agreement with our computational results.
Models in the other two classes have instead a much less trivial behaviour. We consider the CIC and SIS models as representative of each of these two classes and derive a mathematical approach for each. The full development of the theory (and extension to the other dynamical models) is presented in the SM [1]; here, we sketch the main results and the main insights derived from them.
Our analytical approach is based on two successive approximations. The first one is the so-called individual-based meanfield approximation (IBMFA) (see Pastor-Satorras et al. [3] for a review). This analysis starts by describing the evolution of the average value of the state of an individual node in the network, where the average is taken over many realizations of the dynamical process. The approximation consists in neglecting dynamical correlations among variables, so that every node feels only the influence of the average behavior of each of its neighbors. We use the IBMFA for determining how and when the system reaches its long-term dynamical regime. The second approximation consists in deriving a Langevin equation for the overall network activity, written as the sum of the activity variables of all nodes [49]. From this approximation, it is possible to derive the statistics of longterm avalanches based on the equivalence between the resulting Langevin equations and Eq. (3), derived in Ref. [31] to describe the standard branching process and related processes.
Let us first present the derivation of the main results for CIC critical dynamics, in which the only possible change in the state of a node consists in copying the state of a nearest neighbor. In the IBMFA, we focus our attention on the deterministic node variable s i (t) := σ i (t) , defined as the value of the stochastic variable σ i (t), averaged over the realizations of the dynamical process at time t. As we explicitly derive in the SM [1], critical CIC dynamics is described by the IBMFA equation Here, [ s(t)] T = [s 1 (t), s 2 (t), . . . , s N (t)] is the vector describing the average state of the nodes of the network at time t. L = K in − A is the (directed) graph Laplacian of the network, with K in the diagonal matrix whose non-null elements are equal to the in-degree of the nodes, and A is the graph adjacency matrix [50]. In essence, under the IBMFA, the critical CIC coincides with a purely diffusive process. The properties of the solutions of Eq. (5) for arbitrary graphs are described in Refs. [51,52]; we briefly summarize them here. If the underlying network is composed of a single strongly connected component, then, the long-term behavior is such that lim t→∞ s i (t) = s * , for all nodes i. Because the asymptotic limit of the individual variables s i does not depend on i, s * coincides with the asymptotic value of density, r * . The latter is given by the norm of the vector s(t) at large times, i.e., where v (l) 1 and v (r) 1 are the left and right eigenvectors of L T , respectively. These eigenvectors correspond to the eigenvalue ν 1 = 0 of the graph Laplacian. The long-term regime is reached exponentially fast. However, depending on the type of network, we have different diffusion behaviors. For undirected networks, v (l) 1 = v (r) 1 = 1/ √ N; further, the density of active nodes r(t) := 1/N i s i (t) is such that r(t) = r(0) = 1/N for all t. The typical time scale is t * = 1/ν 2 , with ν 2 the second smallest eigenvalue of L (see SM [1]).
On the other hand, if the network is directed, v (r) 1 . This means that, also in this case the vector s has identical components for t → ∞. However, r(t) may increase or decrease depending on the initial condition, so that the steady-state value r * of the density is sensitive to the initial choice of the seed node 5 . Further, in this case, we cannot longer apply the spectral theorem to the corresponding non-symmetric matrix so that the exponential relaxation to the steady-state cannot be easily written in terms of the Laplacian eigenvalues.
To determine the statistical properties of avalanches with duration T t * , we now go back to the stochastic description of the full dynamical system. We take advantage of the previous finding obtained under the IBMFA, and assume that ρ(t) := 1/N i σ i (t) is a quantity that fluctuates around its average value ρ = r * Essentially, we make the hypothesis that the system has reached a stationary state where the number of active nodes is on average constant, but still subjected to demographic fluctuations. In analogy with Ref. [49], we refer to this assumption as the adiabatic approximation. The dynamics of long-term avalanches in critical CIC dynamics thus turns out to obey the following Langevin equation (see where ξ(t) is a zero-mean Gaussian white noise, and k in is the average in-degree of the network. The dependence on ρ of the diffusion coefficient imposes the absence of fluctuations for both ρ = 0 and ρ = 1, corresponding to the two existing absorbing states. Except for higher-order terms, Eq. (6) has the generic form of the representative Langevin equation, Eq.
(3) for avalanches in the class of standard MF branching processes [31]. This implies that long-term avalanches in critical CIC dynamics obey power-law distributions with MF critical exponents, i.e., Eqs. (2). We now briefly illustrate the analytical approach for SIS critical dynamics. We basically repeat the same steps described above for critical CIC dynamics. The IBMFA equation reads as where I is the identity matrix [44]. The solution of the IBMFA equation is a vector whose components are proportional to those of the principal right eigenvector w (r) N of the matrix A T [44]; convergence to the asymptotic solution is exponentially fast. The asymptotic value of density of active nodes is given by r * = ( w (l) N ) T · s(t = 0)] w (r) N , with w (l) N and w (r) N principal left and right eigenvector of the matrix A T respectively 6 . If the network is undirected, the time scale of the exponential relaxation to the steady-state density is given by t * = ω N /(ω N − ω N−1 ), with ω N largest eigenvalue of A, and ω N−1 second largest eigenvalue of A (see SM [1]). If the network is directed, t * is not directly quantifiable in terms of the eigenvalues of the matrix A.
For t t * , the system has reached its long-term dynamical regime. The statistics of long avalanches is described by the Langevin equation 6 In uncorrelated random network models, the components of the vector w (l) N are proportional to the node out-degrees, i.e., w (l) 1,i ∼ k out i (see SM [1]). For undirected configuration models, the previous statement is valid only when the degree exponent γ < 5/2. where w (r) N is the average value of the components of the principal right eigenvector of the matrix A T (see SM [1]). Eq. (8) has the same form as those considered by di Santo et al. [31], valid for avalanche models that are equivalent to standard BP processes. This tells us that long-term avalanches in critical SIS dynamics obey power-law distributions with MF critical exponents, i.e., Eqs. (2).
In summary, the above analytical approach tells us that sufficiently long (large) avalanches in critical CIC and SIS dynamics should follow a standard MF scaling. This conclusion, in principle is true for any network. However, an avalanche is sufficiently long to obey standard BP statistics only if its duration is much longer than the typical time scale that can be deduced from the IBMFA of the process happening on the network. The magnitude of such time scale depends exclusively on the topology of the network, by means of either the Laplacian or the adjacency matrix of the graph. Undirected networks with sufficiently short average distance, for instance, have a relatively small diffusion time scale. There are, however, network topologies where diffusion may be particularly slow to reach its stationary state. Examples are lowdimensional lattices, and networks with long loops. In these cases, the vast majority of observed avalanches may never be long enough to be described by the long-term statistics. We do not have analytical arguments to determine the statistical properties of avalanches in the short-term dynamical regime, but, in principle, one expects that the effective mapping into an ABP should work (for scale-free networks with 2 < γ < 3). Our numerical results seem to indicate that anomalous BP scaling is possible as long as the underlying networks are directed and have power-law out-degree distributions. The cutoff of the out-degree distribution seems also to play an important role for the possible observation of anomalous scaling, at least for the network sizes that we were able to consider in our analysis.

IV. REAL NETWORKS
All considerations, numerical and theoretical, made for synthetic graphs are valid also for real-world networks. In Fig. 4, we summarize the results of numerical simulations performed on three large-scale networks. Additional results are provided in the SM [1]. A pre-asymptotic regime with anomalous scaling for sufficiently small avalanches is seen for example in the Youtube direct social network (Fig. 4c). The distributions for large avalanches are instead very well described by MF critical exponents in all cases. The crossover point may be interpreted as the typical scale that distinguishes local from global avalanches, and it could be employed as a quantitative criterion to tell whether an avalanche is "viral" or not.

V. CONCLUSIONS AND DISCUSSION
In conclusion, we found that any minimal deviation from the assumptions underlying the BP framework brings the system back to the realm of standard MF and its associated super-  Figure 4. Avalanche size in real-world networks. We consider the following networks: (a) undirected graph representing a snapshot of the Internet at the Autonomous system level [53]; (b) directed Twitter network of the Spanish 15M movement [54]; (c) directed graph representing a portion of the Youtube social network [55]. Different symbols and colors refer to different avalanche dynamical models. The red dashed line represent standard BP critical exponents, while the full black line indicates the power-law decay expected for anomalous BP. Note that the out-degree distributions of these networks are all well modeled by power laws with decay exponent γ = 2.1 (see SM [1]). universal exponents, so that anomalous exponents are exceedingly difficult to be observed. We proved this statement to be true for seven well-known avalanche dynamical models, but we believe that it can be extended to many other spreading processes taking place on networks.
Why are is numerical evidence of anomalous scaling so weak, even in the cases when intuition suggests that the dynamical avalanche model could be well mapped to an anomalous branching process? Clearly, if the process is occurring on a directed tree with power-law out-degree distribution, then anomalous scaling occurs. However, avalanche models in networks do not necessarily satisfy such strict conditions. There are many possible ways in which the assumptions of the mapping to an ABP can be violated.
First, both in directed and undirected networks avalanches do not necessarily proceed in a fully feedforward way; already active nodes can be found by a branch of an unfolding avalanche thus breaking the equivalence with a pure branching process. In other words, feedforward loops may exist, meaning that a given node can be reached from a unique seed by following different paths. This is particularly relevant in undirected networks, where activity can attempt to go backwards after any propagation event, following a reversible link. This type of interference reduces the effective number of independent offspring, breaking the BP analogy.
Second, networks in simulations are finite, implying that a finite maximum degree exists, therefore the out-degree variance takes a finite value; this implies that there should be crossovers to the standard exponents for sufficiently large avalanche sizes and durations.
Last but not least, some types of dynamics, even if taking place on top of scale-free networks, do not really involve all neighbors of a single node -as for example in the CP and IP processes-and, thus, have an effective offspring distribution narrowly distributed, implying the emergence of standard MF exponents.
Numerous real-world systems have been investigated in terms of avalanche statistics. Prototypical examples include natural systems, such as neuronal networks [56], γ-ray bursts [57] and earthquakes [58], as well as socio-technical systems, such as power networks [59] and online social media [60][61][62][63]. Among them, some systems display critical avalanche statistics consistent with the MF universality class [56]. However, there are many other systems showing avalanche statistical properties that are not consistent with the MF universality class. Examples can be found especially in the literature studying information avalanches in online social media where measured exponents for the power-law distribution of avalanche size range from τ 4 [60], to τ 2.3 [62] and τ 2 [61].
Our analytical and numerical evidence supports the existence of an extremely robust universality class at the interface between the absorbing and the active phases of many popular models of avalanche dynamics. Such a universality class can be broken only at the expense of making specific assumptions on the shape of the network underlying the spreading model. We believe that it is imperative to understand why there exist real systems that do not conform to such a class, and what alternative hypotheses need to be made to account for their behavior. In other words, a complete analytical theory -extending the approach presented here-and accounting for all types of networks still needs to be constructed.
As a final note, let us stress that our results reveal that constructing a dynamical model characterized by MF avalanche exponents is not a hard task. Consequently, having a model that generates avalanche distributions with MF exponents -being these in agreement with some experimental observation-does not constitute a sufficient evidence that the model is actually right. Other dynamical aspects should be also used to validate the model.