Network processes on clique-networks with high average degree: the limited effect of higher-order structure

One of the problems that has motivated research in network science to a large extent is the assessment of how structural characteristics of real-world networks determine the performance of dynamical processes that take place on them [20]. Most analytical approaches to this problem use networks constructed via configuration models [5] as the substrate for the dynamics. In such models, one specifies the fraction of vertices with $k$ neighbors $p_{k}$. A sequence of vertex degrees $\{k_{1},...,k_{N}\}$ is then drawn independently following $p_{k}$. The network is then assembled by choosing pairs of “half-edges” (or stubs) uniformly at random from this sequence, which are joined to form complete edges. While this method is able to generate networks with any prescribed degree distribution along with offering great analytic tractability, it has the shortcoming that the generated networks are locally tree-like. That is, the density of cycles vanishes asymptotically as the network size increases. This contrasts markedly with the rich topological structure of real-world networks, which often exhibit short cycles, degree correlations and clustering (i.e., the tendency of groups of three vertices to form triangles). Clustering is common to a variety of systems, but it is specially important in social networks, where the average probability that two neighbors of a vertex are also neighbors themselves (also referred to as clustering coefficient) often reaches values of tens of percent [20]. Other classes of systems known to be highly clustered in this sense comprise biological and information networks [20]. Hence, the inclusion of triangles and other types of subgraphs in random network models appears to be a crucial step to model dynamical process on networks accurately.

A practical method to create analytically tractable random networks with a more realistic clustering structure is to extend the standard configuration in order to explicitly include the generation of motifs that yield clustering. The first model of this kind was proposed independently by Newman [21] and Miller [19]. This model sets two degree sequences drawn from a joint degree distribution: the first sequence prescribes how many edges each vertex is incident to, exactly as in the standard configuration model; and the second degree sequence defines the number of triangles to which each vertex is attached. As the model then matches these stubs accordingly into edges and triangles, it generates networks with non-vanishing clustering even in the limit of large sizes [21, 19]. This strategy can be adapted to produce networks not only with triangles, but also with distributions of cliques of larger size [8], different types of subgraphs [16], or edge-multiplicities [30].

A number of previous authors have investigated the impact of added clustering on several types of network dynamics by employing such extensions of the standard configuration model. For instance, using a model that created networks with arbitrary distributions of cliques [8], Gleeson et al. [9] showed that clustered networks exhibit higher bond percolation thresholds in comparison to locally tree-like structures with same degree distributions and correlation properties. Very recently, Mann et al. [18] studied the percolation properties of the model by Karrer and Newman [16] under different combinations of cycles and cliques as building blocks for the networks. The authors confirmed that the increased clustering created by cliques leads to higher percolation thresholds [18]. On the other hand, the dynamics of networks containing only cycles were shown to approach the result obtained for the configuration model when the length of these cycles increases, as the model then becomes more locally tree-like. A different method to add clique structures to standard configuration models is to use household models, where every vertex of the configuration model is exploded into a clique of a specified size [1]. In this model, clustering was found to increase the percolation threshold [2, 6]. However, when including other clustered subgraphs than cliques, the percolation threshold may either increase or decrease compared to a locally tree-like model [15, 28].

Network processes on configuration models with higher-order clustering find their widest application in mathematical epidemiology, because of the natural importance of modeling of outbreaks in real-world scenarios and the close analogy between disease spreading and percolation processes. Indeed, many results uncovered in the context of percolation have counterparts in disease spreading. For instance, the presence of triangles has been found to increase the epidemic threshold while decreasing the outbreak size  [19]. Likewise, networks composed of cycles have been shown to yield epidemic dynamics similar to those of tree-like networks as the length of these cycles increases [26]. Examples of other dynamics investigated with higher-order configuration models include cascade propagation  [13, 12], the Ising model [14], and synchronization of coupled oscillators [24].

