Heterogeneous micro-structure of percolation in sparse networks

TL;DR

This paper investigates the heterogeneous responses of individual nodes in sparse networks to edge removal, revealing variability in node connectivity and cluster sizes, with implications for understanding cascading failures and epidemic spread.

Contribution

It introduces a message-passing framework to analyze node heterogeneity in percolation, linking it to non-backtracking centrality and providing analytical results for large random graphs.

Findings

Node-to-node heterogeneity is significant near the percolation threshold.

Non-backtracking centrality captures heterogeneity in node responses.

Analytical approximations are provided for Erdős-Rényi graphs.

Abstract

We examine the heterogeneous responses of individual nodes in sparse networks to the random removal of a fraction of edges. Using the message-passing formulation of percolation, we discover considerable variation across the network in the probability of a particular node to remain part of the giant component, and in the expected size of small clusters containing that node. In the vicinity of the percolation threshold, weakly non-linear analysis reveals that node-to-node heterogeneity is captured by the recently introduced notion of non-backtracking centrality. We supplement these results for fixed finite networks by a population dynamics approach to analyse random graph models in the infinite system size limit, also providing closed-form approximations for the large mean degree limit of Erd\H{o}s-R\'enyi random graphs. Interpreted in terms of the application of percolation to real-world processes, our results shed light on the heterogeneous exposure of different nodes to cascading failures, epidemic spread, and information flow.

Equations

$$\phi (s)=\frac{1}{N}\sum _i\delta (s-⟨{\sigma }_i⟩).$$

$$\psi (n)=\frac{1}{N}\sum _i\delta (n-⟨n_i⟩).$$

$$H_{i\leftarrow j}=(1-\rho )+\rho \prod _{\ell \in {\mathcal{N}}_j\setminus i}{H}_{j\leftarrow \ell }.$$

$$S=\frac{1}{N}\sum _{i=1}^N\left(1-\prod _{j\in {\mathcal{N}}_i}{H}_{i\leftarrow j}\right).$$

$$⟨{\sigma }_i⟩=1-\prod _{j\in {\mathcal{N}}_i}{H}_{i\leftarrow j}.$$

$$⟨n_i⟩=1+\sum _{j\in {\mathcal{N}}_i}\frac{H_{i\leftarrow j}^{\prime }}{H_{i\leftarrow j}},$$

$$H_{i\leftarrow j}^{\prime }=\rho \left(1+\sum _{\ell \in {\mathcal{N}}_j\setminus i}\frac{H_{j\leftarrow \ell }^{\prime }}{H_{j\leftarrow \ell }}\right)\prod _{\ell \in {\mathcal{N}}_j\setminus i}{H}_{j\leftarrow \ell }.$$

$$1-\mathbf{H}=\epsilon \mathbf{a}+{\epsilon }^2\mathbf{b}+{\epsilon }^3\mathbf{c}+\mathcal{O}({\epsilon }^4),$$

$$\begin{split}\bm{a}&=\rho_{c}B\bm{a}\,,\\ \bm{b}&=\rho_{c}B\bm{b}-\frac{\rho_{c}}{2}\Big{(}(B\bm{a})^{2}-B\bm{a}^{2}\Big{)}+B\bm{a}\,,\\ \bm{c}&=\rho_{c}B\bm{c}-\rho_{c}\Big{(}B\bm{a}B\bm{b}-B(\bm{a}\bm{b})\Big{)}\\ &\hskip 36.98857pt+\frac{\rho_{c}}{6}\Big{(}(B\bm{a})^{3}-3B\bm{a}B\bm{a}^{2}+2B\bm{a}^{3}\Big{)}\\ &\hskip 36.98857pt+B\bm{b}-\frac{1}{2}\Big{(}(B\bm{a})^{2}-B\bm{a}^{2}\Big{)}\,.\end{split}$$

$$\mathbf{v}/{\rho }_c=B\mathbf{v},{\mathbf{u}}^T/{\rho }_c={\mathbf{u}}^{T}B,\Vert \mathbf{v}{\Vert }_1=1.$$

$$\mathbf{b}={\rho }_{c}B\mathbf{b}-{\alpha }^2\frac{{\mathbf{v}}^2}{2{\rho }_c}+{\alpha }^2\frac{{\rho }_c}{2}B{\mathbf{v}}^2+\alpha \frac{\mathbf{v}}{{\rho }_c}.$$

$$⟨{\sigma }_i⟩\approx \alpha (\rho -{\rho }_c)\sum _{j\in {\mathcal{N}}_i}{v}_{i\leftarrow j}.$$

