Network susceptibilities: theory and applications

We introduce the concept of network susceptibilities quantifying the response of the collective dy- namics of a network to small parameter changes. We distinguish two types of susceptibilities: vertex susceptibilities and edge susceptibilities, measuring the responses due to changes in the properties of units and their interactions, respectively. We derive explicit forms of network susceptibilities for oscillator networks close to steady states and offer example applications for Kuramoto-type phase- oscillator models, power grid models and generic flow models. Focusing on the role of the network topology implies that these ideas can be easily generalized to other types of networks, in particular those characterizing flow, transport, or spreading phenomena. The concept of network susceptibil- ities is broadly applicable and may straightforwardly be transferred to all settings where networks responses of the collective dynamics to topological changes are essential.


I. INTRODUCTION
Susceptibility constitutes a key concept in physics, from statistical mechanics to condensed matter theory and experiments. In these fields, susceptibility quantifies the change of a systems' state, typically measured by order parameters, in response to a change in some external field. In simple settings, susceptibility is well approximated by linear response theory and one global order parameter changes in response. Generally, there can be many order parameters, as for instance the site-dependent average spin in the theory of magnetism. While ideal solids are organized in the form of perfectly periodic crystals with, e.g., nearest-neighbor interactions, many natural and engineered complex systems are organized in networks with a rich variety of their underlying interaction topologies [1,2]. The susceptibility of such a networked system, i.e., its response to changes in their parameters, is thus essentially determined by their topology. Furthermore, unlike in periodic systems the response depends crucially on the location of the perturbation. Given that there are different types of local properties that may change, it is not yet clear how to appropriately define susceptibilities in a networked system and consequentially what such susceptibilities would tell us about the collective dynamics.
In this article we introduce two types of susceptibilities in network dynamical systems. Focusing on changes to steadystate operating points, we first systematically study the impact of small local perturbations of single units and effective interactions in networks. As a key class of network dynamics, we analyze the susceptibilities of oscillator networks describing the dynamics of various natural and manmade systems. We define both vertex susceptibilities and edge susceptibilities to qualitatively and quantitatively distinguish the responses to changes of single-unit and single-interaction properties, respectively. In particular, we reveal how the interaction topology of the network jointly with the type and location of the perturbation relative to the response location determine the response strength. These susceptibilities are shown to be related to, but not equal to, established measures of network centrality. Several applications, in particular to Kuramoto phase oscillator and power grid networks, are discussed. We specifically identify certain instances of vertex susceptibilities for electric power grid models as power transfer distribution factors known in electric engineering. Network susceptibilities are readily generalizable to all kinds of supply and transport networks as well as network dynamical systems whose dynamics exhibits a standard flow structure.

II. NETWORK SUSCEPTIBILITIES
A continuous time network dynamical system can be described by the equations of motion of N variables (the "vertices"), where p 1 ,p 2 ,...,p M are tunable. The network interactions (the "edges") are defined by which variables x j appear in the equation of motion of x i . Now we define network susceptibilities in the following. Definition 1. Let x * = (x * 1 ,x * 2 ,...,x * N ) be a steady state for a network dynamical system defined by Eq. (1). Suppose that on applying a small perturbation to one of the network parameters p k , the fixed point changes by a certain amount, Then the network susceptibility due to parameter p k is defined as We note that this definition can easily be extended to dynamics with other invariant sets (e.g., limit cycles) instead of fixed points, and also to stochastic dynamics.

