Deciphering the imprint of topology on nonlinear dynamical network stability

Coupled oscillator networks show a complex interrelations between topological characteristics of the network and the nonlinear stability of single nodes with respect to large but realistic perturbations. We extend previous results on these relations by incorporating sampling-based measures of the transient behaviour of the system, its survivability, as well as its asymptotic behaviour, its basin stability. By combining basin stability and survivability we uncover novel, previously unknown asymptotic states with solitary, desynchronized oscillators which are rotating with a frequency different from their natural one. They occur almost exclusively after perturbations at nodes with specific topological properties. More generally we confirm and significantly refine the results on the distinguished role tree-shaped appendices play for nonlinear stability. We find a topological classification scheme for nodes located in such appendices, that exactly separates them according to their stability properties, thus establishing a strong link between topology and dynamics. Hence, the results can be used for the identification of vulnerable nodes in power grids or other coupled oscillator networks. From this classification we can derive general design principles for resilient power grids. We find that striving for homogeneous network topologies facilitates a better performance in terms of nonlinear dynamical network stability. While the employed second-order Kuramoto-like model is parametrized to be representative for power grids, we expect these insights to transfer to other critical infrastructure systems or complex network dynamics appearing in various other fields.


Introduction
Many critical infrastructure and supply systems (e.g., transportation, health care or power supply) are based on structures which can be described in terms of complex networks [1,2,3,4,5,6]. Such real world systems often evolved for the primary goal of fulfilling a specific function while also subject to certain constraints (e.g., financing or geography). An issue which is increasingly attracting notice in various fields (or from science to policy making) is the resilience of such critical infrastructure against perturbations in the form of external shocks, internal failures or changing environmental conditions [7,8,9,10], i.e. often non-small perturbations. Ultimately, it is highly desirable to gain a deep understanding of how optimal functionality on the one hand and resilience on the other hand can be achieved simultaneously. A still open question in this context is how the stability and resilience of networked systems is interrelated with their topological properties.
There are many concepts for assessing the stability of states within multistable dynamical systems. In classical linear stability analysis (in terms of Lyapunov exponents) the focus is on the local properties of a system's phase space. With the concept of Master Stability Functions this approach has been very successfully transferred to networked dynamical systems [11,12], in particular allowing the prediction of the synchronizability of coupled oscillator networks. The concept of Lyapunov functions [13,14,15] and basin stability in contrast present a nonlinear stability measure which also accounts for non-small perturbations and hence emphasizes a global perspective [16] on the dynamics. These concepts only account for the asymptotic behaviour of perturbed trajectories, however, for real-world systems the transient behaviour can be just as relevant to ensure their proper functioning. The survivability of a deterministic dynamical system is a suitable concept which complements both linear and asymptotic approaches [17] by focussing on transients. A complementary approach, studying the timing aspects of transient behaviour, can be found in [18] and [19,20,21,22]. In the context of control systems, questions of transient stability are explored in Viability theory [23,24] and under the name of robust control.
In this study we investigate the collective dynamics of power grids which can serve as prototypical examples of critical infrastructure systems [25,26,27]. The employed model is the (second order) Kuramoto model of coupled oscillators [28], hence the insights are rather general and can be of relevance to other fields. In particular, we use basin stability and survivability to assess the asymptotic and transient stability of power grids against large nodal perturbations, respectively. As opposed to other approaches, these sampling-based measures allow us to study relatively high dimensional systems, and to localize perturbations on the underlying network. These (nodal) stability measures are then related to purely topological characteristics of the respective nodes in the network. In this way we are able to identify a small number of topological classes of nodes which are characterized by similar stability properties. Remarkably, we also find a novel, previously unknown asymptotic state in the system that can only be accessed by perturbations at a particular topological class of nodes. The findings complement and extend previous results on the relationship between topological motifs and stability within power grids [29,30,31], and demonstrate the power of sampling based methods for studying the properties of dynamical systems on networks.