$$(I-{\rho }_{c}B)\mathbf{b}=-\frac{{\mathbf{a}}^2}{2{\rho }_c}+\frac{{\rho }_{c}B{\mathbf{a}}^2}{2}+\frac{\mathbf{a}}{{\rho }_c},$$

$$\mathbf{b}=(I-{\rho }_{c}B{)}^{+}\left(\frac{{\rho }_{c}B{\mathbf{a}}^2}{2}-\frac{{\mathbf{a}}^2}{2{\rho }_c}+\frac{\mathbf{a}}{{\rho }_c}\right)+\beta \mathbf{v},$$

$$\begin{split}0&=-\bm{u}^{T}\rho_{c}\Big{(}\bm{a}B\bm{w}-\bm{a}\bm{w}\Big{)}+\frac{1}{6\rho_{c}}\bm{u}^{T}\bm{a}^{3}-\frac{\rho_{c}}{2}\bm{u}^{T}(\bm{a}B\bm{a}^{2})\\ &\quad+\frac{\rho_{c}}{3}\bm{u}^{T}\bm{a}^{3}+\bm{u}^{T}\bm{w}-\frac{1}{2\rho_{c}}\bm{u}^{T}\bm{a}^{2}+\frac{1}{2}\bm{u}^{T}\bm{a}^{2}\\ &\quad+\beta\left(\rho_{c}\bm{u}^{T}\bm{a}\bm{v}-\bm{u}^{T}\bm{a}\bm{v}+\bm{u}^{T}\bm{v}\right)\,.\end{split}$$

$$\begin{array}{rl}⟨{\sigma }_i⟩& \approx (\rho -{\rho }_c)\sum _{j\in {\mathcal{N}}_i}{a}_{i\leftarrow j}\\ & +(\rho -{\rho }_c{)}^2\sum _{j\in {\mathcal{N}}_i}\left(b_{i\leftarrow j}-\frac{1}{2}{a}_{i\leftarrow j}\sum _{\ell \in {\mathcal{N}}_i\setminus j}{a}_{i\leftarrow \ell }\right)\end{array}$$

$${\mathbf{H}}^{\prime }=\rho (I-\rho B{)}^{-1}1.$$

$${\epsilon }^{-1}\mathbf{x}+\mathbf{y}+\mathcal{O}(\epsilon )=({\rho }_c-\epsilon )(1+{\epsilon }^{-1}B\mathbf{x}+B\mathbf{y}+\mathcal{O}(\epsilon ))$$

$$\begin{array}{rl}\mathbf{x}& ={\rho }_{c}B\mathbf{x}\\ \mathbf{y}& ={\rho }_{c}B\mathbf{y}+{\rho }_{c}1-B\mathbf{x}.\end{array}$$

$${\xi }^{(-)}={\rho }_c^2\Vert \mathbf{u}{\Vert }_1.$$

$$\begin{split}&\varepsilon^{-1}\bm{x}+\bm{y}+\mathcal{O}(\varepsilon)\\ &=(\rho_{c}+\varepsilon)\Big{(}\bm{1}+B(\varepsilon^{-1}\bm{x}+\bm{y}+\mathcal{O}(\varepsilon))(\bm{1}+\varepsilon\bm{a}+\mathcal{O}(\varepsilon^{2}))\Big{)}\\ &\qquad\times\Big{(}\bm{1}-B(\varepsilon\bm{a}+\mathcal{O}(\varepsilon^{2}))\Big{)}\\ &=\varepsilon^{-1}\rho_{c}B\bm{x}+\rho_{c}\Big{(}\bm{1}+B(\bm{x}\bm{a})+B\bm{y}-B\bm{x}B\bm{a}\Big{)}\\ &\qquad+B\bm{x}+\mathcal{O}(\varepsilon)\,.\end{split}$$

$$\begin{array}{rl}& 0={\rho }_c\Vert \mathbf{u}{\Vert }_1+\alpha {\xi }^{(+)}{\mathbf{u}}^T{\mathbf{v}}^2-\frac{\alpha {\xi }^{(+)}}{{\rho }_c}{\mathbf{u}}^T{\mathbf{v}}^2+\frac{{\xi }^{(+)}}{{\rho }_c}\\ & \Rightarrow {\xi }^{(+)}={\xi }^{(-)}=\xi :={\rho }_c^2\Vert \mathbf{u}{\Vert }_1.\end{array}$$

