Empirical evidence of non-trivial feature distributions

In this section, we report different patterns of feature distributions in real-world complex networks. The chosen examples are arbitrarily selected based on our biases, but they are not, by no means, an exception. Indeed, when collecting real data to construct network structures, most of the times the nodes have some associated properties, or metadata, that individually characterize them beyond their degree.

The first example corresponds to a board interlock network, i.e., a bipartite system of directors and companies. Links between them exist whenever a director holds a seat on the corporate board of a company, and nothing prevents a director to sit on more than one boardroom. We regard the feature as the age of the members of the board. We see that, in this case, the average age of the directors is uncorrelated with respect to their degree, see Fig. 1a. Nevertheless, the spread of the distribution as a function of the degree is not constant. Hence, feature-enriched percolation in this example can help reveal the role played by directors of a certain age on the global connectivity structure of the system.

Fig. 1

Empirical evidence of different feature-degree patterns and sketch of different types of removal protocols in the percolation process.

a Correlation between the degree of a director of a company (number of companies in which she participates) and her age. Data from⁶⁷. b Wikipedia pages show a positive correlation between the hyperlinks present in an article (out-degree) and the number of unique editors that have modified such article. Data from an updated version of the ref. ⁶⁸. c Negative correlation between the number of followers and the number of fake news posted by a user. The subset of data taken from⁶⁹, comprising worldwide users during the first two weeks of April 2020. The black solid lines correspond to the mean value of the data. In the bottom row, we show how different interventions can affect the network. In the random case (d) nodes are chosen equiprobably to be deleted. In the degree case (e), the removed nodes are those with the largest number of connections. In the feature case (f), the selected elements are those with the largest value of the feature (the darkest blue). In this example, the values of the feature have been assigned manually for convenience and the feature vector F has been considered unidimensional.

The second example belongs to the context of crowdsourced creation of cultural content. We take the snapshot of the current version of Wikipedia, in which nodes are articles and edges are the hyperlinks among them. Moreover, for each node we keep track of the number of revisions that it has suffered since its creation, and how many unique users have participated in these edits. We see that the correlation between the degree of a Wikipedia page and its number of unique editors is positive and heterogeneously distributed, see Fig. 1b. With the framework of feature-enriched percolation, we are able to remove nodes with a certain degree-feature pattern, thus, for example, one can be interested to assess how the navigability across this knowledge corpus is modified under the removal of well-connected articles that do not show enough levels of collaborative edition, i.e., that have been written by too few users.

As a final example, we present a system that displays a negative correlation among degree and feature, see Fig. 1c. It corresponds to misinformation propagation in the online social platform Twitter. We take a sample of messages within a two-week time window during the COVID-19 pandemic and consider only messages that share an url in the text. The nodes are users, where the degree is their number of followers, and the feature is the number of fake news that they have posted. We see that the average number of fake news tends to decrease as the visibility of Twitter accounts is higher, varying broadly in this case as well. Feature-enriched percolation can be helpful, for example, to shed light on the problem of how to guarantee the information spreading while cutting off users that systematically share dysfunctional and harmful content.

The model

Let us take the ensemble of networks generated via the configurational model⁹, with a degree sequence obtained from a degree distribution p_k. We assume that each node is characterized by its degree k and by a set of features F = (F₁, F₂, …, F_M), with $$ M\in \mathbb{N}$$. The degree distribution is considered discrete and the feature distribution can be considered either discrete p_F or continuous p(F). Throughout the article, we shall assume that features are continuous, but reproducing the results for a discrete domain is straightforward. The joint probability is indicated by P(k, F). We define the occupation probability ϕ_k,F as the probability that a node with degree k and feature values F within the interval [F, F + dF] have not been removed from the original network. Our goal is to find the position of the critical point and the size of the giant connected component S as a function of the parameters of P(k, F) and ϕ_k,F³².

Let u be the average probability that a node with k connections and feature values F does not belong to the giant component via one of its neighbors. The probability that it does not belong to the giant component is u^k. Thus, the probability that the node belongs to the giant component due to a neighbor state is 1 − u^k and has to be multiplied by ϕ_k,F, which determines whether or not the node itself is present in the network. Averaging this quantity over degrees and features we obtain ∫dF∑_kϕ_k,FP(k, F)(1 − u^k), where ∫dF is an M-dimensional definite integration over the elements of the feature vector. This is identified as the average probability of finding a node in the giant cluster, or equivalently, the fraction of nodes in the giant component. Therefore