In this paper we reveal an effect that seems to have remained unnoticed in previous works; namely, we show that the influence of higher-order subgraphs on network dynamics is negligible when the average degree is large. Specifically, we show that in such a limit the percolation dynamics of clustered networks for large outbreaks as well as the critical percolation value converge to the one expected for locally tree-like networks. We focus on the most clustered subgraphs possible: cliques of different sizes. While our analytical results are for the large average degree limit, our simulations show that this convergence kicks in for average degrees as small as 6 for several degree distributions. We also show that these conclusions hold for the synchronization transition of phase oscillators modelled by the Kuramoto model [27], indicating that the insensitivity to local network structures may hold for a wide range of network processes.

As a random graph model, we employ the random graph model with clustering developed in [21, 16]. This random graph model is a general framework that extends the configuration model to create networks with specified densities of arbitrary specified subgraphs. Including clustered subgraphs in the set of specified subgraphs enables to overcome the locally tree-like property of the standard configuration model. In this manuscript we focus on the most clustered sets of subgraphs, cliques. That is, every vertex has a joint clique degree vector $(s^{(1)},\dots,s^{(m)})$. Here $s^{(1)}_{i}$ denotes the edge-degree of vertex $i$, and $s^{(j)}_{i}$ denotes the clique-degree of size $j+1$ of vertex $i$. The clique-degree of a vertex describes a vertex’ involvement in cliques of a specified size. Thus, a vertex of clique degree $s^{(2)}_{i}=3$ is part of 3 cliques of size 3. We denote the probability that a vertex has clique-degrees $s^{(1)},\dots,s^{(k)}$ by $q_{s^{(1)},\dots,s^{(m)}}$. The degree or the total number of connections of vertex $i$ is then described by $\sum_{j=1}^{m}js^{(j)}_{i}$, because every clique of size $j+1$ adds $j$ connections to the vertex. We denote the degree distribution of a vertex by $p_{k}$, so that

$$p_{k}=\sum_{s^{(1)},\dots,s^{(k)}=1}^{\infty}q_{s^{(1)},\dots,s^{(m)}}\mathbbm{1}_{s^{(1)}+2s^{(2)}+\dots+ms^{(m)}=k}.$$

(1)

After sampling a joint clique-degree for every vertex, the network is then formed by selecting $j$ uniformly chosen clique-edges of size $j$, and pairing the corresponding vertices into a clique for all $j$ until all clique-edges have been paired into a clique. This is an extension of the standard configuration model, where the network is formed by pairing two uniformly chosen half-edges until all half-edges have been paired.

We now investigate the behavior of this network model under bond percolation, where every edge is removed independently with probability $1-\pi$. We first focus on the case where every vertex is part of only $k$-cliques. Let $q_{i}$ denote the probability that a randomly chosen vertex is part of $i$ $k$-cliques. Define the generating functions

$$g(x)=\sum_{i=1}^{\infty}q_{i}x^{i},\quad g_{p}(x)=\frac{1}{\langle s\rangle}\sum_{i=1}^{\infty}iq_{i}x^{i-1}=\frac{g^{\prime}(x)}{\langle s\rangle},$$

(2)

where $\langle s\rangle$ denotes the average number of $k$-cliques a vertex is part of. Let $u$ denote the probability that a randomly chosen clique-edge is not connected to the giant component. We are interested in the fraction of vertices in the largest component after percolation, $S$, which can be obtained by [21],

	$\displaystyle u$	$\displaystyle=g_{p}(\sum_{j=0}^{k-1}h(k,j,\pi)u^{j}),$
	$\displaystyle S$	$\displaystyle=1-g(\sum_{j=0}^{k-1}h(k,j,\pi)u^{j}),$		(3)

where $h(k,j,\pi)$ is the probability that a given vertex of a $k$-clique is still connected to $j$ other vertices of the clique after percolation with probability $\pi$. These implicit equations are in general difficult to solve [16, 17], so that it is difficult to make general observations on the solution of these equations. Therefore, we here focus on an approximation of $S$, first for large outbreaks ($\pi$ large), and then for small outbreaks (approximating the critical value where $S$ becomes larger than zero). In these approximations, we will assume that the number of connections of a vertex is large.