III. DYNAMICS OF OSCILLATOR NETWORKS
As a cornerstone example we analyze the susceptibility of a network of coupled oscillators. The celebrated Kuramoto model [3] characterizes the collective dynamics of a variety of dynamical systems ranging from chemical reactions [4] and neural networks [5] to coupled Josephson junctions [6], laser arrays [7], and optomechanical systems [8]. In the Kuramoto model, N phase oscillators are coupled via their phase differences. The rate of change of each phase φ j is given by where ω j is the intrinsic frequency of the j th oscillator, j ∈ {1,...,N}, and K jℓ = K ℓj denotes the coupling strength of two oscillators j and ℓ. A similar model describes the frequency dynamics of complex power grids and has gained a strong interest recently [9][10][11][12][13][14]. The model describes the dynamics of rotating synchronous generators and motors, representing power plants and consumers, respectively. Each machine is characterized by the power it generates (P j > 0) or consumes (P j < 0) and rotates with a frequency close to the grid's reference frequency of 2π × 50/60 Hz, such that its phase is written as θ j (t) = t + φ j (t). The dynamics of the phases is given by the swing equation [15,16], where M j is proportional to the moment of inertia and D j is proportional to the damping torque of the respective synchronous machine. This "oscillator model" assumes that all consumers can be described as synchronous motors with a nonvanishing inertia M j . (It should be noted that since the oscillator model is valid only for the high-voltage transmission grid, the consumers do not represent individual electrical devices in each household, but rather whole cities or neighborhoods.) In the "structure-preserving model" used in electric power engineering [17] one assumes different consumers. In contrast to a synchronous machine this type of consumer cannot store any kinetic energy, such that the inertia vanishes. Hence, the equations of motion of the structure-preserving model are still given by Eq. (6), but with M j = 0. In the oscillator model as well as the structure-preserving model the power flow from machine k to machine j is given by where K jk is the maximum transmission capacity, which is proportional to the susceptance of the respective transmission line. The relative load of the transmission line is defined as The two models admit different forms of synchrony. The Kuramoto model was initially introduced to study the emergence of partial synchronization when the coupling of the oscillators is increased [3]. A power grid must be operated in a state of perfect synchronization: all phase differences φ k − φ j must be constant in time to enable a steady power flow [Eq. (7)]. In this article we analyze how such a phase-locked state responds to a local change in the network and in particular how this change depends on the topology of the network. Transforming to a corotating frame, the phase-locked states are then just the steady states of Eq. (5)orEq. (6), respectively, which are determined by the algebraic equation such that we can treat the Kuramoto model and the power grid model on the same footing. However, the perspective of a flow network is particularly helpful in understanding the mathematical results introduced below. We note that the steady states do not depend on the mechanical properties of the individual machines, i.e., the moments of inertia M j and the damping coefficients D j .

IV. LINEAR RESPONSE THEORY AND NETWORK SUSCEPTIBILITIES
In a complex network there are two general scenarios for a microscopic change of the dynamical system: (1) the modification of an edge weight (signifying, e.g., an electrical transmission line capacity) or (2) the modification of a vertex property (e.g., the power generation of a power plant) of the system. In the following, we introduce a linear response theory for both scenarios.

A. Perturbation at a single edge
In the first scenario we consider the coupling matrix K ij being perturbed slightly to yield the new perturbed matrix K ′ ij , which differs from K ij only at a single edge (s,t): κ ij = κ for (i,j ) = (s,t) and (i,j ) = (t,s) 0 all other edges. (11) This perturbation causes the steady-state phases of the network to change from φ j to φ ′ j . The new steady-state Eq. (9)n o w reads ∀ j ∈{1,...,N}.
In the following we calculate this perturbation within a linear response theory. We note that the steady state is defined only up to a global phase shift. Throughout this article we fix this phase such that j φ j = 0.
We expand the steady-state condition Eq. (13) to leading order in κ ij and and subtract Eq. (9), to obtain for all j = 1,...,N using the Kronecker symbol δ.I nt h e last step we have used the definition of the flow Eq. (7), the definition of relative load Eq. (8), and the perturbation matrix Eq. (11). Furthermore, we have introduced the matrix In a short-hand vectorial notation, Eq. (15) then reads using the vector q (st) ∈ R N with the components We note that the matrix A is singular such that it cannot be inverted. However, the vector q is orthogonal to the kernel of A, which is spanned by the vector (1,1,...,1) T such that this is no problem. In order to formally solve Eq. (17) we can thus use the Moore-Penrose pseudoinverse of A, which we will call T := A + in the following. Thus, we find The perturbed flow Eq. (7) over an edge (i,j ) is then given by up to first order in κ and ξ .UsingEq. (19), this result reads