To solve Eq. (1) we need to obtain an expression for u. This can be achieved by writing a self-consistent equation that has two contributions. On the one hand, a neighbor might not be in the giant component because it has been deleted, which occurs with probability 1 − ϕ_k,F. On the other hand, if the neighbor has not been removed, it should not belong to the giant component via any of its other k_e neighbors. k_e is called the excess degree distribution, and it is equal to k_e = k − 1. This happens with probability $$ {\varphi }_{k,\mathbf{F}}{u}^{k_e}$$. Averaging over the distributions, we get

Notice that the relation between generating functions held in ordinary percolation is valid here as well, i.e., g₁(z) = ∂_zg₀(z)/〈k〉. It is important to also note that our generating functions, although including the M-dimensional integral in the definition, depend only on one variable since the feature enrichment does not add any new information in terms of connectivity. Therefore, our framework should not be taken equivalent, for instance, to the study of percolation in graphs with colored edges³³, in general multilayers³⁴, in networks with multi-type nodes³⁵ or in interdependent systems³⁶, where multivariable generating functions are common.

The occupation probability ϕ_k,F allows us to understand the role played by certain values of degree and/or features in the connectivity of the network. The classical percolation, where nodes are removed in a uniformly random fashion, is recovered by selecting a constant function ϕ_k,F = ϕ ∈ [0, 1]. The case of removing the most connected nodes is recovered by setting ϕ_k,F = θ( − (k − k₀)), being θ(⋅) the Heaviside step function and k₀ a threshold such that all nodes with degree larger than it are removed. Similarly, one can apply the same arguments in the feature space, and study the case in which all nodes with a feature larger than a threshold are deleted ϕ_k,F = θ( − (F − F₀)). These three examples are sketched in Fig. 1.

Applications

To illustrate and check the validity of the theory, we investigate several examples. For the sake of simplicity we focus on unidimensional feature vectors, i.e., F = F. First we address the case of independent degree and feature, i.e., P(k, F) = p_kP(F). We then move to consider joint distributions which are positively and negatively correlated. These latter cases leave the nature of the feature undetermined. In this section, though, we also address problems in which the features are related to the distance in a geometrical space and to dynamical processes evolving on top of the network.

Independent case

Let us consider a network with degree distribution and feature distribution

Fig. 2

Size of the giant component for the independent degree-feature case.

Solid lines are computed from the theory, the markers come from simulations. In (a), the control parameter a is related to the topology, while in (b) the control parameter α is related to the feature distribution, see Eq. (6). Vertical lines indicate the value of the critical point computed from the theory, see Eq. (8). System size is N = 3000 and each point is computed by averaging over 100 independent realizations. In (c), the size of the giant component with error bars as a function of the system size, together with the fits to the data.

We investigate the universality class of the percolation process for independent feature-degree distributions, in order to figure out whether the introduction of features modifies the critical properties of mean-field percolation. Here we understand mean-field as the behavior of classical percolation in regular lattices of dimension d ≥ d_c = 6. The critical properties of complex networks, which are infinite-dimensional entities, might not be the same as mean-field, if, for instance, the underlying degree distribution is power-law¹⁶ or the removal process is modified³⁷. Note that in feature-enriched percolation we will necessarily have two parameters, one controlling the topology and another the features. Hence, for the size of the giant component, we write $$ S\left(a\right)~{\left(a-a_c\right)}^{{\beta }_a}$$ and $$ S\left(\alpha \right)~{\left(\alpha -{\alpha }_c\right)}^{{\beta }_{\alpha }}$$, being β_a and β_α the corresponding critical exponents. (To be accurate, S(a) also depends on α and S(α) depends on a too, although when studying the critical behavior, they are taken as constants, hence we do not write them for the sake of clearness.) In Supplementary Note 3 we analytically show that β_a = β_α = 1, i.e., it takes the mean-field value. To find other critical exponents, we can employ the finite-size scaling hypothesis

Positively correlated case

Let us consider now the more interesting and realistic case of a joint distribution in which feature and degree are not separable. We take one of the simplest scale-free distributions that are positively correlated,

In practical terms, to assign the values of the degree and the feature in the simulations, we first construct the network from the configurational model with a degree sequence drawn from p_k = ∫dFP(k, F). Then, depending on the degree of each node, we draw a random variable from the conditioned distribution P(F∣k). Randomized versions of the system are possible either by randomly shuffling the feature of the nodes in the correlated network or by constructing a new network from p_k and assign features from P(F).

