Percolation in networks with local homeostatic plasticity.

Percolation is a process that impairs network connectedness by deactivating links or nodes. This process features a phase transition that resembles paradigmatic critical transitions in epidemic spreading, biological networks, traffic and transportation systems. Some biological systems, such as networks of neural cells, actively respond to percolation-like damage, which enables these structures to maintain their function after degradation and aging. Here we study percolation in networks that actively respond to link damage by adopting a mechanism resembling synaptic scaling in neurons. We explain critical transitions in such active networks and show that these structures are more resilient to damage as they are able to maintain a stronger connectedness and ability to spread information. Moreover, we uncover the role of local rescaling strategies in biological networks and indicate a possibility of designing smart infrastructures with improved robustness to perturbations.

Introduction

The resilience of complex networks to withstand failures and attacks is one of their most intriguing properties[1]. When edges have equal strengths, resilience can be studied with the mathematical framework of percolation[2]—a stochastic process akin to the permeation of a liquid through a porous membrane. This process removes edges uniformly at random with some given probability. The large-scale connectivity of the network is then studied by tracking the size of the largest connected component (LCC) as a function of this probability[3,4]. The size of LCC indicates what fraction of the whole network stays connected after damage, and in contexts such as transportation or communication—connected means functional. Studying percolation helps to understand the contribution of the network structure to its resilience, and both first- and second-order phase transitions in the size of LCC have been observed across a wealth of studies[1,4-8]. The most notable example refers to explaining the spreading of a disease driven by the susceptible–infected–recovered contact process. Percolation in a random network is identical to this process when the infectious time is constant[9] and leads to a useful approximation otherwise[10]. Percolation has been also used to understand the early stages of formation of the brain by probing how resilient are interconnected neuronal cultures in vitro[11,12].

When links feature different strengths (also called weights or capacities), both the structure and link weights contribute to network resilience, for example, as it happens with road and airline networks, the Internet, electricity grid, and financial networks. The role of such weights in network resilience can be mathematically studied with filtration[13,14], a process that removes all edges with a weight less than a given threshold. This process can be thought of as a permeation of a liquid with insoluble particles through a porous membrane, wherein the particles may flow through a pore only after they have clogged it up, otherwise, they remain filtrated.

Optimizing the distribution of link weights to secure a more resilient network is an old problem: in epidemiology, it is known as targeted immunization[15], but it also occurs in electrical engineering when designing efficient power grids, and urban planning when improving the traffic flow capacity of road networks. Apart from wide use in the top-down design, several biological complex systems were also observed to employ local self-regulation of link strengths with a positive effect on the global functioning of the whole network[16,17]. Both brain networks and food webs are believed to perform link strengths optimization by virtue of a self-regulatory mechanism that is triggered in response to damage.

In such systems, filtration-like processes are accompanied by an active response that mitigates the damage due to removed links. Neurons are hypothesized to control their activity using synaptic scaling[18,19], a mechanism that allows neurons to adjust their synaptic strengths to conserve the overall neural activity despite external perturbations or damage[19-21]. In a similar fashion, food webs representing the dynamical system for inter-species mass/energy transporting have also been observed to feature a high degree of coherence that may be a consequence of self-regulation and adaptive rewiring[22-25]. Load redistribution and rewiring mechanisms are thought to have shaped the core-periphery structure of the world airline network[26]. At the same time, self-regulation-inspired principles were also proposed to be used in top-down intervention scenarios for preventing species extinction in ecosystems[27].

Results

In the remainder, we present a mathematical theory for damage-response processes in complex networks that explains how maintaining local conservation of the total in- (or out-) weights of a node reflects on the global connectedness of the whole system. To this end, we consider a directed random network model with positive weights on edges. An instance of such a network is defined by a weighted adjacency matrix if node i points to node j, otherwise A = 0.

We model the adaptive network degradation by several iterations of the filtration stage, followed by an update of the weights of the survived edges according to a homeostatic plasticity principle, that is conservation of the sum of all in-weights for each node. This principle is inspired by the conservation of the total pre-synaptic strength observed in neurons[20]. Using this theory, we show that a simple local self-regulatory mechanism actively adjusting link weights in a fashion similar to the synaptic scaling in neurons[20], may significantly improve large-scale connectivity of the whole network and hence maintain network functioning, even if the damage has caused a loss of a large fraction of links. Note that the concept of homeostatic plasticity we adopt here is different from the one of homeostasis that usually appears in the literature of dynamical systems[28,29], which instead is generally related to the properties of stable fixed points in systems of ODEs.

