## Abstract

Identifying key players in coupled individual systems is a fundamental problem in network theory. We investigate synchronizable network-coupled dynamical systems such as high-voltage electric power grids and coupled oscillators on complex networks. We define key players as nodes that, once perturbed, generate the largest excursion away from synchrony. A spectral decomposition of the coupling matrix gives an elegant solution to this identification problem. We show that, when the coupling matrix is Laplacian, key players are peripheral in the sense of a centrality measure defined from effective resistance distances. For linearly coupled systems, the ranking is efficiently obtained through a single Laplacian matrix inversion, regardless of the operational synchronous state. The resulting ranking index is termed LRank. When nonlinearities are present, a weighted Laplacian matrix inversion gives another ranking index, WLRank. LRank provides a faithful ranking even for well-developed nonlinearities, corresponding to oscillator angle differences up to approximately Δθ ≲ 40°.

## INTRODUCTION

Because of growing electric power demand, increasing difficulties with building new lines, and the emergence of intermittent new renewable energy sources, electric power systems are more often operated closer to their maximal capacity (*1*, *2*). Accordingly, their operating state, its robustness against potential disturbances, and its local vulnerabilities need to be assessed more frequently and precisely. Furthermore, because electricity markets become more and more integrated, it is necessary to perform these assessments over geographically larger areas. Grid reliability is commonly assessed against *n* − 1 feasibility, transient stability, and voltage stability, by which one means that a grid is considered reliable if (i) it still has an acceptable operating state after any one of its *n* components fails, (ii) that acceptable state is reached from the original state following the transient dynamics generated by the component failure, and (iii) the new operating state is robust against further changes under operating conditions such as changes in power productions and loads. This *n* − 1 contingency assessment is much harder to implement in real-time for a power grid loaded close to its capacity where the differential equations governing its dynamics become nonlinear—the fast, standardly used linear approximation breaks down as the grid is more and more heavily loaded. Nonlinear assessment algorithms have significantly longer runtimes, which makes them of little use for short-time evaluations. In the worst cases, they sometimes even do not converge. Briefly, heavily loaded grids need more frequent and more precise reliability assessments, which are however harder to obtain, precisely because the loads are closer to the grid capacities.

Developing real-time procedures for *n* − 1 contingency assessment requires new, innovative algorithms. One appealing avenue is to optimize contingency ranking (*3*) to try and identify a subset of *n*_{s} < *n* grid components containing all the potentially critical components. The *n* − 1 contingency assessment may then focus on that subset only, with a substantial gain in runtime if *n*_{s} ≪ *n*. Identifying such a subset requires a ranking algorithm for grid components, following some well-chosen criterion. Procedures of this kind have been developed in network models for social and computer sciences, biology, and other fields, in the context of the historical and fundamental problem of identifying the key players (*4*–*8*). They may be, for instance, the players who, once removed, lead to the biggest changes in the other player’s activity in game theory or to the biggest structural change in a social network. That problem has been addressed with the introduction of graph-theoretic centrality measures (*9*, *10*), which order nodes from the most “central” to the most “peripheral”—in a sense that they themselves define. A plethora of centrality indices has been introduced and discussed in the literature on network theory (*9*, *10*), leading up to PageRank (*11*). The latter ranks nodes in a network according to the stationary probability distribution of a Markov chain on the network; accordingly, it gives a meaningful ranking of websites under the reasonable assumption that web surfing is a random process. Their computational efficiency makes PageRank and other purely graph-theoretic indicators very attractive to identify key players on complex networks. It is thus quite tempting to apply purely graph-theoretic methods to identify fast and reliably key players in network-coupled dynamical systems.

Processes such as web crawling for information retrieval are essentially random diffusive walks on a complex network, with no physical conservation law beyond the conservation of probability. The situation is similar for disease (*12*) or rumor (*13*) spreading and for community formation (*14*) where graph-theoretic concepts of index, centrality, betweenness, coreness, and so forth have been successfully applied to identify tightly bound communities.