Considering again the removal of all nodes with feature value above a threshold F₀, the generating functions are

The theoretical predictions are compared with the simulations in Fig. 3a, finding an excellent agreement. Another curve is also displayed, corresponding to the randomized version of the correlated model Eq. (10), and that serves to single out the role of the degree-feature correlations in the behavior of the size of the giant component. It is immediate to observe that the critical point and the size of the giant component are smaller for the correlated case than for the randomized case. This holds true for all feature thresholds, being the separation among both the critical points and the S more accentuated as F₀ decreases. The rationale behind this is the following: the marginal feature probability P(F) is the same for both cases, so for a given threshold F₀ the same fraction of nodes is removed, although the chances of hitting a low-degree node are much larger in the randomized case than in the correlated case, due to, precisely, the degree-feature correlation. Since targeting hubs it is a very fast way of dismantling a network, it is natural to observe a smaller critical point and a smaller size of the giant component. A quantitative way to measure the deviations caused by degree-feature correlations with respect to the uncorrelated scenario is to compute the area between both curves,

Fig. 3

Percolation in correlated degree-feature networks.

Size of the giant connected component for (a) the positively correlated distribution Eq. (10) and for (b) the negatively correlated distribution Eq. (14), and their randomized counterparts. All nodes with feature F ≥ F₀ = 3, are removed. Points are the results from simulations and the solid line comes from the theory. Network sizes are N = 2¹⁶ and each point is averaged over 100 realizations. In the inset of the panels, we display Δ(F₀), corresponding to the area between the randomized and correlated curves and defined in Eq. (13), as a function of the feature threshold. The integral has been computed in the region α ∈ (1, 5].

Figure 3a also shows a striking peculiarity: the non-monotonicity of the giant component and its approach to the non-percolating regime S = 0 when α → 1. In Supplementary Note 5, we give an explanation for this behavior and show that there are no critical properties in the vicinity of α = 1. This approach to the non-percolating phase is a genuine effect of the degree-feature correlation that would be missed if the correlation is disregarded. In other words, the robustness of this positively correlated system under the attack to its most prominent nodes in the feature space would be overestimated if correlations are ignored. In Supplementary Note 6 (see Supplementary Movies 1 and 2 as well) we also address the critical properties of this positively correlated network model, where we show that it shares the same critical exponents as the mean-field percolation.

Negatively correlated case

Contrary to the last section, there are certain situations in which peripheral, low degree nodes might be the ones carrying the largest feature values. Again, we model this scenario with an ad hoc negatively correlated joint distribution, one of the simplest displaying power-law behavior³⁹:

The generating functions corresponding to Eq. (14), and for the same occupation probability as the other examples, are

The comparison between the theoretical predictions and the simulations for the negatively correlated distribution is given in Fig. 3b, both for the correlated distribution and its randomized version, finding again a very good agreement. We observe that in this case the behavior is reversed with respect to what is found in the positively correlated case. The critical point and the size of the giant component for the anticorrelated networks are slightly higher than the one for the randomized counterparts. However, surprisingly this behavior depends on the feature threshold employed. Depending on the F₀, we might be in regimes in which disregarding the negative degree-feature correlations leads to an underestimation or to an overestimate of the robustness, see the inset of Fig. 3b. Similarly to the positively correlated case, Δ is non-monotonic and there is an intermediate threshold value for which the underestimation is maximum, but increasing F₀ even more we find the point at which the overestimation is maximum. This evinces the non-trivial phenomenology that arises even in these simple cases of ad hoc correlations.

Features related to the network construction

Let us now consider a case where the degree and feature distributions are not imposed exogenously but they emerge naturally from the network construction process. As an example, we consider RGGs, which are especially useful to model situations in which there is some kind of physical contact or proximity between the units of the system, and they find applications in areas as diverse as wireless sensor architectures⁴¹, population dynamics⁴², or consensus formation⁴³, to name but a few. Even though percolation has been widely studied in RGGs in the mathematical literature (see, e.g., Refs. ^44,45), critical points are known up to some bounds⁴⁶, and because of the strong topological correlations of RGGs, tree-like approximations, in general, are not as accurate as in the case of infinite-dimensional networks. On this basis, we will need to quantify the discrepancy between theory and simulations and how it scales with the link density.