Oscillator model for nodal dynamics
Many biological, chemical and technical systems of N coupled oscillators can effectively be described by the Kuramoto model [32,28], which is given by the following temporal evolution law for an oscillator i's phase φ i : where P i /α denotes the oscillators natural frequency and K ij reflects the coupling strength between nodes i and j, satisfying K ij = K ji > 0 if nodes i and j are connected and K ij = 0 otherwise. Hence, the K ij define the topology of a (weighted) network of oscillators. α i denotes a damping coefficient of node i. In many contexts the dynamics of coupled oscillators is described more accurately when inertia is accounted for. This is obtained by the second-order Kuramoto model which additionally describes the temporal evolution of the node's frequency ω i : This model was shown to effectively describe the dynamics of synchronous machines within power grids [33,34,35,36] where the φ i denote the phase angles and the ω i the frequency deviations from the grid's rated frequency. In this context the P i correspond to the net power input at node i which is positive (negative) for generator (consumer) buses and the α i describe the strength of the electro-mechanical damping and droop control. The coupling coefficients K ij correspond to the capacity of the transmission lines of the power grid. For typical parameter values, the dynamical system given by (2.2) and (2.3) features a stable steady state (φ * , ω * ) = (φ * 1 , ...φ * N , ω * 1 , ..., ω * N ) with constant phase angles φ * i − φ * j and vanishing frequency deviations ω * i = 0 at every node i, given by a solution the following system of nonlinear equations: This state corresponds to the desirable synchronous operating mode of the power grid. However, depending on the parameter choices, there might be also undesirable nonsynchronous states which correspond to attracting limit cycles within the phase space of the dynamical system.
Menck et al. showed that for a single-node, connected to an infinite power grid, the limit cycle solution can be approximated as provided that |P | /α 2 1 and |P | 2 /α 2 K [29].

Random growth model for network topologies
In order to generate a representative ensemble of spatially embedded power grid topologies, we used a suitable random growth model [30]. It aims at generating synthetic network topologies that reproduce topological properties of real-world power grids and other spatially embedded infrastructure networks. The two-phase algorithm starts with a minimum spanning tree of size N 0 to which further nodes and (redundant) lines are added iteratively. The network growth is subject to a heuristic redundancy-versuscost optimization, which takes not only the line lengths but also additionally-created redundancy in the form of alternative routes into account. The growth model parameters have been set to (N 0 = 1, p = 1/5, q = 3/10, r = 1/3, s = 1/10), where p, q and s are probabilities related to the creation and splitting of lines, and r specifies the redundancyversus-cost trade-off. For a detailed explanation of all parameters, we refer to [30]. For this choice of parameter values the randomly generated networks match characteristics of real-world power grids, for instance the sparsity with a mean degree of aboutd ≈ 2.7. The ensemble consists of M = 50 random networks of size N = 100, an exemplary topology is shown in Figure 1.

Topological classification scheme for nodes in tree-shaped parts
We further make use of a topological classification scheme for nodes, which particularly distinguishes nodes located in or adjacent to tree-shaped parts of the networks. We start by giving concrete definitions for trees and tree-shaped parts and their roots: is called a tree if it is connected and has no cycles.
Definition 2.2 (Tree-shaped part, root) Let G = (V, E) be an undirected graph that is connected but not a tree. A tree-shaped part is an induced subgraph T = (V , E ) of G which is a tree and is maximal with the property that there is exactly one node r ∈ V that has at least one neighbour in G − T . r is called the root of T and has degree d(r) ≥ 3.
The union of all nodes located in tree-shaped parts is subsequently denoted by T = i T i and called the "forest" part of G, where the T i are the different tree-shaped parts within G. The remaining parts of G are referred to as the bulk, denoted by B = G − T . A finer partition of T is achieved by distinguishing the root nodes R from the non-root nodes N . The non-root nodes can be further subdivided into the leaves L = {l ∈ N | d(l) = 1}  which have degree one, and the inner tree nodes I = N − L which are located between the root and the leaves. We will see that for stability assessment, an even finer partition is useful for whose definition we first need to introduce the following properties of nodes in tree-shaped parts: Definition 2.3 (depth, height) Given a tree-shaped part T = (V , E ) of a graph G = (V, E), the depth δ(x) of a node x ∈ T is the length of the shortest path from x to the root of T . The root r ∈ T has depth δ(r) = 0. The height η(x) of a node x ∈ T is the length of the longest outward path from x to a leaf of T . All leaves l ∈ T have height η(l) = 0.
Note that the presented definitions of height and depth of nodes might appear counter-intuitive when applied to the picture of trees growing upward from the root. However, this terminology originates from the data structure of a "tree" in informatics, which is typically depicted as growing downwards, and became standard in graph theory.
The smallest possible type of tree-shaped part consists of a root and some adjacent leaves. Such leaves are subsequently termed sprouts and form the class S = {x ∈ N | η(x) = 0 ∧ δ(x) = 1} = {x ∈ L | δ(x) = 1}. The leaves of larger tree-shaped parts are called proper leaves and form the class Finally, the group of sprouts can be separated into those which are connected to high-degree roots, called dense sprouts S d = x ∈ S |d N > 5 , and those connected to rather low-degree roots, the sparse sprouts S s = x ∈ S |d N < 6 whered N (x) denotes the (average) degree of the neighbour(s) of x. For this last distinction, we chose the threshold of 5 so that the stability properties of the two groups are separated best.
For an efficient algorithm which provides both the partition of the nodes into the topological groups and their respective height and depth levels see Appendix A.
The nodes in the exemplary network in Figure 1 are coloured according to this classification and a representative node of each group is labelled accordingly. Definitions and total shares of the node categories in the ensemble of randomly generated network topologies are summarized in Table 1. Table 1. Overview of names, symbols and definitions of the hierarchically ordered topological groups of nodes in tree-shaped parts of networks. The last column shows the shares of nodes of each category in the ensemble of the M = 50 randomly generated network topologies. More than half of the nodes belong to tree-shaped structures within the networks and about a quarter is given by leaf nodes. An exemplary network topology with the nodes coloured according to these groups is shown in Figure 1. The simulation results shown in Figure 5 are also coloured according to this scheme.