Coupled dynamical systems such as complex supply networks (*15*), electric power grids (*16*), consensus algorithm networks (*17*), or more generally network-coupled oscillators (*18*, *19*) are, however, fundamentally different. There, the randomness of motion on the network giving, e.g., the Markov chain at the core of PageRank, is replaced by a deterministic dynamics supplemented by physical conservation laws that cannot be neglected. Pure or partially extended graph-theoretic methods have been applied in vulnerability investigations of electric power grids (*20*–*22*) and investigations of cascades of failures in coupled communication and electric power networks (*23*, *24*). They have, however, been partially or totally invalidated by investigations on more precise models of electric power transmission that take fundamental physical laws into account (in this case, Ohm’s and Kirchhoff’s laws) (*25*, *26*). It is therefore doubtful that purely topological graph-theoretic descriptors are able to identify the potentially critical components in deterministic, network-coupled dynamical systems. Purely graph-theoretic approaches need to be extended to account for physical laws (*20*). The influence of the dynamics on transient performance for regular graphs on *d*-dimensional tori has been emphasized in (*27*).

Here, we give an elegant solution to the key player problem for a family of deterministic, network-coupled dynamical systems related to the Kuramoto model (*18*, *19*). While we focus mostly on high-voltage electric power grids whose swing dynamics, under the lossless line approximation, is given by a second-order version of the Kuramoto model (*16*, *28*), we show that our approach also applies to other generic models of network-coupled oscillators. Key players in these systems can be defined in various ways. For instance, they can be identified by an optimal geographical distribution of system parameters such as inertia, damping, or natural frequencies or alternatively as those whose removal leads to the biggest change in operating state. Here, we define key players as those nodes where a local disturbance leads to the largest short-time transient network response. In the context of electric power grids, transient stability is the ability of the grid to maintain synchrony under relatively large disturbances such as loss or fluctuations of power generation or of a large load (*29*). If under such a fault, the system remains in the vicinity of its original state, it has maintained synchrony. There are different measures to quantify the magnitude of the transient excursion, such as nadir and maximal rate of change of the network-averaged frequency (*30*, *31*) or other dynamical quantities such as network susceptibilities (*32*) and the wave dynamics following disturbances (*33*). Here, we quantify the total transient excursion through performance measures that are time-integrated quadratic forms in the system’s degrees of freedom (see Materials and Methods and the Supplementary Materials). Transient excursions typically last 10 to 20 s in large, continental power grids, which sets the time scales we are interested in.

Anticipating on results to come, Fig. 1 illustrates the excellent agreement between analytical theory and numerical calculations for these performance measures. Particularly interesting is that in both asymptotic limits of quickly and slowly decorrelating noisy disturbance, the performance measures are simply expressed in terms of the resistance centrality (*34*, *35*), which is a variation of the closeness centrality (*9*) based on resistance distances (*36*). This is shown in the insets of Fig. 1. Our main finding is that the resistance centrality is the relevant quantity to construct ranking algorithms in network-coupled dynamical systems.

## RESULTS

We consider network-coupled dynamical systems defined by sets of differential equations of the form

The coupled individual systems are oscillators with a compact angle degree of freedom θ* _{i}* ∈ ( − π, π). Their uncoupled dynamics are determined by natural frequencies

*P*(

_{i}*37*), inertia parameters

*m*, and damping parameters

_{i}*d*. Because the degrees of freedom are compact, the coupling between oscillators needs to be a periodic function of angle differences, and here, we only keep its first Fourier term. The coupling between pairs of oscillators is defined on a network whose Laplacian matrix has elements

_{i}*i*≠

*j*and

*m*= 0 ∀

_{i}*i*, Eq. 1 gives the celebrated Kuramoto model on a network with edge weights

*b*> 0, ∀

_{ij}*i*,

*j*(

*18*,

*19*). With inertia on certain nodes, it is an approximate model for the swing dynamics of high-voltage electric power grids in the lossless line approximation (

*16*,

*28*,

*38*). The latter is justified in high-voltage transmission grids, where the resistance is smaller than the reactance typically by a factor of 10 or more. Applied to high-voltage grids, Eq. 1 describes the transient behavior of power grids on time scales of up to roughly 10 to 20 s. Over these time intervals, voltage amplitudes of high-voltage power grids are almost constant; accordingly, it is justified to consider only the dynamics of voltage angles (

*29*). Here, we are interested in that transient time regime and accordingly focus on the voltage angle dynamics given by Eq. 1. When angle differences are small, a linear approximation sin(θ

*− θ*

_{i}*) ≃ θ*

_{j}*− θ*

_{i}*is justified, giving first-order (without) or second-order (with inertia) consensus dynamics (*

_{j}*17*).

