Node-Level Resilience Loss in Dynamic Complex Networks

In an increasingly connected world, the resilience of networked dynamical systems is important in the fields of ecology, economics, critical infrastructures, and organizational behaviour. Whilst we understand small-scale resilience well, our understanding of large-scale networked resilience is limited. Recent research in predicting the effective network-level resilience pattern has advanced our understanding of the coupling relationship between topology and dynamics. However, a method to estimate the resilience of an individual node within an arbitrarily large complex network governed by non-linear dynamics is still lacking. Here, we develop a sequential mean-field approach and show that after 1-3 steps of estimation, the node-level resilience function can be represented with up to 98% accuracy. This new understanding compresses the higher dimensional relationship into a one-dimensional dynamic for tractable understanding, mapping the relationship between local dynamics and the statistical properties of network topology. By applying this framework to case studies in ecology and biology, we are able to not only understand the general resilience pattern of the network, but also identify the nodes at the greatest risk of failure and predict the impact of perturbations. These findings not only shed new light on the causes of resilience loss from cascade effects in networked systems, but the identification capability could also be used to prioritize protection, quantify risk, and inform the design of new system architectures.


Networks -Supplemenatary Information
Giannis Moutsinas, Weisi Guo I. EXPOSITION OF THE METHOD Let G be a weighted directed graph of N vertices and A be its weighted adjacency matrix. We associate each vertex v i with the function of time x i , that we call the state of v i and we assume that x i 's satisfy the dynamical systemẋ where f : R → R and g : R 2 → R, a ji = (A) ji . Our goal is to describe the equilibria of this system and their stability. We start by finding the mean field approximation of the system and we iteratively refine this prediction.
We make the following assumptions: 1) The functions f and g are twice differentiable and there exists R ∈ R such that for all for all |x|, |y| > R and w ∈ R + , (f (x) + wg(x, y)) x < 0.
2) The weights are independent and identically distributed random variables, with positive mean.
3) The graph has low degree correlation.
Assumption 1: This assumption implies that the state of the system cannot escape to infinity. For a given system we may be able to relax this a bit by assuming that it holds for w in a given interval.
In Section I-C we discuss briefly an example of how our method can fail if the first assumption is not satisfied.

Assumption 2:
The second assumption implies that the in-degree and the weighted in-degree of a vertex are correlated, i.e. the weighted degree of a vertex can be used as an indicator of how connected the vertex is. The smaller the mean is in comparison to the standard deviation, the less correlated the two degrees become. This increases the error of our approximation.
The assumption allows us to use the weighted in-degree instead of the in-degree of a vector when we calculate the average effect of the neighbourhood. On the extreme case that the standard deviation of the weights is 0 then the weighted in-degree is a multiple of the in-degree and this substitution does not produce any error. On the other hand, if the average weight is 0, then the two degrees are uncorrelated.
This means that a vertex can have high in-degree and at the same time weighted in-degree 0. In this case our approximation will find that the neighbours have not effect, whereas the neighbours' effect can be significant. Between these two extremes lies a continuum of possibilities. We find that when the standard deviation is of similar size to the mean, our approximation works well, as it is with the cases we discuss in the paper.

Assumption 3:
The last assumption implies the network is not assortative in any type of degree assortativity and that the neighbourhoods of the vertices are in a sense uniform. This is an important assumption that allow us to approximate the effect that the neighbourhood has on the vertex by the effect of the"average neighbourhood". Note that assumption 2 allows us to use the weighted in-degree to calculate the "average neighbourhood" of a vertex. In particular this assumption implies that given a vertex i, the probability of a vertex j being an in-neighbour of i is proportional to the out-degree of j and independent to the in-degree of i.
We find that for graphs with high assortativity our method has a significant error, for example spatial graphs. On the other hand for graphs that are constructed using the WattsStrogatz model that have high clustering coefficient but low assortativity, we find that our method has very good accuracy.

A. Notation
Let w i be the weighted in-degree of vertex v i , i.e.
and we denote by w av the (potentially weighted) average of all weighted in-degrees. We denote by w out j the weighted out-degree of vertex v j , i.e.
Similarly we denote by d i the in-degree of v i and by d out i its out-degree. We define the closed interval and lastly for a vertex, v, we denote by N (v) its neighbourhood.
We denote by 1 the N -vector, whose entries are all 1. We will call X ∈ R N a uniform vector if there exists x ∈ R such that X = x 1. Similarly, if an equilibrium is a uniform vector, we will call it a uniform equilibrium. Let X T = (x 1 , . . . , x N ), then we can rewrite the system (1) in the forṁ with F : R N → R N defined by equation (1). We define F : Note that we have not specified yet how exactly this mean is calculated. We will present two possibilities here, the uniform mean and the out-degree weighted mean. The uniform mean counts every component of the vector F (X) with the same weight and it is the most intuitive one. The out-degree weighted mean we weight the i-th component of the vector F (X) with w out i . This is less intuitive but it captures the fact that given a vertex v i the probability that a vertex v j is its neighbour is proportional w out j , see Section I-D for a brief explanation or [1] for an in depth discussion. We will denote this mean by E out In the latter case the mean field approximation is the same as in [2].