Group
Symbol Definition Share of all nodes

Stability measures
The subsequently introduced measures assess the stability of a dynamical system with respect to large perturbations. They reflect global characteristics of the system's phase space which distinguishes them from the local perspective taken in conventional linear stability analysis in which only small (infinitesimal) perturbations are regarded [11,37]. These measures are particularly suited for assessing the stability of power grids, since in this context large perturbations from the normal operating state are a common threat [38].
The first measure to be introduced, basin stability, is an indicator for the system's likelihood to asymptotically return to a desirable state following a large perturbation [16]. The second, survivability, in turn is sensitive to whether the transient behaviour after a larger perturbation remains within a desirable region of the system's phase space [17].
In complex networks, it is instructive to regard only localized perturbations originating at a single node which makes the stability measures node-wise quantities. However, it should be noted that for the stability assessment the response of the whole system with all nodes is relevant and hence the violation of the stability constraints might happen at other nodes than the perturbed one.
Both measures necessitate the specification of a probability distribution from which large (finite) perturbations are drawn. For the case of power grids modelled by (2.2) and (2.3) these are given by values chosen uniformly at random (δφ, δω) ∈ [−π, π] × [−∆ω, ∆ω] which are added at t 0 = 0 to the state variables of a single node j while all unperturbed nodes are initialized to the desirable steady state: where δ ij denotes the Kronecker delta. Recall that ω * i = 0. Figure 2. Schematic of the single-node model's phase space. The union of the blue, green and yellow areas is the synchronous state's (X ) basin of attraction B while trajectories starting from the remaining parts of the phase space (B LC ; red and white) converge to the (non-synchronous) limit cycle located around ω = P/α (grey dashed). The union of the red, green and yellow areas forms the subset X 0 from which the random perturbations are drawn and (here) coincides with the desirable region X + relevant for the survivability measure. The green coloured region shows the (infinite-time) basin of survival X S . While trajectories starting within the yellow region converge to the synchronous operating state and are thus asymptotically stable their transient leaves the desirable region |ω| ≤ ω + and hence do not "survive" the perturbation.

Basin stability
The concept of basin stability is applicable to multistable dynamical systems with state space X for which there exist attracting states distinct from the set of desirable attracting states, where the latter is denoted by X ⊂ X in the following. The basin of attraction B of X is given by all initial states from which the system asymptotically converges to the desirable attractor: When the perturbations are drawn from a certain subset X 0 ⊆ X, it is instructive to define the basin of attraction restricted to this region: X B = B ∩ X 0 (cf. Figure 2). This is especially the case if the phase space is not compact. In the case of perturbations drawn uniformly at random, the basin stability β is then simply the ratio of the volume of the (restricted) basin of attraction X B to the overall region of perturbations X 0 [16,39]: In other words, this quantity corresponds to the probability for the system to return to the desirable attractor after a perturbation from X 0 . For the synchronous machine power grid model described by (2.2) and (2.3) the desirable attractor X is identical to the set of synchronous states (φ * , ω * = 0), while there exist several non-synchronous attracting states which are undesirable. Since basin stability is determined for each individual node j, the region of perturbations is the subset Hence the single-node basin stability corresponds to the ratio of areas in the phase space cross-section spanned by the dimensions associated to the node j (cf. Figure 2). As the basin of attraction and its geometry are typically not known a priori and difficult to determine, especially in high-dimensional systems [40], β needs to be determined numerically via a Monte-Carlo method. For each node L = 200 independent perturbations have been chosen uniformly at random from X 0 j and the corresponding trajectories simulated for t = 100 time units. The basin stability β of node j can be estimated from the number s of trajectories which asymptotically return to X . More details are given in Appendix B.