$$⟨n_i⟩=\frac{{\rho }_c^2\Vert \mathbf{u}{\Vert }_1}{|\rho -{\rho }_c|}\sum _{j\in {\mathcal{N}}_i}{v}_{i\leftarrow j}+\mathcal{O}(1).$$

$$\frac{{\rho }_c^2\Vert \mathbf{u}{\Vert }_1{v}_{\min }}{|\rho -{\rho }_c|}{k}_i\lesssim ⟨n_i⟩\lesssim \frac{{\rho }_c^2\Vert \mathbf{u}{\Vert }_1{v}_{\max }}{|\rho -{\rho }_c|}{k}_i.$$

$$\hat{\phi }(h)=\sum _{k\ge 1}\frac{kp(k)}{c}\int \prod _{\ell =1}^{k-1}\text{d}\hat{\phi }(h_{\ell })\delta \left(h-1+\rho -\rho \prod _{\ell =1}^{k-1}{h}_{\ell }\right)$$

$$\phi (s)=\sum _{k\ge 0}p(k)\int \prod _{\ell =1}^k\text{d}\hat{\phi }(h_{\ell })\delta \left(s-1+\prod _{\ell =1}^k{h}_{\ell }\right).$$

$$\begin{array}{rl}{\phi }_k(s)& =\int \prod _{\ell =1}^k\text{d}\hat{\phi }(h_{\ell })\delta \left(s-1+\prod _{\ell =1}^k{h}_{\ell }\right).\end{array}$$

$$\begin{array}{rl}& \hat{\psi }(h,h^{\prime })=\\ & \sum _{k\ge 1}\frac{kp(k)}{c}\int \prod _{\ell =1}^{k-1}\text{d}\hat{\psi }(h_{\ell },h_{\ell }^{\prime })\delta \left(h-1+\rho -\rho \prod _{\ell =1}^{k-1}{h}_{\ell }\right)\\ & \delta \left(h^{\prime }-\rho \left(1+\sum _{\ell =1}^{k-1}\frac{h_{\ell }^{\prime }}{h_{\ell }}\right)\prod _{\ell =1}^{k-1}{h}_{\ell }\right),\end{array}$$

$$\begin{array}{rl}{\psi }_k(n)& =\int \prod _{\ell =1}^k\text{d}\hat{\psi }(h_{\ell },h_{\ell }^{\prime })\delta \left(n-\left(1+\sum _{\ell =1}^k\frac{h_{\ell }^{\prime }}{h_{\ell }}\right)\right)\end{array}$$

$$\quad h_{\star}=\begin{cases}1-\rho-\frac{1}{c}W\big{(}-c\rho e^{-c\rho}\big{)}\quad&\rho\geq\rho_{c}\\ 1&\rho<\rho_{c}\,,\end{cases}$$

Full text

I Introduction

Since the very begining of the modern fascination with networked systems, researchers have been interested in questions of propagation. Across many applications, the analysis of bond percolation provides a simple framework with which to analyze the capability of a network to transmit information, disease, influence, or failure Grassberger (1983); Callaway et al. (2000); Newman (2002); Watts (2002). Early work in this area mainly concentrated on understanding the global properties of percolation in the ensemble average of randomly generated model networks. Surprisingly, detailed results for single instances of fixed networks have only been available relatively recently Karrer et al. (2014); Hamilton and Pryadko (2014), and very little is known exactly about the responses of individual nodes Rogers (2015).

We consider bond percolation for fixed networks defined as follows: starting from an arbitrary large (connected) network, we evaluate each edge independently, keeping it with probability $\rho$ and deleting it with probability $1-\rho$. The largest connected component remaining after this random edge removal process is referred to as the percolating cluster or giant component; write $S$ for its size measured as a fraction of the total number $N$ of nodes in the network. For large sparse networks it was shown in Karrer and Newman (2010); Karrer et al. (2014) that this quantity can be computed to close approximation using a message-passing protocol.

Knowledge of $S$ gives global information about the robustness of a network to attack or infection; in particular, there is a critical value of $\rho$ below which no percolating cluster survives in the thermodynamic limit, and $S=0$. One of the main results of Karrer et al. (2014); Hamilton and Pryadko (2014) was to identify the percolation threshold $\rho_{c}$ as the reciprocal of the largest eigenvalue of the non-backtracking (or Hashimoto) matrix that encodes the relationship between variables in the message-passing equations. The eigenstructure of this matrix has received considerable attention recently, having been proposed as an efficient tool for both network clustering Krzakala et al. (2013) and centrality analysis Martin et al. (2014).