When the degree of a vertex is large, the probability that a randomly chosen clique-edge is not connected to the giant component becomes small. Thus, we expand (III) with a first-order Taylor expansion around $u=0$. This yields

$$u=g_{p}(h(k,0,\pi)+h(k,1,\pi)u).$$

(4)

Filling in the expressions $h(k,0,\pi)=(1-\pi)^{k-1}$ and $h(k,1,\pi)=(k-1)\pi(1-\pi)^{2(k-2)}$ yields

	$\displaystyle u$	$\displaystyle\approx g_{p}((1-\pi)^{k-1}+u\pi(1-\pi)^{2(k-2)})$
		$\displaystyle\approx g_{p}((1-\pi)^{k-1})$
		$\displaystyle\quad+g_{p}^{\prime}((1-\pi)^{k-1})(k-1)\pi(1-\pi)^{2(k-2)}u$		(5)

This results in

$$u\approx\frac{g_{p}((1-\pi)^{k-1})}{1-g_{p}^{\prime}((1-\pi)^{k-1})(k-1)\pi(1-\pi)^{2(k-2)}}.$$

(6)

Using a first order Taylor expansion, (III) then yields for $S$,

	$\displaystyle S$	$\displaystyle\approx 1-g(h(k,0,\pi)+h(k,1,\pi)u)$
		$\displaystyle\approx 1-g(h(k,0,\pi))+g^{\prime}(h(k,0,\pi))h(k,1,\pi)u$
		$\displaystyle=1-g((1-\pi)^{k-1})$
		$\displaystyle\quad-\frac{g_{p}((1-\pi)^{k-1})(k-1)(1-\pi)^{2(k-2)}\pi g^{\prime}((1-\pi)^{k-1})}{1-g_{p}^{\prime}((1-\pi)^{k-1})(k-1)\pi(1-\pi)^{2(k-2)}}$
		$\displaystyle=1-g((1-\pi)^{k-1})$
		$\displaystyle\quad-\frac{\langle s\rangle g_{p}((1-\pi)^{k-1})^{2}(k-1)(1-\pi)^{2(k-2)}\pi}{1-g_{p}^{\prime}((1-\pi)^{k-1})(k-1)\pi(1-\pi)^{2(k-2)}},$		(7)

where $\langle s\rangle$ again denotes the average number of cliques a vertex is part of.

Now $g((1-\pi)^{k-1})=g_{D}(1-\pi)$, where $g_{D}(x)=\sum_{k}p_{k}x^{k}$ is the generating function of the vertex degrees from (1). This means that for a given degree distribution $D$, the leading order term of the approximation of the largest component size does not depend on the clique size in which the vertex degrees are split. Furthermore, we show in Appendix C that the numerator of the second term also only depends on the degree distribution, not on the clique structure. Furthermore, $g_{p}^{\prime}((1-\pi)^{k-1})$ decreases when the network degrees increase. Thus, large outbreaks become asymptotically independent of the clique structures in the networks.

We now investigate networks where cliques of different sizes are present. By introducing different clique sizes, it is possible to create degree-degree correlations that have often been said to influence the largest component size after percolation [4, 11, 3]. Thus, we now investigate to what extent the introduction of mixed clique sizes influences the size of the largest component.

Under bond percolation of networks where every vertex is part of $s_{1}$ cliques of size $k_{1}$, and $s_{2}$ cliques of size $s_{2}$ with probability $q_{s_{1},s_{2}}$ the generating-function methodology gives the following results. Let $g(x,y)=\sum_{s_{1},s_{2}>0}p_{s_{1},s_{2}}x^{s_{1}}y^{s_{2}}$ be the generating function of the clique degrees. Furthermore, let

	$\displaystyle g_{p}(x,y)$	$\displaystyle=\frac{1}{\langle s_{1}\rangle}\sum_{s_{1},s_{2}>0}{s_{1}}q_{s_{1},s_{2}}x^{s_{1}-1}y^{s_{2}},$		(18)
	$\displaystyle g_{q}(x,y)$	$\displaystyle=\frac{1}{\langle s_{2}\rangle}\sum_{s_{1},s_{2}>0}{s_{2}}q_{s_{1},s_{2}}x^{s_{1}}y^{s_{2}-1},$		(19)