Survivability
In contrast to basin stability the survivability concept [17] presumes a desirable region X + ⊆ X which must not be left by a trajectory for a time t following a large perturbation in order to call the system "survived". The finite-time basin of survival X S t is given by the fraction of those initial states of the system which give rise to evolutions that stay within X + until time t: The infinite-time basin of survival is obtained by taking the limit X S = lim t→∞ X S t . The survivability σ of a system with respect to uniformly drawn random perturbations from the region X 0 is then analogously given as the ratio (cf. Figure 2) Hence it can be regarded as the probability for the system to remain within a desirable region X + after being hit by a perturbation from X 0 . For the operation of power grids it is necessary to keep the frequency deviations of all generators and consumers below a certain level [41]. Hence for the synchronous machine model the desirable region of the state space is given as with ω + > 0 being the maximally tolerable frequency deviation. Again for each node j the perturbations are chosen at random from X 0 j (cf. (2.10)) and the basin of survival is determined node-wise as X S j ⊆ X 0 j ( Figure 2). Note that it is not of relevance at which node the frequency constraint |ω| ≤ ω + is violated for the perturbation to count as not survived. Subsequently, the desirable region boundaries are chosen identical to the maximal perturbation level of the frequency deviations, ω + ≡ ∆ω.
In order to estimate the single-node survivability σ j for each node, L = 200 trajectories with independent perturbations to the synchronous operating state, drawn uniformly from X 0 j (cf. 2.10), have been simulated. The fraction of trajectories which did not leave the desirable region X + has been used as a statistical estimator for σ j (cf. Appendix B for details).

Simulations
Single-node basin stabilities and survivabilities have been estimated for all nodes in an ensemble of M = 50 randomly generated networks with N = 100 nodes each. Within each network half of the nodes act as net generators (P = +1), while the remaining nodes are net consumers (P = −1) such that there is an overall power balance, N i=1 P i = 0. Even though the synthetic power grid topologies generated by the random growth model are spatially embedded, the coupling strengths of the transmissions lines have been chosen uniformly to K = 6 for simplicity and to highlight purely topological effects. The damping coefficients have been set uniformly to α = 0.1. Note that for this choice of parameters the preconditions for the approximation of the limit cycle position, equation (2.5), are met.

Interrelations between a node's stability and its topological properties
Firstly, we investigated interrelations between the node-wise stability measures and topological properties of the perturbed node. We highlight topological effects in the example of a node's degree d and shortest-path betweenness b [3]. They turned out to reveal the most prominent insights for transient and asymptotic stability, respectively, and are basic established local/mesoscale network characteristics. Furthermore, this choice facilitates a comparison with previous findings [29,30].

Basin Stability
We regard different maximal perturbation levels ∆ω = 2.5 to 12.5, and observe that as expected the mean basin stability of all N · M = 5000 nodes decreases with growing perturbation levels ∆ω.  In order to detect which topological node characteristics influence asymptotic stability, the basin stability scores are regarded in dependence of the degree d of the node which is defined as the number of neighbouring nodes (Figure 3 a). While there is no significant dependency for lower perturbations levels (∆ω ≤ 5.0) there is a slight increase of β with d for larger perturbation levels, for which there are, however, relatively large standard deviations in the stability estimates. Hence, degree alone is a weak predictor for basin stability which is in line with previous findings [29,30].
There is, however, a characteristic dependence of basin stability on the betweenness b of the node, which is defined as the number of shortest paths between all pairs of nodes within the network which are passing through the regarded node (Figure 3 b). While there is no general trend for any perturbation level, distinctive down-peaks of basin stability are observable at certain betweenness values (as illustrated by Menck et al. [29]). These particular values of b correspond to nodes which lie within tree-shaped parts of the network (see the green-coloured nodes of the "inner tree node" category shown in Figure 1). A similar dependency was found in [29] for a maximal perturbation level of ∆ω = 100. Hence we were able to qualitatively reproduce these findings for considerably smaller perturbation levels which appear more realistic by comparison to real-world cases.
3.1.2. Survivability Next, we studied the transient stability against large perturbations as measured by single-node survivability, again for different values of the maximal perturbation strength ∆ω (Figure 4). The generally larger survivability scores achieved with increasing perturbations can be explained by the fact that the desirable region boundaries ω + are increased simultaneously with ∆ω (cf. Section 2). Another convention, which is not followed in this study, would be to hold the desirable region (ω + ) constant when increasing the perturbations (∆ω). In this case the average survivability of all nodes would decrease since for ∆ω > ω + a certain fraction of trajectories would start outside X + and hence could not "survive". In contrast to basin stability, the single-node survivability features a significant negative correlation with the degree which is found for all levels of perturbations (Figure 4 a). This means that the frequency deviations within the network following a perturbation tend to be larger when nodes with a high degree are hit, thereby making the dynamics leave the desirable region. For example, for ∆ω = 12.5 the probability of the trajectory to survive a perturbation is close to 1 if a node with degree d = 1 is hit, while it lies below 0.5 for nodes with degree d = 13. For smaller perturbation levels the same trend is observable. Hence the degree may serve as a suitable predictor for the node's survivability.
For the betweenness measure the same but less striking trend is found (Figure 4  b). Survivability is generally decreasing with the node's betweenness. Inner tree nodes which are characterized by specific values of b do not stick out as dominantly as it was the case for basin stability. For lower perturbations (∆ω ≤ 5.0) we observe small uppeaks in σ at the specific b-values while for ∆ω = 10.0 the same down-peaks as with basin stability occur. This reveals that nodes adjacent to tree-like structures are also crucial for predicting survivability, however, in a more subtle way compared to basin stability. Overall, betweenness alone is only a weak predictor for survivability, showing different features at different perturbation levels.
As for basin stability comparing survivability to other nodal network measures did not lead to more insights. Instead, a direct comparison of basin stability and survivability of the nodes was found to reveal non-trivial aspects of the system's dynamics and helped in identifying the node classification scheme introduced in 2.3.