In this article we will be concerned with more detailed questions about the typical outcomes for individual nodes in the network, when averaged over many instances of the percolation process. For a given random instantiation of percolation, write $\sigma_{i}=1$ if node $i$ appears in the largest connected component, and $\sigma_{i}=0$ if not. Taking the ensemble average of this variable yields the probability $\langle\sigma_{i}\rangle$ for node $i$ to appear in the percolating cluster. Heterogeneity in the responses of individual nodes to percolation is captured by the empirical distribution of $\langle\sigma_{i}\rangle$, defined as

[TABLE]

Notice that the total fractional size of the percolating cluster is given by the mean of $\varphi$, that is, $S=\int s\,\varphi(s)\,\text{d}s$. When node $i$ does not appear in the percolating cluster, write $n_{i}$ for the size of the component it belongs to, and $\langle n_{i}\rangle$ for the average over many instances. The node average $\frac{1}{N}\sum_{i}\langle n_{i}\rangle$ was again analysed for finte networks in Karrer et al. (2014), and previously results for the distribution of finite cluster sizes in large random graphs was presented in Callaway et al. (2000). The empirical distribution small clusters is defined analogously to $\varphi$ in Eq. (1);

[TABLE]

It turns out that many networks exhibit extreme differences between nodes in both $\langle\sigma_{i}\rangle$ and $\langle n_{i}\rangle$, which are not well-represented by the average value of percolation probability or small cluster size. Figure 1 showsn an illustraive example for the probability to appear in the percolating cluster of nodes in a real-world dataset. This behaviour was previously observed in the particular context of epidemic spreading on networks, and exploited in Moreno et al. (2002) to formulate and analyse a heterogeneous dynamic mean-field approximation of epidemic spreading. Within that approximation the infection probability of nodes is postulated to depend only on the degrees of individual nodes, an assumption which allows of a self-consistent solution. This degree-based approximation has become a mainstay in the analysis of epidemic spreading on networks, including in particular also in the search for optimal vaccination strategies; see Pastor-Satorras et al. (2015) for a recent review. The heterogeneity was again observed within a cavity formulation of the problem Rogers (2015) which for the analysis of SIR dynamics allows an exact mapping on bond-percolation. In that work it was found that node degrees play a dominant role in the behaviour of $\langle\sigma_{i}\rangle$ near $\rho=1$, but also that the picture becomes much more complex near $\rho_{c}$.

In this paper we explore these issues in detail. After recapping the message passing formulation in the next section we move on in section III to consider node variability in fixed finite networks in the neighbourhood of the percolation transition. In this regime we apply a weakly non-linear analysis to compute the probability of a node to appear in the percolating cluster, and the expected size of non-percolating clusters containing that node. In fact, we will show that the measure of non-backtracking centrality proposed in Martin et al. (2014) determines the leading-order behaviour of both $\langle\sigma_{i}\rangle$ and $\langle n_{i}\rangle$ near percolation. This gives and interesting physical interpretation to a quantity that was originally devised for purely practical reasons. In section IV we move on to use a population dynamics approach compute the distributions of node percolation probability and expected small component size in the ensemble average for large random graphs with specified degree distributions. This analysis is more generally applicable to the whole range of $\rho$. Numerical simulations reveal the fine structure of the distributions $\varphi$ and $\psi$. We are also able to compute a closed-form approximation for the large mean degree limit of Erdős-Rényi random graphs.

II Message passing

As detailed previously in Karrer et al. (2014), analysis of the probability generating function of component sizes yields a set of self-consistency equations which can be solved efficiently by iteration. For a network with $M$ edges we define the $2M$-vector $\bm{H}$ to be the smallest solution in $[0,1]$ of the system

[TABLE]

Here we write $\mathcal{N}_{j}$ for the neighbourhood of node $j$, and the entries of the vector $\bm{H}$ are indexed by ordered pairs of nodes attached by an edge. The expected size of the percolating cluster is then given by

[TABLE]

Unpacking the sum above yields additional information about the likely outcomes of the percolation process for individual nodes, specifically

[TABLE]

Similarly, for the expected sizes of small clusters one obtains

[TABLE]

where $\bm{H^{\prime}}$ solves

[TABLE]

For fixed finite graphs, both (3) and (7) can be solved quickly by simple iteration, making analysis of percolation computationally tractable even in very large networks.

Notice that equation (3) admits a trivial fixed point in which $H_{i\leftarrow j}\equiv 1$, yielding $S=0$. This solution may or may not be stable, and the transition boundary exactly corresponds to the percolation threshold $\rho_{c}$. In the following section we examine in detail the behaviour of node-specific probability $\langle\sigma_{i}\rangle$ to belong to the giant cluster and the expected size $\langle n_{i}\rangle$ of finite clusters in the vicinity of the percolation phase transition.