When the natural frequencies *P _{i}* are not too large, synchronous solutions exist that satisfy Eq. 1 with

*i*. Without loss of generality, one may consider Eq. 1 in a frame rotating with the synchronous angular frequency ω

_{0}, in which case, these states correspond to stable fixed points with

*39*). Linearizing the dynamics about that solution, Eq. 1 becomes

*i*≠

*j*and

We assess the nodal vulnerability of the system defined in Eq. 1 via the magnitude of the transient dynamics determined by Eq. 2 under a time-dependent disturbance δ*P _{i}*(

*t*). We take the latter as an Ornstein-Uhlenbeck noise on the natural frequency of a single node, with vanishing average,

_{0},

*k*= 1, …,

*n*nodes. This noisy test disturbance is designed to investigate network properties on different time scales by varying τ

_{0}and identify the set of most vulnerable nodes, i.e., the key players, as the nodes where the system’s response to δ

*P*(

_{k}*t*) is largest. Besides being a probe to test nodal vulnerabilities, these noisy disturbances alternatively model fluctuating renewable energy sources in electric power grids. In this latter case, however, the correlation time τ

_{0}is no longer a free parameter and is typically of the order of a minute or more, i.e., larger than any dynamical time scale in the system, as we discuss below. We quantify the magnitude of the response to the disturbance with the following two performance measures (

*40*)

They are similar to quadratic performance measures based on L_{2} or ℋ_{2} norms previously considered in the context of electric power grids, networks of coupled oscillators, or consensus algorithms (*30*, *40*–*46*) but differ from them in two respects. First, here, we subtract the averages Δ(*t*) = *n*^{−1}∑* _{j}* δθ

*(*

_{j}*t*) and

_{1,2}by

*T*before taking

*T*→ ∞ because we consider a noisy disturbance that is not limited in time and that would otherwise lead to diverging values of 𝒫

_{1,2}.

Here, we calculate 𝒫_{1,2} for the network-coupled dynamical system defined in Eq. 1 when (i) both inertia and damping parameters are constant, *m _{i}* ≡

*m*

_{0}and

*d*≡

_{i}*d*

_{0}, (ii) the inertia vanishes,

*m*≡ 0, (iii) the ratio γ ≡

_{i}*d*/

_{i}*m*is constant, and (iv) both inertia and damping vary independently. In cases (i) to (iii), 𝒫

_{i}_{1,2}can be analytically expressed in terms of resistance centralities that will be introduced in the next section (see Materials and Methods and the Supplementary Materials). The next paragraphs focus on case (i), following which we present numerical data for case (iv), which illustrate the general applicability of these results for not too short noise correlation time.

The performance measures 𝒫_{1,2} can be computed analytically from Eq. 2 via Laplace transforms (see Materials and Methods and the Supplementary Materials), for homogeneous damping and inertia, i.e., *d _{i}* =

*d*= γ

*m*, ∀

_{i}*i*. In the two limits of long and short noise correlation time τ

_{0}, they can be expressed in terms of the resistance centrality of the node

*k*on which the noisy disturbance acts and of graph topological indices called generalized Kirchhoff indices (

*36*,

*40*). Both quantities are based on the resistance distance, which gives the effective resistance Ω

*between any two nodes*

_{ij}*i*and

*j*on a fictitious electrical network where each edge is a resistor of magnitude given by the inverse edge weight in the network defined by the weighted Laplacian matrix. One obtains

*36*). The resistance centrality of the

*k*th node is then defined as

*C*

_{1}(

*k*) = [

*n*

^{−1}∑

*Ω*

_{j}*]*

_{jk}^{−1}. It measures how central is the node

*k*th in the electrical network, in terms of its average resistance distance to all other nodes—more central nodes have smaller

*C*

_{1}(

*k*). A network descriptor, the Kirchhoff index, is further defined as (

*36*)

Generalized Kirchhoff indices Kf* _{p}* and resistance centralities

*C*(

_{p}*k*) can be defined analogously from the

*p*th power of the weighted Laplacian matrix, which is also a Laplacian matrix (see Materials and Methods and the Supplementary Materials). In terms of these quantities, the performance measures defined in Eq. 3 depend on the value of the noise correlation time τ

_{0}relative to the different time scales in the system. The latter are the ratios

*d*/λ

_{α}of the damping coefficient

*d*with the nonzero eigenvalues λ

_{α}, α = 2, …,

*n*, of

^{−1}=