B. Edge susceptibilities
Depending on the application we want to measure different effects caused by the perturbation at the edge (s,t). First, we quantify how much the phase of a single oscillator j is affected by the edge-to-vertex susceptibility, using Eq. (19), To measure the change of the oscillator state on a global scale in response to perturbation at a single edge, we define the global edge susceptibility as the norm of the local susceptibilities For applications to flow networks, such as the power grid model Eq. (6), we are especially interested in how the flows change as this determines the stability of the grid. In particular, stability can be lost when a single edge becomes overloaded. Thus, we define the edge-to-edge susceptibility as the change of flow at another edge, Using Eq. (21) this relation reads We conclude that the effects of a perturbation at a single edge (s,t) as measured by the susceptibilities defined above are proportional to the load of edge L st . We note that this edge susceptibility formalism can be used to detect [18]t h e phenomenon of Braess paradox, where the flow at the most loaded edge (i,j ) increases on increasing the coupling strength of one edge (s,t). To be precise, Braess paradox will occur on increasing the coupling strength at edge (s,t) if the edge-toedge susceptibility of the most loaded edge (i,j ) and the flow at that edge has the same sign, Furthermore, the susceptibilities are essentially given by the matrix T , the pseudoinverse of A. The properties of these matrices will be analyzed in detail in the following sections.

C. Perturbation at a single vertex
The above calculations can be readily generalized to analyze the change of the steady state in response to a local perturbation of a single vertex property. To this end we consider a change of the power injected at a single vertex s. However, a steady state of Eq. (6) exists only if the power is balanced 012319-3 such that we consider a small perturbation of the power vector of the form Expanding the definition of a steady state to leading order in p and ξ j := φ ′ j − φ j then yields Solving this equation for the changes ξ yields with the vector r s ∈ R N whose components are given by

D. Vertex susceptibilities
In analogy to the case of a perturbed edge discussed in Sec. (IV B) we define the vertex-to-vertex susceptibility as the global vertex susceptibility as and the vertex-to-edge susceptibility as We note that measures similar to the vertex-to-vertex susceptibility χ s→j are used in electric power engineering where they are called power transfer distribution factors [19,20]. In this context one generally uses a fixed reference or slack node, which absorbs the power change p, such that Eq. (27)i s modified to

E. Properties of the matrix A
We have shown that the response of a network to a local perturbation is essentially given by the matrix T , which is the Moore-Penrose pseudoinverse of the matrix A defined in Eq. (16). Before we discuss the potential applications of the network susceptibilities we thus have a closer look at the properties of the matrix A.
The matrix A encodes the dynamical stability and synchrony of steady states [21,22]. A steady state of the Kuramoto model or the power grid model defined by Eq. (6)i s dynamically stable if and only if A is positive semidefinite, i.e., all its eigenvalues a j ,1 j N are nonnegative. For the sake of simplicity we fix the ordering of the eigenvalues such that 0 = a 1 a 2 a 3 ...a N . We have to take into account that A always has one eigenvalue a 1 = 0. The corresponding eigenvector is (1,1,...,1), signifying that a small perturbation that is exactly the same in all phase angles is neutrally stable. However, this is merely due to the steady state itself being arbitrary up to a constant global phase shift. Stable steady states can emerge or disappear when a system parameter is varied through an (inverse) saddle node bifurcation at which one eigenvalue vanishes, a 2 → 0.
In particular, A is positive semidefinite if the relation cos(φ i − φ j ) > 0 holds for all edges (i,j ) of the network and the network is globally connnected. Stable steady states that do not satisfy this relation typically exist only at the edge of the stable parameter region [21]. We can thus assume that during normal operation we always have cos(φ i − φ j ) 0for all edges such that we can use the following relations: The expression K ij can be understood as the free capacity of an edge (ij ), which can be used to respond to the perturbation and is thus referred to as the responsive capacity.
For normal operation, cos(φ i − φ j ) 0 for all edges (i,j ), the nondiagonal entries of the matrix A are all nonpositive such that A is a Laplacian matrix for which many properties are known [2]. In particular, the eigenvalues of a Laplacian matrix satisfy 0 = a 1 a 2 ··· a N , where a 2 is an algebraic measure for the connectivity of the underlying network [23,24].