III Weakly non-linear analysis

III.1 Percolation probability

The value of the percolation threshold for a fixed finte network is revealed by linear stability analysis of the message passing equations. For a given ordered pair $i\leftarrow j$, the pair $k\leftarrow\ell$ will appear on the right hand side of (3) if and only $j=k$ and $\ell\neq i$. The matrix $B$ encoding this relationship is exactly the non-backtracking matrix mentioned in the introduction. The instability of the $H_{i\leftarrow j}\equiv 1$ solution is thus seen to occur at $\rho_{c}$ satisfying $\rho_{c}\lambda_{\rm max}(B)=1$.

Here we will consider $\rho=\rho_{c}+\varepsilon$ for small positive $\varepsilon$, and postulate a Taylor expansion for the $2M$-vector of messages

[TABLE]

where $\bm{a}$, $\bm{b}$ and $\bm{c}$ are constant vectors, and $\bm{1}$ is the vector of ones, i.e. $\bm{1}=(1,1,\dots,1)^{T}$. Inserting into both sides of (3) and matching powers of $\varepsilon$ to third order we obtain the equations

[TABLE]

Here we use the notational shorthand of applying multiplication and exponentiation element-wise so that, for example, $\bm{a}\bm{b}$ denotes the vector with entries $a_{i\leftarrow j}b_{i\leftarrow j}$ and $\bm{a}^{2}$ that with entries $a_{i\leftarrow j}^{2}$.

The first equation in (9) simply states that $\bm{a}$ is an eigenvector of $B$ associated to the Frobenius eigenvalue $\lambda_{\max}(B)=1/\rho_{c}$. To obtain the proper normalisation for $\bm{a}$ it is necessary to check the solvability of the second order equation. Let us write $\bm{a}=\alpha\bm{v}$, where $\bm{v}$ is the right Frobenius eigenvector of $B$ with positive entries, summing to one. We also introduce the corresponding left eigenvector $\bm{u}$, again with positive entries. That is,

[TABLE]

Following standard conventions by normalising left and right eigenvectors to form a bi-orthonormal system and, having fixed $\|\bm{v}\|=1$, we need to choose $\|\bm{u}\|_{1}$ such as to achieve $\bm{u}^{T}\bm{v}=1$. Then, replacing $\bm{a}$ by $\alpha\bm{v}$, the second order equation from (9) reads as:

[TABLE]

Multiplying through on the left by $\bm{u}^{T}2\rho_{c}/\alpha$ cancels some terms, yeilding $\alpha=2/\bm{u}^{T}\bm{v}^{2}(1-\rho_{c})$. Returning to (5), we find the first order approximation

[TABLE]

In fact the right hand side here is exactly the so-called “non-backtracking centrality” that was proposed in Martin et al. (2014). If the original system was simply connected to begin with, our choice of normalisation for $\bm{v}$ implies that $0<v_{\rm min}\leq v_{i\leftarrow j}\leq v_{\rm max}<1$, entailing that close to the percolation threshold we have lower and upper bounds for $\langle\sigma_{i}\rangle$ which are proportional to degree and of the form $\alpha(\rho-\rho_{c})v_{\rm min}k_{i}\lesssim\langle\sigma_{i}\rangle\lesssim\alpha(\rho-\rho_{c})v_{\rm max}k_{i}$. If the system was not simply connected to begin with, we have $v_{\rm min}=0$, and the lower bound becomes trivial.

By the same methodology we are also able to obtain the curvature of $\langle\sigma_{i}\rangle$ near criticality. The $\bm{b}$ equation from (9) can be rewritten as

[TABLE]

and thus

[TABLE]

where $\,{}^{+}$ denotes the Moore-Penrose pseudoinverse and $\beta\in\mathbb{R}$ is a constant yet to be established. Introducing $\bm{w}=\bm{b}-\beta\bm{v}$, we multiply third order equation in (9) through on the left by $\bm{u}^{T}\rho_{c}$ to obtain

[TABLE]

The constant $\beta$ is easily obtained by rearranging.

Having solved for $\bm{a}$ and $\bm{b}$, a complete second order expansion for the probability of appearing in the percolating cluster is then given by

[TABLE]

Figure 2 shows some numerical examples.

III.2 Small clusters

Turning attention now to the expected size of finite clusters containing a particular node, we first note that below percolation we have $\bm{H}=\bm{1}$ and hence (7) admits the exact solution