To consider a general setting, our directed network model is defined by an arbitrary degree-weight distribution, f(x), the probability that a uniformly chosen directed edge has weight in the interval [x, x + dx] and terminates at the node of in-degree k. The joint distribution can be factorized, f(x) = lw(x), where the excess-degree distribution l, ∑l = 1 is the probability that a randomly chosen edge terminates at a node of in-degree k, and w(x) satisfying is the probability density function for edges that terminate at a node of in-degree k. One step of the damage-response cycle is introduced as an operator acting on the degree-weight distribution, where is the distribution before the edge removal and is the distribution after the damage-response cycle. This operator can be further decoupled as a convolution product of the damage and response operatorsHere, the contribution to the probability density of degree-k nodes comes from the nodes with higher degrees due to the fact that some edges are removed but none are added, and the surviving edges increase their weights as captured by the convolution operation (“*”) along the weight dimension x. See the “Methods” section for more details.

Our operator provides a general framework for studying the phenomenon of adaptive degradation. Let us now consider a particular instance of degradation/response mechanism wherein the former is represented by filtration, removing all edges with the weight below a given threshold y, and the latter is given by the following redistribution rule. For a node of degree k, the weights of the in edges x, i = 1, . . . , k are updated according towhere m is the total number of edges for which x > y, and is the sum of all removed weights. A simple check shows that similarly to homeostatic response in neurons, ∑x is conserved. Methods section provides the explicit forms of and for this mechanism.

In the first panel of Fig. 1, we present a sketchy representation of the process described by rule (2), while the second panel illustrates multiple successive steps of the filtration process with and without homeostatic response on an example of a small synthetic network. One can see that the homeostatic response mechanism mitigates the degradation process by enabling the network to withstand larger filtration thresholds. In Fig. 2 we use master equation (1) to quantify the evolution of the degree-weight distribution caused by the multi-step damage/response process in a larger synthetic random network of N = 5 × 104 nodes, and normally distributed weights. Validation of these predictions with stochastic simulations shows that our equation is highly capable to predict complex patterns revealed by the simulations.

Fig. 1

Simple filtration vs. filtration with homeostatic response.

a Pictorial representation of the damage–response process for a given node i. Here, in the damaged state, one of the three initial in-edges is removed, hence its weight, w2, is equally redistributed among the surviving edges, in order to conserve the local in-strength of the initial state. b Several snapshots illustrating an example of multi-step degradation on a small network, governed by: filtration, where every edge is removed if its weight is below a given threshold (first row), and filtration with the homeostatic response that maintains the total weight of in-edges at a constant level (second row). The bars below the panels indicate the distribution of edge weights (black) and the value of the threshold red.

Fig. 2

Impact of the homeostatic response on the weight distribution.

Evolution of the total weight distribution w(x) on four successive instances of filtration without response (upper panel), or with the homeostatic response (lower panel), on a large random regular graph of N = 5 × 104 nodes, degree z = 4 and normally distributed edge weights, at increasing values for the threshold y. The network has been generated according to the directed configuration model[42], while the edge weights are distributed according to Gaussian distribution with mean μ = 0.5 and standard deviation σ = 0.1. Each plot shows the empirical weight distribution from stochastic simulations (blue histogram) compared with the prediction derived from our model (orange solid line). The vertical dashed red line in each plot represents the value of y.

From the structural point of view, the connectivity of the whole network is characterized by the largest component in which nodes are path-connected regardless of the direction of the links—the largest weakly connected component. We estimate the size of this component, S, from the degree distributions for different increasing values of the threshold y, using the theory introduced by one of us[30]. In the multi-step setting, when the degradation occurs in a series of filtration-response steps, we observe that the disruption of S is significantly postponed when compared to the case of no active response.

Figure 3 illustrates the percolation transition being significantly delayed as a result of the homeostatic response in a random graph and an empirical brain network[31,32]. The values obtained from direct numerical simulations of S (colored markers) are well captured by the estimates (black× and +markers) obtained by the directed configuration model[30] supplied with the degree-distributions derived from the master equation (1).

Fig. 3

Homeostatic response delays the disruption of the giant component.