A. Scaling properties of network susceptibilities
The susceptibilities are especially large in the limit of a weakly connected network. For a power grid this corresponds to the scenario of high loads when the responsive capacities K ij become small. In the following we analyze this case in detail for a perturbation at a single edge (s,t) [cf. Eq. (11)]. The case of a vertex perturbation is discussed briefly at the end of this section.
Throughout this section we assume the case of "normal" operation; i.e., we assume that K ij 0 for all edges (i,j ). Then the matrix A is a Laplacian matrix with eigenvalues 0 = a 1 a 2 ··· a N and the associated eigenvectors v n . We can then formally solve Eq. (17)forξ with the result The term n = 1 does not contribute since we have fixed the global phase such that j ξ j = 0. This expression shows four important properties of the network susceptibility: (1) The response ξ and thus also the edge susceptibilities scale with the load of the perturbed edge L st = F st /K st .F or a complete breakdown of an edge (s,t), we have κ =−K st , such that ξ scales with the flow F st of the defective edge. The scenario of a complete breakdown is further discussed in Sec. VII.
(2) The prefactors 1/a n decrease with n. In particular for a weakly connected network the algebraic connectivity a 2 becomes very small [2,23,24], such that the term n = 2 dominates the sum. Then the susceptibility of all edges in the network scale inversely with the algebraic connectivity a 2 .This proves our claim that the susceptibility is large if the network defined by the responsive capacities K ij is weakly connected.
(3) For a weakly connected network, the edge susceptibility scales with the overlap |v 2 · q st |, where v 2 is the so-called Fiedler vector. This overlap can be interpreted as a measure of the local algebraic connectivity of the nodes s and t.T o see this, note that the Fiedler vector can be used to partition a graph into two weakly connected parts [2,24]. The overlap with the vector q st is largest if the two nodes s and t are in different parts and thus weakly connected.
(4) In the limit of a disconnected network the response ξ to a perturbation at the edge (s,t) diverges if the edge links the weakly connected components. If the perturbation occurs within one component, then the response remains finite. This will be shown in detail in the following section.
(5) The global edge susceptibility defined in Eq. (23) can be expressed as w h e r ew eh a v eu s e dE q . ( 36) for the phase response. This quantity measures the average phase response to the perturbation of a single edge (s,t). An example is shown in Fig. 1 for a synthetic power grid model based on the topology of the British high-voltage grid. One observes that the global susceptibility of an edge (s,t) is essentially determined by the load L st ,the connectivity of the network, and the location of the edge within the network. Edges are highly susceptible if they are heavily loaded or connect two components of the grid. In the shown example we observe two highly susceptible edges connecting the northern part to the rest of the grid. Averaging the global susceptibilities χ st over all edges (s,t) in the network, we find an almost perfect proportionality with the inverse algebraic connectivity 1/a 2 . If the transmission capacity K of the edges increases, the algebraic connectivity a 2 also increases and the grid becomes less susceptible to perturbations.

