Critical Network Cascades with Re-excitable nodes: Why tree-like approximations usually work, when they breakdown, and how to correct them

TL;DR

This paper rigorously analyzes the accuracy of tree-like approximations in network cascade models with re-excitable nodes, identifying when they work well and how motifs cause deviations, with implications for real-world network analysis.

Contribution

It provides a theoretical framework to quantify deviations from tree-like approximations in network cascades and derives correction terms for non-tree-like motifs, especially the bi-parallel motif.

Findings

Tree-like approximations are accurate for networks with higher average degree.

Motifs like the bi-parallel cause deviations from tree-based theories.

Corrections improve predictions of network criticality in complex networks.

Abstract

Network science is a rapidly expanding field, with a large and growing body of work on network-based dynamical processes. Most theoretical results in this area rely on the so-called \emph{locally tree-like approximation}. This is, however, usually an `uncontrolled' approximation, in the sense that the magnitudes of the error are typically unknown, although numerical results show that this error is often surprisingly small. In this paper, we place this approximation on more rigorous footing by calculating the magnitude of deviations away from tree-based theories in the context of discrete-time critical network cascades with re-excitable nodes. We discuss the conditions under which tree-like approximations give good results for calculating network criticality, and also explain the reasons for deviation from this approximation, in terms of the density of certain kinds of network motifs. Using this understanding, we derive results for network criticality that apply to general networks that explicitly do not satisfy the locally tree-like approximation. In particular, we focus on the bi-parallel motif, the smallest motif relevant to the failure of a tree-based theory in this context, and we derive the corrections due to such motifs on the conditions for criticality. We verify our claims on computer-generated networks, and we confirm that our theory accurately predicts the observed deviations from criticality. Using our theory, we explain why numerical simulations often show that deviations from a tree-based theory are surprisingly small. More specifically, we show that these deviations are negligible for networks whose average degree is even modestly large compared to one, justifying why tree-based theories appear to work well for most real-world networks.

Equations

$$A=p\times \tilde{A},$$

$$x_i(t)=H\left[\sum _{j=1}^N{x}_j(t-1)s_{ji}(t)\right],$$

$$\lambda =p\times \tilde{\lambda }.$$

$${\rho }_n=\underset{t\ge 0}{\min }\left\{t|x_k(t)=0\forall k\right\}.$$

$$s_n=\sum _{t=0}^{\infty }\sum _{k=1}^N{x}_k(t).$$

$$c_n(t)=\mathrm{Pr}({\rho }_n\le t),$$

$$c_n(t)=\prod _m[(1-A_{nm})+A_{nm}{c}_m(t-1)].$$

$${\lambda }_c=1\text{ or }{p}_c=(\tilde{\lambda }{)}^{-1},$$

$$c_{n_0}(t)=\prod _{(n_0,n_1)\in \mathcal{E}}[(1-p)+p{c}_{n_1}(t-1)],$$

$$0=c_{n_0}(t)-{\left[(1-p)+p\frac{1}{d(n_0)}\sum _{(n_0,n_1)\in \mathcal{E}}{c}_{n_1}(t-1)\right]}^{d(n_0)}.$$

$$\prod _{i=1}^m{X}_i$$

$$\approx ⟨X{⟩}^m,$$

$$0$$

$$=\frac{\partial c_{n_0}(t)}{\partial b}-\sum _{(n_0,n_1)\in \mathcal{E}}\sum _{(n_1,n_2)\in \mathcal{E}}{p}^2\times \frac{\partial c_{n_2}(t-2)}{\partial b}\times \xi (t|1)\times \xi (t|2),$$

$$⋮$$

$$=\frac{\partial c_{n_0}(t)}{\partial b}-\sum _{(n_0,n_1)\in \mathcal{E}}\dots \sum _{(n_{(k-1)},n_k)\in \mathcal{E}}{p}^k\times \frac{\partial c_{n_k}(t-k)}{\partial b}\times \left\{\xi (t|1)\times \dots \times \xi (t|k)\right\},$$

$$\xi (t|j)={\left[(1-p)+p\frac{1}{d(n_{(j-1)})}\sum _{(n_{(j-1)},n_j)\in \mathcal{E}}{c}_{n_j}(t-j)\right]}^{d(n_{(j-1)})-1}.$$

$$0=1-\underset{k\to \infty }{\lim }\sum _{(n_0,n_1)\in \mathcal{E}}\dots \sum _{(n_{(k-1)},n_k)\in \mathcal{E}}(p_c{)}^k,$$

$$p_c=\underset{k\to \infty }{\lim }{\left(\prod _{l=0}^{k-1}⟨d(n_l)⟩\right)}^{-1/k},$$

$$p_c=d_{\infty }^{-1}.$$

$$c_{n_0}(t)=\sum _{k=0}^{d(n_0)}\sum _{{\mathcal{E}}_k^0}(1-p{)}^{d_{n_0}-k}{p}^k\prod _{v\in {\mathcal{V}}_k^0}{c}_v(t-1),$$

$$\prod _{v\in {\mathcal{V}}_k^0}{c}_v(t-1)$$

$$\prod _{v\in {\mathcal{V}}_k^0}{c}_v(t-1)=\sum _{r=0}^{D[{\mathcal{V}}_k^0]}\sum _{{\mathcal{E}}_r^1}(1-p{)}^{D[{\mathcal{V}}_k^0]-r}{p}^r\prod _{w\in {\mathcal{V}}_r^1}{c}_w(t-2),$$

$$\begin{gathered} c_{n_0}(t)=\sum _{k=0}^{d(n_0)}\left[\sum _{{\mathcal{E}}_k^0}(1-p{)}^{d(n_0)-k}{p}^k \\ \times \sum _{r=0}^{D[{\mathcal{V}}_k^0]}\left(\sum _{{\mathcal{E}}_r^1}(1-p{)}^{D[{\mathcal{V}}_k^0]-r}{p}^r\prod _{w\in {\mathcal{V}}_r^1}{c}_w(t-2)\right)\right], \end{gathered}$$

$$\begin{gathered} \delta r=\sum _{\alpha }{Q}_{\alpha }(n_0){\left(\frac{D[{\mathcal{V}}_k^0]}{d(n_0)⟨d(n_1)⟩}\right)}^{\alpha } \\ \times \left[\alpha \left(\frac{r}{D[{\mathcal{V}}_k^0]}\right)-1+{\left(1-\frac{r}{D[{\mathcal{V}}_k^0]}\right)}^{\alpha }\right]. \end{gathered}$$

$$\delta r=\left\{\frac{Q_2(n_0)}{[d(n_0)⟨d(n_1)⟩{]}^2}+\frac{3Q_3(n_0)}{[d(n_0)⟨d(n_1)⟩{]}^3}D[{\mathcal{V}}_k^0]\right\}{r}^2-\frac{Q_3(n_0)}{[d(n_0)⟨d(n_1)⟩{]}^3}{r}^3.$$

$$\frac{\partial c_{n_0}(t)}{\partial b}=\sum _{k=0}^{d(n_0)}\left\{\sum _{{\mathcal{E}}_k^0}(1-p{)}^{d(n_0)-k}{p}^k\sum _{r=0}^{D[{\mathcal{V}}_k^0]}\left[\sum _{{\mathcal{E}}_r^1}(1-p{)}^{D[{\mathcal{V}}_k^0]-r}{p}^r\left(\left\{\prod _{w\in {\mathcal{V}}_r^1}{c}_w(t-2)\right\}\sum _{w\in {\mathcal{V}}_r^1}\frac{1}{c_w(t-2)}\frac{\partial c_w(t-2)}{\partial b}\right)\right]\right\}.$$

$$\begin{gathered} 1=\sum _{k=0}^{d(n_0)}\left\{\sum _{{\mathcal{E}}_k^0}(1-p{)}^{d(n_0)-k}{p}^k\times \\ \sum _{r=0}^{D[{\mathcal{V}}_k^0]}\left[\left(\genfrac{}{}{0ex}{}{D[{\mathcal{V}}_k^0]}{r}\right)(1-p{)}^{D[{\mathcal{V}}_k^0]-r}{p}^r\left(r-\delta r\right)\right]\right\}, \end{gathered}$$

$$1=\left(\frac{p^2}{d(n_0{)}^2⟨d(n_1){⟩}^3}\right)\times \left\{d(n_0{)}^3⟨d(n_1){⟩}^4-d(n_0)⟨d(n_1)⟩Q_2(n_0)\left[⟨d(n_1)⟩-p⟨d(n_1)⟩+p⟨d(n_1{)}^2⟩+p^2⟨d(n_1){⟩}^2(d(n_0)-1)\right]+Q_3(n_0)\left[-3⟨d(n_1{)}^2⟩(p-1{)}^2+p⟨d(n_1{)}^3⟩(p-3)-3p(p-1{)}^2⟨d(n_1){⟩}^2(d(n_0)-1)+p^3(p-3)⟨d(n_1){⟩}^3(d(n_0)-1)(d(n_0)-2)+⟨d(n_1)⟩(1-3p+2p^2+3p^2(p-3)(d(n_0)-1)⟨d(n_1{)}^2⟩)\right]\right\}.$$

$$1=\left(\frac{p_c^2}{d_{\infty }^5}\right)\times \left\{d_{\infty }^7-d_{\infty }^2\overline{Q_2}\left[d_{\infty }-p_c{d}_{\infty }+p_c⟨d^2⟩+p_c^2{d}_{\infty }^2(d_{\infty }-1)\right]+\overline{Q_2}\left[-3⟨d^2⟩(p_c-1{)}^2+p_c⟨d^3⟩(p_c-3)-3p_c(p_c-1{)}^2{d}_{\infty }^2(d_{\infty }-1)+p_c^3(p_c-3)d_{\infty }^3(d_{\infty }-1)(d_{\infty }-2)+d_{\infty }(1-3p_c+2p_c^2+3p_c^2(p_c-3)(d_{\infty }-1)⟨d^2⟩)\right]\right\},$$

$$\sigma =\frac{{\sum }_{\alpha }[\overline{Q_{\alpha }}\times (\alpha -1)]}{d_{\infty }^2}.$$

$$\mathrm{Pr}(|{\mathcal{V}}_r^1|)=\left(\genfrac{}{}{0ex}{}{d}{|{\mathcal{V}}_r^1|})\times (\genfrac{}{}{0ex}{}{r-1}{|{\mathcal{V}}_r^1|})/(\genfrac{}{}{0ex}{}{r+d-1}{d-1}\right),$$

Full text

I Introduction

The study of dynamical processes on networks is among the most important areas of research in network scienceBarrat et al. (2008); Boccaletti et al. (2006); Porter and Gleeson (2016). Theoretical understanding of these processes on networks found in real-world studies has wide potential impact, owing to the applicability of such systems to fields of study as diverse as epidemiologyMiller (2009); Valdano et al. (2015); Prakash et al. (2012), neuroscienceShew et al. (2009); Kinouchi and Copelli (2006); Tanaka et al. (2009), ecologyCohen et al. (2012); Sole and Montoya (2001), electrical engineering Motter and Lai (2002), social sciencesDraief and Massouli (2010); Iribarren and Moro (2011) and several others. Obtaining rigorous analytical results for such systems on networks that are found in nature is generally very hard. At the present time a majority of central results in this area rely on the so-called ‘locally tree-like approximation’, which neglects the effect of small loops and cycles in the network structureNewman (2018); Strogatz (2001); Dorogovtsev et al. (2008). However, despite being commonly used, this approximation is uncontrolled, i.e., there do not exist clear indications of the extent of validity of the approximation, nor are estimates of the scaling of the order of error to be expected typically available. This issue is additionally compounded due to two observations: first, real-world networks tend to have a high clustering with a significant density of small loops and are hence far from tree-likeNewman (2018); Strogatz (2001); and, second, theoretical results on the basis of locally tree-like approximations appear to be in close agreement with numerically obtained results for dynamical processes on such far from tree-like real-world networksMelnik et al. (2011). This leads to the central foundational questions of our paper: Why do locally tree-like approximations appear to work well on real-world networks, when can we safely use these approximations, and what are the sizes of expected errors in the approximation? In this paper we work in the context of cascade processes on networks with re-excitable nodes, and aim to put the locally tree-like approximation on more rigorous footing by calculating the expected order of error in using this approximation. We also discuss the application of our results to other types of network dynamics in our discussion in Sec. VI

We restrict our analysis to the important dynamical process of network cascades. Cascade processes on networks, also referred to as avalanches on networks, have been widely studied due in part to their wide range of applicability, including problems relating to epidemiologyMiller (2009); Valdano et al. (2015); Prakash et al. (2012), neuroscienceShew et al. (2009); Kinouchi and Copelli (2006); Larremore et al. (2011a); Tanaka et al. (2009), genealogyWatson and Galton (1875), social sciencesDraief and Massouli (2010); Iribarren and Moro (2011), and network securityAcemoglu et al. (2016). In a cascade process, if the average number of nodes excited by a single node is sufficiently large, then small initial perturbations can give rise to activity that may persist indefinitely. In contrast, if the average number of single-event induced excitations is too small, then cascades die out and the network can suppress the future activity resulting from even large initial perturbations. At the boundary between these two phases, the network is said to be critical. Network criticality has been studied in relation to a wide range of phenomena, such as percolation thresholds, epidemic thresholds and phase transitions in cooperative models, among othersDorogovtsev et al. (2008). In the particular context of neuronal networks, several studies suggest that networks of neurons tend to operate in this critical regimeShew et al. (2009); Poil et al. (2008); Beggs and Plenz (2003), which admits an increased dynamical rangeShew et al. (2009) and enlarged information capacityShew et al. (2011). Such critical phenomena are generally characterized by the presence of power-law statistics in various relevant distributions Beggs and Plenz (2003).

Watson and GaltonWatson and Galton (1875) first studied the problem of a branching processes on a network having a tree topology in the context of the extinction of family names, and examined the case for which the ‘cascade’ of family names would die out. Here, we consider a related problem of cascades on networks with general topologies. Motivated by the case of neuronal cascades, we specifically consider the situation in which nodes can be re-excited multiple times during the same cascade, similar to discrete SIS (susceptible-infected-susceptible) models in epidemiology (see for example Refs.Shi et al. (2008); Parshani et al. (2010); Boguná and Pastor-Satorras (2002)).

The problem of criticality in branching processes on networks with general network topologies has been studied over the last several years, resulting in some analytical results regarding the conditions for criticalityGómez et al. (2010); Larremore et al. (2012); Pastor-Satorras and Vespignani (2001); Larremore et al. (2011a, b); Assis and Copelli (2008); Kinouchi and Copelli (2006). However, as discussed above, these results generally rely on ‘locally tree-like approximations’Gómez et al. (2010); Larremore et al. (2012, 2011a); Kinouchi and Copelli (2006), or pertain to specific classes of networks having specific topologiesMoore and Newman (2000); Pastor-Satorras and Vespignani (2001); Assis and Copelli (2008). The ‘locally tree-like approximation’Newman (2018) is the assumption that subgraphs that extend a short distance from a given node are typically trees. For a wide range of dynamical processes on networks, the results predicted by a tree-based theory (i.e., for networks satisfying the locally tree-like approximation) have been shown to often be close to results obtained for real networksMelnik et al. (2011) that are not tree-like. In this paper, we provide analytic justifications for the success of the locally tree-like approximation in predicting when to expect critical cascades, analyze the reasons for the breakdown of such a tree-based theory, and develop an approach for capturing deviations from tree-like behavior.

We use the framework previously employed by Larremore et al.Larremore et al. (2012) to analyze the conditions for criticality in a general network in the thermodynamic limit of large network sizes. We demonstrate that network motifs that we call $k$-parallel motifs (see Fig.1) are especially relevant to the failure of tree-like approximations, and we study the corrections introduced due to such motifs. In particular, we discuss how the bi-parallel motif is the most relevant motif for deviations away from tree-like behavior, and we derive the condition for criticality for networks containing such motifs. Our results demonstrate why the tree-based theory always gives good results for real-world networks that are not locally tree-like and may have high clustering — this is consistent with previous numerical studiesMelnik et al. (2011). Notably, we explain why critical dynamics will always be consistent with tree-like approximations in networks with even modestly large average degree. We show that when the fourth power of the average degree is large compared to 1, critical dynamics will be consistent with tree-like approximations, even when the networks is far from tree-like. Empirically, we see strong agreement with tree-like approximations in networks with average degree as small as 10. We also consider networks constructed to have a large number of bi-parallel motifs so as to demonstrate an observable deviation from the tree-based theory for criticality and we find a close agreement of our derived results with numerical experiments.

The remainder of this paper is organized as follows. In Sec. II we detail the model of network cascades that we consider. In Sec. II.1, we discuss the locally tree-like approximation and briefly describe results obtained under this approximation. We then discuss the reasons for deviations away from such a theory in terms of relevant network motifs (Sec. II.2), and introduce the ‘second-level approximation’ that we develop to estimate criticality in networks that are explicitly allowed to break the locally tree-like approximation via the presence of bi-parallel motifs (Sec. III). We then present numerical results demonstrating the validity of our claims (Sec. IV), followed by details of the derivations for the results under both the locally tree-like approximation and the second-level approximation (Sec. V). We conclude in Sec. VI with a summary of our results and a discussion of limitations of our techniques, as well as extensions to other network-based dynamical processes.

II Network cascade dynamics

Our basic setup is motivated by previous work of Larremore et al.Larremore et al. (2012). We consider excitable dynamics in a network of nodes, using the following definitions:

• •

Network setup

• $\circ$

$N$ : Number of nodes in the network

• $\circ$

$\widetilde{A}$ : Adjacency matrix of the network; $\widetilde{A}_{mn}=1$ if there is a directed edge from node $m$ to node $n$, and $\widetilde{A}_{mn}=0$ otherwise.

• $\circ$

$\widetilde{\lambda}$ : Perron-Frobenius eigenvalue of $\widetilde{A}$

• $\circ$

$x_{n}(t)$ : The state (either [math] or $1$) of vertex $n$ at time $t$

• •

Network dynamics

• $\circ$

$p$ : Transmission probability of excitation across any edge; constant for all edges

• $\circ$

$A$ : Probability weight matrix; $A=p\times\widetilde{A}$

• $\circ$

$\lambda$ : Perron-Frobenius eigenvalue of $A$; $\lambda=p\times\widetilde{\lambda}$

• $\circ$

$\rho_{n}$ : Duration of an avalanche starting at node $n$

• $\circ$

$c_{n}(t)$ : $\Pr(\rho_{n}\leq t)$

• $\circ$

$b_{n}$ : Probability that $\rho_{n}$ is finite; $b_{n}=\lim_{t\to\infty}c_{n}(t)$

• $\circ$

$s_{n}$ : Size of an avalanche starting at node $n$; $s_{n}=\sum_{t=0}^{\infty}\sum_{k=1}^{N}x_{k}(t)$

• •

Critical network cascades

• $\circ$

$p_{c}$ : Transmission probability such that the network exhibits criticality

• $\circ$

$\lambda_{c}$ : Perron-Frobenius eigenvalue of $A$ such that the network exhibits criticality; $\lambda_{c}=p_{c}\times\widetilde{\lambda}$

We consider a directed, unweighted, strongly-connected (i.e., every node is reachable from every other node) network of $N\gg 1$ nodes, labeled by the integers from $1$ to $N$, connected according to an adjacency matrix $\widetilde{A}$. We also assume our network to have no self-loops or multiple edges. We further assume discrete time network dynamics, in which the state of the $n^{\text{th}}$ node at time $t$ is represented by $x_{n}(t)$, which can take a value of either [math] or $1$. If $x_{n}(t)=1$ ($x_{n}(t)=0$) then the node $n$ is said to be in the active or excited (inactive or resting) state at time $t$. The dynamics of activation and deactivation on the network are governed by a probability weight matrix,

[TABLE]

where $p$ is a transmission probability, $0<p<1$. If a node $m$ is active at time $t-1$, it sends an activating signal at time $t$ to node $n$ with probability $A_{mn}$. If, at time $t$, node $n$ receives an activating signal from any node it is activated, and all edges that have sent activating signals are said to be active at that time. An active node will relax back to the inactive state at the next time step unless it is activated again. To summarize, we can write the dynamics of cascade propagation that we consider as

[TABLE]

where for each $i$, $j$ and $t$, the quantity $s_{ji}(t)\in\{0,1\}$ is an independent Bernoulli random variable such that $\Pr[s_{ji}(t)=1]=A_{ji}$ and $\Pr[s_{ji}(t)=0]=1-A_{ji}$, and $H(x)=1$ for $x>0$ and $H(x)=0$ for $x=0$. Note that in our model the state of each node is updated in parallel at each timestep.

We denote the Perron-Frobenius eigenvalue (i.e., the eigenvalue having the largest magnitude, which is unique, real, and also referred to as the leading or dominant eigenvalue) of $\widetilde{A}$ as $\widetilde{\lambda}$, and similarly denote the Perron-Frobenius eigenvalue of $A=p\times\widetilde{A}$ as

[TABLE]

In practice, we use $\lambda$ as the tunable parameter to vary the edge-weights $A_{mn}$, rather than directly tuning $p$.

We start from an initial condition in which all the nodes in the network are inactive, except for a single randomly chosen node, $n$ which is set to the excited state at time $t=0$: $x_{n}(0)=1$. The network is then allowed to evolve under the aforementioned dynamics. We refer to the sequence of resulting excited nodes as the avalanche or cascade starting at $n$. The duration, $\rho_{n}$, of an avalanche starting at node $n$ is defined to be the minimum number of time steps after which all nodes are in the resting state, i.e.,

[TABLE]

If no such minimum $t$ exists, the avalanche is said to last for an infinite duration. The size of an avalanche starting at node $n$ is defined as

[TABLE]

Note this is not the number of nodes involved in the avalanche, but rather the number of node activations, and a single node may be activated multiple times during one avalanche. The size of avalanches lasting for infinite duration is therefore infinite. Due to the probabilistic nature of activity propagation across edges, $\rho_{n}$ and $s_{n}$ are both random variables whose distributions characterize criticality in networks.

A cumulative distribution function of avalanche durations can be defined as

[TABLE]

which is the probability that an avalanche that begins at node $n$ has a duration that is less than or equal to $t$. Note that from the initial conditions in the problem, $c_{n}(0)=0$ for all $n$. Further, since $c_{n}(t)$ is necessarily a non-decreasing function that is bounded above by $1$, it must converge to some limit $b_{n}=\lim_{t\to\infty}c_{n}(t)$, which is the probability that an avalanche starting from node $n$ has a finite duration. In particular, for small values of $p$, all avalanches eventually subside, and hence have a finite duration, giving $b_{n}=1$. Such networks are called subcritical. For sufficiently large values of $p$, there is a positive fraction of avalanches that last for an infinite duration and hence $b_{n}<1$. Such networks are called supercritical networks. As we increase $p$ to go from a subcritical network to a supercritical network, there is a transition between the two phases, corresponding to the largest value of $p$ such that $b_{n}=1$. At this ‘critical’ transition probability, denoted by $p_{c}$, the network is said to be in the critical state, or equivalently, the network is said to exhibit criticality.

This phase transition is demonstrated in Fig. 2 for a variety of different networks, each with approximately $N=4\times 10^{4}$ nodes (for the algorithms used to construct these networks, and to demonstrate that these networks are sufficiently large to study criticality, see Appendix A). For each network we consider $10^{6}$ avalanches initialized at random nodes, and measure the fraction of avalanches that last for durations longer than $N/15$ timesteps. Since we consider the avalanches initialized in this random fashion, this fraction can be interpreted as the average value of $(1-b_{n})$ over all nodes, which we represent as $(1-b)$. Note how $(1-b)$ remains close to zero for low values of Perron-Frobenius eigenvalue, $\lambda$, until a critical phase transition, after which it rapidly increases. We denote the Perron-Frobenius eigenvalue of the critical network $A=p_{c}\times\widetilde{A}$ as $\lambda_{c}$. We are interested in estimating this quantity, and in Sec. IV we compare the estimates according a tree-based theory with our ‘second-level approximation’ to empirically determined values for $p_{c}$ and $\lambda_{c}$. Note that our work does not address the distributions of cascades away from criticality. While Ref. Melnik et al. (2011) discusses numerical results for network cascades both at criticality and away from criticality, our work will only address theoretical results for estimating $\lambda_{c}$, corresponding to networks at criticality.

We use the $\kappa$-metric developed in Ref.Shew et al. (2009) to characterize whether or not a given network is displaying criticality (For completeness, Appendix B gives the definition of $\kappa$). Since criticality is often characterized by power-law distributions of the avalanche sizes and the avalanche durationsShew et al. (2009); Beggs and Plenz (2003), this $\kappa$-metric was introduced to characterize the deviations of a given distribution of avalanche sizes away from the power-law fit to the data. $\kappa=1$ indicates a network at criticality, corresponding to an observed distribution of avalanche sizes close to a power-law distribution. $\kappa<1$ indicates subcritical network avalanches, with the distribution of avalanche sizes exponentially decaying at large duration values, and $\kappa>1$ indicates supercritical network avalanches, with a measurable fraction of avalanches that are infinitely long. Examples of such distributions are shown in Fig. 5, where we present log-log plots of the Complementary Cumulative Distribution Functions for cascade sizes in each of the three regimes (i.e., the probability that a cascade starting with a random initial perturbation has a size less than or equal to $s$, as a function of $s$) in the subcritical, critical and supercritical regimes.

II.1 Locally tree-like approximation

We say that a network satisfies the locally tree-like approximation if, for any two nodes $m$ and $n$ that are separated by a path with a small number of edges $l_{1}\ll N$, the probability that there exists another path from $m$ to $n$ of length $l_{2}\leq l_{1}$ is negligible.

Under this approximation we assume that events occurring along different edges that lead away from the same node are independent from each other, since the approximation implies that the corresponding activation paths do not share common nodes. Using this approximation, the following recursion relation for $c_{n}(t)$ holds:

[TABLE]

This recursion relation was used by Larremore et al.Larremore et al. (2012) to demonstrate that the network exhibits criticality when the Perron-Frobenius eigenvalue of the matrix $A$ is one, i.e.,

[TABLE]

where $\widetilde{\lambda}$ is the Perron-Frobenius eigenvalue of the adjacency matrix $\widetilde{A}$. This can also be seen in Fig. 2, where the phase transition to criticality for locally tree-like networks, such as a large random Erdős-Rényi network (shown in red squares) appears to occur at $\lambda=1$.

II.2 Reasons for success and conditions for failure of the tree-based theory

In order to derive the results in Eq. (8), the main assumption made about the network structure is the locally tree-like approximation, which is manifested in Eq. (7). In order to calculate $c_{n}(t)$, Eq. (7), the probability distributions $c_{m}(t-1)$ are assumed to be independent distributions for each node $m$. This assumption is equivalent to assuming that the cascades of excitation propagating from different edges starting from the same node do not interact. Let us assume that this were not the case, and that cascades resulting from different edge excitations do interact. Note that under the dynamics that have been defined, two excitation cascades starting at a node $n$ can only interact at time $k$ if some node $m$ is being excited by activation of two different edges simultaneously. This would mean that there exist two $k$-length paths between the nodes $n$ and $m$. We shall call the motif generated by the two distinct paths of length $k$ from the same initial node to the same final node as shown in Fig. 1 to be $k$-parallel. (This nomenclature has been chosen to be similar to earlier nomenclature of the bi-parallel motif, such as in Ref.Milo et al. (2002), which is then equivalent to what we term the $2$-parallel motif.) Thus we see that in the described dynamics, the motifs that are primarily responsible for deviations from the tree-based theory are $k$-parallel motifs, the smallest such motif being the $2$-parallel motif.

For the case of epidemiological models involving SIR-like dynamics (for a description of SIR dynamics see Refs.Volz (2008); Miller (2011); May and Lloyd (2001)), i.e., dynamics in which each node is activated exactly once and is then removed from the network, the smallest motif that allows interaction between cascades is the feed-forward triangular motif. Since our consideration is restricted to SIS-like dynamics, i.e., dynamics in which each node may be activated repeatedly, such triangular motifs are not relevant to the breakdown of the tree-based theory.

Thus, we note that while SIR and SIS dynamics are identical on directed networks near the critical point under the locally tree-like approximation, they can be distinguished from each other once this approximation is invalidated. (For recent work analyzing SIR-like dynamics on networks beyond the locally tree-like approximation see Ref.Radicchi and Castellano (2016).) For additional discussion regarding SIR-like dynamics on networks see Sec. VI

Note that while these statements are strictly true in our model described in Sec. II, there are some caveats that must be taken into account to make a similar statement applying to real-world scenarios. The model that we have described is a discrete state and discrete time model. For our model to be a reasonable representation of a real-world continuous time model, we are effectively making the assumption that the time-scale of propagation of node activity through each edge is approximately the same for all edges, and similarly the time-scale for which a node remains active before relaxing back into an inactive state is also approximately the same for all nodes. Further, in our discrete time model we have assumed that nodal state updates occur approximately synchronously for all nodes in the network.

As discussed above, triangular motifs are not directly relevant to the deviations from the tree-based theory for the dynamics we study, which is consistent with the observations Melnik et al.Melnik et al. (2011), where the authors note that clustering coefficients are not highly relevant to the deviations from the tree-based theory. To test the hypothesis that $k$-parallel motifs are essential as opposed to triangular motifs, for deviations from the locally tree-like approximation, we construct networks by randomly choosing triples of nodes that are connected as triangular feed-forward motifs. These networks have a low number of $k$-parallel motifs. We then set $p$ such that $\lambda=1$ and numerically simulate $10^{7}$ avalanches with random initially activated nodes. We obtain a value of $\kappa=1-1.3\times 10^{-3}\approx 1$, indicating that the prediction according to the locally tree-like approximation, i.e., $\lambda=\lambda_{c}=1$, accurately predicts criticality in the network. This can also be seen in Fig. 2, where the phase transition to criticality for such a network (shown in blue triangular markers) appears to occur at $\lambda=1$. We observe a similar behavior for triplets of nodes connected in cyclic motifs instead of feed-forward motifs.

Importantly, we also note that the presence of a large number of $k$-parallel motifs does not guarantee large deviations from the tree-based theory. To demonstrate this, we use the example of a network with all-to-all connections. By construction, all-to-all networks have the maximum possible number number of $k$-parallel motifs for each $k$ at each node. Numerical simulations of cascades on all-to-all networks with $N=4\times 10^{4}$ nodes with $p$ such that $\lambda=1$ yielded $\kappa=1-1.4\times 10^{-3}\approx 1$, indicating criticality and hence demonstrating no significant deviation from the prediction according to the locally tree-like approximation. This can also be seen in Fig. 2, where the phase transition to criticality for an all-to-all network (shown in green circular markers) appears to occur at $\lambda=1$. To see why this is the case, observe that in an all-to-all network, the number of paths of length $k$ between two given nodes scales as $N^{k-1}$, whereas the total number of paths of length $k$ starting at a given node scales as $N^{k}$. Thus the probability for two excitation cascades to meet at any given node scales as $N\times(N^{k-1}/N^{k})^{2}\sim 1/N$ and hence tends to zero as $N$ goes to $\infty$.

In general, the deviations from a tree-based theory are suppressed by large average degree in the network for similar reasons, since it is increasingly unlikely for cascades beginning at a given node to interact at a later time. This is demonstrated in Fig. 3, which shows a comparison between our theoretical results for $\lambda_{c}$ (shown in the black solid curve) with numerical estimates for $\lambda_{c}$ using the $\kappa$ metric ($\kappa$ shown in color; the numerical estimate for $\lambda_{c}$ corresponds to the green (light gray) region representing $\kappa\approx 1$) for networks with varying average degree with a constant density of 2-parallel motifs. Note that both the numerical estimate as well as our theoretical result rapidly approach 1 with increasing average degree. Particularly, average degrees as small as 10 are sufficient for the deviations away from a tree-based theory to be negligible. We further quantify this intuition in Sec. V.2.2. In Fig. 2 we see that all-to-all networks (which have the maximum possible density of 2-parallel motifs with the maximum possible average degree) continue to show a phase transition at approximately $\lambda=1$, whereas a high density of 2-parallel motifs in a network with a low average degree of 3.3 (magenta diamond markers) shows a phase transition at a significantly larger value of $\lambda$ (to demonstrate that the networks we studied are sufficiently large to study criticality, see Appendix A).

In addition, we note that the effect of interaction between two excitation paths is always to reduce the number of active nodes at any given time with respect to the number of active nodes expected according to a tree-based theory. This is because interaction effectively makes the transmission of an excitation through an edge immaterial if the corresponding node is excited by another node. The probability of an avalanche to last for any given duration is less due to the presence of such an interaction, and the addition of interaction effectively suppresses the overall number of transmitted excitations. Hence, a tree-based theory always under-predicts the value $p_{c}$ with respect to the actual critical transmission probability, i.e., for networks that may not satisfy the locally tree-like approximation, $(\widetilde{\lambda})^{-1}\leq p_{c}$ or $\lambda_{c}\geq 1$. We verify this numerically in Sec. IV.

Previous workMelnik et al. (2011) demonstrated that shorter mean inter-vertex length correlates well with smaller deviations from the results of a tree-based theory (in particular, the difference between the mean inter-vertex length for the given network and the mean inter-vertex length for a random network with the same size and joint degree distribution as the given network). Our interpretation of this result is as follows: For networks with a low mean inter-vertex length, $l$, the number of nodes that can be reached from any given node in a short number of steps rises very rapidly with number of steps. Assume that in $\mathcal{O}(l)$ number of steps, the number of nodes reached is $\mathcal{O}(N)$. Thus if we consider two random paths of the same length $k$ from the same starting node, the probability that they both end at the same node is very small. In particular, if $k\sim\mathcal{O}(l)$, then the probability that a $k$-parallel motif might be relevant to the interaction of two excitation paths scales as $\mathcal{O}(1/N)$. Thus, a network with a small mean inter-vertex length would be expected to have a low density of $k$-parallel motifs for small $k$. Following our earlier discussion on the importance of $k$-parallel motifs to deviations away from the locally tree-like approximation, we note that our result is hence consistent with the observations in Ref. Melnik et al. (2011).

In any iteration, the probability of a given $k$-parallel motif being relevant to such interactions scales as $p^{2k}$, and hence only $k$-parallel motifs for small $k$ are relevant to the deviations from a tree-based theory (since $0<p<1$). Under the assumption of no double edges in the network, the smallest $k$-parallel motif is the $2$-parallel motif. Hence in what follows we neglect the effect of $k$-parallel motifs with $k>2$ and focus on the correction to an estimate for $p_{c}$ due to the presence of $2$-parallel motifs. We note that in a similar fashion to the discussion above, ignoring the interaction between excitation cascades due to $k$-parallel motifs for $k>2$ leads to an underestimate of the transmission probability $p$ for criticality.

III Corrections due to 2-parallel motifs

According to the tree-based theory, in Eq. (7), when calculating the cumulative distribution at a node, we assume that the cumulative distributions at time $(t-1)$ are all independent of one another. We now consider what happens if we break this assumption, and assume instead that while distributions at time $(t-1)$ may not be independent of one another, distributions at time $(t-2)$ are independent of one another. We refer to results obtained assuming this condition as the ‘second-level approximation’. This approximation only takes into account the effect of $2$-parallel motifs. While we present the details and derivations of our second-level approximation in Sec. V, a short summary of our argument is as follows:

Note that the product in Eq. (7) explicitly assumes the locally tree-like approximation. This product can be rewritten as a summation (Eq.(10)) where an important quantity in the expression is the number of nodes reached from paths of lengths 1 and 2 starting at any node. In this summation form of the expression, the locally tree-like approximation is now encoded in writing this number of nodes reached in terms of the degree distributions of the network. 2. 2.

For networks that explicitly do not satisfy the locally tree-like approximation, we write a variant of the summation form Eq. (7) where the relevant quantity in the expression is again the number of nodes reached (Eq. (17)) — however, when the locally tree-like approximation does not hold, the number of nodes reached does not depend only on the degree distribution, but also on the density of small motifs. As discussed earlier, the smallest relevant motif is the $2$-parallel motif, hence we calculate the average number of nodes reached by paths of length 1 and 2 in terms of the degree distribution as well as the density of $2$-parallel motifs. 3. 3.

As discussed in Sec. II, for critical cascades the probability that cascades will not persist indefinitely, $b$, is equal to one. Following the analysis in Ref. Larremore et al. (2012), it can be shown that at criticality $b=1$ is in fact a double root of Eq. (7) in the limit of $t\to\infty$ (In this limit note that $b=1$ is always a solution, however at criticality an additional root at $b=1$ appears, resulting in the double root). We use a similar double root criteria on Eq. (17) to obtain an expression relating the probabilities of activity transmission across edges in the network to quantities depending on the degree distribution and the density of $2$-parallel motifs (Eq. (22), and in Appendix G).

IV Numerical verification of results

Our second level prediction for $p_{c}$ is derived in Sec. V and given in Eq. (22). To test this prediction we generate networks designed to exhibit a large deviation from a tree-based theory by growing networks with a high density of $2$-parallel motifs. In Appendix C we describe the algorithms used to construct such networks. In particular, Appendix C.2 describes the construction of a network with a scale-free degree distribution, and Appendix C.1 describes the construction of a network with a sharper, approximately Poisson degree distribution. Fig. 4 shows an example of a small network generated via the algorithm described in Appendix C.1.

Unless otherwise specified, in all cases we consider network with approximately $5\times 10^{4}$ nodes with algorithm parameters chosen as described in the Appendices to result in an average degree of about $3.5$. We numerically simulate $10^{6}$ iterations of the branching process dynamics as described in Sec. II to obtain avalanche distributions for varying values of $\lambda$. As discussed earlier in Sec. II, critical avalanches are characterized by power-law statistics, and we use the $\kappa$ metric to evaluate whether the obtained avalanche distributions correspond to criticality. We use the Python package ‘powerlaw’ created by Alstott et al.Alstott et al. (2014), which uses the tools developed by Clauset et al.Clauset et al. (2009) and Klaus et al.Klaus et al. (2011) to calculate the value of $\kappa$ for the obtained distributions for each value of $\lambda$.

To demonstrate the shift in criticality due to the presence of 2-parallel motifs, we present representative distributions of avalanche sizes on one such generated network. We calculate the complementary cumulative distribution function (CCDF) for the distribution of avalanche sizes, i.e., the probability that an avalanche has a size greater than or equal to a given size $s$, as a function of $s$, $[1-\langle\Pr(s_{n}\leq s)\rangle_{n}]$, where $\langle\ldots\rangle_{n}$ denotes an average over $n$. We present log-log plots of this CCDF in Fig. 5, for three values of $p$, corresponding to Perron-Frobenius eigenvalues of $\lambda=1.0$ (the critical value of $\lambda$ predicted by the locally tree-like approximation), $\lambda=1.0094$ (the critical value predicted by our second level approximation Eq. (22)) and $1.019$. Consistent with our theoretical analysis, networks with $\lambda=1$ are seen to be subcritical, while our second level approximation, $\lambda=1.0094$, yields the closest correspondence to criticality, and $\lambda=1.019$ is supercritical. The $\kappa$ values corresponding to these three curves are $0.992$, $0.998$, and $1.009$, respectively.

We note that the $\lambda=1.0094$ curve is almost perfectly linear out to a size of $10^{4}$ and only noticeably begins to deviate away from linearity as the size nears $10^{5}$ where the curve starts to have discernible downward curvature. This indicates a slight degree of subcriticality (reflected by the value of $\kappa$, $\kappa=0.998<1$). This behavior is to be expected, since, as discussed subsequently, our “second level” theory only takes into account 2-parallel motifs, while neglecting the effect of $k$-parallel motifs (Fig. 1) for $k>1$, and thus is expected to slightly underestimate the critical value of $\lambda$. For reference, numerically stepping through values of $\lambda$ indicates that $\kappa=1$ at an eigenvalue of $\lambda\approx 1.011$.

To test the effect of 2-parallel motifs on deviations from a tree-based theory we start with a network with a relatively high density of 2-parallel motifs, and then generate a set of networks having the same joint in-out degree distribution by swapping the destination nodes of a collection of randomly chosen edge pairs. By increasing the number of swapped edge pairs, we generate networks that have a decreasing density of $2$-parallel motifs, ranging from a maximum in the initial network, to close to zero for a network in which a very large number of edges have been swapped. For each such network, we vary the excitation transmission probability across edges by varying the Perron-Frobenius eigenvalue $\lambda$ of the weighted adjacency matrices of the networks. At each eigenvalue, we numerically simulate $10^{6}$ cascades following the dynamics described in Sec. II and calculate $\kappa$ based on the distribution of avalanche sizes. Values of $\kappa$ closer to $1$ indicating criticality are represented by the green regions in Fig. 6. We compare this to the value of $\lambda_{c}$ as predicted by a tree-based theory (described in Sec. II.1, results proven in Sec. V.1) and the prediction made by the analysis using the second-level approximation equations (described in Sec. III and derived in Sec. V.2). We present the results of this comparison in Fig. 6, where we plot the results as a function of the density of 2-parallel motifs, $\sigma$ (defined in Eq. (23)). For a locally tree-like network, without the presence of any 2-parallel motifs, $\sigma=0$, and for an all-to-all network, with each node containing the maximum possible number of two parallel motifs, $\sigma=1$. We see that the prediction for $\lambda_{c}$ according to the locally tree-like approximation (dashed black line) is always less than or equal to the observed values of $\lambda_{c}$, shown in green. Further, we see that in all cases, the second-level approximation (solid black line) is significantly better than the tree-based theory at predicting the observed values of $\lambda_{c}$.

V Derivations

Our proof is structured according to the summary presented in Sec. III. In this section and in the Appendices, apart from the definitions made earlier in Sec. II, we also make use of the definitions presented below:

• •

$\mathcal{V}$ : Set of vertices in the network

• •

$\mathcal{E}$ : The set of the pairs of vertices corresponding to directed edges in the network; $(m,n)\in\mathcal{E}\iff\widetilde{A}_{mn}=1$

• •

$d(n)$ : Out degree of node $n$

• •

$v^{PF}$ : Perron-Frobenius eigenvector of $\widetilde{A}^{T}$

• •

$d_{\infty}$ : Average degree of nodes weighted by $v^{PF}$

• •

$n_{0}$ : Initial node that we consider

• •

$n_{j}$ : Variable that spans over nodes at a distance $j$ from node $n_{0}$

• •

$\mathcal{E}^{0}_{k}$ : A $k$ element subset of $\mathcal{E}$ such that the edges in the set start at node $n_{0}$

• •

$\mathcal{V}^{0}_{k}$ : Set of vertices that the edges in $\mathcal{E}^{0}_{k}$ point to; $|\mathcal{V}^{0}_{k}|=|\mathcal{E}^{0}_{k}|$

• •

$\mathcal{E}^{1}_{r}$ : An $r$ element subset of $\mathcal{E}$ such that the edges in the set start at a node in $\mathcal{V}^{0}_{k}$

• •

$\mathcal{V}^{1}_{r}$ : Set of vertices that the edges in $\mathcal{E}^{1}_{r}$ point to; $|\mathcal{V}^{1}_{r}|\leq|\mathcal{E}^{1}_{r}|$

• •

$\delta r$ : $|\mathcal{E}^{1}_{r}|-|\mathcal{V}^{1}_{r}|=r-|\mathcal{V}^{1}_{r}|$

• •

$D[\mathcal{V}^{0}_{k}]$ : Sum of degrees of nodes in set $\mathcal{V}^{0}_{k}$

• •

$Q_{\alpha}(n)$ : The number of nodes that are reachable from node $n$ via exactly $\alpha$ edge-independent paths of length $2$

• •

$\overline{Q_{\alpha}}$ : Average of $Q_{\alpha}(n)$ over all nodes weighted by $v^{PF}$

• •

$\sigma$ : A single parameter of the network defined to count the normalized effect of all $\overline{Q_{\alpha}}$s; defined in Eq. (23)

V.1 Estimating the critical transition probability using the locally tree-like approximation

To derive an explicit expression for the critical transition probability, $p_{c}$, we restrict ourselves to networks that are near criticality, i.e., networks for which $p$ is approximately the same as $p_{c}$. This condition implies that $b_{n}$ (the probability that an avalanche starting at node $n$ is finite) and correspondingly, $c_{n}(t)$ (the probability that an avalanche has duration less than $t$) for large $t$, are both close to $1$, and hence also close to each other.

Thus we ignore differences between $c_{n}(t)$ for various $n$ when compared with quantities comparable to $1$. We rewrite Eq. (7) as

[TABLE]

where $\mathcal{E}$ is the edge set of the adjacency matrix $\widetilde{A}$ (i.e., $(m,n)\in\mathcal{E}$ if and only if $\widetilde{A}_{mn}=1$) and hence $n_{1}$ spans over all nodes that have a directed edge from $n_{0}$ to $n_{1}$. More generally, we will use the index $n_{k}$ to span over all nodes that are at a distance of exactly $k$ steps away from $n_{0}$. Writing the out-degree of a node $n_{0}$ as $d(n_{0})$, we can approximate the above expression and rewrite it as

[TABLE]

The approximation is valid since for any set of arbitrary quantities $X_{i}$ that are close to each other

[TABLE]

when $m\times\delta X_{i}=m\times(X_{i}-\langle X\rangle)\ll\langle X\rangle$.

We can then Eq. (10) as a recursion relation to write $0=c_{n_{0}}(t)-f(c_{n_{k}}(t-k))=0$ in terms of a known function $f$, for some $t\gg k\gg 1$. We then take the limit $t\to\infty$, and the limit $k\to\infty$. Note that since Eq. (10) is true for all times and at all nodes, the two limits can be interchanged. In taking the limit of $t\to\infty$, we have an equation of the form $0=F(\{b_{i}\})$. As discussed in Sec. II, for networks that are critical or subcritical, $b_{n}=1$ for each $n$, and for networks that are supercritical, $b_{n}<1$ for each $n$. It should be noted that $b_{n}=1$ for each $n$ is always be a solution for the equation $0=F(\{b_{i}\})$, since setting $c_{n}(t)=1$ for each $n$ is always a solution to Eq. (7). In general, from arguments given in Ref. Larremore et al. (2012), we note that there are always two solutions for $b_{n}$. One solution of $0=F({b_{i}})$ is at $b_{n}=1$. For supercritical networks, the other solution is less than 1, and is the value of $b_{n}$ to which $c_{n}(t)$ converges. For subcritical networks, the other solution is greater than 1, and hence $c_{n}(t)$ converges to the solution at $1$ in this case instead. Since we are interested in networks at criticality, we are interested in calculating the conditions for a degeneracy of the two roots of $b_{n}$ at $b_{n}=1$ at each $n$. We do this by looking at the subspace of $b_{n}=b$ for all $n$, and by solving for the vanishing derivative of $F$ with respect to $b$ at $b=1$, which results from the existence of the double root. We first take the derivative on Eq. (10), and then, after taking the appropriate limits, use the result to derive the condition on $p$ for the vanishing derivative for $F$ at $b=1$:

[TABLE]

where

[TABLE]

In the appropriate limits, Eq. (13) reduces to the condition for the vanishing derivative of $F(b)$, i.e., $0=\partial F/\partial b$. Since $c_{n}(t)\to b=1$ as $t\to\infty$, we set $\partial c_{m}(t)/\partial b=1$ for each $m$. Since the derivatives are being evaluated at $b=1$, in the limit of $t\to\infty$, each $\xi(t|j-1,j)$ evaluates to $1$, and can hence be ignored. We take the limit of $k\to\infty$ to give

[TABLE]

where $n_{0}$ is the node considered initially, and for each $k$, $n_{k}$ spans over all nodes reachable after traversing $k$ edges from $n_{0}$. The summand is now independent of the variables of summation and hence the entire series of summations is reduced to the total number of terms that are summed over. Thus,

[TABLE]

where $\langle d(n_{l})\rangle$ is the average degree of nodes connected to the node $n_{0}$ after traversing exactly $l$ edges (under this definition $\langle d(n_{0})\rangle=d(n_{0})$). We define $d_{\infty}=\lim_{l\to\infty}\langle d(n_{l})\rangle$. By our initial assumption of the network being strongly connected, $d(n)>0$ for each node $n$. Thus, $\langle d(n_{l})\rangle>0$ for each $l$, and hence we obtain

[TABLE]

It can be shown that $d_{\infty}$ is the average of the degrees at each node of the network when weighted by the component of $v^{PF}$ at that node, where $v^{PF}$ is the Perron-Frobenius eigenvector of $\widetilde{A}^{T}$. We use this idea of weighing quantities by $v^{PF}$ later when the assumption of the locally tree-like approximation is broken, since it arises naturally as a consequence of taking limits of $t\to\infty$ in the recursion relation, which in effect includes terms spanning the nodes of the network according to the paths that connect them to $n_{0}$. We can then show that this weighted average of the degrees, $d_{\infty}$ is just the Perron Frobenius eigenvalue of $\widetilde{A}$, giving the result shown earlier in Eq. (8). The details of the derivations can be found in Appendix D.

Thus from Eq. (14), this method of determining criticality under the tree-like approximation gives $p_{c}=(\widetilde{\lambda})^{-1}$ and hence at criticality the Perron-Frobenius eigenvalue for $A$, i.e., $\lambda=p_{c}\times\widetilde{\lambda}=1$, replicating the condition for criticality from Ref.Larremore et al. (2012).

V.2 Estimating the critical transmission accounting for corrections due to $2$-parallel motifs

In Sec. V.2.1 we first setup a recursion relation in analogy to Eq. (10) that takes into account the effect of $2$-parallel motifs by assuming independence of distributions at times $(t-2)$ (as opposed to independence at times $(t-1)$, as is assumed for the locally tree-like approximation). Then, in Sec. V.2.2, we use the same idea of evaluating derivatives to find the condition for a double root at $b=1$ to estimate $p_{c}$.

V.2.1 Recursion relation

We rewrite the recursion relation in Eq. (10) as

[TABLE]

where $\sum_{\mathcal{E}^{0}_{k}}$ denotes a sum over all possible sets $\mathcal{E}^{0}_{k}$ which are $k$-element subsets of the set of edges that begin at $n_{0}$; and $\mathcal{V}^{0}_{k}$ is the set of nodes to which the edges of the set $\mathcal{E}^{0}_{k}$ point.

We can interpret the terms in Eq. 15 as follows: the index $k$ counts the number of edges connected to the node $n_{0}$ that are activated due it and ranges from [math] to $d(n_{0})$; the set $\mathcal{E}^{0}_{k}$ is the $k$-element set of activated edges; $\mathcal{V}^{0}_{k}$ is the corresponding set of activated nodes, which also has $k$ elements, since we have assumed that the network has no double edges; $p^{k}$ is the probability that the $k$ edges were activated; $(1-p)^{d(n_{0})-k}$ is the probability that the remaining $(d(n_{0})-k)$ edges remained unactivated; and finally, the product term,

[TABLE]

, is the probability that after each of the $k$ activated nodes, all avalanches had a duration of less than or equal to $(t-1)$ units of time. This final product term from Eq. 15 can be rewritten using the same equation as a recursion relation to obtain an expression dependent on probabilities as a function of $(t-2)$. Rather than using the recursion relation to write the expression directly, we rewrite the product term in an equivalent form as

[TABLE]

where we define $D[\mathcal{V}^{0}_{k}]=\sum_{v\in\mathcal{V}^{0}_{k}}d(v)$ as the sum of the degrees of the nodes in the set $\mathcal{V}^{0}_{k}$, which is the total number of edges that begin from the $k$ nodes activated by $n_{0}$; as earlier, $\sum_{\mathcal{E}^{1}_{r}}$ is the sum over all possible sets $\mathcal{E}^{1}_{r}$ which are $r$ element subsets of the set of edges, $\mathcal{E}$, that begin anywhere in the set $\mathcal{V}^{0}_{k}$; and, $\mathcal{V}^{1}_{r}$ is the set of vertices to which the edges of the set $\mathcal{E}^{1}_{r}$ point.

Here, $r$ is the variable that counts the number of edges activated due to any of the $k$ activated nodes at the previous time step, and analogous to the earlier equation, $\mathcal{V}^{1}_{r}$ is the set of activated nodes due to the $r$ activated edges.

In the case of the locally tree-like approximation: all edges present at one edge away from the initial node $n_{0}$ are independent of each other, and hence the set $\mathcal{V}^{1}_{r}$ has exactly $r$ elements; the final product term then represents the probability that the avalanches beginning from these nodes have a duration of no longer than $(t-2)$; and then Eq. (16) is equivalent to the original recursion relation in Eq. (7).

In the more general case (where the locally tree-like approximation may not be valid), it is possible for some of the edges that are present one edge away from the initial node, $n_{0}$, to end at the same node. This is due to the presence of $2$-parallel motifs in the network, and in this case the set $\mathcal{V}^{1}_{r}$, which contains the nodes activated due to the $r$ activated edges, contains $\leq r$ elements. Thus, we have,

[TABLE]

with $|\mathcal{V}^{1}_{r}|\leq r$. We henceforth refer to this equation as the second-level approximation equation.

We now estimate the quantity $|\mathcal{V}^{1}_{r}|$, which can then be combined with the the recursion relation of Eq. (17 to estimate $p_{c}$ in a similar fashion to the technique used in Sec. II.1 to calculate $p_{c}$ from Eq. (10). Since $|\mathcal{V}^{1}_{r}|$ is necessarily less than or equal to $r$, we write it as $|\mathcal{V}^{1}_{r}|=r-\delta r$.

We define the coefficients $Q_{\alpha}(n_{0})$ at a given node $n_{0}$ to count the number of 2-parallel structures originating at the node $n_{0}$. In particular, $Q_{\alpha}(n_{0})$ is the number of nodes that are reachable from node $n_{0}$ via exactly $\alpha$ edge-independent paths of length $2$. This is a quantity that is dependent on the topology of the network which can be measured independent of the dynamics on the network. For a locally tree-like network, $Q_{1}=d(n_{0})\langle d(n_{1})\rangle$, where as earlier $\langle d(n_{1})\rangle$ is the average degree of nodes connected to the node $n_{0}$ after traversing exactly one edge, and $Q_{\alpha}=0$ for $\alpha>1$. In terms of these coefficients, we show in Appendix E, that by accounting for 2-parallel motifs we can approximate $\delta r$ as

[TABLE]

Equations (17) and (18) contain all the required information to treat the second-level approximation equation in a similar fashion to Eq. (10) and derive the conditions and equations for $p_{c}$.

While the sum over $\alpha$ in principle goes up to a maximum value of $\alpha=D[\mathcal{V}^{0}_{k}]$, for further simplification we can apply a cut-off on $\alpha$, by only considering terms for which $\alpha<\alpha_{thr}$. We apply this cut-off on both the $Q_{\alpha}(n_{0})$ coefficients, as well as in the expansion of $(1-r/D[\mathcal{V}^{0}_{k}])^{\alpha}$. The cutoff is justified because when averaging across nodes $Q_{\alpha}(n_{0})$ falls off very rapidly with $\alpha$, i.e., there are few nodes with significant values of $Q_{\alpha}(n_{0})$ for large $\alpha$, while most nodes only have nonzero values of $Q_{\alpha}(n_{0})$ for small values of $\alpha$. Further, $r/D[\mathcal{V}^{0}_{k}]$ can be assumed to be much smaller than $1$, since terms for larger $r$ are exponentially suppressed with a $p^{r}$ term in the second level recursion relation, and would correspond to a large fraction of the edges connected from a single node being activated simultaneously in the same time step. Hence we can apply a cut-off on the expansion of $(1-r/D[\mathcal{V}^{0}_{k}])^{\alpha}$ as well.

Note that $Q_{\alpha}(n_{0})$ is the number of nodes that are reachable from $n_{0}$ via exactly $\alpha$ edge-independent paths of length two. Since our algorithm for network generation (described in Sec. IV) involves adding new paths of length two between nodes that already are separated by a path of length two, many of our generated networks tend to have $Q_{\alpha}>0$ for larger values of $\alpha$ than might be normally expected in real-world networks. We find that for the networks that we have constructed, and for other networks we have observed, $Q_{\alpha}\approx 0$ for $\alpha>6$, and hence for our final results we use a cutoff of $\alpha_{thr}=6$. For the remainder of Sec. III however, we use the an approximation to $\delta r$ assuming a cut-off of $\alpha_{thr}=3$ as a representative cut-off to demonstrate the subsequent algebra. This gives us

[TABLE]

It should be noted that this cut-off can be made higher without any significant change to the method of analysis presented below. In general, for a cut-off at $\alpha_{thr}$, the resulting approximation for $\delta r$ is an $(\alpha_{thr}-2)$ degree polynomial in $D[\mathcal{V}^{0}_{k}]$. We present relevant expressions of our final results for larger values of the cut-off in Appendix G.

V.2.2 Finding the critical transmission probability

Treating the addition of $2$-parallel motifs to the network structure as a small change to the overall network dynamics around the new point of criticality in the network, we use the same method for estimating $p_{c}$ as discussed previously in Sec. II.1. The derivatives are evaluated of the second-level approximation equation, Eq. (17), after which appropriate limits are taken.

[TABLE]

Once again, since the derivatives are evaluated at $b=1$, several terms evaluate to $1$ in the appropriate limits as performed earlier in Sec. II.1, yielding

[TABLE]

which results from using $|\mathcal{V}^{1}_{r}|=r-\delta r$, and observing that $\mathcal{E}^{1}_{r}$ is an $r$ element subset of a $D[\mathcal{V}^{0}_{k}]$ element superset of edges. Thus, there are $\binom{D[\mathcal{V}^{0}_{k}]}{r}$ such subsets, and the summation over $\mathcal{E}^{1}_{r}$ contains $\binom{D[\mathcal{V}^{0}_{k}]}{r}$ terms in the summation. We can now use the previously derived approximation for $\delta r$, Eq. (19), to simplify the above expression.

The term in the square brackets when summed over $r$ from [math] to $D[\mathcal{V}^{0}_{k}]$ is equivalent to the expectation value of the $3^{\text{rd}}$ degree polynomial in $r$, i.e., $r-\delta r$ (with $\delta r$ given by Eq. (19)) over a binomial distribution with probability $p$ over $D[\mathcal{V}^{0}_{k}]$ trials. This simplifies to give a $3^{\text{rd}}$ degree polynomial in $D[\mathcal{V}^{0}_{k}]$. Further, since the only term dependent on the set $\mathcal{E}^{0}_{k}$ is $D[\mathcal{V}^{0}_{k}]$, to simplify the sum over all such sets it suffices to evaluate this sum for powers of $D[\mathcal{V}^{0}_{k}]$. We demonstrate how to evaluate these sums in Appendix F. Using this, we simplify the final binomial summation over $k$ to obtain a polynomial equation for $p$ as

[TABLE]

Under the locally tree-like approximation, if we evaluate derivatives for only a single step in the recursion relation of Eq. (10), we obtain $1=p\times d(n_{0})$. Then, in the process of taking the appropriate limits, this form can be reduced to $1=p_{c}\times d_{\infty}$. We analogously posit that if we look at the right-hand side of Eq. (21) and evaluate each term not at $n_{0}$ or $n_{1}$, but rather in terms of an average over $v^{PF}$ as earlier, then the resulting equation allows us to determine $p_{c}$ by taking appropriate limits in Eq. (20). This gives the following equation to be solved for $p_{c}$

[TABLE]

where $\langle d^{q}\rangle$ is the average of the $q^{\text{th}}$ power of the degrees in the network when weighted by $v^{PF}$, and $\overline{Q_{q}}$ is the quantity $Q_{q}(n_{0})$ averaged over all nodes, weighted by $v^{PF}$. For the case of a network satisfying the locally tree-like approximation, i.e., when there are no $2$-parallel motifs in the network, $\overline{Q_{q}}=0$ for all $q>1$, and hence the above equation reduces to $p_{c}=1/d_{\infty}$, which is the same result as obtained in Eq. (14). Similar to the expression derived for the case of $\alpha_{thr}=3$, we can derive polynomial expressions in $p_{c}$ for higher values of $\alpha_{thr}$. For the results presented in Sec. IV we have used a cut-off of $\alpha_{thr}=6$. The polynomial expression solved in this case has been shown in Appendix G. Solving these high-degree polynomial expressions generally results in multiple roots for $p_{c}$. However, we find numerically that only one root satisfies $0<p_{c}<1$, which is the solution that we use.

We can also use the coefficients $\overline{Q_{\alpha}}$ to construct a single network parameter that we can use to compare the number of $2$-parallel motifs in different networks. In Eq. (18) if we set $r$ to its maximum possible value of $D[\mathcal{V}^{0}_{k}]$, and set $D[\mathcal{V}^{0}_{k}]$ to its maximum possible value of $d(n_{0})\langle d(n_{1})\rangle$, then $\delta r=\sum_{\alpha}[Q_{\alpha}(n_{0})\times(\alpha-1)]$. For this $\delta r$, the quantity $r-\delta r$ represents the total number of unique nodes that can be reached by traversing two edges. To normalize $\delta r$, we consider the ratio $\delta r/r=\sum_{\alpha}[Q_{\alpha}(n_{0})\times(\alpha-1)]/d(n_{0})\langle d(n_{1})\rangle$. Accordingly, we define the parameter $\sigma$, which we use in Sec. IV to compare the networks of different sizes, as

[TABLE]

Note that $0\leq\sigma\leq 1$, and for a network with all-to-all connections, $\sigma=1$.

We now go back to the earlier claim that for graphs with large average degree, a tree-based theory is valid, independent of the number of $2$-parallel motifs in the network (see Fig. 3). Let us assume that the nodes in $\mathcal{V}^{0}_{k}$ have some large average degree $d\gg 1$. Let us consider the case of the highest density of 2-parallel motifs where we assume that each of the nodes in $\mathcal{V}^{0}_{k}$ connect to the same $d$ nodes. If we consider the set $\mathcal{V}^{1}_{r}$ (set of vertices reached from an $r$ element subset of edges starting at $\mathcal{V}^{0}_{k}$), the probability of a certain value of the cardinality $|\mathcal{V}^{1}_{r}|$ over the ensemble of all possible activated $r$ edges is given by

[TABLE]

which is obtained by considering the number of ways to first choose the $|\mathcal{V}^{1}_{r}|$ out of the $d$ vertices to get activated, and then count the number of ways to be able to write $r$ as the sum of $|\mathcal{V}^{1}_{r}|$ positive integers. We then divide by the total number of possibilities, which are the number of ways to write $r$ as the sum of $d$ non-negative integers Under this distribution, the expected value of $|\mathcal{V}^{1}_{r}|$ is $rd/(r+d-1)$, which in the limit $d\to\infty$ gives $|\mathcal{V}^{1}_{r}|\to r$. Thus $\delta r\to 0$ and hence networks with large average degree can be treated directly under the locally tree-like approximation.

We further examine this heuristic argument in the context of Eq. (22). Since $\widetilde{\lambda}=d_{\infty}$ (as derived in Appendix D) we write $\lambda_{c}=p_{c}\widetilde{\lambda}=p_{c}d_{\infty}$ and note that $\lambda_{c}\sim\mathcal{O}(1)$ to imply that $p_{c}\sim\mathcal{O}(d_{\infty}^{-1})$. Thus, if $\mathcal{O}(\overline{Q_{\alpha}})<\mathcal{O}(\overline{Q_{2}})$ for $\alpha>2$, we see from Eq. (22) that

[TABLE]

We first consider the above equation under the approximation of a tightly-peaked degree distributions, such that $\langle d^{q}\rangle\sim\mathcal{O}(d_{\infty}^{q})$. The expression in the square brackets is dominated by the leading term of $1/d_{\infty}^{4}$, and hence

[TABLE]

Thus in the limit of large degrees $\lambda_{c}=1$ and $p_{c}=1/d_{\infty}=(\widetilde{\lambda})^{-1}$. However, when the fourth power of the average degree is large compared to 1, deviations from a tree-based theory are suppressed.

What happens in the case that networks do not satisfy the an approximation of a tightly-peaked degree distribution? In this case, where $\langle d^{q}\rangle\sim\mathcal{O}(d_{\infty}^{q})$ is no longer valid, the magnitude of the terms in the terms in the square brackets cannot be compared in a straight-forward fashion and it is unclear which of the positive or negative terms dominates over the others. In fact, for the case of scale-free degree distributions, averages over higher powers of degrees may diverge and not be well defined in the limit of infinitely large networks. However, as verified in Fig. 6(b), for large-but-finite network sizes our derived result Eq. (22) continues to hold, and only a small deviations from a tree-based theory are observed for a network with average degree as low as $3.6$. Further, as can be seen in Fig. 3, deviations away from a tree-based theory decay at a similar rate to the case for sharply peaked distributions. Thus despite a straight-forward comparison for increasing average degree not being possible for scale-free networks from Eq. (24), numerical experiments on large-but-finite networks suggest a similarly rapid decay of $\lambda_{c}$ towards 1.

In Fig. 6 the deviation away from a tree-based theory for a scale-free degree distribution is marginally greater than the deviation in the case of a Poisson degree distribution when compared at the same $\sigma$ with despite having a slightly larger average degree (3.6 for the scale-free network versus 3.3 for the poissonian network). See Fig. 8 for an explicit comparison. As noted above however, this is not a systematic comparison apparent from Eq. (24); In Appendix C.3 we present an example of a network with a degree distribution that is intermediate between a scale-free degree distribution and a tightly-peaked degree distribution which exhibits a larger deviation at the same $\sigma$ and average degree. Thus our results do not give a clear interpretation to the expected magnitude of deviation away from a tree-based theory dependent on the nature of the degree distribution. However, in each case we observe from Eq. (22) that the deviations from a tree-based theory appear to be of small magnitudes even at low average degree, and numerical results in each case demonstrate that our theory continues to give good results in all cases.

In practice (as demonstrated in Fig. 3), we see that for large networks with average degrees $\gtrsim 5$, often give deviations of less than $0.1\%$ between the empirically determined $\lambda_{c}$, and $\lambda_{c}=1$ as predicted by the tree-based theory. Since most networks that are encountered in physical systems tend to have large average degrees (for example almost all networks considered by Melnik et. alMelnik et al. (2011)), the locally tree-like approximation is usually a useful tool.

VI Discussions and Conclusions

In studying criticality of network cascades, tree-based approximations have been shown to often be accurate for a wide variety of real world networksMelnik et al. (2011). We develop analytic reasoning to justify the effectiveness of this approximation in predicting criticality for network dynamics that are a variant of discrete time SIS epidemiological models on a network, i.e., dynamics that allow for nodes to be re-excited an arbitrary number of times (model description in Sec.II).

Under the locally tree-like approximation, the condition for criticality of network avalanches is that the largest eigenvalue $\lambda$ of the probability weight matrix is one. Examining the assumptions made in deriving conditions of criticality according to a tree-based theory, we study the factors that contribute to deviations introduced when the locally tree-like approximation is no longer valid (Sec.II.2). We show that the network motifs contributing to these deviations are the $k$-parallel motifs (see Fig. 1), the smallest of which is the $2$-parallel motif (also known as the bi-parallel motif). In particular, for the simple discrete time SIS model that we study, the presence of a high density of triangular motifs does not imply large deviations away from a tree-based theory. In addition, large densities of 2-parallel motifs do not necessarily imply large deviations from a tree-based theory — large average degree can suppress these deviations independent of the density of 2-parallel motifs. For example, critical points for large all-to-all networks are well predicted by tree-based theories. We also explain how tree-based theories necessarily under-predict the required largest eigenvalue for criticality, i.e., $\lambda_{c}\geq 1$.

To study critical network cascades in networks that may not satisfy the locally tree-like approximation, we derive an expression for $\lambda_{c}$ on network topologies that allow for the presence of 2-parallel motifs (Sec.V). We obtain a polynomial equation for the transmission probability, $p_{c}$, corresponding to criticality in the network, which depends on the density of 2-parallel motifs, $\sigma$, in the network. Given $p_{c}$ we get $\lambda_{c}$ from the direct relationship: $\lambda_{c}=p_{c}\tilde{\lambda}$, where $\tilde{\lambda}$ is the Perron-Frobenius eigenvalue of the network’s adjacency matrix $\tilde{A}$. We use the obtained expressions to demonstrate that for networks with large average degree, the relevance of the $2$-parallel motifs is reduced, and results according to a locally tree-like approximation agree well with empirically determined quantities.

We then verify the derived equation, by performing numerical simulations of the described dynamics on a class of networks generated to have a large number of $2$-parallel motifs, and hence exhibit a critical transmission probability that is different as compared with the predictions of the locally tree-like approximation. We observe that the predictions made by taking into account the effect of $2$-parallel motifs agree well with empirical observations made on these generated networks, demonstrating the effectiveness of our derived results (Sec. IV, Fig. 6).

For networks with tightly-peaked degree distributions we show that deviations from a tree-based theory will be suppressed by a factor of the fourth power of the degree (see Eq. (25)). For networks with skewed degree distributions, although a straight-forward relationship with average degree is not apparent, application of our results to generated scale-free networks continues to show a similarly rapid decay of deviations from a tree-based theory with increasing degree as can be seen in Fig. 3. Thus, we note that when the average degree is even modestly large – more specifically, when the fourth power of the average degree is large compared to 1 – deviations away from a tree-based theory are strongly suppressed. Empirically, from Fig. 3, we observe that an average degree as small as 10 is sufficient to almost entirely supress any deviations. Since most real-world networks (including a large majority of networks considered in Ref. Melnik et al. (2011)) have large average degrees that are at least modestly large, we expect tree-based theories to provide sufficiently accurate descriptions of the network dynamics.

It is important to note that we have restricted our analysis to the limit of large network sizes. In general, small networks also promote the interaction between network cascades and contribute to deviations away from a tree-based theory. In particular, small networks also promote larger than expected values of $\lambda_{c}$. We expect that in general with increasing network sizes the predictions made via our ‘second-level approximation’ will agree with empirical estimates for criticality (as measured via the $\kappa$ metricShew et al. (2009)) with increasing accuracy. For skewed degree-distributions, in the limit of an infinite network size, the averages over higher powers of degrees may diverge. However, for large but finite network sizes we expect our derived expressions to be valid, and our broader conclusions of suppression of deviations from a tree-based theory with increasing average degree to continue to hold.

We note that our analysis assumes the system is near criticality. While empirical evidence suggests that tree-based approximations appear to hold away from criticalityMelnik et al. (2011), our analysis does not address this case.

Note that we assume the independence of distributions at two time steps before the current time, as opposed to independence at a single time step before the current time as is done in a tree-based theory. This idea may in principle be extended towards a greater number of time steps to observe higher-order effects relating to the deviations from a tree-based theory.

While we have focused on the case of SIS-like dynamics on networks (i.e., where nodes are re-excitable), we note that in the case of a directed network that does satisfy the locally tree-like approximation, SIS-like dynamics are identical to SIR-like dynamics (i.e., percolation problems, wherein after activation, nodes are removed from the network) near the critical point, since if excitation paths of short length do not interact, then the re-excitability of a node is immaterial to the system dynamics. In our analysis of SIS dynamics corresponding to Eq. (2), we argued that the 2-parallel motif was the smallest motif relevant to deviations away from a tree-based theory, and studied the correction introduced to Eq. (16) as a result of 2-parallel motifs in the network [Eqs. (17, 18)]. Analogously, we expect these ideas to apply to SIR-like dynamics as well, by including the correction to Eq. (16) as a result of feed-forward triangular motifs, which is the smallest relevant motif in this case (in a discrete time SIR model, if a node $x$ activates nodes $y$ and $z$, and node $y$ then attempts to activate $z$ again, $z$ will remain inactive, having been effectively removed from the network following the first time it was activated. Thus such motifs result in behavior that is not locally tree-like, and would manifest as fewer activated nodes as compared to a locally tree-like network). Based on the number of triangular motifs and 2-parallel motifs, SIR-like dynamics result in a different expression for $\delta r$ in Eq. (18), which can be used in the remainder of our analysis to construct an analogous expression for $\lambda_{c}$. We leave the details of this study for future research work.

We expect the results for $\lambda_{c}$ for SIR-like dynamics to similarly demonstrate that in the limit of large average degree, tree-based theories will result in close agreement to observed dynamics on networks that are far from trees. While we do not show this analysis here, we now consider a motivating example of a large-$N$, all-to-all network, and demonstrate that networks with a high density of triangular motifs will continue to behave tree-like in the limit of large average degrees. Starting from any node in the network, consider the ensemble of pairs of paths of length 2 starting at that node. The total number of such pairs of paths scales as $N^{4}$. We compare this with the number of pairs of paths that result in dynamics that are not locally tree-like, i.e., pairs of paths that meet each other at any point. There are $N^{3}$ pairs of paths that meet each other at the end points, corresponding to deviations as a result of 2-parallel motifs. This is relevant to both SIS-like dynamics (as discussed earlier in Sec. II.2) and SIR-like dynamics. However, for SIR-like dynamics, we also need to consider contributions due to feed-forward triangular motifs formed by these pairs of paths. Again, there are $N^{3}$ pairs of paths that result in such a motif. Since we are considering an all-to-all network, each of the pairs of paths will be equally likely, and thus the probability that excitation paths will interact, leading to dynamics that are not locally tree-like scales as $N^{3}/N^{4}\sim 1/N$, which tends to zero as $N$ goes to $\infty$. Thus, in a large all-to-all network the probability that the network will display deviations away from a tree-based theory is negligible for both SIS-like and SIR-like dynamics.

Thus, we expect the type of analysis we have developed here to not only apply to SIS-like dynamics, but also to SIR-like dynamics, and hence also to percolation problems on complex networks. Consequently, we expect our ideas and results to provide insight towards a wide range of problems in network science relating to critical dynamics.

VII Acknowledgments

This work was supported by the National Science Foundation through Award # 1632976 and by the Army Research Office under Grant No. W911NF-12-1-0101.

Appendix A Effects of network size

In our paper we generally perform our numerical simulations on networks with $N\lesssim N_{0}=5\times 10^{5}$ nodes. To demonstrate that this network size is sufficiently large, we present a series of phase transition curves similar to Fig. 2 for varying network sizes. In Fig. 7 we show a series of 4 curves for each type of network considered:

1. (a)

Networks with a high density of feed-forward triangular motifs, but a negligible density of $k$-parallel motifs. These networks were generated by initially considering $N_{0}$ nodes, then considering $N_{0}$ random triples which are then connected in a feed-forward fashion., and finally by taking the strongly connected component of the resulting network. The network obtained will have $N$ nodes for $N\lesssim N_{0}$ 2. (b)

All-to-all networks on $N_{0}$ nodes. 3. (c)

Networks with a high density of $2$-parallel motifs, generated using the algorithm shown in Appendix C.1 4. (d)

The strongly connected component of Erdős-Rényi networks corresponding to the networks in panel (c).

In each case, $N_{0}$ was varied from $5\times 10^{3}$ to $5\times 10^{4}$. Since only the strongly connected component of the network was always considered, the final network sizes, $N$, are slightly smaller than $N_{0}$ (except for the case of all-to-all networks, where $N=N_{0}$). The solid curves with markers represent this largest network size, and have also been shown in Fig. 2. Since it appears that increasing network sizes results in phase transition curves approaching the curve corresponding to $N_{0}=5\times 10^{4}$, we use this largest size for all of our numerical simulations.

Appendix B Definition of the $\kappa$ metric

The $\kappa$ metric was first introduced in Ref. Shew et al. (2009) in the context of cascades of neuronal excitation in the cortex. In this paper we use the python package powerlaw written by Alstott et alAlstott et al. (2014), which implements the following definition of $\kappa$

[TABLE]

where $F^{fit}$ is the cumulative distribution of the best-fit power-law distribution through the data (fit using the techniques of Refs. Clauset et al. (2009); Klaus et al. (2011)), $F^{obs}$ is the cumulative distribution of the observed data, and the $\langle\cdots\rangle_{s}$ denotes an average over logarithmically spaced bins over the data. This results in a sensitive metric for measuring deviations away from critical avalanchesShew et al. (2009), and is such that $\kappa<1$, $\kappa=1$ and $\kappa>1$ corresponds to subcritical, critical and supercritical network avalanches.

Appendix C Network construction algorithms to generate high densities of 2-parallel motifs

To construct networks with high densities of 2-parallel motifs, we devise algorithms that connect random nodes to triplets of connected nodes in a fashion that generates a 2-parallel motif. The precise algorithm used to make these random choice governs the degree distribution of the resultant network. Below we present three algorithms to generate networks with high density of 2-parallel motifs, the first having a Poisson degree distribution, the second having a scale-free degree distribution, and the third having a degree distribution intermediate between the two. The out-degree distributions for networks generated following each of the three algorithms are shown in Fig. 8(a). A comparison between the observed deviations from the locally-tree like approximation in each of the three networks is shown in Fig. 8(b).

For use in each of the algorithms below, we define the set $triplets$ such that $(x,y,z)\in triplets$ if and only if there is a directed edge from $x$ to $y$ and from $y$ to $z$, and the set $tnodes$ such that $x\in tnodes$ if and only if $(x,y,z)\in triplets$ for some node $y$ and some node $z$.

C.1 Poisson degree distribution

• •

Consider a set of $N_{0}$ nodes with no connections made initially between the nodes. Since there are no edges in the network yet, we initialize the sets $triplets$ and $tnodes$ as two empty sets.

• •

Randomly choose $M_{1}\ll N_{0}^{2}$ pairs of nodes and connect them via directed edges to generate an initial seed network. As each edge is connected update the set $triplets$ and $nodes$

• •

For $M_{2}$ iterations do the following:

• $\circ$

Choose a random node $w$ uniformly from the set of $N_{0}$ nodes

• $\circ$

Choose a random node $x$ uniformly from the set $tnodes$. For this node $x$, choose a random triplet $(x,y,z)$ from $triplets$ such that the first element of the triplet is $x$.

• $\circ$

Use the four nodes $x$, $y$, $z$ and $w$ to make a 2-parallel motif by making directed edges from $x$ to $w$ and from $w$ to $z$

• $\circ$

Update the set $triplets$ and then the set $tnodes$

• •

Finally, to ensure that the networks that we use to test our predictions are strongly connected, we then take the strongly-connected component of the network generated.

This results in a network with a high density of 2-parallel motifs with a degree distribution that is approximately Poisson. Unless otherwise specified, we choose $N_{0}=5\times 10^{4}$, $M_{1}=N_{0}$, and $M_{2}=N_{0}$. To vary the average degree of this network the parameter $M_{2}$ is varied.

C.2 Scale-free degree distribution

• •

Consider a set of $N_{0}$ nodes with no connections made initially between the nodes. Since there are no edges in the network yet, we initialize the sets $triplets$ and $tnodes$ as two empty sets.

• •

Consider the nodes $1,2,\ldots,M_{1}\ll N_{0}^{2}$. Choose $M_{2}$ pairs of nodes among these $M_{1}$ nodes and connect them via directed edges to generate an initial seed network. As each edge is connected update the set $triplets$.

• •

For each node $w$ among the remainder of the nodes $(M_{1}+1),(M_{1}+2),\ldots,N_{0}$ do the following $m$ times:

• $\circ$

Choose a random triplet $(x,y,z)$ uniformly from the set $triplets$.

• $\circ$

Use the four nodes $x$, $y$, $z$ and $w$ to make a 2-parallel motif by making directed edges from $x$ to $w$ and from $w$ to $z$

• $\circ$

Update the set $triplets$

• •

Finally, to ensure that the networks that we use to test our predictions are strongly connected, we then take the strongly-connected component of the network generated.

Since choosing the triplet uniformly from the set of triplets inherently biases the choice of nodes $x$ and $z$ to be proportional to their out-degree and in-degree respectively, thus the edges are created in a preferential attachment fashion, similar to the Barabási-Albert model used to generate scale-free networks. This algorithm results in a network with a high density of 2-parallel motifs with a degree distribution that is approximately scale free. Unless otherwise specified, we choose $N_{0}=5\times 10^{4}$, $M_{1}=N_{0}/5$, $M_{2}=2M_{1}$, and $m=2$. To vary the average degree of this network the parameter $m$ is varied.

C.3 Intermediate degree distribution

• •

Consider a set of $N_{0}$ nodes with no connections made initially between the nodes. Since there are no edges in the network yet, we initialize the sets $triplets$ and $tnodes$ as two empty sets.

• •

Randomly choose $M_{1}\ll N_{0}^{2}$ pairs of nodes and connect them via directed edges to generate an initial seed network. As each edge is connected update the set $triplets$ and $nodes$

• •

For $M_{2}$ iterations do the following with a $95\%$ probability:

• $\circ$

Choose a random node $w$ uniformly from the set of $N_{0}$ nodes

• $\circ$

Choose a random triplet $(x,y,z)$ uniformly from the set $triplets$.

• $\circ$

Use the four nodes $x$, $y$, $z$ and $w$ to make a 2-parallel motif by making directed edges from $x$ to $w$ and from $w$ to $z$

• $\circ$

Update the set $triplets$

• •

and with $5\%$ probability

• $\circ$

Choose 4 random nodes $x,y,z$ and $w$ uniformly from the set of $N_{0}$ nodes

• $\circ$

Use the four nodes $x$, $y$, $z$ and $w$ to make a 2-parallel motif by making directed edges from $x$ to $y$, $x$ to $w$, $y$ to $z$, and from $w$ to $z$

• $\circ$

Update the set $triplets$

• •

Finally, to ensure that the networks that we use to test our predictions are strongly connected, we then take the strongly-connected component of the network generated.

This results in a network with a high density of 2-parallel motifs with a degree distribution that is approximately intermediate between the sharply peaked Poisson distribution and the scale-free distribution for the above two algorithms (see Fig. 8(a) ). We choose $N_{0}=5\times 10^{4}$, $M_{1}=N_{0}$, and $M_{2}=N_{0}$.

Appendix D Proof that $d_{\infty}=\widetilde{\lambda}$

Let $v^{PF}$ be the Perron-Frobenius eigenvector of $\widetilde{A}^{T}$, and let $\widetilde{\lambda}$ be the corresponding eigenvalue. We first show that $d_{\infty}$ is the average of the degrees at each node of the network when weighted by the component of $v^{PF}$ at that node. Then we show that as a consequence of this $d_{\infty}$ is equal to $\widetilde{\lambda}$.

We define $e_{n}=[0\ldots 0\,1\,0\ldots 0]^{T}$ where the $1$ is at the $n^{\text{th}}$ position. The nodes reachable on traversing $k$ edges from the initial node $n_{0}$ are given by the vector $v_{k}=(\widetilde{A}^{T})^{k}e_{n_{0}}$, in which each entry of the vector is interpreted as the number of times the corresponding node is reached on traversing $k$ edges. Since $\widetilde{\lambda}$ is the Perron-Frobenius eigenvalue, and hence the eigenvalue with the largest magnitude, in the limit of $k\to\infty$ the vector $v_{k}$ asymptotically approaches $(\widetilde{\lambda})^{k}v^{PF}\propto v^{PF}$. Hence $d_{\infty}$ is just the average degree calculated using $v^{PF}$ as the set of weights.

Now, since the degree of node $i$ is given by $\sum_{j}{\widetilde{A}^{T}}_{ji}$, we can write

[TABLE]

Appendix E Approximation of $\delta r$ via the network coefficients $Q_{\alpha}$

Recall that $Q_{\alpha}(n_{0})$ is the number of nodes that are reachable from node $n_{0}$ via exactly $\alpha$ edge-independent paths of length $2$. We further define $Q_{\alpha}(n_{0}|\mathcal{V}^{0}_{k})$ to be the number of nodes that are reachable from node $n_{0}$ via exactly $\alpha$ edge-independent paths of length $2$, such that the first edge leads to a node in $\mathcal{V}^{0}_{k}$. Assuming that for each $\alpha$ the paths leading to the nodes being counted are distributed uniformly across all paths of length $2$,

[TABLE]

We use these $Q_{\alpha}(n_{0}|\mathcal{V}^{0}_{k})$ coefficients to estimate the quantity $\delta r$.

Recall that $r$ edges are activated out of a maximum possible of $D[\mathcal{V}^{0}_{k}]$ edges starting from nodes in $\mathcal{V}^{0}_{k}$. Since each of the $D[\mathcal{V}^{0}_{k}]$ edges are equally likely to be chosen, we assume the probability of choosing any one edge is given by $(r/D[\mathcal{V}^{0}_{k}])$. Now, for any given $\alpha$, let us consider a set of $\alpha$ edges that start from nodes in $\mathcal{V}^{0}_{k}$ and end at a common node. Out of the $r$ activated edges, if $\beta\geq 1$ edges are chosen to be activated from this set of $\alpha$ edges, then the contribution to $\delta r$ is $\beta-1$, since only the number of distinct activated nodes are relevant to $|\mathcal{V}^{1}_{r}|=r-\delta r$, and each of the $\alpha$ edges lead to the same node. For $\beta=0$, the contribution to $\delta r$ is [math]. We can find the expected contribution to $\delta r$ by summing over all possible $\beta$’s weighted with their respective probabilities as

[TABLE]

Since there are $Q_{\alpha}(n_{0}|\mathcal{V}^{0}_{k})$ such sets of $\alpha$ edges, we can add their contributions together, and sum over all possible $\alpha$ to write $\delta r$ as

[TABLE]

If we further approximate the binomial expressions in Eq. (28) to write

[TABLE]

we obtain the approximation described in Eq. (18) of the text. The motivation for this final step of approximation is to allow each term in the summation of the expression for $\delta r$ to be expanded to powers of $D[\mathcal{V}^{0}_{k}]$, which we later use in simplifying various terms.

Appendix F Expressions for sum of powers of $D[\mathcal{V}^{0}_{k}]$ over all possible sets $\mathcal{E}^{0}_{k}$

Recall that $D[\mathcal{V}^{0}_{k}]=\sum_{v\in\mathcal{V}^{0}_{k}}d(v)$, where $\mathcal{V}^{0}_{k}$ was a $k$ element subset of the nodes connected to $n_{0}$. We show that for arbitrary $q\in\mathbb{N}$,

[TABLE]

where $\langle d(n_{1})^{q}\rangle$ is to be interpreted as the average value of the $q^{\text{th}}$ power of the degrees of nodes reachable from $n_{0}$ after traversing $1$ edge, as opposed to $\langle d(n_{1})\rangle^{q}$, which is the $q^{\text{th}}$ power of $\langle d(n_{1})\rangle$, which is the average value of the degrees of nodes reachable from $n_{0}$ after traversing $1$ edge. We present a sketch of the argument for $q=1$ and $q=2$, which may be extended to arbitrary $q$.

Recall that $\mathcal{E}^{0}_{k}$ is the $k$ element set of activated edges that start at node $n_{0}$. Since we have assumed that the network has no double edges, hence each of these edges ends at a unique node, and we can treat the sum over all possible sets of edges $\mathcal{E}^{0}_{k}$ as a sum over all possible $k$ element sets of nodes $\mathcal{V}^{0}_{k}$ that are connected to $n_{0}$ via a single edge. For the case of $q=1$, we now have the summation

[TABLE]

There are a total of $\binom{d(n_{0})}{k}$ ways of choosing the $k$ element set of vertices connected to $n_{0}$ out of the maximum possible number of such vertices, $d(n_{0})$. Further, since the summation is performed over possible subsets $\mathcal{V}_{k}$, each vertex connected to $n_{0}$ appears the same number of times in the overall sum. Hence, we can replace the summand with the expected value of the summand over all possible sets, which is then just $k$ times the average degree of nodes connected to $n_{0}$. This gives

[TABLE]

For the case of $q=2$, we have

[TABLE]

The square of the sum on the right-hand side of the above expression results in terms that are either of the form $d(v_{i})^{2}$, or of the form $d(v_{i})d(v_{j})$. There are $k$ terms of the first type, with the degree of each node connected to $n_{0}$ being represented uniformly. Similarly, there are $k(k-1)$ nodes of the second form, which also appear uniformly for all nodes across the summation over all sets $\mathcal{V}_{k}$. As earlier, we can replace the summand with it’s expectation value before evaluating the summand to obtain

[TABLE]

We can use a similar reasoning to argue that for an arbitrary $q\in\mathbb{N}$ we have

[TABLE]

which then reduces to the final expression in Eq. (29).

Appendix G Polynomial expressions for $p_{c}$ for different values of $\alpha_{thr}$

The smallest value of the cut-off that leads to a result distinct from the locally tree-like approximation is $\alpha_{thr}=2$. For this value of $\alpha_{thr}$ we obtain

[TABLE]

where we follow the same notation as discussed earlier at the end of Sec. V.2.2 and write $p_{c}$ as $p$ for simplicity. At a cut-off of $\alpha_{thr}=3$ we obtain the result shown in Eq. (31)

[TABLE]

For $\alpha_{thr}=6$ we obtain

[TABLE]

Such expressions can be easily generated for large values of $\alpha_{thr}$ by analytically evaluating the required sums using Eqs. (20), (18) and (29) via computational tools such as Wolfram Mathematica(v 10.2, Student Edition), which we used to generate the above expressions. We note however that our requirement for a large value of $\alpha_{thr}$ stemmed from our algorithm used for network generation. In practice, one can easily measure $\overline{Q_{\alpha}}$ for large values of $\alpha$ for a given network and choose $\alpha_{thr}$ to be sufficiently large to capture most nonzero values of $\overline{Q_{\alpha}}$; we do not expect large values of $\alpha$ to be required for real-world networks.