Joint distributions of basin stability and survivability
While this section focuses on the general distribution patterns of basin stability and survivability, the stability characteristics of the different node classes is described in the next section (3.3). In order to reveal more insights about the interdependencies of asymptotic and transient dynamics as well as the relation to the network topology, we plotted the joint distribution of single-node survivability σ and single-node basin stability β for different perturbation levels ∆ω ( Figure 5). Each panel shows a total of N × M = 5000 individual estimates.
Let us focus on the two marginal distributions first. For a rather low maximal perturbation level (∆ω = 5.0) the distribution of basin stability is extremely skewed with 73.0 % of nodes featuring β ≥ 0.95. The survivabilities in turn show a widely spread bimodal distribution. This fact shows that survivability is generally much more influenced by topology than basin stability at lower perturbation levels.
At larger perturbation levels (∆ω ≥ 7.5) the distributions of both basin stability and survivability are unimodal and skewed. For ∆ω = 7.5 still 62.4 % of nodes feature a high asymptotic stability of β ≥ 0.9. The distribution of transient stability is still wide but shifted towards larger values. While some groups of nodes show a strong correlation between β and σ, the overall Pearson correlation coefficient is close to zero (ρ = 0.07).
This picture is reversed when looking at the highest perturbation level of ∆ω = 12.5. Here 62.9 % of nodes have survivabilities of σ ≥ 0.9 and the distribution of basin stability is rather widespread. It is known that at considerably larger perturbation levels (∆ω = 100) also basin stabilities show a widespread multimodal distribution [29,42]. Hence basin stability might be a more useful measure for the system's stability if very large perturbations are expected.
The intermediate case with ∆ω = 10.0 is particularly interesting. Here the distributions of β and σ are very similarly shaped and numerous nodes show a strong correlation. However, there are also many nodes for which there is no such correlation, with the overall Pearson correlation coming out at ρ = 0.41.