*m*/

*d*of damping to inertia parameters. In high-voltage power grids, they are approximately given by

*d*/λ

_{α}< 1 s and

*m*/

*d*≅ 2.5 s. Performance measures Eq. 3 can be obtained for any correlation time τ

_{0}(see Materials and Methods and the Supplementary Materials). However, it is interesting to consider the specific cases where τ

_{0}is the smallest (τ

_{0}≪

*d*/λ

_{α}, γ

^{−1}) or the largest (τ

_{0}≫

*d*/λ

_{α}, γ

^{−1}, appropriate for noisy power injections from new renewables) time scale in the probed system. The performance measures particularly take the asymptotic values

_{0}is the smallest or the largest time scale in the system. After averaging over the location

*k*of the disturbed node,

*40*,

*42*,

*43*) for the global robustness of the system.

These results are remarkable: They show that the magnitude of the transient excursion under a local noisy disturbance is given by either of the generalized resistance centralities *C*_{1}(*k*) or *C*_{2}(*k*) of the perturbed node and the generalized Kirchhoff indices Kf_{1,2}. The latter are global network descriptors and are therefore fixed in a given network with fixed operational state. One concludes that perturbing the less central nodes—those with largest inverse centralities _{0} (see Materials and Methods and the Supplementary Materials) is further confirmed in the main panel of Fig. 1 and by further numerical results obtained for different networks shown in Materials and Methods and the Supplementary Materials.

The generalized resistance centralities and Kirchhoff indices appearing in Eq. 6 depend on the operational state via the weighted Laplacian *P _{i}* ≪ ∑

*, ∀*

_{j}b_{ij}*i*, angle differences between coupled nodes remain small, and the weighted Laplacian is close to the network Laplacian,

*47*) and includes particularly many, but not all dead ends, which have been numerically found to undermine grid stability (

*48*).

The asymptotic results of Eq. 6, together with the numerical results of Fig. 1, make a strong point that nodal sensitivity to fast or slowly decorrelating noise disturbances can be predicted by generalized resistance centralities. One may wonder at this point how generalized resistance centralities differ in that prediction from other, more common centralities such as geodesic centrality, nodal degree, or PageRank. Table 1 compares these centralities to each other and to the performance measures corresponding to slowly decorrelating noisy disturbances acting on the 10 nodes shown in Fig. 2A. As expected from Eq. 6, 𝒫_{1} and 𝒫_{2} are almost perfectly correlated with the inverse resistance centralities

## DISCUSSION

Once a one-to-one relation between the generalized resistance centralities *C*_{1}(*k*) and *C*_{2}(*k*) of the disturbed node *k* and the magnitude of the induced transient response is established, ranking of nodes from most to least critical is tantamount to ranking them from smallest to largest *C*_{1} or *C*_{2}. From Eq. 6, which of these two centralities is relevant depends on whether one is interested (i) in the transient response under fast or slowly decorrelating noise or (ii) in investigating transient behaviors for angles (using the performance measure 𝒫_{1}) or frequencies (𝒫_{2}). While this gives a priori four different rankings, Eq. 6 leads to only two rankings, based on either _{1} only, in asymptotic limit of either very fast (shortest time scale τ_{0}) or very slowly (largest τ_{0}) decorrelating noise. From here on, we therefore focus on the angle performance measure 𝒫_{1} of Eq. 3a and consider the two asymptotic limits in Eq. 6a.

We therefore define WLRank_{1} and WLRank_{2} (*49*) as two rankings, which order nodes from smallest to largest *C*_{1} and *C*_{2}, respectively. Smallest WLRank_{1,2} therefore identify the most vulnerable nodes in a given network. Figure 3 shows that they differ very significantly. In particular, a number of nodes are among the most critical according to WLRank_{1} but not to WLRank_{2} and vice versa. This discrepancy means that nodes are not central in an absolute sense, instead, their centrality and hence how critical they are depend on details of the disturbance—in the present case, the correlation time τ_{0}—and the performance measure of interest. One should therefore choose to use one or the other centrality measure, according to the network sensitivity one wants to check.

The resistance centralities in Eq. 6 correspond to the network defined by the weighted Laplacian **θ**^{(0)}, and consequently, WLRank depends not only on the network topology but also, as expected, on the natural frequencies and the coupling between the nodal degrees of freedom. As mentioned above, in the strong coupling limit, angle differences between coupled nodes remain small and _{1,2} as the rankings using resistance centralities *P _{i}* ≲ ∑