[TABLE]

The matrix inverse here implies an order one pole at $\rho_{c}$, around which we develop a first order perturbation theory. If $\rho=\rho_{c}-\varepsilon$ and $\bm{H^{\prime}}=\varepsilon^{-1}\bm{x}+\bm{y}+\mathcal{O}(\varepsilon)$ then

[TABLE]

Separating orders yields

[TABLE]

Evidently $\bm{x}=\xi^{(-)}\bm{v}$ for some constant $\xi^{(-)}$. To determine the constant $\xi^{(-)}$, we consult the second order equation, as usual multiplying on the left by $\bm{u}^{T}$, to find

[TABLE]

To examine behaviour on the other side of the critical point, we set $\rho=\rho_{c}+\varepsilon$ and do the expansion again, this time using the result for $\alpha$ from the previous calculation.

[TABLE]

At leading order we again find a multiple of the Frobenius eigenvector, $\bm{x}=\xi^{(+)}\bm{v}$. The same trick of multiplying the second order term by the left eigenvector determines the constant:

[TABLE]

We conclude that near percolation, typical size of finite clusters involving node $i$ is symmetric around $\rho_{c}$, with the asymptotic form

[TABLE]

Note that this expression is again proportional to the non-backtracking centrality of node $i$. Figure 3 shows some numerical examples. Note also that Eq. (22) also allows to obtain an upper bound on the $\langle n_{i}\rangle$ in the vicinity of the percolation transition. As before, we can exploit the bounds $0<v_{\rm min}\leq v_{i\leftarrow j}\leq v_{\rm max}<1$ for systems which are origninally simply connected; they entail that in the vicinity of the percolation transition we have diverging lower and upper bounds for $\langle n_{i}\rangle$ of the form

[TABLE]

Upper and lower bounds for the $\langle n_{i}\rangle$ therefore proportional to the degree $k_{i}$ of the site in question, with the constant of proportionality diverging as the transition is approached. As for the $\langle\sigma_{i}\rangle$ the lower bound becomes trivial for systems which are not simply connected to begin with.

IV Population dynamics for $N\to\infty$

IV.1 Percolation probability

We now proceed to analyse the bond percolation problem in the thermodynamic limit $N\to\infty$ for random networks in the configuration model class. In this limit Eqs. (3) can be interpreted as a stochastic recursion for the values of $H_{i\leftarrow j}$ of randomly chosen neighbouring nodes $i$ and $j$. A probability density function $\widehat{\varphi}$ for the values of $H_{i\leftarrow j}$ will then have to satisfy a self-consistency equation of the form

[TABLE]

in which $p(k)$ is the degree distribution, and $\text{d}\widehat{\varphi}(h_{\ell})$ is shorthand for $\text{d}h_{\ell}\widehat{\varphi}(h_{\ell})$. It is obtained by averaging the r.h.s of (3) over all realisations for which $H_{i\leftarrow j}\in[h,h+\text{d}h]$. The factor of $kp(k)/c$ above is the branching distribution of the configuration model, giving the distribution of the number of additional edges incident upon node $j$ from the randomly chosen pair $i\leftarrow j$. Equation (24) simply states that if the $H_{j\leftarrow\ell}$ appearing on the right hand side of (3) are random variables with distribution $\widehat{\varphi}$, then across neighbourhoods $\mathcal{N}_{j}$ the resulting $H_{i\leftarrow j}$ should have the same distribution.

The distribution $\varphi$ of the probabilities to remain in the percolating cluster is expressed in terms of the solution $\widehat{\varphi}$ of Eq. (24) via

[TABLE]

Equation (24) does not offer much hope for exact solution, however, highly accurate numerical solutions can be obtained by iterating a large sample population $\{h_{\omega}\}_{\omega=1}^{\Omega}$, a technique known in this context as population dynamics Mézard and Parisi (2001). In Figure 4 we show the results of the population dynamics algorithm for random graphs with a power-law degree distributions of the form $p(k)\propto k^{-3}$, for $k\geq k_{\rm min}=2$, and a edge occupation probability of $\rho=1/2$. Also shown is a partial unfolding of $\varphi(s)$ according to degree $k$. The distributions show a significant amount of structure, including some sharply peaked ‘band-edges’ of a type that also appear in the Gnutella data shown in Fig. 1 for sufficiently large $\rho$.

More insight into the structure of the distribution $\varphi$ can be gained by disentangling the contributions coming from nodes of different degrees. We decompose (25) into $\varphi(s)=\sum_{k\geq 0}p(k)\varphi_{k}(s)$, where