We consider RGGs composed of N nodes, placed uniformly at random on [0,1)² with periodic boundary conditions. Two nodes are connected if they are within a distance r. We take the feature as the distance between a node and its closest neighbor d_min, i.e., F ≡ d_min, and we delete those nodes that have the nearest neighbor at a distance smaller than a certain threshold r₀ > 0, i.e., ϕ_F = θ(F − r₀) (see Fig. 4a). This scenario can be relevant, for instance, in spatial ecological models in which there is a competition for resources^47,48. Note that this particular choice of the occupation probability is made for illustrative purposes, and without much effort, one could envision other relevant scenarios, for example, with ϕ_F = θ( − (F − r₀)) or that take the distance to the furthest neighbor d_max as the feature. In all these mentioned cases, the generating functions can be calculated.

To compute the generating functions we first need the joint probability function. Its calculation is detailed in Supplementary Note 7. Setting a maximum bound k_max = N − 1 in the degree distributions, the generating functions are

The analytical calculations are compared with simulations in Fig. 4b. We keep the radius of interaction r constant and discuss the results as a function of the number of nodes in the network. First of all, we see that the position of the percolation point is inversely proportional to N. This behavior is somehow expected: the larger the system size, the denser the network, so the probability of nodes to have neighbors at distance smaller than r₀ becomes high, hence the rapid network dismantling. Moreover, we see that when the system size is low, the theoretical results systematically overestimate the value of the giant component. Note that this is not a finite-size effect as it occurs in other systems displaying phase transitions, but inherent limitations of the theory due to the topological correlations present in RGGs and not captured in their degree-feature distribution. We can quantify this discrepancy following the rationale of Eq. (13) but taking the absolute value in the integrand. Thus,

Fig. 4

Feature-based percolation in random geometric graphs.

In (a), a sketch showing how the removal of nodes proposed in the main text works: if a node i (central orange marker in this case) has its closest neighbor at a distance smaller than r₀ (shaded area), then i is removed. This removal condition is checked synchronously for all nodes at once, hence node j in the sketch is also removed because, since the network is undirected, i is a neighbor at a distance smaller than r₀. In (b), the size of the giant connected component is a function of r₀. The interaction radius is set to r = 0.1, above which links cannot be drawn, and kept the same for all simulations. Different system size N are explored, indicated in the legend. Solid lines correspond to the theoretical approximation from Eq. (17) and markers are results from simulations. In (c), a discrepancy between theory and simulations, Eq. (18), as a function of the system size. All simulation points are averaged over 100 independent realizations.

Features related to a dynamical process

As a final application of the feature-enriched percolation, we investigate the case in which the features are coupled to dynamics running on top of the network. There are many situations in which the dynamical evolution of some processes on the network generates some quantity, or attribute, that might not be evenly distributed across nodes. For instance, when studying the problem of synchronization, it is known that in the desynchronized phase, the synchronization error—the time-averaged distance in the phase space between the state of a node and the average state of the system—displays a power-law decrease with the degree⁵⁰. In such a scenario, one might be interested to study what is the surviving network after the removal of the most (or least) synchronized nodes with respect to the mean activity of the system. Another example is found in communication networks when modeling traffic, where nodes with a higher degree tend to be more congested on average⁵¹.

Here we focus on the SIS model, well-known in the context of epidemiology⁵². It consists of nodes that can be in either of two states, susceptible or infected. Infected nodes transmit the disease to susceptible ones at a certain rate upon encounter, and infected nodes recover spontaneously at a different rate. The dynamical evolution of the model is written as

For certain dynamical processes on certain types of topologies, one can calculate the exact joint degree-feature distribution and plug it into Eqs. (1)– (5). However, these will be marginal cases over all the ensemble of dynamical models and network architectures, and for many applications, the joint distribution will not be analytically available. To illustrate how one can proceed in this latter case, here we take an agnostic approach and use only information about the node degrees and their feature value to infer an approximate P(k, F). We first collect all the pairs (k_i, F_i) and compute a non-normalized two-dimensional histogram, see Fig. 5a. The collapse in the k − F plane tells us the type of correlation between degree and feature, that for the case of the SIS turns out to be positive. Note that for each value of the degree (recall that k is considered discrete), the distribution P(F∣k) has a bell-shape curve, of different height, different mean value, and different width, see the ridgeline plot of Fig. 5b for a better appreciation. The first strong approximation is to consider that each P(F∣k) is proportional to a normal form

Fig. 5

Feature-enriched percolation of the SIS model.