*, and angle differences become very large. This case has been investigated for an inertialess coupled oscillator system on a random rewired network with constant couplings (see Materials and Methods and the Supplementary Materials) (*

_{j}b_{ij}*50*). It is shown in green in Fig. 4D and corresponds to max(Δθ) = 106°.

In Fig. 5, we investigate more closely when the approximate ranking LRank starts to differ from the true ranking WLRank. To that end, we used the randomly rewired model of inertialess coupled oscillators of Fig. 4D and calculated the percentage of nodes with highest LRank necessary to give the top 15% ranked nodes with WLRank_{2}. The results are plotted as a function of the maximal angle difference between directly coupled nodes. Each of the 12,000 red crosses in Fig. 5 corresponds to one of 1000 natural frequency vectors **P**^{(0)}, with components randomly distributed in (−0.5,0.5) and summing to zero, multiplied by a prefactor β = 0.4,0.6,…2.4,2.6. The blue crosses correspond to running averages more than 500 red crosses with consecutive values of max(Δθ). One sees that, up to almost max(Δθ) ≃ 40°, the set of the 18% of nodes with highest LRank_{2} always includes the top 15% ranked nodes with WLRank_{2}. Similar results for obtaining the top 10 and 20% ranked nodes with WLRank_{2} and for rankings using *C*_{1} instead of *C*_{2} are shown in Materials and Methods and the Supplementary Materials.

It is quite unexpected that nodal ranking remains almost the same up to angle differences of about 40°, since coupling nonlinearities are already well developed there. This is illustrated in the inset of Fig. 5, which plots the Frobenius distance *n*_{s} most critical nodes for any configuration with max(Δθ) ≲ 40°, including cases with nonegligible nonlinearities, is achieved with a single matrix inversion of the network Laplacian *n*_{s} + δ*n*_{s} nodes with highest LRank, δ*n*_{s}/*n*_{s} ≪ 1. This is a moderate price to pay, compared to the price of calculating WLRank for each configuration, which each time requires inverting the weighted Laplacian matrix

So far, we have assumed constant inertia and damping parameters, which led us to the analytical expressions given in Eq. 6 for the performance measures. Analytical results can further be obtained for inertialess systems with *m _{i}* = 0 as well as in the case of homogeneous damping to inertia ratio,

*d*/

_{i}*m*≡ γ. In this latter case, the ranking is again given by a resistance centrality, but this time, related to the inertia-weighted matrix

_{i}**M**being the diagonal matrix whose

*i*th diagonal entry is given by

*m*(see Materials and Methods and the Supplementary Materials) but not in the case of independently varying

_{i}*m*and

_{i}*d*. We therefore lastly address this more general case using a purely numerical approach. This question is especially important for electric power grids where only nodes connected to rotating machines (such as conventional power plants) have inertia and consumer nodes have significantly smaller damping parameters (

_{i}*38*). Time scales in electric power grids have typical values

*m*/

_{i}*d*∈ [1,3]s and

_{i}*d*/λ

_{i}*≲ 1s, and accordingly, we focus on the regime of large noise correlation time τ*

_{α}_{0}≫

*m*/

_{i}*d*,

_{i}*d*/λ

_{i}_{α}, which is appropriate for persisting power fluctuations such as those arising from renewable energy sources. Figure 6 shows results corresponding to inertia and damping parameters fluctuating randomly from node to node by up to 40%. The ranking obtained from a full numerical calculation is compared to the ranking obtained from a direct calculation of the centrality of the weighted Laplacian

_{0}in a much wider range of parameters than their derivation would suggest.

## CONCLUSION

We have formulated a key player problem in deterministic, network-coupled dynamical systems. The formulation is based on the dynamical response to a nodal additive disturbance of the initial problem, and the most critical nodes—the key players—are defined as those where the response to the disturbance is largest. While this manuscript focused on (i) noisy Ornstein-Uhlenbeck disturbances, (ii) network-coupled systems on undirected graphs, particularly with symmetric couplings *b _{ij}* =

*b*in Eq. 1, and (iii) performance measures of the transient response that are quadratic forms in the system’s degrees of freedom, the method is not restricted to these cases. First, it can be used to deal with different disturbances, and in Materials and Methods and the Supplementary Materials, we calculate