[TABLE]

This unfolding according to degree is also shown in Figure 4. It reveals that – to the resolution shown – the peaked band-edges mentioned above are associated with degree 2 and degree 3 vertices, with the exception of the peak at $s=1$ which originates from $\varphi_{k\geq 10}(s)$, i.e. from large-degree vertices. Moreover, for each degree $k$, there is an upper cutoff $s_{k}<1$ of the probability, beyond which $\varphi_{k}(s)=0$. A bound for this cutoff is easily read off from the single instance equations. Indeed Eq. (3) implies that $H_{i\leftarrow j}\geq 1-\rho$, which via Eq. 5 in turn entails that $\langle\sigma_{i}\rangle\leq 1-(1-\rho)^{k_{i}}$. Hence we have the bounds $s_{k}\leq 1-(1-\rho)^{k}$ for the $k$-dependent cutoffs. Numerical evidence suggests that these bounds are rather sharp. Their exact form is very likely amenable to a more detailed probabilistic analysis.

IV.2 Small clusters

The population dynamics approach is also available for the distribution of expected finite cluster sizes in the thermodynamic limit. From Eqs. (3) and (7) we derive a recursion equation for the joint distribution $\widehat{\psi}(h,h^{\prime})$ of messages $H$ and $H^{\prime}$,

[TABLE]

in which we use an analogous shorthand for integration measures, $\text{d}\widehat{\psi}(h_{\ell},h^{\prime}_{\ell})\equiv\text{d}h_{\ell}\text{d}h^{\prime}_{\ell}\,\widehat{\psi}(h_{\ell},h^{\prime}_{\ell})$. This equation generalises (24), as $\widehat{\varphi}$ is recovered as the marginal $\widehat{\varphi}(h)=\int\text{d}h^{\prime}\,\widehat{\psi}(h,h^{\prime})$. The distribution of finite cluster sizes is then given by $\psi(n)=\sum_{k\geq 0}p(k)\psi_{k}(n)$, where

[TABLE]

Figure. 5 shows the results of numerical solution of the population dynamics equations for $\psi(n)$ in the case of a Poisson random graph of mean degree $c=4$ at bond occupation probability $\rho=0.3$, together with a deconvolution according to degree $k$.

The figure exhibits two $\delta$-peaks at $n=1$ and $n=1+\rho=1.3$, the first corresponding to isolated sites with $k_{i}=0$, which remain isolated after deleting bonds, the second giving the average cluster size of pairs of vertices, which formed isolated dimers in the original graph (before deleting bonds), and consequently appears in $\psi_{1}(n)$. There are in principle many more such $\delta$-peaks related to the contribution of originally isolated clusters to the $\langle n_{i}\rangle$-statistics. For the present system their weight is, however, too small to be detected in the continuum at the resolution chosen. The continuous part of $\psi(n)$ reveals several discernible peaks, which according to the deconvolution can be associated with vertices of degrees $k=1,2,3$, and 4, whereas the peaks of $\psi_{k}(n)$ for larger $k$ are too narrowly separated compared to their widths to give rise to discernible features in the overall distribution $\psi(n)$ at larger $n$. Each of the degree dependent $\psi_{k}(n)$ also exhibits an upper cutoff $n_{k}$ beyond which $\psi_{k}(n)=0$.

IV.3 Large mean degree

Further analytical progress can be made in the limit of large mean degree and small edge occupancy; here we develop a single defect approximation Biroli and Monasson (1999) for Poisson random graphs. In this approximation scheme we assume that $\widehat{\varphi}(h)$ is for large mean degree $c$ well approximated by a a Dirac delta distribution, $\widehat{\varphi}(h)=\delta(h-h_{\star})$. Inserting this ansatz into (24) and integrating over $h$ we obtain $h_{\star}=1-\rho+\rho e^{-c(1-h_{\star})}$, from which it follows that

[TABLE]

where $W$ is the Lambert W function, and the percolation transition occurs at $\rho_{c}=1/c$. To obtain non-trivial behaviour in the limit $c\to\infty$, it is therefore necessary to introduce the scaling $\rho=\varrho/c$ where $\varrho>1$. Then $h_{\star}=1-\eta/c$, where $\eta=\varrho+W(-\varrho e^{-\varrho})$.

Imagining a single “defect” node $i$ with degree $k_{i}$ attached to the otherwise homogeneous network, we find the probability of this node to appear in the percolating cluster to be

[TABLE]