B. The weakly connected limit
To obtain a more quantitative understanding of the susceptibility in a weakly connected network we assume that the network is decomposed into two components of size N 1 and N 2 = N − N 1 , respectively. In the limit of complete disconnection, the Laplacian matrix also decomposes with As usual for a Laplacian matrix the lowest eigenvalue vanishes, a 1 = 0, and the associated eigenvector is given by In the disconnected limit also the second eigenvalue (the algebraic connectivity) vanishes, a (0) 2 = 0. The associated eigenvector, the Fiedler vector, is given by Here and in the following the superscript (0) denotes the limiting case of a complete disconnection of the network. For 012319-5 simplicity we assume that the two components are not further disconnected, such that a (0) 3 > 0. To analyze the case of a weakly connected network, we consider a single weak link at position (c,d) between the two components. The Laplacian is then given by A = A (0) + A ′ , with and A ′ ij = 0 otherwise. The connection strength k of the edge (c,d) is assumed to be small, such that we can calculate the eigenvalues and eigenvectors using Rayleigh-Schrödinger perturbation theory (see, e.g., Ref. [26]). We then find the algebraic connectivity, and the Fiedler vector, where q cd is defined as in Eq. (18).
To calculate the response of the network ξ we need the overlap of the vector q st [see Eq. (36)] with the eigenvectors of A, in particular the overlap with the Fiedler vector. The result depends crucially on the location of the perturbed edge (s,t). If this edge connects the two components, i.e., (s,t) = (c,d), we find such that the response diverges as k −1 : If the edge (s,t) lies within one component, then such that the response remains finite in the limit k → 0: For a perturbation at a single vertex as defined in Eq. (27) the response will always diverge in the limit k → 0. Assuming without loss of generality that the perturbed vertex s is an element of the component 1 we find that to leading order.

A. The relation to centralities
Various centrality measures have been defined to quantify the importance of single vertices and edges in complex networks [27]. Centrality measures based on current flows [28] are heavily used in different areas of network science and are directly related to susceptibility measures as defined in the present article. To illustrate this, consider a network of ohmic resistors with conductances G ij . An electrical current flows through the network with I source j being the current in-or outflow at vertex j . The current through a particular edge (i,j )ofthe network is given by the voltage drop across the edge, such that At each vertex the current is conserved such that Kirchhoff's law, is satisfied for all j = 1,...,N. Defining the Laplacian matrix of the conductances (the so-called nodal conductance matrix), and its Moore-Penrose pseudoinverse T := A + , the voltages are given by For the definition of centrality measures [27] one considers the situation that a unit current flows into the network at a single vertex s and out at a different vertex t. Then we have the voltages, and the current flowing over the edge (i,j ) is given by The current flow betweenness centrality of an edge (i,j )i s then defined as the absolute current flowing through the edge averaging over all scenarios of the in-and outflow, i.e., all pairs (s,t) [27]: Correspondingly, the betweenness centrality of a vertex j is defined as We directly see the analogies to the definition of the network susceptibilities if we identify the conductance G ij with the responsive capacity K ij . In particular, the edge betweenness centrality defined in Eq. (54) coincides with the average of the normalized edge-to-edge susceptibility η st→ij /|L s,t | except for a slight difference in the term (s,t) = (i,j ) that vanishes as 1/N 2 . However, in this article we generalize the idea of centralities based on current flow in several ways. First of all, we consider two different scenarios for the in-and outflow of the network: First, for the edge susceptibilities we consider an inflow at vertex s and outflow at vertex t with strength L st as in Ref. [27], Second, for the vertex susceptibilities we assume a unit inflow at vertex s and equal outflow at all other vertices, such that Third, we analyze not only the change of the flows, as in the edge-to-edge susceptibilities, but also the change of the state variables ξ j , which correspond to the voltages in the resistor networks. In this sense, the global susceptibilities χ 2 (st) and χ 2 s are given by the variance of the voltages in the network. Therefore, they quantify the global response of the network to a local in-or outflow in terms of the average variation of all voltages.

B. Relation to resistance distances
In a manner similar to Sec. VI A, the concept of susceptibilities can be understood in terms of resistance distance, which is defined as follows. As in the previous section we consider a network of Ohmic resistors with conductances G ij and suppose a unit current enters the node s and exits through node t. Then the resistance distance R st is given by the voltage drop between the nodes s and t. Using the relation Eq. (52), this yields using the symmetry of the matrix T . This relation can be inverted with the result [29] T ij =− where we have defined R tot i = j R ij . Substituting Eq. (59) into Eqs. (19) and (31), we can express all susceptibilities equivalently in terms of the matrix T or the resistance distances. For the vertex-to-vertex susceptibility we find and subsequently the global average of susceptibilities take the simple form t =s This relation clearly demonstrates that nodes that are on an average "close" to the rest of the network (i.e., with high centrality values), tend to have higher global susceptibility. In a similar manner, the vertex-to-edge susceptibilities can be expressed as and the edge-to-vertex susceptibility follows an almost identical form, apart from the prefactor: The global edge susceptibilities are given by (derivation in the Appendix) (64)