Evolution of the weak giant component’s size, S, as a function of the threshold y, for a a regular random graph of N = 2 × 104 nodes and degree z = 4 with Gaussian distributed weights (the same distribution of Fig. 2), and b a brain network of N = 983 nodes representing a volume of the mouse neocortex[31, 32] with uniformly distributed weights in the interval [0,1]. Colored markers represent the outcome of Monte Carlo simulations, obtained by averaging over 100 realizations of the weight distribution for the regular graph, and 1000 realizations for the real dataset. In both cases, we considered 20 successive steps of filtration without response (red circles) and with the homeostatic response (blue squares). Black markers represent the corresponding estimates from the model for each value of y, while gray dashed and dash-dotted lines represent linear interpolations between theoretical values, which serve only as visual guidelines.

Successive steps of the damage–response process may only occur in a system that quickly responds to an external perturbation. In some systems, synaptic scaling is reported to have a typical time of hours/days after external damage[20,33]. Therefore, we also analyze a single step of the process, mimicking a relatively slow response. As explained in the “Methods” section, in the single-step case, our model can be solved analytically provided the initial network does not have weight-degree correlations. Then even a single step introduces nontrivial correlations between weights and degrees (see “Methods” section for details), by which we lose the analytical tractability. Therefore letting w(x, y = 0) = w(x, y = 0) be the initial weight distribution, the average weight after a single step of the damage–response process at threshold y, , readswhere is the initial cumulative weight distributionand G(z) = ∑lz—the generating function of l.

The second term of Eq. (3) represents the contribution from the homeostatic response. As a matter of fact, one may verify that β1(y) coincides with the average weight in the case of no response (see the “Methods” section). On the other hand, the average degree after the damage is simply given bywhere stands for the initial average degree.

Figure 4 shows a very good agreement between our analytical predictions and stochastic simulations for the values of the average weight and the average degree , in both synthetic and real networks.

Fig. 4

Average degree and average weight for single instances of the damage–response process.

Colored markers represent the outcome of a single stochastic simulation. Solid and dashed lines represent the analytical results for average weight and degree (derived in the Methods section). Red dots and continuous lines represent the average degree (left scale), while blue triangles and dashed lines represent the average weight (right scale). a Poisson network of N = 25 × 103 nodes and uniform weights. b Random regular graph with N = 25 × 103 nodes and normally distributed weights (same distribution of Fig. 2). c Scale-free network of N = 25 × 103 nodes with power-law exponent γ = 2.5 and uniform weights. d Empirical network with uniform distributed weights (same network as in Fig. 3b).

Note that, by construction, our framework naturally takes into account the correlation between weights and degrees. In the “Methods” section we prove that, for a single instance of the damage–response process, the sign of this correlation is strongly affected by the underlying network topology, in particular, it is positive for scale-free networks, e.g., edges pointing to higher degree nodes have higher weights, while is negative for both random regular and Erdös–Rényi networks, e.g., edges pointing to higher degree nodes have lower weights.

From the results of Fig. 4, it is evident that is a decreasing function in y, while is an increasing function in y. One can verify from equations (3) and (6) that this is true in general. In the Methods section, we show that the product of the two, that is , coincides with the network average strength . This quantity tends to be quasi-conserved for small values of y in the absence of low-degree nodes (see Methods section), which is easily understood when we look back at the local rule (2): for each node i the local strength is conserved only if at least one edge survives the filtration process, otherwise, the lost weight is not redistributed.

We can use the analytical estimate for in order to approximate the leading eigenvalue of the network λmax: in the case of weighted directed networks, λmax is approximated by [34,35], which approaches the average strength as we neglect in/out strength correlations.

Hence, by combining Eqs. (3) and (6) we can assess the behavior of λmax as well as other spectral properties. To this end we consider, as an example, the trace of the communicability matrix, i.e., the Estrada index (EE)[36], motivated by the recent applications in the fields of both brain networks[37,38] and traffic flows in cities[39]. Note that in this case, since both λmax and EE depend on the weights, the effect of the damage/response process is observable also in the single-step case, while when we considered S, because the giant component is determined by the degree distribution only, the effect of the homeostatic response can only be seen in the multi-step case, for which we are sure that weight-degree correlations have formed.

Figure 5 shows that, similarly to the case of S, the homeostatic response has a clear effect of sustaining high values of both λmax and EE, compared to the case where no response is present. Moreover, in the case of a random regular network with uniformly distributed weights, the estimation of λmax by means of is quite accurate, both in the single-step case (panel a) and in the successive-steps case (panel b). Hence, we approximate the value of [36], by keeping the contribution of the largest eigenvalue plus the 0th order of the remaining part of the spectrum. Therefore, for this particular case we have . From Fig. 5, we see that this approximation is less accurate, but does not fail to capture the general behavior emerging from stochastic simulations.