Taking the Gaussian limit of the Poisson distribution at large mean degree $c$, we note that the distribution of the ratio $k_{i}/c$ of a randomly selected node $i$ is well-approximated by a normal random variable $x$ with mean one and variance $1/c$. Under these assumptions we change probability variables from $x$ to $s$, computing

[TABLE]

where $f(x)=1-e^{-\eta x}.$ The inverse and its derivative are

[TABLE]

Putting all this together we arrive at

[TABLE]

Although it is exact only in the limit of large $c$, this approximation holds remarkably well for smaller values; see Fig. 6 for an example with $c=8$. As $c\to\infty$, the distribution $\varphi$ approaches a Gaussian with mean $1-e^{-\eta}$ and variance $\eta^{2}e^{-2\eta}/c$.

Applying the same approach to $\psi$, the distribution of finite clusters, we obtain a simple Gaussian law

[TABLE]

where $\gamma=\varrho/(e^{\varrho+W(-\varrho e^{-\varrho})}-\varrho)$.

V Summary and Discussion

To summarise, we have taken a new look at the message passing approach Karrer and Newman (2010); Karrer et al. (2014) to bond percolation on complex networks, revealing a considerable degree of heterogeneity. We have seen that this approach allows one to determine the distribution of probabilities of individual nodes to belong to the percolating cluster, as well as the distribution of the average sizes of non-percolating clusters to which individual nodes may belong. We found both distributions to be typically broad, so that the average percolation probabilities and average sizes of finite clusters that are typically reported in analyses of percolation on random networks must be regarded as poor representations of the actual heterogeneity that is present in this problem. We have also obtained deconvolution of both distributions according to degree, in analogy to the deconvolution of sparse random matrix spectra advocated in Kühn (2008), and found that the component distributions $\varphi_{k}$ and $\psi_{k}$ are themselves non-degenerate (except at $k=0$) and indeed typically broad as well. A fairly detailed analysis of node-specific percolation probabilities $\langle\sigma_{i}\rangle$ and average cluster-sizes $\langle n_{i}\rangle$ was provided in the vicinity of the percolation transition, and formulated in terms of spectral properties of the Hashimoto matrix $B$, notably the right and left eigenvectors corresponding to its largest eigenvalue.

In the present paper we have used the message passing approach to analyse a specific finite large real world instance, percolation on the network of the Gnutella file sharing platform Leskovec and Krevl (2014), and presented methods to analyse the problem for random graphs in the configuration model class in the thermodynamic limit of infinite system size, $N\to\infty$. We have presented examples for the distribution of percolation probabilities on a scale free graph with power-law degree distribution, and for the distribution of finite cluster sizes in the case of an Erdős-Rényi network. These two examples can only scratch the surface of the variability of phenomena that might be observed. A few general trends may be noted though. If the original random network ensemble contains finite clusters to begin with, as is indeed the case for Erdős-Rényi networks of finite mean degree, then the resulting distribution $\psi(n)$ of the average finite cluster sizes will contain a family of $\delta$-peaks originating from the distribution of the resulting finite cluster sizes generated by bond removal, whereas the broad continuum is generated by clusters obtained as a result of disconnecting a finite connected set of vertices from the original percolating cluster. Only if the non-percolating fraction in the original graph is sufficiently large will the delta-peaks carry sufficient weight to be detectable in the population dynamics. This is in particular the case for low mean degree Erdős-Rényi graphs. Cluster size-distributions typically get broader as the percolation transition is approached from above or below; the same is true after deconvolution meaning that peaks contributed by individual $\psi_{k}$ will cease to be discernible in the sum. The same broad trends are observed in $\varphi(s)$ and it’s deconvolutions.

We have been able to obtain closed form analytic approximations for the distribution of percolation probabilities and mean cluster sizes in the large mean connectivity limit of Erdő-Renyi graphs, which produces excellent results already for fairly moderate values of $c$. The same methodology could be adapted to other degree distributions, and should be efficient whenever these distributions become narrow in the large mean degree limit.

Returning to the heterogeneity of percolation probabilities across a network, and the practical relevance of this phenomenon in the context of cascading failures, epidemic spreading, or probabilistic information spreading, we note that the considerable detail which our methods allow to unearth might be useful for instance in the design of optimal vaccination strategies that exploit information beyond degree. Indeed the fact that percolation probabilities conditioned on degrees are themselves broadly distributed is a clear indicator of the fact that the degree based approximation on which the majority of attempts to design optimal vaccination strategies has been based misses important information which could be exploited to improve upon such strategies. We believe that this point is worth investigating.

*Acknowledgements

*TR is supported by the Royal Society.