with

$$\langle s_{1}\rangle:=\sum_{s_{1},s_{2}>0}s_{1}q_{s_{1},s_{2}},\quad\langle s_{2}\rangle:=\sum_{s_{1},s_{2}>0}lq_{s_{1},s_{2}},$$

be the generating functions of the number of cliques that are reached by following a randomly chosen clique-edge. Let $u$ denote the probability that a randomly chosen $k_{1}$-clique-edge is not connected to the giant component. Similarly, let $v$ denote the probability that following a randomly chosen $k_{2}$-clique edge does not lead to the largest component.

We show in Appendix A that $u$ and $v$ can be approximated by

	$\displaystyle u$	$\displaystyle\approx\frac{g_{p}((1-\pi)^{k_{1}-1},(1-\pi)^{k_{2}-1})}{A(k_{1},k_{2},\pi)}$
	$\displaystyle v$	$\displaystyle\approx\frac{g_{q}((1-\pi)^{k_{1}-1},(1-\pi)^{k_{2}-1})}{A(k_{1},k_{2},\pi)}.$		(20)

Furthermore, the giant component size $S$ is then approximated as

	$\displaystyle S=$	$\displaystyle 1-g_{D}(1-\pi)-\pi(1-\pi)^{2(k_{1}-2)}(k_{1}-1)$
		$\displaystyle\times\frac{g_{p}((1-\pi)^{k_{1}-1},(1-\pi)^{k_{2}-1})^{2}\langle s_{1}\rangle}{A(k_{1},k_{2},\pi)}$
		$\displaystyle-\pi(1-\pi)^{2(k_{2}-2)}(k_{2}-1)$
		$\displaystyle\times\frac{g_{q}((1-\pi)^{k_{1}-1},(1-\pi)^{k_{2}-1})^{2}\langle s_{2}\rangle}{A(k_{1},k_{2},\pi)},$		(21)

where

	$\displaystyle A(k_{1},k_{2},\pi)=$
	$\displaystyle 1-(k_{1}-1)\pi(1-\pi)^{2(k_{1}-1)}\frac{\partial g_{p}((1-\pi)^{k_{1}-1},(1-\pi)^{k_{2}-1})}{\partial x}$
	$\displaystyle-(k_{2}-1)\pi(1-\pi)^{2(k_{2}-1)}\frac{\partial g_{q}((1-\pi)^{k_{1}-1},(1-\pi)^{k_{2}-1})}{\partial y},$

and where $g_{D}(x)$ is the generating function of the total vertex degrees. Thus, the leading order term of the giant component size does not depend on the distribution of the clique degrees $k_{1}$ and $k_{2}$, but only on the total vertex degree. Furthermore, the numerators of the second order terms also only depend on the degree distribution, and not on the clique sizes, similarly to the one clique-size case.

It is not difficult to extend this analysis to include more than two different clique sizes, where (IV) contains terms for all size biased generating functions of the clique sizes $g_{k_{i}}$, instead of only $g_{p}$ and $g_{q}$ in (IV). Therefore, even in the presence of multiple clique sizes that can generate degree-degree correlations, large outbreaks are clique-structure independent for large average degrees.

In this section we illustrate the limited effect of higher-order structure on more complex dynamic network processes than bond percolation. In particular, we focus on the dynamics of coupled oscillators. For this purpose, we employ the paradigmatic Kuramoto model [27] that can describe synchronization phenomena on complex networks. In the Kuramoto model, the oscillator of vertex $i$ is characterized by a phase variable $\theta_{i}$, and the dynamics on a heterogeneous network is dictated by the following equations [27]:

$$\frac{d\theta_{i}}{dt}=\omega_{i}+K\sum_{j=1}^{N}A_{ij}\sin(\theta_{j}-\theta_{i}),\;i=1,...,N,$$