B. Mean Field Approximation
For the mean field approximation, we try to approximate the equilibrium of the system by a uniform vector. If the equilibrium is uniform, then there exists e ∈ R such that F (e1) = 0. Of course this is not true in general, so we are looking for e ∈ R such that F(e1) = 0.
As we noted in the previous section there are two somewhat obvious choices for the mean. We find that the uniform mean tends to underestimate value of the equilibrium and the out-degree weight mean tends to overestimate it. However when we used the latter, the next approximations have the same or sometimes better accuracy than when we used the former. Nevertheless, the uniform mean may still be proven useful, especially as a second guess about the dynamics. This opens up the possibility of defining a continuum of ways to take the mean using as weights (d out ) α with α ∈ [0, 1]. We defer the discussion of finding the optimal value for α as a function of the topological properties of the graph to future work.
We will be using the out-degree weighted average, so w av = E out [w i ]. This implies that This implies that every e that satisfies equation (3) potentially approximates an equilibrium of the system.
If e1 approximates well the equilibrium, we can look at the stability of the "average vertex" by assuming that the state of every vertex is fixed at e. More precisely we look at the slope s = f (e) + w av |g (1,0) (e, e)|.
If s > 0 the equilibrium that e approximates is unstable. However the condition s < 0 is not sufficient for stability of the system.
We will call the mean field approximation the zeroth order approximation of the equilibrium. is where Id N is the N × N identity matrix. For the eigenvalues, µ i , of J e it holds that µ i = h(e, e) λ i + f (e).

C. First Order Approximation
Once we have the mean field approximation, we can look at the dynamics of a single vertex. We choose a vertex v i and we assume that the state of all its neighbours is e {0} and that the state of v i does not affect them. Under these assumptions the dynamics of v i is given bẏ The equilibrium of the above dynamics will approximate the value of the equilibrium for v i . Let us define g {1} (x) = g x, e {0} and let χ {1} (w) denote the root of the function that is closer to e {0} . If χ {1} is not defined everywhere on W or is discontinuous, then this is an indication that the approximation is not accurate and the equilibrium we try to approximate might not exist. Similarly to the mean field approximation, in order to determine the stability of the system, we define Since e {1} approximates the equilibrium, F (e {1} ) approximates the Jacobian. Then by the Gershgorin circle theorem [3] we get that max s

Unboundedness of solutions
It is interesting at this point to discuss briefly the the case where Assumption 1 does not hold. As an example let us use We see that function f has an unstable equilibrium close to 5 and if a trajectory is on the right of this, then it escapes to infinity.
This means that we may have the case where the mean field approximation of the systeṁ has a stable equilibrium. However, when we look at the vertex level dynamicṡ of a node with high w i we may see that there is actually no equilibrium at all and that the trajectory 6 escapes to infinity. This will have a domino effect on the whole graph and every vertex will escape to infinity. This means that in such system there exists the possibility of a catastrophic failure that could in principle be undetectable, thus extra care is required before drawing any conclusions.

D. Average Effect of the Neighbours
For the first order approximation we looked at a vertex and we assumed that all its neighbours have the same state. Now that we have a better guess about the state of each vertex, we can update our guess about the states of the neighbours.
We pick a vertex v i and an in-edge of that vertex and we ask what is the probability that a given vertex v j is on the other end of the edge. This question is not solved in general for any distribution of random graphs. However, in the configuration model for random graphs, the probability for v j is approximately weighted by its out-degree d out j . The error of this approximation for large graphs constructed by the configuration model is of the order of 1 N . Under our assumptions it is a reasonable approximation for the graph G.
Let q i be a scalar quantity associated with vertex v i and Q = {q 1 , . . . , q N }, since the out-degree and the weighted out degree of a node are correlated, we define the weighted average We can define higher moments the usual way, for example the variance is defined by