C. Scaling with distance
The effect of a linear perturbation generally decays with distance. To obtain a better understanding of this decay, we consider a continuum version of the linear response theory, concentrating on the vertex-to-vertex susceptibility. We consider a two-dimensional square lattice with equal weights, as power grids are naturally embedded into a two-dimensional plane and most grids can be assumed to be approximately planar. In the continuum limit the Laplacian matrix tends to the two-dimensional Laplace operator and Eq. (28) becomes a Poisson equation, where ξ (x) is the local response at position x (e.g., the local phase angle), p is the power injection, which occurs at position x 0 , and is the 2D Laplace operator. The solutions to this equation are well known. On an infinite two-dimensional domain it is where b is a constant of integration. Generally, no unique notion of Euclidean distance between nodes exists for networks. The closest analog is the shortest path distance, denoted by d(s,t) in the following, which is related to the Euclidean distance, for instance, in regular grids. Figures 2(a) and 2(b) show the decay behavior in a uniform square grid, compatible with the continuum results. Realistic network topologies are more complicated as shown in Figs. 2(c) and 2(d). We computed the susceptibility of the Continental European Transmission Network [30]t o perturbing one vertex for two cases of free capacitiesK. First, we obtained realistic valuesK ij,real from Ref. [30], then we considered a uniform model in which all free capacities are replaced by the averageK ij,unif = K real . In the vicinity of the perturbation, monotonic decay can be seen in both cases. However, there exist several vertices in the periphery of the network that are much more susceptible than the rest for realistic free capacities [dark blue in Fig. 2(c)]. These vertices are highly susceptible independent of the perturbed vertex.
Analogously to the case of vertex perturbation, the effects of edge perturbations can also be solved in the continuum limit, the result being the same as the potential due to an electrical dipole: The shaded region represents a 95% confidence interval. We clearly see logarithmic decay. (c) Color coded plot of χ s→t in the Continental European Transmission Network topology with realistic free capacities (taken from Ref. [30]). The vertex positions are not realistic. One central vertex was perturbed. There are several highly susceptible vertices in the network periphery (dark blue). (d) Decay behavior of the mean χ s→t in the same topology as in (c) as a function of d(s,t). We show both realistic free capacities K ij,real as well as uniform free capacitiesK ij,unif = K real set to the mean realistic value. The same vertex as in (c) was perturbed. In the realistic case, few highly susceptible vertices in the network periphery lead to a high variance at large distances, in contrast to the uniform case.
where q is the unit vector in the direction along which the perturbed edge lie.
This equation shows that unlike the response to vertex perturbation, the response to edge perturbation in a network will be highly directional. The susceptibilities should decay the fastest in the direction along the edge perturbed, according to the power law d −2 , consistent with the results presented in Ref. [31], but much slower in the orthogonal direction. In Figs. 3(a) and 3(c), we see this direction dependence in a regular square lattice. The lower envelope of the distancesusceptibility plot decays approximately as d −2 as expected.
We repeat the same analysis on the Continental European Transmission Network. We see that the susceptibilities are spread even wider for constant distance, indicating a stronger dependence on the orientation of the edge. We notice that the upper envelope in Fig. 3(d) decays very slowly: ≈d −0.4 , i.e., there exists a small but nonzero number of nodes that are heavily affected by the perturbation, despite being very far away from the perturbed edge.