(25)

where $\omega_{i}$ is the natural frequency of oscillation of oscillator $i$, and $A_{ij}$ is the network adjacency matrix. If there is an edge connecting $i$ and $j$, $A_{ij}=1$ (0 otherwise), and the interaction between the vertices is weighted by the coupling $K$. If $K$ is lower than a certain $K_{c}$, the oscillators rotate incoherently, each one at its own rhythm set by the natural frequency $\omega_{i}$. For $K>K_{c}$, the incoherent state loses stability: a cluster of oscillators is formed around an average phase value, and these units begin to rotate locked in the same frequency  [29, 7]. This transition from asynchrony to a partially synchronized state is measured by the Kuramoto order parameter $R$ given by [29, 7]

$$Re^{i\psi}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}}\quad(0\leq R\leq 1),$$

(26)

where $R$ quantifies the level of synchrony achieved by the oscillators, and $\psi$ is their average phase. While one can monitor the synchronization transition of a heterogeneous network with Eq. (26), it is not possible to decouple Eqs. (25) in terms of a global order parameter. Instead, in order to perform a self-consistent analysis and characterize the onset of synchronization analytically, we need to employ heterogeneous degree mean-field approximations  [27]. This is equivalent to replacing the terms of the adjacency matrix $A_{ij}$ by their ensemble averages in the configuration model, which in the single-edge version is $\langle A_{ij}\rangle=d_{i}d_{j}/N\langle d\rangle$. In the model that generates networks with a single clique type the expression is analogous, namely, $\langle A_{ij}^{(c)}\rangle=(k-1)c_{i}c_{j}/N\langle c\rangle$, where $c_{i}$ is the number of cliques attached to vertex $i$. Replacing $A_{ij}$ by $\langle A_{ij}^{(c)}\rangle$ in Eq. (25) we obtain

$$\frac{d\theta_{i}}{dt}=\omega_{i}+\frac{K(k-1)c_{i}}{N\langle q\rangle}\sum_{j=1}^{N}c_{j}\sin(\theta_{j}-\theta_{i}),\;i=1,...,N,$$

(27)

which motivates the definition of the following order parameter

$$re^{i\phi}=\frac{1}{N\langle c\rangle}\sum_{j=1}^{N}c_{j}e^{i\theta_{j}},$$

(28)

which in turn allows us to rewrite Eq. (27) as

$$\frac{d\theta_{i}}{dt}=\omega_{i}+Kr(k-1)c_{i}\sin(\phi-\theta_{i}).$$

(29)

In the limit of $N\rightarrow\infty$, we assume that the assignment of cliques and natural frequencies is well described by distributions $q_{c}$ and $g(\omega)$; we further assume that the collections of vertices with clique number $c$ and frequency $\omega$ form a phase density $\rho(\theta,t|c,\omega)$. Thus, we rewrite Eq. (28) in the continuum limit as

$$re^{i\phi}=\frac{1}{\langle c\rangle}\sum_{c}cq_{c}\int\int d\omega d\theta\rho(\theta,t|c,\omega)g(\omega)e^{i\theta}$$

(30)

By choosing $g(\omega)=(\sqrt{2\pi})^{-1}e^{-\omega^{2}/2}$, we can set $\phi=0$ without loss of generality. Substituting the stationary synchronous solution of Eq. (29)–i.e. $\rho(\theta|\omega)=\delta\{\theta-\arcsin[\omega/Kr(k-1)c]\}$, for $|\omega|\leq Kr(k-1)c$–into Eq. (30), we arrive at the following implicit equation

	$\displaystyle\frac{\langle c\rangle}{K}\sqrt{\frac{8}{\pi}}$	$\displaystyle=(k-1)\sum_{c}c^{2}q_{c}e^{-K^{2}(k-1)^{2}c^{2}r^{2}/4}$
	$\displaystyle\times$	$\displaystyle\left\{I_{0}\left[\frac{K^{2}(k-1)^{2}c^{2}r^{2}}{4}\right]+I_{1}\left[\frac{K^{2}(k-1)^{2}c^{2}r^{2}}{4}\right]\right\},$		(31)