_{ji}*P*(

_{i}*t*) = δ

*δ*

_{ik}*P*

_{0}Θ(

*t*)Θ(τ

_{0}−

*t*) with the Heaviside function Θ(

*t*). This disturbance gives the same ranking as the Ornstein-Uhlenbeck noise disturbance considered above. Second, asymmetric couplings occurring, e.g., in directed graphs (

*51*), in Kuramoto models with frustration (

*19*), or in electric power grids with Ohmic dissipation (

*16*) can also be considered. In this case, the internodal coupling is given by asymmetric real matrices instead of symmetric Laplacian matrices. However, the definition of the resistance distance (Eq. 4) remains valid even if L is replaced by an asymmetric matrix A, in that it still gives Ω

*= 0, Ω*

_{ii}*≥ 0, and Ω*

_{ij}*≤ Ω*

_{ij}*+ Ω*

_{ik}*, ∀*

_{ki}*i*,

*j*,

*k*as long as the synchronous fixed point considered remains stable. Third, nonquadratic performance measures can, in principle, be considered within the spectral decomposition used in this article. One may think of average frequency nadir and rate of change of frequency, which are linear performance measures (

*30*,

*31*). It is, at present, unclear whether these quantities can be analytically related to the location of disturbances via resistance or other centralities.

We gave an elegant answer to this key player problem: Ranking nodes from most to least critical is tantamount to ranking nodes from least to most central in the sense of resistance centralities. Depending on how the problem is formulated—mostly on details of the disturbance and on how the magnitude of the transient response is measured—different centralities have to be considered, giving different rankings. The key player problem in deterministic systems is therefore not uniquely defined, and its formulation must be tailored to reflect the most relevant dynamical properties one wants to evaluate. Averaged rankings, reflecting several such properties, simultaneously could also be considered. Last, we found numerically that resistance centralities are still accurate to identify the most critical nodes even when nodal dynamical parameters (damping and inertia) are not homogeneous.

The results shown in Fig. 6 are rather unexpected, and further inspection of our analytical results (Eq. 6 and eq. S14B) suggests that an inertia dependence could emerge in the opposite limit of short correlation time τ_{0} ≪ *m _{i}*/

*d*,

_{i}*d*/λ

_{i}_{α}. This point deserves further investigations. It would be furthermore interesting to extend our investigations to cases of distributions of inertia and damping parameters corresponding to realistic electric power grids. Work along those lines is in progress.

## MATERIALS AND METHODS

Four different networks were considered. Data for implementing the IEEE 57 bus test case and the MATPOWER Pegase 2869 test case have been obtained from (*52*, *53*). The random networks were obtained through the rewiring procedure of initially regular networks discussed in (*50*). We constructed the European power grid model from publicly available data on power plants and bus locations. More details on the procedure were given in the Supplementary Materials and in (*54*).

Several operational states of the European power grid were obtained from an optimal power flow. The latter minimizes a cost function including production costs for all available power plants, under constraints that power flows do not exceed the thermal line capacity of each power line and that productions do not exceed rated powers for each power plant.

For each network considered, Eq. 1 was numerically integrated using a fourth-order Runge-Kutta algorithm. Values for the performance measures were then calculated by numerical integration of the obtained time-dependent voltage angle and frequencies. Resistance distances, centralities, and Kirchhoff indices were calculated through a direct inversion of the Laplacian matrix.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/11/eaaw8359/DC1

Section S1. Calculation of the performance measures

Section S2. Resistance distances, centralities, and Kirchhoff indices

Section S3. Numerical models

Section S4. Numerical comparison of LRank with WLRank.

Fig. S1. Comparison between theoretical predictions and numerical results for both performance measures 𝒫_{1} and 𝒫_{2}.

Fig. S2. Comparison of the performance measures 𝒫_{1}, 𝒫_{2} obtained numerically and in eq. S14.

Fig. S3. Percentage of the nodes with highest LRank necessary to include the nodes with 10% and 20% highest WLRank.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank R. Delabays and T. Coletta for interesting discussions.

**Funding:**This work has been supported by the Swiss National Science Foundation under grants PYAPP2_154275 and 200020_182050.

**Author contributions:**M.T. and L.P. wrote the numerical codes. L.P. constructed the model of the European power grid. M.T. performed the analytical and numerical calculations. P.J. formulated the ideas and wrote the initial draft. All authors participated in discussions of the results and writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).