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.

T he 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 firstand 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 selfregulatory 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][20][21] . In a similar fashion, food webs representing the dynamical system for interspecies 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][23][24][25] . Load redistribution and rewiring mechanisms are thought to have shaped the coreperiphery 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 damageresponse 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 A ij 2 R þ if node i points to node j, otherwise A ij = 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 k (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 k (x) = l k w k (x), where the excess-degree distribution l k , ∑ k≥0 l k = 1 is the probability that a randomly chosen edge terminates at a node of in-degree k, and w k (x) satisfying R 1 0 w k ðxÞdx ¼ 1 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 A acting on the degree-weight distribution, f 1 k ðxÞ ¼ Af 0 k ðxÞ; where f 0 k is the distribution before the edge removal and f 1 k is the distribution after the damage-response cycle. This operator can be further decoupled as a convolution product of the damage D k;y and response R k;y operators Here, 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 A 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 , i = 1, . . . , k are updated according to where m is the total number of edges for which x i > y, and Δ ¼ ∑ fi;x i < yg x i is the sum of all removed weights. A simple check shows that similarly to homeostatic response in neurons, ∑ i x i is conserved. Methods section provides the explicit forms of D k;y and R k;y 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 × 10 4 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.
From the structural point of view, the connectivity of the whole network is characterized by the largest component in which 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, w 2 , 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.  30 . In the multi-step setting, when the degradation occurs in a series of filtration-response steps, we observe that the disruption of S W 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 W (colored markers) are well captured by the estimates (black× and +markers) obtained by the directed configuration model 30 supplied with the degreedistributions derived from the master equation (1).
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 k (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, where FðyÞ ¼ R y 0 wðxÞdx is the initial cumulative weight distribution and G(z) = ∑ k≥0 l k z k -the generating function of l k .
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 by where kð0Þ 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 wðyÞ and the average degree kðyÞ, in both synthetic and real networks.
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 kðyÞ is a decreasing function in y, while w 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 w k, coincides with the network average strength S. 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 S in order to approximate the leading eigenvalue of the network λ max : in the case of weighted directed networks, λ max is approximated by S in S out =S 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 W , 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 W , 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 S ¼ w k 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 EE ¼ ∑ n e λ n 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 EE $ N À 1 þ e λ max $ N À 1 þ e w k . From Fig. 5, we see that this approximation is less accurate, but does not fail to capture the general behavior emerging from stochastic simulations.

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 humandesigned 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 outedges and adopt the following notation to analyze the evolution of the network structure:  p k is the in-degree distribution. Related to this quantity is k :¼ ∑ k ≥ 0 kp kthe average degree, and l k :¼ kp k = k-the excess degree distribution, which denotes the probability of uniformly at random picking an edge sitting on a node of in-degree k.
• w k (x) is the probability density function such that Pðw>xjdeg ¼ kÞ ¼ R 1 x w k ðτÞdτ 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 Pðw > xjdeg ¼ kÞ ≠ Pðw > xÞ.
• f k (x) ≔ l k w k (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 k (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 f 0 k ðxÞ to filtration with threshold y > 0 is given by In Eq. (7), the contribution to the kth-degree nodes after filtration comes from nodes with degree n ≥ k. For these nodes, Damage operator D k;y conceptualizes the effect of filtration on the node degrees and edge weights, whereas the Response operator R k;y represents the response to filtration by increasing weights of surviving edges. The convolution operation Ã is defined by ðf Ã gÞðxÞ :¼ 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 D k;y f n ðxÞ we compute the damaged excess degree distribution l k (y) and the damaged weight distributions w D k ðx; yÞ. The latter is simply given by the original weight distribution w k (x, 0) cut and renormalized where F k (y) indicates the cumulative distribution function of w k (x, 0), defined as On the other hand by definition l k ðyÞ ¼ kp k ðyÞ= kðyÞ, hence we just need to compute the damaged degree distribution. For every node of degree n, each of the n inedges may be removed with probability F n (y), therefore we have Bðk; n; 1 À F n ðyÞÞp n ð0Þ; where B(k, n, 1 − F n (y)) denotes a binomial distribution with parameters (k, n, 1 − F n (y)). The average degree after damage is then given by were we defined hf k i y :¼ ∑ k ≥ 0 f k l k ðyÞ. With this notation one can rewrite Eq. (10) in terms of the excess degree distribution l k , which, combined with Eq. (8), finally yields the explicit form of the Damage operator D k;y f n ðxÞ ¼ k n Bðk; n; 1 À F n ðyÞÞθðx À yÞ ð1 À hF n ðyÞi 0 Þð1 À F n ðyÞÞ f n ðx; 0Þ: ð12Þ The Response operator. Here, we present the explicit form for the response operator R k;y 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 w c k ðx; yÞ which is given by From Eq. (13) we get the explicit form of R k by considering the distribution associated with the sum of n − k deleted edges scaled by a factor k R k;y f n ðxÞ ¼ ½kw c n ðkx; yÞ ÃnÀk ¼ w n ðkx; 0Þθðy À kxÞ F n ðyÞ=k Overall we can now write the explicit form of the master equation (7) for this process f k ðx; yÞ ¼ ∑ n ≥ k k n Bðk; n; 1 À F n ðyÞÞθðx À yÞ ð1 À hF n ðyÞi 0 Þð1 À F n ðyÞÞ f n ðxÞ ! Ã w n ðkx; 0Þθðy À kxÞ F n ðyÞ=k By repetitively applying Eq. (15) at increasing values for y we get the evolution of f k (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 k (x, y) at each step. In particular, the in-degree distribution is derived from l k ¼ R 1 0 f k ðxÞdx, combined with Eq. (11) and the definition p k ¼ kl k =k. On the other hand, the weight distribution is given by w(x, y) = ∑ k f k (x, y). Finally, the out-degree distribution is simply given by directly applying Eq. (10) to the initial p out k 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 W , 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 k (x, 0) = w(x, 0), so that the initial weight-degree distribution is given by f k (x, 0) = l k (0)w(x, 0) and the only cumulative weight distribution is given by where we defined After having introduced the two functions above, we can conveniently write wðyÞ as wðyÞ ¼ β 1 ðyÞ þ β 2 ðyÞ FðyÞ À GðFðyÞÞ 1 À FðyÞ ; ð19Þ where GðzÞ ¼ hz k i 0 ¼ ∑ k l k ð0Þz k is the generating function of l k (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 S is quasi-conserved. We start by showing that S ¼ w k. By definition, the in/out strength of a node S i 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 as Note 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, S, therefore reads Let 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 k (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 by wðyÞ ' wð0Þ 1 À FðyÞ ¼ wð0Þ kðyÞ= kð0Þ ; where we used Eq. (11) with F n (y) = F(y). The above equation then implies wðyÞ kðyÞ ' wð0Þ kð0Þ; ð20Þ meaning that S 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 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 × 10 3 nodes random regular graph with uniform weight distribution and b a N = 25 × 10 3 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.
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 cov ðk; wÞ ¼ E½kw À E½kE½w, which in our notation reads hk w k i y À hki y wðyÞ. 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 obtain hk w k i y ¼ β 1 ðyÞhki y þ β 2 ðyÞðhki 0 À hki y Þ; ð21Þ which combined with the previous results yields cov ðw; k; yÞ ¼ β 2 ðyÞ hki 0 G FðyÞ À À FðyÞð1 À GðFðyÞÞÞ 1 À FðyÞ : ð22Þ Note 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 condition cov ðw; k; yÞ ≥ 0 () GðFðyÞÞ ≥ FðyÞ hki y : One 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.

Data availability
The authors declare that all data supporting the findings of this study are available within the article.

Code availability
The code for reproducing the main findings of the article is available at the GitHub repository: https://github.com/grapisar/PercHomeostatic.