where $I_{0}(\cdot)$ and $I_{1}(\cdot)$ are the modified Bessel functions of the first kind. Thus, with Eq. (31) we can find the dependence of the order parameter $r$ on the coupling strength $K$, and thereby assess the impact of different clique sizes on the onset of synchronization. Letting $r\rightarrow 0^{+}$, we also obtain the expression for the critical coupling

$$K_{c}=\frac{1}{(k-1)}\frac{\langle c\rangle}{\langle c^{2}\rangle}\sqrt{\frac{8}{\pi}}.$$

(32)

Thus, again, plugging in $d=(k-1)c$ shows that the critical coupling is independent of the clustering structure. However, an immediate problem we face is the fact that the mean field approximations behind Eq. (31) are accurate only for sufficiently dense networks, typically when the average degree is at least of order of a few dozen [27, 10]. This limits the analytical verification of the effect of cliques on the dynamics of networks as sparse as the ones considered in the previous sections. Nonetheless, in the appropriate regime in which the mean field approach is valid, Eq. (31) suggests that the conclusions drawn for bond percolation may be similar for synchronization processes: Notice that in Eq. (31) $(k-1)c$ is the actual degree of a vertex; substituting $d=(k-1)c$ in the implicit equation for $r$ and rewriting it in terms of the new variable, we find that the emergence of a synchronous component depends only on the final degree sequence of the network and not on the sizes of the cliques. Therefore, clustered and unclustered networks are expected to exhibit similar dynamics also in the synchronization of coupled oscillators.

In order to confirm the above result, let us first investigate how the critical point $K_{c}$ changes according to the clique structure. In Fig. 6 we compare the predictions of $K_{c}$ by Eq. (32) with the corresponding quantities obtained via numerical integration of the system (25) for several average degrees $\langle d\rangle$. We numerically detect the transition point between incoherence and partial synchronization by identifying $K_{c}$ as the position of the divergent peak of the susceptibility $\chi=N(\langle r^{2}\rangle_{t}-\langle r\rangle_{t})/\langle r\rangle_{t}$ [22], where $\langle\cdot\rangle_{t}$ is a long temporal average. As can been in Fig. 6, the agreement between simulation and theoretical values is satisfactorily good for low $\langle d\rangle$, but it is progressively improved as the networks get denser. Furthermore, Fig. 6 indicates that the transition to synchrony tends to occur sooner as the clique size, and hence the clustering, increases. The numerical value of $K_{c}$ for different clique sizes becomes statistically equivalent at high $\langle d\rangle$. Yet, the solutions obtained from Eq. (32) in Fig. 6 suggest that clustering always ameliorates the network synchrony, as seen in Fig. 6, an effect that asymptotically vanishes as $\langle d\rangle$ increases. This is in apparent contradiction with our analysis of Eq. (31), in that networks with the same degree distribution, regardless of their clustering structure, ought to have identical dependence $r=r(K)$ and critical couplings $K_{c}$. However, similarly to the experiments of Fig. 2, these networks of different clique sizes do not have the same degree distributions.

To verify whether the differences in the onset of synchronization shown in Fig. 6 are due to discrepancies in the degree sequences or are, on the other hand, a true effect of the clustered structure, we repeat the methodology of the experiments depicted in Fig. 3. That is, we simulate the oscillators on networks with different clique sizes but adjusted to have the exact same degree sequence. The result is seen Fig. 7. Again, as in the percolation experiment in Fig. 3, the synchronization values match almost perfectly despite the difference in the clustering levels: in the examples in Fig. 7, the single-edge networks ($K_{2}$) have transitivity coefficient [21] equal to $C\approx 0$, while networks with $K_{4}$ cliques exhibit $C\approx 0.18$–a significant structural difference that is not reflected in the dynamics. It is also noteworthy that this phenomenon occurs at a low average degree ($\langle d\rangle=6$), i.e., the regime depicted in Fig. 6 with the most prominent discrepancies between clustered and unclustered networks. As discussed for the percolation case in Sec. III, actually those discrepancies are due the fact that the degree distributions are not identical for different clique sizes; as a consequence, the mismatches in the degree sequence end up generating different solutions for Eq. (31). The results in Figs. 3 and  7 therefore show that networks with similar degree distributions and correlations may exhibit equivalent dynamical behavior regardless of their subgraph structure and clustering levels.