Basin stability and survivability for different classes of nodes
The findings presented so far are in line with the reasoning of Hellmann et al. who state that basin stability and survivability are generally not correlated and hence the asymptotic behaviour does not allow conclusions about the transient one [17]. We now want to achieve a more precise statement by putting a stronger focus on the patterns observed in the joint distribution of the stability measures in Figure 5.
In order to decipher how these patterns relate to topological characteristics, it turned out to be helpful to focus on nodes located inside or adjacent to tree-shaped parts of the network and to distinguish several types of nodes by how far inside the tree they are located (cf. Section 2.3 and Table 1). This is also suggested by the findings presented in Section 3.1.
Each of the previously defined classes features typical characteristics regarding asymptotic and transient stability at different perturbation levels ( Figure 5). We begin with discussing ∆ω = 5.0. In this case, for all nodes the survivability is lower than the basin stability. This indicates that there are no undesirable attractors within the survival region X + . Hence basin stability sets an upper limit to survivability. All trajectories which stay within the desirable region X + have to converge to the synchronous operating state of the grid, X [17]. This is expected as a desynchronization is expected to lead to oscillators rotating with their natural frequency, which here is ω LC = P/α = 10.0 (cf. (2.5)).
At rather low perturbation levels (∆ω = 5.0) the nodes forming the upper mode of the survivability distribution are the leaves. Dense sprouts show lower survivabilities than sparse sprouts, indicating that also a low neighbour degree is beneficial for survivability. The lower mode of the distribution is formed by bulk, root and inner tree nodes. As expected from the degree dependence (cf. Figure 4) inner nodes show higher survivabilities than root nodes which are most critical. Hence, the transient dynamics following a perturbation of a root node tend to exhibit large frequency deviations, leading to a transgression of the boundaries of the desirable region.
For ∆ω ≥ 7.5 we observe surprising new behaviour. While perturbations at most nodes still behave as if there was no undesirable attractor in the survival region, perturbations originating at a dense sprout almost all have σ > β. This means a novel asymptotic state or a very long transient, is reached, in which the system is not synchronized, but as the system does not leave the survival region X + , no node is swinging at its natural frequency ω LC either. Exemplary trajectories for perturbations at a dense sprout are shown in Figure C2 in Appendix C and reveal that we indeed observe a novel asymptotic state, with a solitary desynchronized oscillator not swinging at its natural frequency. To our knowledge, this state has not previously been observed in these systems, for example in the bifurcation studies of [43,44]. Indeed we expect that this state would be very hard to observe, if the initial conditions are drawn fully randomly, and not localized at individual nodes.
For sparse sprouts, proper leaves and inner tree nodes β and σ are strongly positively correlated. This means that for these nodes trajectories which leave the desired region tend to converge to non-synchronous states. In other words, the basin of attraction of the synchronized state is entirely contained in X + . Root nodes show a pattern contrary to dense sprouts with a wide range of survivabilities and rather high basin stabilities.
The patterns yielded at ∆ω = 10.0 are very similar to those at ∆ω = 7.5. As the frequency of the limit cycle fluctuates around 10.0, drawing the boundary exactly there means that still all nodes that fully desynchronize, and go to their natural frequency, will hit the survivability threshold. Now a few nodes besides the dense sprouts feature σ > β. A further anomaly that can be observed here is that the mean survivabilities of the leaves are smaller for ∆ω = 10.0 than for ∆ω = 7.5, opposed to the general observation of increasing survivabilities for larger values of ∆ω.
Finally, at ∆ω = 12.5, the natural frequency of the oscillators is fully within the survival region, and a desynchronized node is characterized as having survived. Consequently, the majority of nodes has larger survivabilities than basin stabilities. There are numerous non-synchronous states whose trajectories lie completely within the desirable region X + . The limit cycle trajectory is hence also prominent in the multinode network model. While the leaves of the tree-like parts show a similar pattern in the σ-β-space, the inner tree nodes are clearly separated at slightly lower survivabilities and fairly lower basin stabilities. Nodes from the bulk feature the largest basin stabilities.
Independent of the perturbation level ∆ω, dense sprouts tend to feature lower basin stabilities than sparse sprouts. Menck et al. studied the dependence of basin stability on the (average) degree of the neighbour(s) for all leaves ("dead ends") but did not find a significant correlation [29]. It is our finer partition of leaves into proper leaves and sprouts, the neccesity of which was recognized by studying the joint plots, that allows more detailed insights.
Complementary to Figure 5, the partition achieved using the defined topological classes of nodes is highlighted in three-dimensional representations of the joint distributions of basin stability and survivability at different perturbation levels ( Figure  C1 in Appendix C).

Relevance for designing stable power grids
The insights gained from this study are particularly relevant for both the stability assessment and the design of stable power grids. The employed nonlinear stability measures, basin stability and survivability, provide useful information on the asymptotic and transient stability of a grid against large perturbations, thereby surpassing the insights gained from local, linear stability assessments. The probabilistic measures have the benefits of being intuitive to understand and efficient to estimate for the overall system, irrespective of the complexity of the dynamics. They thus allow, for the first time, a systemic understanding of dynamical effects in the overall system.
Of particular relevance to the systemic stability of the power system are the novel asymptotic states we discovered, which primarily arise from perturbations at dense sprouts. The perturbed nodes are desynchronized and oscillate at a frequency much smaller than their natural. They are thus a pure network effect that can not be understood by studying individual machines. If such desynchronized states exist within sufficiently narrow frequency bounds, they would be a severe systemic risk to the power grid. They might be related to observed phenomena like Inter Area Oscillations, in which a deviation from perfect synchrony is observed. While an extensive analysis of these states is not part of this work, the topological characterization of dynamic systems already suggests a mean to prevent them preemptively in the design of the power grid, by avoiding dense sprouts.
Beyond this novel dynamical phenomenon, our study revealed various interrelations between the pure topological properties of a node and its stability. Tree-like structures, which contain about half of the nodes in the network ensemble, were found to be characterized by stability properties which significantly distinguish them from the remaining bulk of nodes. We were able to identify a small set of topological groups whose nodes feature similar patterns of asymptotic and transient stability at various perturbation levels. The knowledge about these patterns can in turn be utilized to predict a node's relevance for the whole grid's stability knowing only its topological embedding within the network. Given a concrete real-world power grid topology, our classification scheme hence enables the identification of potentially vulnerable nodes. While outside the scope of this work we expect that it will be possible to extend this identification to more nodes by means of statistic regressions taking network measures as input [42].
Menck et al. found that the weak basin stability of inner nodes can be cured by reconnecting leaf nodes (there termed "dead ends") to the grid [29]. This agrees well with our observation above to avoid sprouts and that "detour" nodes are favoured by our stability measures [42]. Another strategy suggested by our findings would be to avoid high-degree nodes as these feature the worst survivabilities, independent of the perturbation level. By combining the above-mentioned rules, reconnection of leaves and avoidance of hubs with high centralities, a possible design principle for stable power grids would be to strive for rather homogeneous network topologies, characterized by narrow degree distributions. Using linear stability techniques such homogeneous topologies have also be found to be generally easier to synchronize [45], while tree-shaped structures apparently show rather bad synchronizability [46]. Interestingly, while tree-shaped structures feature both poor linear (synchronizability) and nonlinear (basin stability, survivability) stability properties, in small-world topologies synchronizability and basin stability were found to behave contrarily [16].
It should be pointed out that for our analysis we employed conceptual models of the transmission grid. A direct transfer of the findings to lower grid levels might not be valid. Another simplifying assumption which does not hold for real-world grids is the uniform distribution of both the power injections and loads. Whether rather heterogeneous power distributions would affect the findings has not been investigated. We are also aware that the maximal perturbation levels assumed for our simulations are rather high. Instead of assuming uniformly distributed perturbations it would be more realistic to assume a unimodal distribution in which small perturbations in both the phase and frequency deviations are more frequent and large perturbations are rare events. Realistic distributions of perturbations could, for instance, be derived from the research on intermittency in power fluctuations from wind or solar power systems [47,48]. This work should thus be seen as an exploration of phase space properties like the structure of a state's basin of attraction, rather than a concrete study of realistic power grid models.