Fig. 5

Leading eigenvalue and Estrada index for both single and multiple instances of the damage-response process.

Continuous and dashed black lines represent the values from our model in the case of homeostatic response and no response respectively, while colored markers and colored shaded areas represent mean values and interquartile ranges derived from stochastic simulations. We averaged the damage-response process over 200 realizations of a random regular network of N = 103 nodes and average degree z = 10 with uniform weights. a, c Single instances of the damage–response process. In this case is simply given by combining Eqs. (3) and (6). b, d Successive instances of the damage–response process. Here, both and are computed at each step from the master equation (1). In every panel, β1 is computed from Eq. (4).

Discussion

To summarize, percolation-like processes that conceptualize random damage in networks have been since long viewed as prototypical models for complex systems resilience. However, systems that actively maintain their homeostatic response tend to have the means to respond and actively adapt to such damage. We have presented a theory that naturally extends the classical percolation framework to a more complex one that incorporates the principle of homeostatic plasticity. Being a natural extension of simple percolation, our framework is still a theoretical abstraction but is not size-limited. Therefore, it is able to project the effects of local homeostatic mechanisms, similar to the ones observed in real biological systems at small sizes[16,40], at arbitrarily large scales, thus overcoming the practical limitations of laboratory experiments. By means of our model, it is possible to show that a simple local self-regulatory mechanism may significantly improve the large-scale functionality of the whole network compared to the case in which such a mechanism is absent. Our results reproduce the evolution of the joint weight-degree distribution of the network, which allows predicting the behavior of several global indicators of network structure and dynamics, such as the size of the largest connected component, the largest eigenvalue, and the Estrada index.

Overall our results provide a first mathematical framework for studying the link between local homeostatic plasticity rule in complex networks and its effect on the global functionality, and may also shed light on how the self-regulatory mechanisms observed in biological systems might be transferred to improve the resilience of human-designed infrastructures, for example, communication or transport networks, wherein it is reasonable to assume that homeostatic response might be adopted to mitigate external damage.

Methods

Notation and master equation

We consider a general model for homeostatic plasticity in a random directed network wherein all edges have a weight x > 0. According to this process, the network updates the edges' weights in a response to filtration, which is a deterministic and simultaneous removal of all edges with a weight strictly below a certain threshold y. Such a response can be adopted, for example, to mitigate further damage due to filtration in the future[20,41]. We consider networks wherein each node has at most one in-edge and possibly many out-edges and adopt the following notation to analyze the evolution of the network structure:

p is the in-degree distribution. Related to this quantity is —the average degree, and —the excess degree distribution, which denotes the probability of uniformly at random picking an edge sitting on a node of in-degree k.

w(x) is the probability density function such that is the probability of picking an in-edge with a weight greater than x, given that it points to a node with in-degree k. Note that plasticity may result in weight-degree correlation (see the last section), therefore .

f(x) ≔ lw(x) is the joint probability density function, corresponding to the event of picking an edge with weight w ∈ [x, x + dx] pointing to a node of in-degree k. Let us denote the marginal distributions of f(x) as

Since the response mechanism involves the weights of in-edges only, we refer to the in-degree as k. Response of degree-weight density function to filtration with threshold y > 0 is given bywhereIn Eq. (7), the contribution to the kth-degree nodes after filtration comes from nodes with degree n ≥ k. For these nodes, Damage operator conceptualizes the effect of filtration on the node degrees and edge weights, whereas the Response operator represents the response to filtration by increasing weights of surviving edges. The convolution operation is defined by .

The Damage operator

The Damage on the network is represented by a filtration process in which all the edges with weight below a fixed threshold y are removed. In order to get the explicit form of we compute the damaged excess degree distribution l(y) and the damaged weight distributions . The latter is simply given by the original weight distribution w(x, 0) cut and renormalizedwhere F(y) indicates the cumulative distribution function of w(x, 0), defined asOn the other hand by definition , hence we just need to compute the damaged degree distribution. For every node of degree n, each of the n in-edges may be removed with probability F(y), therefore we havewhere B(k, n, 1 − F(y)) denotes a binomial distribution with parameters (k, n, 1 − F(y)).

The average degree after damage is then given bywere we defined .

With this notation one can rewrite Eq. (10) in terms of the excess degree distribution l, which, combined with Eq. (8), finally yields the explicit form of the Damage operator

The Response operator