E. Second and Higher Order Approximations
It is clear that we can continue the above process indefinitely. Let us assume that we know the n-th order approximation of the equilibrium, e {n} . Then just like before we define the function g {n+1} (x) = E out g x, e {n} and the function We find the roots of G {n+1} to get the function χ {n+1} and define e We can assess how well each step approximates the equilibrium by looking at the vector F e {n} . We know that if e * is an equilibrium, then F e * = 0. This means that the closer F e {n} is to 0, the better e {n} approximates e * . Then at each step we calculate F e {n} and compare with the previous. We choose to stop when the improvement becomes marginal.
Since we can continue this process indefinitely, we can define a sequence of approximations {e {n} } n≥0 .
Even though each approximation seems to be better than the previous and the sequence appears to converge to a limit in all examples we tried, we do not currently know whether this is actually the case and we will make no claim in this direction.
In order to determine the stability of the system, we define Similarly to what we had in the case of the first order approximation, max s  The bounds are not very useful when the system is away from a bifurcation point, because the approximations tend to be quite accurate. However it can be used to detect the range in which the bifurcation may happen. We will see their usage in two examples: mutualistic interactions in ecological networks and gene regulatory networks.

G. Assessing the accuracy of the approximation
Our approximation may not be accurate if the system is close to a bifurcation point. By using the bounds we can detect that. If the upper bound gives a resilient state and the lower bound gives a non-resilient state, then we know that the system is close to a bifurcation point.
Outside this region, the approximation is close to the equilibrium. Since for the equilibrium e * we have

II. APPLICATIONS
Similarly to [2], we will apply our technique in the case of mutualistic interactions in ecological networks and in gene regulatory networks. We have tried our method using graphs created by different models. We found that it works equally well in the case of Erdős-Rényi, Barabási-Albert and Watts-Strogatz models.
We decide to use the Erdős-Rényi model in our exposition.

A. Mutualistic Interactions in Ecological Networks
The dynamics in this model is given by We follow [2] and we set B = 0.1, C = 1, D = 5, E = 0.9, H = 0.1 and K = 5. As we saw earlier, the serves as the mean field approximation. In this case, as the parameter w varies the system undergoes two saddle-node bifurcations, see [4]. The first bifurcation happens at w − * ≈ −2.37 and the second at w + * ≈ 6.97. The bifurcation diagram of the dynamical system (6) is shown in Figure 1. For w < w − * the system has only one stable equilibrium x 0 * (w) that is close to 0, 0 < x 0 * (w) < 0.3. Similarly when w > w + * the system has only one stable equilibrium x + * (w) that does not approach 0, x + * (w) > 2.3. In the case when w − * < w < w + * , the system has both the above mentioned stable equilibria and between them another unstable one. Since we model mutualistic interactions, we do not have a physical interpretation of the case w < 0.
In this case x describes the population of a plant species. The stable equilibrium x 0 * (w) corresponds to the extinction of the species and we consider it undesirable. The fact that this equilibrium is not exactly 0 is the effect of migration. A system will be called resilient if the only stable equilibrium that exists is x + * (w). In this case the migration rate is enough bring the population back to abundance. In the non-resilient case, when x 0 * (w) is present, an eradicated species cannot return to abundance due to migration.
When the system is away from the region where the bifurcation can happen, our method predicts the equilibrium with very good accuracy. Figure 3a show the typical picture for Erdös-Rényi graphs.
Because the second order approximation is indistinguishable from the first in this scale, only the first order approximation is drawn.

1) Critical Weight:
The first order approximation of the vertex dynamics is given bẏ A saddle-node bifurcation happens also in this dynamics, so it is natural to ask at which is the critical value of w i , i.e. w crit ∈ R such that when w i > w crit , then there is only one stable equilibrium away from 0 and when w i < w crit , there are two stable equilibria, one of which is close to 0.
The critical weight, w crit , is a function of w av since it is a function of e {0} and e {0} is a function of w av . In Figure 2 we see the graph of w crit versus w av for the mutualistic interactions dynamics. Notice that because e {0} is discontinuous, w crit is also discontinuous.
The critical weigh can reveal some basic properties for the dynamics on a nodal level. For example we see in Figure 2 that when when w av > w + * , w crit is almost 0. This implies that if the system on average is in the resilient region, a vertex will also be in the resilient region even if it is very weakly connected to the rest of the network.
Notice that by definition the critical weight is graph-agnostic, since the mean field approximation is graph agnostic. We can use e {n} , n > 0, instead of e {0} to calculate the critical weight to get a better estimate. In this case we need to know the distribution of weighted in-degrees of the graph.  Finally in the case of weight reduction, we multiply each weight by a factor r ∈ [0, 1], starting with r = 1 and moving to r = 0 in 100 steps. Since this is a deterministic perturbation we do this only once.
Once more we plot the values of the equilibrium at each vertex. We see again that the uncertainty interval is useful to determine the region were the loss of resilience will happen.