Relevance for complex dynamical networks in general
Beyond the direct implications for the resilience of power grids, our findings show the power of sampling-based dynamical methods, like survivability and basin stability. These are evaluated by sampling the system's dynamical reaction to non-small perturbations, here in phase space, but more generally also in parameter space. Combined with network measures and topological classifications, they provide general insights for complex dynamical networks. We expect that the crucial role of tree-shaped parts for the dynamics found here, is not specific to power grids but might rather be a general phenomenon in oscillator networks, infrastructure networks or even biological networks, e.g., in neural or genomic dynamics. The suggested topological classification scheme and the terminology for nodes in trees are independent of the nodal dynamics and are thus easily applicable to networks of different types.
As the employed synchronous machine model (or "Swing Equation") is equivalent to the Kuramoto model with inertia [28], the findings are particularly relevant for the general study of oscillator networks. The combined analysis of transient and asymptotic behaviour via basin stability and survivability allows indirect insights into the geometry of the system's phase space. The identification of nodes for which σ > β allowed us to infer the existence and position of new types of limit cycles. Particularly, for the group of dense sprouts we found that non-synchronous states exist besides those given by the approximated solution to the single-node system (2.5).
In a first order approximation the coupling strength K in (2.5) is proportional to a node's degree d, while P and α are independent of topological characteristics. Hence the approximation (2.5) suggests that the amplitudes A LC = P K/α of a limit cycle are proportional to the degree of the respective non-synchronously rotating nodes. This is a contributing mechanism to the low survivability of high degree nodes for very large perturbations. We suspect that it will be possible to gain further understanding of the relationship between topology and dynamics through more sophisticated approximations. Understanding the dependencies of the asymptotic spectrum of networked dynamical systems on the ambient topology in more detail will however require significant further work, both numerical and analytical.