We also note that the critical coupling $K_{c}$ can be estimated via “quenched” mean-field approximations and, for the parameters considered here, expressed in terms of the largest eigenvalue $\lambda_{1}$ of the adjacency matrix as $K_{c}=\lambda_{1}^{-1}\sqrt{8/\pi}$ [25, 22]. We can complement the latter expression with the recent results of Ref. [23] in which the largest eigenvalue for Poisson random networks constructed with $K_{3}$ cliques has been estimated as $\lambda_{1}=2\langle c\rangle+1+1/\langle c\rangle$. In the limit of high average degrees, $\langle c\rangle\rightarrow\infty$, the third term of $\lambda_{1}$ vanishes, and the corresponding result for tree-like Poisson random networks is recovered ($\lambda_{1}=\langle d\rangle+1$). Therefore, also in the quenched mean-field formulation, the value of $K_{c}$ of clustered networks is expected to asymptotically approach the calculations for tree-like networks, in agreement with the results in Fig. 6.

In this paper, we have investigated the influence of the presence of clustered structures in random graphs in the form of cliques on two network processes: bond percolation and synchronization. Percolation on such clustered networks has been investigated frequently, but as the equations for the giant component size under percolation are given by several implicit equations that are difficult to analyze mathematically, the factors dominating the behavior of percolation processes on such networks are largely unknown. By approximating the size of the giant component under large outbreaks as well as the critical percolation value where a giant component starts to form, we have found that the degree distribution is the dominant factor in these approximations, especially when the average degree of the network is large. In particular, our approximations are independent of the amount of clustering in the network. This means that introducing clustering by locally inserting cliques or other types of subgraphs in the frequently used locally tree-like random graph models, barely influences the size of the largest component.

We also show that differences in percolation behavior due to the introduction of cliques in the configuration model that were found in previous works can be ascribed to the fact that the degree distribution changed in those experiments as well. When keeping the degree distribution fixed while introducing more clustering, this difference disappears.

While our approximations show that the dominant factor for large outbreaks as well as the critical percolation value is the degree distribution, and not the clustering in the network, our simulations show that actually the entire percolation curve seems to become independent of the clustering in the network once the average degree becomes large. Showing this analytically would be an interesting point for further research.

Furthermore, while we have primarily focused on the process of bond percolation, we also showed that for a different network process of oscillator synchronization, the same independence of higher-order structures is present when the average degree is large. We therefore believe that other processes such as opinion dynamics or the contact process could be independent of the clustering structure of this model as well. Investigating which types of dynamics are independent of the clique structures is therefore an interesting avenue for further research.

Another interesting line of research following from these results is in higher-order dynamics. We showed that singe-edge dynamics on networks where a clique structure is imposed behave similarly as in networks without the clique structure. However, when studying the network model for example as a simplicial complex instead, it is possible to impose simplicial dynamics on top of it, where the dynamics involve all clique vertices in the interactions. It would be interesting to see under which conditions on the dynamic process on such a simplical extension of this model depends on the clique structure, and under which conditions it does not.

Finally, while this work shows that inserting clustering in a locally tree-like model barely affects the behavior of an epidemic process under bond percolation, we believe that in different models, where clustering is introduced by the presence of geometry, bond percolation can behave very differently under two models of the same degree distribution. Showing general conditions on the network structure under which clustering does or does not affect the size of a giant component compared to tree-like network models would therefore also be an interesting avenue for further research.

T.P. acknowledges FAPESP (grant No. 2016/23827-6).