D. Explaining the vulnerability of dead ends
The topology of a supply network determines its local [12,32] as well as global stability [9,10,33]. Recently, Menck et al. have shown that dead ends are particularly prone to instabilities [34]. They have measured the robustness of  d(s,j). The same edge as in (c) was perturbed. The susceptibilities for a single distance are even more widely distributed than in a regular lattice . The straight lines in (b) and (d) are algebraic fits to the upper and lower envelopes of the data set to obtain the exponent of the power-law decay. As the power-law decay breaks down near the boundary due to finite-size effects, we have to choose a cutoff, restricting the fit to the shaded region.
a power-grid model to large perturbations at a single node in terms of the so-called basin stability. To this end, the dynamics is simulated after a random perturbation to the steady state at a single node of the network. The basin stability is then defined as the probability that the network relaxes back to the steady state. Extensive Monte Carlo studies show that nodes adjacent to a dead end or dead tree have a particularly small basin stability.
The particular sensitivity of dead ends is directly related to the vertex-to-vertex susceptibility introduced in Sec. IV D. The main mechanism causing desynchronization at a dead end is shown in Fig. 4. The generation or power injection P s at a vertex s adjacent to a dead end is increased for a short period of time. This perturbation has a strong influence on the vertex s itself but also at the dead end t, causing a transient loss of synchrony. For longer times, the vertex s relaxes and resynchronizes with the rest of the network, whereas the dead end t does not. In summary, a perturbation at the vertex s has a large influence on the dynamics of the dead end, while its influence on the bulk of the network is small.
This property if directly mirrored by the vertex-to-vertex susceptibility χ s→j . Generally, the susceptibility is largest locally, i.e., for j = s, while it decreases with the distance as discussed above. Only if s is adjacent to a dead end t, the nonlocal susceptibility χ s→t is comparably large. One particular example is shown in Fig. 4(a). This shows that the nonlocal impact is strongest at dead ends and thus provides an explanation for their low basin stability.

A. From small to large changes
Linear response theory readily predicts how the flow in a network changes after a small perturbation of the network topology. But can it be used to estimate the effects of major changes such as the complete outage of an edge? This is especially important for electric power grids, where transmission line failures repeatedly induce large-scale outages (see, e.g., Refs. [35][36][37][38][39][40][41]). Thus, any method that helps to predict the stability of a grid after the failure of a single edge is extremely valuable. For an ad-hoc analysis of network stability in practical applications such methods should be only based on the topological and load properties of the original network and avoid time-consuming direct numerical simulations.
We can treat macroscopic changes within a linear response approach if we slightly modify the derivation of edge susceptibilities introduced in Sec. IV A. As before we keep only terms linear in ξ but we drop the assumption that the perturbations κ ij are small. Then Eq. (15) has to be modified as This set of linear equations is rewritten in matrix form as with the matrix where the superscript T denotes the transpose of a vector or matrix. The change of the local phases is then obtained by formally solving Eq. (69), In particular, we will need the phase differences between two nodes, which is given as This expression suggests that we need to calculate the inverse separately for each edge (s,t) if we want to assess the impact of all possible edge failures. However, we can greatly simplify the problem using the Woodbury matrix identity [42], which yields The network flows after the perturbation are now given by for all edges (i,j ) = (s,t). This expression differs from Eq. (21) only by the denominator, which tends to one in the limit of small perturbations κ → 0. For a macroscopic perturbation the denominator is essential to predict the magnitude of the flow changes correctly. The complete failure of an edge is described by κ =−K st , such that we obtain for all edges (i,j ) = (s,t) and F ′′ st = 0 for the failed edge. Similar formulas are used in power engineering, where the fraction is referred to as a line outage distribution factor (LODF) [19,43].
An example of how the damage of a single transmission line affects the flows in a power grid is shown in Fig. 5. We plot the change of the flow magnitude |F ij | predicted by the simple linear response approach, Eq.