Conclusions
Our results form another step towards a better understanding of the interrelations between topology and stability in complex dynamical networks. Tree-shaped topologies which are particularly prominent in infrastructure networks, have been found to feature stability properties which considerably deviate from those of the remaining bulk of nodes. A topological classification scheme for nodes adjacent to those tree-shaped parts has been suggested which enables a prediction of a node's transient and asymptotic stability against large perturbations. This classification of nodes can hence aid both the stability assessment and the design of stable infrastructure systems.
The sampling based stability measures we employed were shown to enable surprising novel insights into the asymptotic dynamics of networked dynamical systems, revealing both, previously unknown asymptotic states and surprisingly precise relationships between the topology and these novel states.
Due to the parametrization of the model equations, the results are particularly relevant in the context of power grid research. If both high asymptotic stability (reflected by single-node basin stability) and transient stability (reflected by survivability) of power grids are desired, avoiding both sparsely connected tree-shaped structures and highdegree hub nodes appears to be a promising design principle.
Independent of the particular application the presented study shows how the nonlinear stability concepts of basin stability and survivability can be combined to gain a better understanding of a -not necessarily networked -dynamical system.
Appendix A. Algorithm for the identification and classification of nodes in tree-shaped parts of networks Let G = (V, E) be an undirected graph that is connected and not itself a tree (i.e., contains at least one cycle). Our goal is to identify all nodes x ∈ V that are in a treeshaped part of G and classify them using height and depth. As defined in the main text, a tree-shaped part of G is an induced subgraph T = (V , E ) = G| V of G (i.e., a subset of nodes V ⊆ V together with the set E of all edges in E between nodes in V ) that is maximal with the property that there is exactly one node r ∈ V that has at least one neighbour in G − T . r is then called the root of T , and one can see easily that it must have degree at least three. For any graph G = (V , E ) and node x ∈ V , we denote by A simple algorithm to identify all tree-shaped parts of G, their roots, and the parents, children, heigths, depths, and branches of all their members is the following.
In the first part, we iteratively define • a decreasing sequence of node sets V 0 ⊃ V 1 ⊃ V 2 . . ., • the respective induced subgraphs • a sequence of disjoint height level sets H i , • parents π(x), • sets of children C(x), • branches B(x), • and height labellings η(x), by successively removing leaves from the remaining graph as follows. Put V 0 := V and initially C(x) := ∅ for all x ∈ V . Given V i and G i : 1} be the set of leaves of G i . For each x ∈ H i , let the parent of x, π(x), be the unique neighbour of x in G i ; add x to its set of children, C(π(x)). Note that π(x) ∈ V i − H i . The branch of x is B(x) := {x} ∪ y∈C(x) B(y), and the height is η(x) = i. As long as To finishing the first part after these iterations, let N := i H i be the set of all thus identified non-root nodes, let R := {π(x) : x ∈ N } − N , and call each r ∈ R a root. Put B(r) := {r} ∪ y∈C(r) B(y) and η(r) := 1 + max{η(x) : x ∈ N, π(x) = r} for all r ∈ R. The tree-shaped parts T of G are now exactly the subgraphs T = G| B(r) induced by the branches of any roots r ∈ R.
In the second part, we define a depth δ(x) for each x ∈ N ∪ R, counted outwards starting from the roots, in addition to the height, which is counted inwards starting from the leaves. This is again done iteratively by defining a sequence of disjoint depth level sets D i . Put D 0 := W 0 := R, and put η(x) := 0 for each x ∈ D 0 . Having defined D i−1 and W i−1 , define D i := x∈D i−1 C(x) − W i−1 and W i := D i−1 ∪ D i , and put δ(x) := i for each x ∈ D i , iterating this until D i = ∅. Note that δ(x) is the distance from x to the root of its tree-shaped part.

Appendix B. Estimation of basin stability and survivability from simulations
Both employed stability measures, basin stability β (2.9) and survivability σ (2.12), have been estimated using Monte-Carlo sampling. For each of the M × N = 5000 nodes L = 200 trajectories with perturbed initial conditions (2.6) and (2.7) have been simulated.
If s of these trajectories return (sufficiently close) to X after the simulation time of t = 100, their fraction is used as an estimator of β: Since the perturbed trajectories either converge or not, the sampling of initial conditions can be regarded as a Bernoulli experiment. Thus the standard error of the probability estimator is given by If either all or none of the perturbed trajectories converge the estimated standard error becomes zero which is troublesome for further statistical calculations. Hence we used the Agresti-Croull method for a better estimation of β and its standard error e β [49]. For the desired confidence this corresponds to adding one trial which is "half success" and "half failure". Definings = s + 1 2 andL = L + 1, the corrected estimatorβ is given byβ Analogously, σ and its standard error e σ have been estimated using the fraction of trajectories which did not leave the desirable region X + given by (2.13) within the simulation time.  ω(t) Figure C2. Exemplary trajectories of the frequency deviations ω(t) for perturbations at a dense sprout which are representative for the different observed asymptotic states. The trajectory of the perturbed node is coloured blue while those of all other nodes in the network are shown in grey. For low perturbations (upper panel) the frequency deviation of the perturbed node declines over a long transient phase until the network settles to the synchronous state after about 100 seconds. At medium perturbations levels (lower left panel) we observe the novel asymptotic state for which the perturbed node oscillates around a frequency deviation of about half of its natural frequency, while relatively strong oscillations around ω = 0 are observed in the remaining network. Finally, for rather high perturbations (lower right panel) the dense sprout asymptotically oscillates around its natural frequency ω LC = P/α = 10.0, as suggested by the approximation 2.5, while the remaining network shows oscillations around ω = 0 but with a smaller amplitude than in the previous case. Note that the impression of accelerating oscillations is due to the logarithmic time axis.