We integrate the dynamical Eq. (19) for a reshuffled Barabási–Albert⁷⁰, with m = 3 and system size N = 2000. A non-normalized histogram of the pairs (k_i, F_i), where the feature value is F_i = x_i(t → ∞), is shown in (a). The input data for the histogram comes from 100 independent realizations. In (b), ridge plot of the 25 first curves of the histogram, i.e., those with the lowest degree and highest peaks. The bell shape behavior with varying mean, standard deviation, and peak height can clearly be appreciated. In (c), projection of the k − F plain, together with the expressions μ_F(k) (see text) given by the Bayesian machine scientist method. The standard deviation σ_F(k) is not shown because it visually overlaps μ_F(k) for large k. The expressions are $$ {\mu }_F\left(k\right)=\exp \left(a_1/\left(b_1-k\right)\right)$$ and $$ {\sigma }_F\left(k\right)=\exp \left(a_2/\left(b_2-k\right)\right)-{\mu }_F\left(k\right)$$, and have been obtained imposing only one free parameter per function, namely v = − 1.1641 and w = 0.4317. The constants in the exponential are given by a₁ = − v, $$ b_1=-{\cosh }^v\left(v^2\right)$$, a₂ = 2w and $$ b_2=-{\cos }^2\left(w\right)$$. d To find the height of the probability peaks we use the BMS on the pairs of variables $$ \left(k,\max \left(P\left(F\mid k\right)\right)\right)$$, indicated by points, obtaining $$ h\left(k\right)=a_3{\left(x\left(b_3+x\right)\right)}^{-1}$$ where the free parameter is t = 26.1442 and the constants in the function are $$ a_3={\tan }^t\left(t\right)/2$$ and b₃ = t/2, indicated by the solid line. In (e), feature-enriched percolation of the SIS model, with the simulation (points) and the theoretical curve obtained from Eq. (21). Each point in (d) is averaged over 100 independent realizations of the dynamics.

In Fig. 5e we compare the results of the simulations with the curve obtained from the theory. We have employed the joint degree-feature distribution (21) to compute the generating functions, but we proceed numerically because no closed expressions can be found. We see that the agreement is quite good, despite the strong approximations used during the process. A small discrepancy around the critical region can be appreciated, rooted in finite-size effects and probably a systematic deviation that could potentially be reduced by employing more complicated functions in the process of constructing P(k, F), such as skewed Gaussian distributions. Anyway, we have shown that with little information one can find a quite satisfactory description of the percolation properties of systems in which the features are related to a dynamical process. This works reasonably well for synthetic networks, and we proceed next to test the accuracy of the feature-enriched percolation in real networks, which are characterized by topological correlations not taken into account in the theory.

To this goal, we use three different dynamics on three different real networks and proceed similarly as before, i.e., we will find the approximate degree-feature joint distribution by means of the BMS. The first example is a mutualistic dynamics arisen in symbiotic ecosystems, where the time dependence of the abundance x_i of the species i is given by the equation^54,55

Fig. 6

Feature-enriched percolation of dynamical processes on top of real-world networks.

In (a) mutualistic dynamics are given by Eq. (22), in (b) population dynamics from Eq. (23) and in (c) the biochemical dynamics (24) The main plots display the size of the largest connected component from simulations (markers) and its approximation given by the theory (lines). Insets show the data used to feed the Bayesian machine scientist (markers) and the resulting curves (dashed lines), used to calculate the degree-feature joint distribution (21). The output functions of the BMS are given in Supplementary Note 8. The feature is the value of the dynamical variable x_i at the steady-state, normalized by its maximum value, i.e., $$ F=x_i\left(t\to \infty \right)/\underset{i}{\max }\left[x_i\left(t\to \infty \right)\right]$$. The range of F in all the insets is the unit interval. The range of $$ \max \left[h\left(k\right)\right]$$ is arbitrary since its values depend on the binning to compute the histograms. All constants {τ_i} are set equal to 1. The initial condition is a uniform distribution in the unit interval. Each point of the main plot is computed as the average of 100 independent realizations of the dynamics.

The second example corresponds to birth-death processes⁵⁷. In the context of population dynamics, the temporal evolution of the population x_i in a site i can be described by

As a final dynamics, we use the mass-action kinetics model often employed in biochemistry⁵⁸. The equation

The dynamical processes presented here are merely illustrative examples of the potential and flexibility of the theory. For all the chosen examples the microscopic rules and the time evolution of the variable of interest were known, but we could have proceeded even without that information. There are two minimal ingredients to apply the feature-enriched percolation: the degree and the feature value for every node. If for whatever reason, a relation between the feature and the degree cannot be obtained, one can always proceed by employing the BMS technique, or other methods whose goal is to provide closed-form equations from data^60–62.