B. Gene Regulatory Networks
The dynamics in this model is given by We follow [2] and set B = 1, α = 1, and β = 2. By looking at the bifurcation diagram, Figure 4, of we see that a saddle-node bifurcation happens at w * = 4. If w < w * , the only stable equilibrium of the system is x = 0. If w > w * a pair of a stable and an unstable equilibrium appears in the system. The equilibrium x = 0 implies cell death and it will be considered undesirable. The positive stable equilibrium corresponds to cell survival. We will call the system resilient if a positive stable equilibrium exists. x Fig. 4: The bifurcation diagram of (7). The stable equilibria are marked with a solid blue line and the unstable equilibrium is marked with a dashed red line.
The peculiarity of this dynamics is that the vertex dynamics is trivial, The second term of the right-hand side is constant so the equilibrium is just x = w i e β /(1 + e β ). Because of this the idea of critical weight that we discussed in the previous section cannot be applied here, as no bifurcation happens in (8). For the same reason the first order approximation is basically the best we can do in this dynamics. However, as Figure 3b shows, the first order approximation is already very accurate.
A final implication of the triviality of (8) is that once the mean field approximation becomes 0, then the next approximations and bounds all become 0. At each step we compute the 1st order approximation of the system along with the bounds, as explained above in this dynamics this is the can do. Because of the peculiarity of this system, the lower bound drops gradually to 0 and it never becomes 0 before the whole network collapses. The actual critical value below which activity is assumed to have ceased depends on the biological system that is being considered. Here we have chosen 0.1 as this critical value and we mark the region where the lower bound has dropped below it. The upper bound performs a more dramatic drop so we can use it as before. We also perform a numerical simulation of the system and compare this with our approximation. The results are shown in Figure 4 of the main text.
We consider three cases of perturbation, removal of vertices, removal of edges and uniform reduction of weights. These are performed the same way they were performed in the case of ecological networks. We see that in this system the loss of resilience happens with a less dramatic drop compared to the mutualistic interaction, yet it happens in the uncertainty interval we compute.
In the case of uniform weight reduction, because of the peculiarity of gene regulatory dynamics, once the mean field approximation becomes 0 everything collapses to 0. This is the reason that the uncertainty interval becomes very narrow and just misses the actual loss of resilience. We could get a wider uncertainty interval by repeating the process using the uniform mean for the mean field approximation. This would make the interval large enough that it would include the point where the loss of resilience happens.
However, in order to keep the exposition clean and consistent we decided against this idea.

C. Counter-Intuitive Dynamics
The dynamics studied in the previous sections follow the conventional wisdom, where a more connected network is always more resilient. This happens because the function G(x, y) was in both cases strictly increasing with respect to y. We will see that we can define dynamics where even though we can make a vertex resilient by adding more edges to it, we can also have the same result by removing edges from the network. This phenomenon is presented in Figure 6 of the main text. We will use the dynamics defined by g(x, y) = 15 x 2 y y 2 + 1 .
For this dynamics we used the same f as in mutualistic interaction in ecological networks but we changed the function g. Figure 9 shows the bifurcation diagram of the dynamics (9). Notice the similarities with the bifurcation diagram of the dynamics of ecological networks in Figure 1. We will keep the convention of section II-A and we will say that a vertex is resilient if there is only one stable equilibrium which is away from 0.
As an example we use an undirected graph with 30 vertices, 121 edges and edge weights equal to 1.
For this graph, the mean weighted degree is 8.9 and the average degree is 8. The weighted degree of vertex 11 is 2, which is the minimum. In Figure 6 we see the graph of the critical degree function for the dynamics chosen here. When w av = 8 the critical degree is approximately 2.31 and w av = 8.8 the critical degree is approximately 2.32. Because of this we expect that vertex 11 is not resilient. Numerical simulation proves that our prediction was correct.
If we add a new edge at vertex 11, we will bring the its degree to 3 which is in the resilient region.
Again we can check and see that numerical simulation proves us correct.
However, by looking at Figure 6 we can see that there is an alternative way to make vertex 11 resilient.
We can prune edges from our network that connect to high degree vertices and bring the average degree down. We remove 46 edges and bring the average degree down to 5 and the weighted average degree to 5.26, these correspond to critical degrees of 1.9 and 1.95 and it predicts that vertex 11 is not resilient.
Numerical simulation proves this prediction correct.

D. Approximation of equilibria of different graph models
In