Here, we present the explicit form for the response operator introduced in Eq. (7) for our particular case in which the response is locally specified by the rule (2). In terms of a probability distribution, similarly to Eq. (8) we first compute the distributions of the deleted weights which is given byFrom Eq. (13) we get the explicit form of by considering the distribution associated with the sum of n − k deleted edges scaled by a factor kOverall we can now write the explicit form of the master equation (7) for this processBy repetitively applying Eq. (15) at increasing values for y we get the evolution of f(x) for a multi-step process. The plots presented in Figs. 2 and 3 of the Results section are derived by computing the marginals of f(x, y) at each step. In particular, the in-degree distribution is derived from , combined with Eq. (11) and the definition . On the other hand, the weight distribution is given by w(x, y) = ∑f(x, y). Finally, the out-degree distribution is simply given by directly applying Eq. (10) to the initial at any value for y, since the response mechanism involves weights of in-edges only, hence no additional correlation between out-degrees and edge weights is introduced in the process. With both in and out-degree distributions, we then estimate the size of the weakly largest connected component, S, for a directed configuration model[30], as threshold y is increased.

Single-step from an initial state with no weight-degree correlation

Here, we show that if the initial configuration is given by a network with no weight-degree correlation, the outcome of a single damage-response step is analytically tractable.

Let w(x, 0) = w(x, 0), so that the initial weight-degree distribution is given by f(x, 0) = l(0)w(x, 0) and the only cumulative weight distribution is given by .

Under these assumptions, Eq. (15) can be rewritten aswhere we highlighted in the square brackets the damaged weight distribution and the response operator. Let the convolution of the two be W(x). Its Laplace transform, defined as , readsFrom the properties of the Laplace transform, it is then straightforward to obtainwhere we definedAfter having introduced the two functions above, we can conveniently write aswhere is the generating function of l(0).

The results were obtained by testing the predictions of Eqs. (19) and (11) on several networks (both synthetic and real) are shown in Fig. 4 of the main article.

Average strength quasi-conservation

From Eq. (2) it is easily verified that the local strength is conserved if at least one edge survives. Here we show that from Eq. (19), one can verify that under certain assumptions the average strength of the whole network is quasi-conserved. We start by showing that . By definition, the in/out strength of a node S is defined by the sum of all the in/out weights. The strength distribution P(S = x) (either in or out) can therefore be written asNote that we allow the presence of weight-degree correlation for the most general case. However, we do not allow the correlation between weights belonging to the same neighborhood.

The average strength, , therefore readsLet us now consider Eq. (19). We note that since F(y) ∈ [0, 1], the contribution given by G(F(y)) can be neglected if low degree nodes are not present and F(y) is small enough. Under these assumptions, F(y) ≪ 1, i.e., the probability that for a node of degree k there are no surviving edges after damage, is very small. Therefore Eq. (19) can be approximated bywhere we used Eq. (11) with F(y) = F(y).

The above equation then impliesmeaning that is quasi-conserved by the response mechanism.

The emergence of weight-degree correlation

In this last section, we show that a single instance of the damage response process generates either positive or negative weight-degree correlations. Since the sign of the correlation coefficient is determined by the one of covariance, we can focus ourselves on the covariance only. By definition , which in our notation reads . Note that in our notation 〈k〉 is not the average degree, but it’s the first moment of the excess-degree distribution. By using again Eq. (17), one can obtainwhich combined with the previous results yieldsNote that β2(y) = 0 ⇒ cov(w, k; y) = 0, meaning that in the case of no response no additional weight-degree correlation is introduced, as expected.

From the last equation, we derive the conditionOne can verify the condition above is never satisfied for either Erdös–Rényi or Random Regular networks (for the case of Random Regular networks although the strict inequality is never satisfied, the equality holds for extreme case of constant degree z = 1), while for Scale-Free networks the opposite result holds.

In Fig. 6, we show the comparison between the analytic theory and stochastic simulations on a random regular network and a scale-free one, both with uniform weight distribution.

Fig. 6

Weight-degree covariance as a function of threshold y for a single instance of the damage–response process.

The initial configurations are given by: a a N = 25 × 103 nodes random regular graph with uniform weight distribution and b a N = 25 × 103 nodes scale-free network with exponent γ = 2.5, with uniform weight distribution (right panel). In both plots red markers represent the results from a single stochastic simulation, while continuous black lines show the analytical values from the theory. As we can see in the regular graph case the covariance is always negative, while the opposite happens in the scale-free network case.

Peer Review File