B. Identification of critical edges
The modified formula, Eq. (74), can be used to predict impeding overloads and large-scale outages in complex supply networks [44]. Figure 6 shows the effect of the breakdown of a single transmission line for two examples. In the first example, Eq. (74) predicts that no overload occurs in agreement with the direct solution of the steady-state condition, Eq. (9). Thus, we expect that the grid relaxes to a new steady state after the failure of the respective edge. This prediction is confirmed by a direct numerical simulation of the equations of motion, Eq. (6). In the second example, Eq. (74) predicts that further overloads occur, i.e., that |F ′′ ij /K ij | > 1 for at least one edge (i,j ), after the transmission line (s,t) failed. Indeed, numerical simulations show that no steady state solution of Eq. (9)exists and that the grid becomes unstable and looses synchrony. In the following, we call an edge "critical" if its breakdown induces a desynchronization of the grid. If the grid relaxes back to a steady operation, i.e., an attractively stable synchronized state withφ j = 0 for all j , we call the edge "stable." Based on these results we propose to use the maximum load max (i,j ) |F ′′ ij /K ij | predicted by the modified linear response formula, Eq. (74), as a criterion to infer network stability. An edge (s,t) is predicted to be "critical" or "stable" according to the following classification system: where h is a threshold value. Bridges, i.e., edges whose removal disconnects the grid are always predicted to be critical.
To test this method, we perform direct numerical simulations of the equations of motion, Eq. (6), for a large number of test grids, each starting from a stationary state of normal operation and study the influence of the breakdown of a single edge. Examples for both scenarios are shown in Fig. 6. We analyze the coarse-grained structure of the British high-voltage transmission grid [9,25], which has 165 edges. We consider 100 random realizations with random generator positions, thus testing 16 500 edges in total. For each out of 100 random realizations, we fix the network topology by randomly selecting half of the nodes to be generators (P j =+1 P 0 ) and the others to be consumers (P j =−P 0 ), with P 0 = 1s −2 . The transmission capacity of all edges is fixed as K ij = K 0 = 4s −2 . One example of such a network is depicted in Fig. 1. Networks not supporting a steady state before any edge breakdown were discarded.
To evaluate the performance of the proposed classification scheme, Eq. (75), we must first define the possible outcomes of a prediction, where we distinguish between two different kinds of prediction errors: True positive: edge is predicted critical and is critical; False positive: edge is predicted critical but is stable; False negative: edge is predicted stable but is critical; True negative: edge is predicted stable and is stable. Generally, it is impossible to rule out both false-negative and false-positive predictions such that a compromise must be achieved. In the current setting, the number of false-positive predictions can be minimized by choosing a high value of h, while the number of false-negative predictions can be minimized by choosing a small value of h.
A quantitative assessment of the performance of a classifiers is then provided by a receiver operating characteristics (ROC) curve (Fig. 7) [45]. Here, the true-positive rate of the test, also called the sensitivity, SEN := no. of true-positive predictions no. of true-positive pred. + no. of false-negative pred.
, 012319-11 approaches the perfect limit (FPR,SEN) = (0,1) and clearly outperforms a classifier based on the load of the edges. Therefore, linear response theory provides a very promising approach to identifying critical infrastructures in complex supply networks.

VIII. CONCLUSION AND OUTLOOK
In summary, we introduced the concepts of vertex susceptibilities and edge susceptibilities as measures of responses to parametric changes in network dynamical systems. They qualitatively distinguish-and quantify-the responses due to changes in the properties of units and their interactions, respectively. Focusing on steady-state responses of oscillator network characterized by phases or phases and their velocities, we derived explicit forms of such network susceptibilities.W ein particular analyzed the role of irregular interaction topologies as those are the least investigated compared to the susceptibilities that are standard in physics. Specifically, we have analyzed how the responses of a network in some given phase-locked state depend on the relative location of perturbation and response sites and how the network topologies enter. We linked susceptibilities to established measures, for instance, special cases are known as line outage distribution factors in power-grid engineering and susceptibilities are closely related to centrality measures. We explicated an accurate prediction of network responses not only to small perturbation but also after the full breakdown of edges. In power grids, this may be applied, for instance, for an ad hoc security assessment. Furthermore, network susceptibilities directly reveal weak points of flow networks and may thus be used in the planning and design of future grid extensions and establishing other supply network infrastructures. Finally, the two types of network susceptibilities are generic measures of responses to parameter changes and as such may be straightforwardly generalized across flow, transport, and supply networks as well as other network dynamical systems where responses are nonlocal due to genuine collective dynamics.