Targeted suppression of failure spreading in multistable oscillator networks

Fluctuations and damages crucially determine the operation and stability of networked systems across disciplines, from electrical powergrids, to vascular networks or neuronal networks. Local changes in the underlying dynamics may affect the whole network and, in the worst case, cause a total collapse of the system through a cascading failure. It has been demonstrated that certain subgraphs can reduce failure spreading drastically, or even inhibit it completely. However, this shielding effect is poorly understood for non-linear dynamical models. Here, we study the effect of perturbations in networks of oscillators coupled via the Kuramoto model. We demonstrate how the network structure can be optimised for suppressing specific, targeted fluctuations at a desired operational state while letting others pass. We illustrate our approach by demonstrating that a significant reduction in time-dependent fluctuations may be achieved by optimising the edge weights. Finally, we demonstrate how to apply the developed method to real-world supply networks such as power grids. Our findings reveal that a targeted shielding of specific solutions in multistable systems is possible which may be applied to make supply networks more robust.


Introduction
Complex supply networks in nature are subject to numerous external stimuli and perturbations. Whereas many complex systems are surprisingly resilient to most random failures [1], single, targeted attacks or failures of critical elements may trigger cascades of failures that lead to an entire collapse of the network [2][3][4]. Counteracting strategies against these failures are heavily sought after due to the increasing dependence of modern society on infrastructural networks such as power grids and are an important part of complex networks research [5][6][7][8][9]. Whereas traditional approaches rely on weakening connectivity which has been demonstrated to successfully limit the spreading of perturbations [10][11][12], recent work has analysed counteracting strategies that are based on an increase in connectivity [6,13,14]. Remarkably, it has been demonstrated that certain subgraphs-termed 'network isolators'-can completely suppress failure spreading at any degree of connectivity in linear flow networks [6]. Nevertheless, many supply networks are in fact described by non-linear relationships and evolve dynamically. The main question addressed here is the following: to what extend can we make use of network isolators in non-linear systems?
To model non-linearly coupled networks, the Kuramoto model stands out as the prime model. It has found numerous applications [15,16]-ranging from powergrids [17,18] to biological systems [19][20][21] and any dynamical system of weakly coupled oscillators near a limit-cycle [22,23]. Thus, a targeted suppression of failure spreading in the Kuramoto model may find immediate application throughout disciplines and models.
However, non-linearly coupled oscillator systems such as the Kuramoto model introduce an additional complexity to the problem: they display multiple stable fixed points [24][25][26][27]. This introduces an additional challenge to control strategies, since they should aim at not only suppressing cascade propagation, but also keeping the system in the desired state of operation, i.e. control the multistability [28][29][30][31]. Here, we demonstrate that both goals-suppressing failure spreading and controlling multistability-can be achieved at the same time in the Kuramoto models.
In this manuscript we focus on network isolators to shield parts of the network against perturbations and fluctuations in non-linearly coupled oscillator networks. To this end, we first discuss the extension of network isolators from linear flow networks to non-linearly coupled oscillators in the Kuramoto model. We then demonstrate a possible application of network isolators in the Kuramoto model: shielding targeted, desirable fixed points while allowing failure propagation for others. Finally, we illustrate how to suppress time-dependent fluctuations in oscillatory networks by means of network isolators. Our results allow to significantly reduce the impact of perturbations at a desired fixed point of a multistable system and may be very useful for making real-world networks such as power grids more resilient to failures.

Theoretical background: network isolators in the Kuramoto model 2.1. Kuramoto dynamics on networks
Many natural systems such as fireflies or synapses in the brain, and many man-made systems such as the electrical power grid exhibit collective synchronous behaviour [17,32,33]. Commonly, this behaviour is described using the Kuramoto model which models the time evolution of N non-linearly coupled limit cycle oscillators [34]. It is defined by the following set of differential equations: Here, ϑ i ∈ R is the phase angle of oscillator i, ω i ∈ R is its natural frequency, K ij ∈ R denotes the interaction strength of oscillators i and j andθ i is the differentiation of the phase angles with respect to time, i.e. the angular frequency at node i. The population of coupled oscillator may be described in graph theoretical terms by assigning to each oscillator i a vertex i ∈ V of a graph G = (E, V) and the mutual coupling between two oscillators i and j to an edge e = (i, j) ∈ E of that graph, with a total number of edges L = |E|. Here, we consider graphs that are weighted, undirected and connected and do not contain self loops. We label all edges in the graph by an index l = 1, . . . , L and define the graph's coupling strength matrix as the diagonal matrix K ∈ R L×L K = diag(K 1 , K 2 , . . . , K L ), where we use the notation K l of the index representing edge = (m, n) and K mn interchangeably. In this paper we assume that this coupling is positive, K ij > 0, ∀(i, j) ∈ E, i.e. the oscillators attract each others. We further introduce the graph's edge-vertex incidence matrix I ∈ R N×L . For this purpose, we arbitrarily fix an orientation of each edge and say that the edge e = (i, j) is oriented from node i to node j. Thus, I is defined as The edge incidence matrix I as well as the coupling matrix K ∈ R L×L may be used to express the Kuramoto model (1) in a compact, vectorised formθ = ω − IK sin(I ϑ), (2) with ϑ = (ϑ 1 , ϑ 2 , . . . , ϑ N ) ∈ R N and where the sine is taken elementwise. The graph theoretical approach allows to intuitively connect the dynamics to networked objects such as power grids. When approximating the dynamics of power grids using the Kuramoto model, vertex i represents a synchronous machine or generator with voltage phase angle ϑ i [15,17,[35][36][37]. In case of the synchronous motor model (SM model) [17], both, consumers and generators are modelled as second-order Kuramoto oscillators. Each oscillator is assigned an inertia constant M i , i ∈ V and a damping coefficient D i , i ∈ V, that both depend on the rotating generator or motor attached to node i. The network dynamics of a lossless power grid consisting of synchronous motors and generators is then compactly summarised as Here, the diagonal matrices M = diag (k 1 , . . . , M N ) ∈ R N×N and D = diag (D 1 , . . . , D N ) ∈ R N×N summarise the inertia and damping, respectively, of each synchronous machine. The natural frequencies Figure 1. Multistability in the Kuramoto occurs in simple networks. Dynamically stable states in the Kuramoto model for the cycle graph of six nodes, C 6 , may be characterised uniquely by their winding numbers. We show the phase angles for the phase synchronised fixed point with ϑ i = π, ∀i ∈ V(G), i.e. winding number q = 0 (a), and the two other stable solutions given by (8)). ω = (ω 1 , . . . , ω N ) may be interpreted as the power each oscillator produces or consumes depending on the sign [27]. In the context of power grids, the vector denotes the flow of real power along the edges in the graph in the lossless line approximation [17,18,35] such that the flow along an edge e = (i, j) is given by

Fixed points and cycle-flow induced multistability in the Kuramoto model
Consider the Kuramoto model (2) on a network and assume that the system is in a stable fixed point which is characterised by a vanishing time derivative,θ i = 0, such that equation (1) reduces to where ϑ * i denotes the phase angles at the fixed point. We assume that the sum over all natural frequencies ω i is zero, i ω i = 0. Otherwise the system does not relax to a stable fixed point with vanishing angular frequencies, ϑ i = 0, but is constantly drifting, since the coupling in the Kuramoto model is diffusive such that From this steady state condition equation (6) we can observe straightforwardly that a possible fixed point in the case where the natural frequencies are zero ω i = 0 is the fully synchronised state where all phase angles are equal ϑ i = c, c ∈ R. There are, however, more complex fixed points in the Kuramoto model as demonstrated in recent publications [24][25][26][27]. More complex fixed points may occur due to a combination of cycle flows on different cycles of the underlying graph [26,27] or due to phase lags in the coupling [24,25]. Here, we do not consider phase lags and thus focus on the former path to multistability. Furthermore, we focus on fixed points where angular differences between neighbouring vertices are small such that which is known as normal operation and guarantees the fixed points to be linearly stable [27]. As a simple example, consider the cycle graph C 6 consisting of N = 6 nodes V = {1, 2, 3, 4, 5, 6} that are connected in a circle via the edges E = {(1, 2), (2,3), (3,4), (4,5), (5,6), (6, 1)}. The graph is illustrated in figure 1. For simplicity, assume further that all natural frequencies vanish, i.e. ω i = 0, ∀i ∈ V, and all couplings strengths are unity, i.e. K ij = 1, ∀(i, j) ∈ E. In this case, equation (6) characterising the fixed points reduces to Set ϑ * 1 = 0 without loss of generality. Due to the periodicity of the sine function, this gives rise to the set of solutions to the above equation given by for any node i ∈ {1, 2, 3, 4, 5, 6} of which the solutions characterised by m * ∈ {0, 1, 5} can be identified as stable solutions and the remaining ones as unstable [38]. Thus, multistability occurs due to the non-linearity already in this simple network. By virtue of equation (4), the phase differences induce flows which are referred to as cycle flows [27] and do not originate from a source or sink as ω i = 0, ∀i ∈ V. As a result, different fixed points may be characterised by the resulting cycle flows in the network (see figure 1).
We will now illustrate how to generalise this reasoning to networks of arbitrary size. The set of all cycles in a graph defines a vector space known as the graph's cycle space. The basis of this space is known as the cycle basis B C which has the dimensionality L − N + 1 [39]. Thus, consider a graph G whose cycle basis is given by . . , C L−N+1 } and assume that the edges in each cycle C i = {e i1 , e i2 , . . .} are oriented in the same direction as the edges in the graph. Then the cycle flow along each cycle C i can be quantified by the winding number The winding number q i can take integer values in q i ∈ {− n i /4 , . . . , n i /4 } ⊂ Z for each cycle C i in the network, where n i is the number of nodes in the cycle, such that a fixed point in a planar network can be fully characterised by a unique winding vector q = (q 1 , . . . , q L−N+1 ) [26,27] which collects the winding number of each individual cycle in the network.

Linear flow networks: linear approximation of flows in the Kuramoto model
In many applications, the relevant-and most stable [40]-fixed points in the Kuramoto dynamics turn out to be characterised by small phase differences between neighbouring nodes. For example, for power grids at the transmission level, phase differences along transmission lines are typically small, i.e. |ϑ i − ϑ j | 1, ∀e = (i, j), if the power grid is not very heavily loaded [41]. Thus, the sine function in equation (5) can be approximated by its argument by using its Taylor series resulting in the linear flow which is known as the DC approximation of power flow [42]. If we further assume that the nodal phase angles do not vary over time, i.e. the system is at a fixed point such thatθ i =θ i = 0, ∀ i ∈ V, equation (3) reduces to the following set of equations, known as Kirchhoff's current law Together, this set of equations fully determines the nodal potential ϑ i and the flows F ij in the steady state up to a constant phase shift ϑ i + c, c ∈ R, once the natural frequencies ω i and coupling constants K ij are known for all nodes and edges, respectively. This is in contrast to the Kuramoto model, where multiple solutions to the corresponding set of equations exist [27]. Inserting the linear flow equation (9) and rewriting the equation in vectorised notation yields ω = Lϑ, (11) where we defined the graph's Laplacian matrix L = IKI ∈ R N×N . Componentwise we may write the Laplacian in the following way This matrix plays a crucial role in network science since it fully characterises the network topology and its eigenvalues characterise the systems stability [15].
In the following sections, we will also make use of another matrix which allows to characterise the connectivity structure of a graph. We define the graph's weighted adjacency matrix A ∈ R N×N as [43] which thus fully describes all neighbourhood relations of the graph.

A linear response theory approach to structural perturbations
Real-world dynamical networks are subject to external influences which may induce damages to edges or nodal perturbations. This causes phase shifts of the affected nodes to which adjacent nodes respond with phase shifts and so forth. The impact of structural disturbances or nodal perturbations may be analysed within the framework of linear response theory [44]. Before we proceed to illustrate the impact of perturbations in detail, we note that structural perturbations, i.e. modifications of edge weights, result in a similar set of equations.
We here restrict ourselves to perturbations of the natural frequencies which leave the sum i ω i -the power balance in the case of electric power grids-invariant. As a prototypical example we consider dipolelike frequency perturbations here, where the frequency at one node i is increased by a value of Δω and the frequency at another node j is decreased by the same value (see appendix A). Assume that after the system relaxed to a fixed point ϑ * as described in equation (6), we apply a dipole-like perturbation Δω to nodes i and j such that the system relaxes to a new fixed point ϑ , Thus the steady state equation for the system after we applied the perturbation reads where δ in is the Kronecker symbol. Expanding this to leading order with respect to the new fixed point and subtracting the steady state equation of the original fixed point, equation (6), leads to (see appendix A) which describes how the system evolves from a fixed point ϑ * into a different fixed point ϑ due to a small perturbation Δω i . The sum in the expression may be used to define an effective LaplacianL with elements Since we focus on fixed points in normal operation, all effective edge weights K ij cos(ϑ * j − ϑ * i ) are positive and the above matrix is Laplacian. Using this effective Laplacian matrix, we can write the above equation in vectorial form asL where e i ∈ {0, 1} N is a vector with entry one at position i and zero otherwise. Since the effective Laplacian matrix depends directly on the initial fixed point ϑ * , we will also denote it byL(q) in the following, where q is the winding vector characterising the initial fixed point to include this dependence. We are thus left with a discrete Poisson equation with a dipole source term that has been analysed thoroughly in the context of linear flow networks [6,42,45]. In particular, several rigorous results on failure spreading were derived for this type of equation and may be transferred to the non-linear Kuramoto model as we will see in the next subsection. Importantly, the Laplacian is not purely topological anymore, but depends on the fixed point ϑ * in which the system was before the perturbation. Based on this, we can calculate the change in real power flows along the lines by making use of equation (4) as Figure 2. Illustration of node labelling in two subgraphs connected by a network isolator (thick black lines). The subgraphs G(V 1 , E 1 ) (red nodes) and G(V 2 , E 2 ) (blue nodes) characterised by their vertex sets V 1 and V 2 , respectively, of magnitudes |V 1 | = n 1 and |V 2 | = n 2 are connected via subsets of their nodes (light red nodes, light blue nodes) with magnitudes k 1 and k 2 , respectively. As discussed in section 2.5, we label the vertices that are connected to the other subgraph and contained in V 1 from 1 to k 1 and the vertices contained in V 2 from n 1 + 1 to k 2 .

Network isolators in linear systems
In principle, perturbations at single nodes will spread throughout the whole network via the coupling between the oscillators. However, it has recently been demonstrated that certain subgraphs, termed 'network isolators', inhibit failure spreading in linear flow networks completely [6,13]. Before we proceed, we fix some notation. We consider a graph G consisting of two subgraphs, i.e. its vertex set V is written as V = V 1 ∪ V 2 . We now label the nodes in V as follows without loss of generality Then we can write the weighted adjacency matrix A of the network as where A 12 is the submatrix of the adjacency matrix representing the mutual connections between the two subgraphs G 1 and G 2 induced by the vertex set V 1 and V 2 , respectively. For this setup, the following theorem characterises network isolators in the context of linear flow networks [6]. (9) and (10). Assume a dipole-like perturbation at two nodes in the induced subgraph G(V 1 ), r, s ∈ V 1 , or equivalently the failure of a (r, s) with r, s ∈ V 1 . The failure induces a change in the nodal potentials Δϑ which can be calculated as [6,42]

Theorem 1 (Network isolators in linear flow networks). Consider a linear flow network consisting of two modules with vertex sets V 1 and V 2 characterised by equations
where L is the Laplacian matrix of the network, η rs is a source term and e r is a vector with entry one at position r and zero otherwise. If the subgraph of the adjacency matrix A 12 of the mutual connections (cf equation (19)) has unit rank, i.e. rank(A 12 ) = 1, then the change in nodal potentials is given by where Δϑ 1 ∈ R n 1 is the change in potentials in the first induced subgraph G(V 1 ) and Δϑ 2 = c1 2 ∈ R n 2 the change in potentials in the other induced subgraph G(V 2 ). As a consequence, the flows on all links in the induced subgraph G(V 2 ) are not affected by the failure, that is their flow changes ΔF vanish where E 2 is the edge set of the induced subgraph G(V 2 ). The subgraph G corresponding to the mutual connections is referred to as network isolator.
A proof can be found in reference [6].The theorem of network isolators thus claims that, if the adjacency matrix of mutual connections has unit rank, i.e. rank(A 12 ) = 1, then the flow on all edges in G(V 2 ) stays constant upon a link failure in G(V 1 ). Strikingly, this implies that subgraph G(V 2 ) is completely shielded from perturbations in subgraph G(V 1 ) and vice versa. To quantify to which extend a graph fulfils theorem 1 on network isolators in linear flow networks, a measure based on it coherence statistics has been suggested in reference [6] that reads as follows Here, a i are the rows-or columns-of the adjacency matrix A 12 and , is the standard scalar product. Note that if all rows-or columns-are parallel vectors which indicates that the matrix rank is unity, rank (A 12 ) = 1, the latter expression approaches unity such that ξ(A 12 ) = 0 in this case. Thus, the coherence statistics quantifies to what degree the subgraph under consideration deviates from being a perfect network isolator.

Network isolators for targeted suppression of failure spreading in the Kuramoto model
Theorem 1 is only valid for linear flow networks, but many systems such as the Kuramoto model are in fact non-linear. In this section, we expand this approach in two ways: (1) we study the applicability to non-linear systems with the Kuramoto model serving as the prime example. (2) We study to what extent network isolators prevent failure spreading or the propagation of disturbances also for perturbations that fluctuate dynamically in time. In particular, we demonstrate that network isolators can be used to greatly reduce failure spreading in trivial fixed points. Furthermore, we illustrate how to tune the coupling weights in the network isolator in such a way that it shields specific fixed points.

Network isolators in the Kuramoto model
As a first step, we now illustrate how to construct network isolators in the Kuramoto model. As we demonstrate in figure 2, part of the results obtained for linear flow networks may be translated directly: we analyse the potential of a network isolator to shield perturbations in a Kuramoto network and show that it strongly decreases the flow changes |ΔF| (cf equation (18)) induced by frequency perturbations. To this end, we consider a small network which consists of two hexagonal subgraphs that are connected by three edges (cf panel (a)). We then add a single link to build a network isolator as characterised by the condition rank(A 12 ) = 1 (panel (b)) and compare the effect of a dipole-like perturbation Δω (black arrows, cf equation (A.2)) in both networks that are initially at a stable fixed point: building a network isolator strongly decreases the flow changes |ΔF| as visible in the colour code. Next, we analyse the mean of the flow changes |ΔF| with respect to the edge-distance to the perturbed nodes (panel (c)) and find that that it decreases rapidly beyond the network isolator. We then analyse the impact of the perturbation strength |Δω i | in a quantitative manner (panel (d)): we compute the mean of all absolute flow changes in the passive subgraph ΔF passive , i.e. the subgraph shielded by the network isolator, at different perturbation strengths. We find that the isolator effect is lost at a certain threshold which is where the fixed point ϑ loses its stability. Now consider again a Kuramoto network consisting of two subgraphs connected via a network isolator G. Assume once again that the nodes of the network are sorted in the way discussed in section 2.5 and pictured in figure 3. In this case, the matrixÃ 12 (cf equation (19)) in the upper right of the effective Laplacian matrixL (equation (16)) is given byÃ 12 Here,ã ∈ R k 1 ×k 2 is the adjacency matrix of connections in the isolator graph G and has the entries a ij = K (i,n 1 +j) cos(ϑ * i − ϑ * n 1 +j ), such that-due to the labelling of the nodes-it contains the cosine of angular differences measured along the connections between the two subgraphs. Thus, we can extend network isolators to the Kuramoto model in the following way: The effective adjacency matrixã does not in general have unit rank, even if the original subgraph G is a network isolator in the linear sense, i.e. rank(a 12 ) = 1. This is due to the fact that it explicitly dependents on the original fixed pointã =ã(q), where q is the winding number characterising the fixed point ϑ * . However, we may adjust the weights in the network isolator K ij , (i, j) ∈ E(G) in an iterative manner to achieve rank(ã) = 1 using the following algorithm:

Algorithm 1 (Network isolators in the Kuramoto model)
(a) Tune the edge weights in the network isolator such that K [n] ij is replaced by (b) Relax the system to a new fixed point ϑ [n] → ϑ [n+1] .
(c) Compare the updated to the old coupling matrix. If stop.
Here, ε > 0 is a predefined threshold value and K [n] denotes the edge weights at the nth iteration of the algorithm. In appendix B, we illustrate that this procedure converges in most cases, and-as a result-we may build a network isolator in the linear response regime characterised by rank(ã) = 1. Thus, dipole-like perturbations are then shielded in the Kuramoto model by virtue of theorem 1. In appendix C, we give more details on the theoretical background. Next, we will discuss an important consequence of this finding.

Targeted shielding of fixed points to static perturbations in the Kuramoto model
As discussed above, the Kuramoto model admits multiple stable fixed point solutions. However, in many applications, one fixed point may be the desired state of operation for the system such that it is supposed to be shielded against fluctuations whereas this is not the case for other fixed points [29]. This targeted shielding may be achieved for the Kuramoto model by making use of network isolators as we discuss in the following paragraph.
In section 3.1, we demonstrated how to extend network isolators to the Kuramoto model. Here, we discuss an important implication of this result: network isolators in the Kuramoto model allow to shield targeted fixed points. Thus assume that the network of Kuramoto oscillators admits at least two different stable fixed point solutions ϑ (0) , ϑ (1) ∈ R N as defined by equation (6) that are characterised by two different winding vectors q (0) = q (1) (cf equation (8)). Using the algorithm outlined in the last section, we may then tune the edge weights K ij , (i, j) ∈ E(G) in the network isolator such that rank (ã(q (0) ) = 1, i.e. failure spreading is suppressed in the fixed point ϑ (0) . In most cases, this suppression of failure spreading in one fixed point implies that there is no perfect network isolation in the other fixed point, i.e. rank (ã(q (1) ) = 1, such that we observe a targeted suppression of failure spreading only in the first fixed point. Note that in order for this statement to hold, we need to exclude trivial network isolators, e.g. consisting only of a single line, where rank(ã) = 1 is trivially fulfiled for any fixed point, such that a targeted shielding is no longer possible. In appendix C, we discuss the theoretical background in detail.
We illustrate this procedure in figure 4 where we tune the edge weights such that a specific fixed point is shielded. We consider two different fixed points ϑ (0) (panel (a)) and ϑ (1) (panel (b)) characterised by their winding vectors q (0) = (0, 0) and q (1) = (1, −1) , respectively, on a Kuramoto network consisting of two hexagonal subgraphs that are connected by a network isolator. Here, we focus on the winding numbers of the two large facets to the left and the right of the network and assume all other facets to have a zero winding number. We then analyse how both fixed points react to frequency perturbations in the initial setup where all edges have uniform weights K = diag(1, 1, . . . , 1): we apply a dipole-like perturbation to the nodes marked by black arrows ((c) and (d)). Consequently, the isolator reduces the flow changes |ΔF| in the passive subgraph, but the effect is substantially more pronounced for the first fixed point ϑ (0) than for the other one ϑ (1) . To specifically shield the second fixed point we then tune the edge weights in the network isolator according to algorithm 1. Consequently, the first fixed point ϑ (0) is no longer shielded completely to the perturbation (panel (f)), whereas the shielding is significantly improved for the other one ϑ (1) (panel (g)). To confirm the quality of this result, we compare the flow changes with respect to all possible dipole perturbations in the left subgraph for both, the initial weights (panel (e)) and the modified weights (panel (h)). Thus, we conclude that network isolators can shield flow changes due to dipole perturbations in the Kuramoto model: by tuning the edge weights according to algorithm 1 we can shield a specific fixed point while at the same time decreasing the effect for other fixed points. A caveat of this technique is that due to the modification of the weights the fixed point may change, therefore one has to use an iterative approach to find the correct compensation weights. We discuss this iterative approach in detail in appendix B.

Targeted suppression of dynamical perturbations in the Kuramoto model
Whereas static perturbations are adequate to describe permanent changes to a network, e.g. by modifying its topology, other perturbations are dynamically fluctuating over time. This is for example the case in modern power grids were an increasing penetration of renewable energy sources leads to a fluctuation in the generation of electrical power. However, as we will demonstrate here, the effect of network isolators in reducing the impact of perturbations is not limited to static perturbation, but they can also be used to greatly reduce the impact of fluctuations.  ) and (e) We monitor the angular frequenciesθ i for both topologies for the driven nodes (top), a node located in the same subgraph as the driven nodes (blue, centre) and another node at the same distance that is located in the other subgraph (yellow, bottom) and show both, the signal and a low-pass filtered signal (thick lines). Introducing a network isolator significantly reduces the fluctuations at the shielded node. (c) and (f) The control effort per node P i (equation (23)) calculated over the whole period of time shown here is reduced drastically in the subgraph shielded by the network isolator. Note that both graphs shown here are six-regular and taken from [6].
To illustrate the reduction in fluctuations, we make use of a measure that is commonly used in power grid research to quantify the degree of angular frequency excursions that represents the primary control effort necessary to keep the system at its nominal frequency [46][47][48]. This control effort for a single node i ∈ V can be obtained as where is the averaged angular frequency at time t weighted by the damping constants D i (cf equation (3)). The overall primary control effort P(T) can then be obtained by summing over this control effort per node Note that this performance metric may also be interpreted as a measure of stability since it quantifies the degree of frequency excursions that might result in the network ending up in a different fixed point. It may thus also be interpreted as weighted mean squared frequency deviation. In figure 5 we analyse the impact that a network isolator has on this metric for the subgraph of the network shielded by the isolator. To this end, we consider the second order Kuramoto equation (3) subject to time dependent fluctuations [50,51] Mθ + Dθ = ω − IK sin(I ϑ) + F(t). (25) Here, F(t) is a vector of time dependent-fluctuations that is generated by a Wiener process with zero mean and unit variance. We first consider a six-regular network that does not contain a network isolator (a) on which we analyse the second-order Kuramoto model (3) (see section appendix B for details on the simulation). After the system has relaxed to a stable fixed point, we apply the time-dependent, dipole-like driving F(t), i.e. opposite driving at both nodes (shaded brown). In panel (b), we monitor the response at the driven node (top, brown), a node in the same subgraph as the driven node (blue, centre) and a node in the other subgraph (yellow, bottom). Finally, in panel (c), we colour-code the control effort per node. We find that the control effort is slightly smaller in the passive module of the grid, but differences are limited. After rewiring two links to build a network isolator (d), fluctuations are significantly reduced in the subgraph shielded by the isolator ((e), bottom) and the control effort in this subgraph is reduced (f). Thus, we may conclude that isolators reduce fluctuations and therefore the control effort. We monitor the mean flow changes |ΔF| due to a perturbation (arrow, (c)) as they decay with distance to the perturbation. Tuning the edge weights in a network isolator as described in algorithm 1 allows to significantly reduce the flow changes (dark grey, colour code in (e)) even compared to the case of a non-tuned isolator (grey, colour code in (d)) and the case without an isolator (light grey, colour code in (c)). Dashed lines represent the perturbed (right) subgraph, solid lines represent the passive (left) subgraph and shaded regions display the variance |ΔF| at a given distance.

Case study: IEEE RTS-96 power grid
In this section, we discuss potential applications of our findings to increase the resilience of real-world power grids. As a case study, we focus on the three-area IEEE reliability test system 1996 (RTS-96) in our analysis [49]. This test system has been used frequently in recent publications for stability assessment [51][52][53].
We first analyse the impact of frequency perturbations (equation (14)) in figure 6. The initial setup of the RTS-96 grid is shown in panel (a) where we indicate the generators (triangles pointing up, ω i > 0), loads (triangles pointing down, ω i < 0) and neutral nodes (circles, ω i = 0) of the power grid (see appendix B for details). Both perturbed nodes (panel (c), black arrow) are generators of power, such that a dipole-like shift in the natural frequencies, Δω, applied to both nodes corresponds to reducing the generated power at one of the nodes and having the other generator compensate for the reduction, i.e. a power transfer between the generators. For this setup, we then analyse the absolute changes in real power flow (equation (18)) as a result of the transfer in generation: we demonstrate that flow changes |ΔF| can be significantly reduced by building a network isolator according to theorem 1 (panel (d)), and decreased even further by tuning the isolator's edge weights (panel (e)). Even without tuning the edge weights, flow changes in the passive module of the grid are reduced significantly when comparing to the initial setup without a network isolator (panel (c)). This finding is summarised in panel (b) where we show the mean flow changes |ΔF| against the edge distance for the three scenarios and illustrate the reduction of flow changes in the passive subgraph with respect to the respective isolator-type.
As a next step, we consider how the introduction of a network isolator changes the control effort in the grid due to a stochastically fluctuating generation source. To this end, we add a fluctuating power generation time series F(t) to a single generator node and simultaneously subtract the same time series from another generator to make sure the generation is balanced at any point in time such that N i=1 ω i = 0. The stochastic time series F(t) is generated from the procedure outlined in reference [51], such that it reproduces the main features of a wind power generation time series, such as the Kolmogorov exponent of its power spectrum and the non-Gaussianity of the increment statistics. Adding the time series F(t) can thus be understood as adding a wind farm to the respective node and having a device, e.g. a battery storage, at another node taking up the fluctuating power.
In figure 7, we illustrate how the RTS-96 grid reacts to the introduction of such a fluctuating energy source for different scenarios. In particular, we make use of the control effort per node P i (equation (23)) to quantify Figure 7. Network isolator reduces control effort due to stochastically fluctuating power feed-in F(t) in the power system test case 'RTS-96' [49]. (a) We consider a network of second order Kuramoto oscillators on the power grid 'RTS-96' that is initially at a stable fixed point and add a stochastic, time-dependent signal F(t) (blue) to two nodes (black arrows, cf equation (25)). We then monitor the control effort per node P(T) (colourcode, cf equation (23)) with and without a network isolator with weights adjusted to the fixed point (see figure 6). To model the dynamics of the generators and loads in the power grid, we consider four different scenarios characterised by heterogeneous coupling constants K ε (panel (b)) and homogeneous coupling constants K 0 (panel (c)) and heterogeneous inertia constants M ε (left boxplots) and homogeneous inertia constants M 0 (right boxplots). For each scenario, we quantify the reduction of the control effort achieved by the network isolator by dividing the control effort of every node with the isolator P w iso by the control effort without the isolator P no iso . We notice that in all scenarios, a significant reduction is achieved in the passive subgraph shielded from the fluctuations by the network isolator with the effect being strongest for homogeneous coupling constants. The scenario colourcoded in panel (a) is indicated by a darker colouring of the boxes.
to what extend a network isolator that is adjusted to an initial, operational fixed point ϑ is able to suppress fluctuations from this fixed point. We analyse four different scenarios (panels (b) and (c)) in which we vary whether coupling strengths K ij are homogeneous (K 0 , (c)) or not (K ε , (b)) and whether inertia constants are homogeneous (M 0 ) or not (M ε ). We discuss these different choices in detail in appendix B. For each combination of these attributes, we then compare the control effort before introducing the isolator, P no iso , to the control effort after introducing the isolator, P w iso , for both the 'driven' subgraph where the wind farm is introduced and the 'passive' subgraph that is shielded by the isolator. Remarkably, the control effort in the shielded subgraph is reduced drastically in most scenarios. The effect is most strongly pronounced when edge weights are homogeneous where we achieve almost a tenfold reduction of the control effort on average. Notably even though we increased the connectivity between the two subgraphs, thus potentially allowing for fluctuations to propagate more easily, we reduced the control effort necessary to keep the system at its nominal frequency.
The fact that the reduction works best for homogenised networks can be understood as follows; for a network with homogeneously distributed coupling constants K ij and inertia constants M i , the two driven nodes react similarly to the stochastic driving signal and thus the perturbation is comparable to the constant load shift studied in figure 6 since both nodes are stochastically perturbed by signals with opposite signs. As a consequence, algorithm 1 may be applied at each point in time and the resulting signal is effectively shielded. This is in contrast to the situation were all nodes-and edges-have very heterogeneous parameters such that the nodal frequenciesθ i and nodal phase angles ϑ i behave very differently.

Discussion and conclusion
In this article, we outlined possible applications of network isolators in networks of coupled Kuramoto oscillators. We demonstrated that network isolators allow to control multistability by suppressing targeted fluctuations at a specific, desirable fixed point. Our approach may be applied to stochastic, time-dependent fluctuations as well as to structural perturbations applied to both, nodes and edges. Finally, we illustrate that network isolators may be readily used to achieve these goals in real-world networks such as power grids.
Our results offer a new framework to deal with multistability in the Kuramoto model. Static perturbations may be almost completely suppressed for a specific, desired state of the network by tuning the edge weights in a network containing a network isolator-or constructing new edges that are equipped with the necessary edge weights. Recent analyses have demonstrated the possibility of multistability in real-world networks such as power grids [24,25,27] that result in undesired phenomena such as cycle flows that do not result in any power supply [27] or power flows resulting in particularly high losses [25]. Our approach allows to suppress fluctuations specifically for the states with low losses and a minimum amount of cycle flows, thus increasing the efficiency of the power grid.
While we focussed on the Kuramoto model with positive coupling in this manuscript, we expect our results to be valid also in case some edges have negative coupling. Kuramoto oscillators that are joined by a negative coupling repel one another instead of attracting and eventually synchronising. As a result, phase differences between neighbouring vertices in a fixed point tend to increase when negative coupling is present. Thus, we expect our approach of targeted shielding to yield even better results since larger phase differences correspond to larger weight modifications in algorithm 1. This may affect the effective adjacency matrix in a similar way as multistability induced by cycle flows does. In appendix D we examine potential applications of the weight tuning algorithm 1 in a network with negative coupling and demonstrate that, for this example, a significant improvement in terms of the shielding can be achieved by tuning the edge weights.
Besides increasing the resilience to static, permanent perturbations, we also illustrate that the isolator suppresses stochastic fluctuations and thus the control effort. This largely increases the applicability of network isolators to real-world power grids where an increase in fluctuating renewable power generation may result in the propagation of fluctuations through the power grid [50,54]. In a cellular grid where individual cells are coupled via network isolators, the propagation of failures and fluctuations is strongly suppressed as demonstrated here. In our article, we analysed two different types of noise-Wiener noise and a noise signal with properties similar to the statistics of wind power generation time series. Since both signals are effectively shielded by the isolator, we expect this to hold also for noise signals with stochastic properties that were not explicitly incorporated into the present analysis.

Acknowledgments
We thank Vito Latora for valuable discussions, Matthias Gries for help with processing the data of the RTS-96 test case and Oliver Kamps for help creating the synthetic wind power generation time series. We gratefully acknowledge support from the German Federal Ministry of Education and Research (BMBF Grant No. 03EK3055B) and the Helmholtz Association (via the joint initiative 'Energy System 2050-a contribution of the research field energy' and the Grant No. VH-NG-1025 to DW).

Data availability
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A. Linear response theory
In this section, we briefly discuss how to analyse structural perturbations, i.e. modifications of edge weights, and nodal perturbations, i.e. perturbations to the nodal frequencies ω, within the framework of linear response theory. We will follow reference [44] in the following discussion.

A.1. Frequency perturbation
We consider a Kuramoto-network in a fixed point, described by the steady-state equation Next we apply a dipole-like perturbation Δω i of the form to nodes (n, m) on the network, where Δω ∈ R is the perturbation strength. We assume that the perturbation is small enough such that the system relaxes to a new fixed point which can be described by ϑ i = ϑ * i + Δϑ i . Thus, we can write down the steady-state equation of the new fixed point: where δ nm is the Kronecker delta. Expanding this equation at ϑ * i and subtracting the original fixed point condition (A.1) leads to Dropping the higher order terms, and re-writing the equation in a vectorised form leads tõ where e i ∈ {0, 1} n is a vector with entry one at position i and zero otherwise and Δϑ = (Δϑ 1 , Δϑ 2 , . . . , Δϑ n ) ∈ R N is the vector of phase changes due to the perturbation. We observe that matrixL N×N has the form of a Laplacian matrix and thus refer to it as the effective Laplacian of the system which is of the formL

A.2. Structural perturbations to edges
To describe frequency perturbations we again assume the system to be in a steady state according to equation (A.1). Now, assume a perturbation of edge e = (n, m) in the graph G(V, E), such that we can write Again, assume that the edge perturbation is small enough such that the system relaxes to a new fixed point described by: We expand the steady-state condition of the new fixed point at ϑ * i leading to Next, we subtract the steady-state condition of the original fixed point equation (A.1) to obtain Expanding this in K ij at K ij leads to which describes the response of the system to the edge perturbation. Using equation (A.5) we can write the sum including ΔA ij in the expression above in terms of: which has precisely two entries with opposing signs and thus is a dipole-like perturbation. Next, we drop the higher order terms yielding In a vectorised form this can be written as s ≈LΔϑ, which is exactly the same form as for frequency perturbations, equation (A.3) and thus yields the same effective LaplacianL (see equation (A.4)).  (21)) shows the deviation from unit rank. Thus, as ξ(Ã 12 ) = 0 is reached, the fixed point is perfectly shielded.

Appendix B. Numerical methods
In this section, we discuss the numerical methods used in this article. First we present the algorithm to tune the edge weights of the network isolator in order to shield a specific fixed point and then explain in detail the numerical methods used to simulate the Kuramoto model for every single figure.

B.1. An algorithm for fixed point shielding
In section 3.1 we showed that we can tune the edge weights K ij , ∀ (i, j) ∈ E(G) of the isolator nodes such that rank(ã) = 1 and thus improve the shielding behaviour of the network isolator for a certain fixed point ϑ * . As a modification of the edge weights K ij → K mod ij slightly alters the fixed point ϑ * → ϑ mod , we repeatedly applied the following procedure until the absolute changes reached below a predefined threshold value ε = 1e − 10. We start with the initial weights K [0] ij = K ij and denote the edge weights and angles at iteration step n by K [n] ij and ϑ [n] , respectively.
(a) Tune edge weights in the network isolator such that K [n] ij is replaced by (b) Relax the system to a new fixed point ϑ [n] → ϑ [n+1] . (c) Compare the updated to the old coupling matrix. If stop.
Note that in rare-2 out of 10 000-cases if choosing the same parameters as in figure A.1-the weight modification may lead to a loss of stability of the fixed point ϑ * . In these cases we reject the fixed point and chose a new one. Also note that there are several possibilities of tuning the weights such that rank(ã) = 1 is fulfiled and that comparing the updated to the old weight matrix is not the only possible termination condition.
In figure A.1, we illustrate the convergence of this algorithm. To this end, we considered a hexagonal Kuramoto network with N = 60 nodes which is subdivided into two subgraphs with n 1 = n 2 = 30 nodes connected by a network isolator G which consists of k 1 = 3 nodes that are fully connected with k 2 = 3 nodes in the other subgraph. We performed the algorithm for 50 different fixed points ϑ and displayed the mean (lightblue, circles) and the maximum to minimum (lightblue, shaded) values for three different measures. We also show the mean values for the IEEE RTS96 test grid. In panel (a), we show that the absolute difference in the edge weights between two subsequent steps n and n + 1, |K [n+1] − K [n] |, decreases monotonously with the number of iteration steps n. After n = 10 steps all 50 runs converged s.t. ε < 10 −10 ∀ (i, j) ∈ E(G). In panel (b) we illustrate the convergence of the phase shift |ϑ [n+1] − ϑ [n] | per iteration step. The deviation of the rank of matrix a from unit rank is given by its coherence statistics ξ(ã) (see equation (21)) which is shown in panel (c). The condition ξ(ã) = 0 characterising a perfect network isolator in the linear response regime is reached before the capacity updates reach the threshold value ε.

B.2. Simulating the Kuramoto dynamics
In this section we give a detailed analysis on how we built figures 3-7. In figures 3, 4 and 6, we have indicated in light grey any flow changes ΔF that have fall below a certain threshold value, which is indicated with the lowest value on the colourbar of each figure, respectively.
To create figure 3, we set up a network of first order Kuramoto oscillators, cf equation (2) on two connected hexagon graphs as shown in panels (a) and (b). We draw the natural frequencies from a Gaussian distribution, ω ∈ N (0, 0.005), and distribute the coupling constants uniformly K ij = 1, ∀(i, j) ∈ E(G). After the system is relaxed to a stable fixed point, we perturb two nodes in the left subgraph (black arrows) with a dipole-like frequency perturbation (cf equation (14)) with perturbation strength Δω = ±0.1. Subsequently, the system relaxes to a new fixed point, leading to flow changes |ΔF| on the edges. In (c), we show the mean flow changes |ΔF| with respect to the edge distance, i.e. the shortest path between two edges, to the perturbation location for both networks. In (d) we plot the mean of all flow changes of the passive subgraph |ΔF passive | against the perturbation strength applied to the same nodes as marked in (a) and (b). The dashed line at Δω = 0.1 represents the value of the perturbation strength which is applied in panels (a), (b) and (c).
To create figure 4, we set up a network of first order Kuramoto oscillators (cf equation (2)) on two hexagonal graphs connected via a network isolator. Again, the natural frequencies are drawn from a Gaussian distribution, ω ∈ N (0, 0.005). The two systems were relaxed to two different fixed points with winding numbers q = (0, 0) and q = (1, −1) , respectively. These fixed points were obtained with uniformly distributed coupling constants K ij = 1, ∀(i, j) ∈ E(G). Subsequently, we applied dipole-like perturbations in the driven module with strength Δω = ±0.01 in both fixed points that result in flow changes |ΔF|. We then tuned the weights as described above and relaxed the system in an iterative manner until the overall, absolute change in the edge weights between two subsequent iterations was below a threshold value ε = 1e − 10. To create panels (e) and (h) we perturbed every possible node pair in the driven module and computed the subsequent flow changes |ΔF|.
To create figure 5, a network of second order Kuramoto oscillators as in equation (3) was set up on the sixregular graph containing a network isolator shown in figure 5(a). The natural frequencies were drawn from a Gaussian distribution, ω ∈ N (0, 25), the coupling constants were distributed uniformly K ij = 10, ∀(i, j) ∈ E(G) and the ratio between inertia and damping constant was fixed to α i = M i D i = 1, ∀ i ∈ V(G). The network was then first relaxed to a stable fixed point. Then the edge weights in the isolators were tuned by repeatedly-adjusting the edge weights and relaxing the system again until no further changes occur. Finally, the two nodes shaded brownish were driven stochastically by a Wiener process with zero mean and unit variance and the response at the indicated nodes was monitored.
For both figures 6 and 7, we used the dataset of the test grid IEEE 'RTS-96' published in reference [49]. For each node, we aggregated the load P (l) and generation P (g) data of all machines and loads attached to the respective node resulting in an effective load P = P (g) + P (l) which is identified with the natural frequency ω = P of the resulting node in the Kuramoto model, i.e. we model each node as a second order Kuramoto oscillator in correspondence with the SM model [17]. We converted the loads to the p.u. system with a base value of P 0 = 100 MW. Since we consider the Kuramoto model as described in equation (3), we effectively model a lossless power grid and neglect voltage fluctuations, thus assuming a value of V i = 1 p.u. for all nodes i. In this setup, the load and generation values do not add up to zero, but there is an overproduction of power resulting in an unbalanced grid. We therefore increase the load at each node by 5% with respect to the original dataset and subtract the remainder from the slack bus. The coupling constants are simply taken to be the inverse line reactances K ij = x −1 ij given in the dataset. To build a network isolator in the RTS-96 we add edges to the node pairs (121, 318) with weight K 121,318 = K 121,325 and (325, 223) with weight K 325,223 = K 223,318 such that rank(A 12 ) = 1. We chose the weights to be similar to the weights of their respective edge neighbours in the isolator to create a realistic scenario. We then tuned the network isolator weights as described in the above algorithm and illustrated in figure A.1.
To create figure 6, we relaxed the RTS-96 network of Kuramoto oscillators to a stable fixed point ϑ * . Afterwards we applied a dipole perturbation at the node pair (318, 321) with Δω 318 = +1p.u. and Δω 321 = −1p.u. which lead to a new fixed point ϑ . We computed the subsequent flow changes |ΔF| (Equation (18)) for the three cases: The subgraphs are connected with (c) no isolator, (d) an isolator, (e) a tuned isolator. Then, we colour-coded the subsequent flow changes |ΔF i | on their respective edges. To create panel (b) we compute the mean flow changes |ΔF| and their variance Var(|ΔF|) at all distances to the perturbation. We also distinguish between the respective subgraphs, as the same edge distance may be found in both subgraphs.
To create Figure 7, we relaxed the network of second order Kuramoto oscillators to a stable fixed point and then added and subtracted the stochastic signal F(t) to the upper and the lower generator node, respectively, that are marked by arrows in panel (a). The stochastic signal is produced using the procedure described in reference [51] with the resulting time series having part of the stochastic properties of the power generation time series of a wind farm. For all panels, we used a penetration of p = 0.5 to add the signal to the nodes. When considering heterogeneous inertia constants M ε as in figure 7, we estimated each node's inertia constant M i based on the node's natural frequency ω i = P i using the estimate discussed in references [17,52]. When considering homogeneous inertia constants M 0 , we set all constants to M i ≈ 0.13s 2 . Damping constants were set homogeneously to D i ≈ 0.13s for both cases as discussed in references [17,52] as well. When considering homogeneous coupling constants K 0 , we used the mean coupling constant of the original dataset for all edges.

Appendix C. Extended theoretical analysis of network isolators in the Kuramoto model C.1. Network isolators in the Kuramoto model
In this section we extend on the results on network isolators in the Kuramoto model presented in section 3.1.
We consider the vectorised Kuramoto model on a weighted graph G(V, E) with coupling matrix K ∈ R L×L . We furthermore assume that the graph can be divided into two subgraphs G 1 and G 2 that are connected by a network isolator G such that we can write the graphs adjacency matrix A ∈ R N×N as where A 1 and A 2 are the adjacency matrices of the subgraphs G 1 and G 2 , respectively. The submatrices A 12 and A 12 represent the adjacency of the mutual connections between G 1 and G 2 . Without loss of generality A 12 can be written as A 12 = a 0 0 0 , with a ∈ R k 1 ×k 2 , i.e. dimension of the isolator is k 1 × k 2 . Now, we assume that the system is relaxed into a stable fixed point ϑ * , described by the fixed point condition ω = IK sin(I ϑ * ).
Next, we apply a sufficiently small dipole-like perturbation of the form equation (A.2), i.e. a frequency shift Δω, to two nodes in the subgraph G 1 and let the system relaxes into a new fixed point ϑ . We can compute the response to the dipole perturbation by making use of linear response theory as in appendix A yielding the effective LaplacianL (cf equation (A.4)) which depends on the fixed point ϑ * before the perturbation (cf equation (A.3))L Δϑ = Δω.
We are thus left with a discrete Poisson equation as in theorem 1. However, the effective adjacency matrix of the system now readsÃ
We thus now longer have rank(ã) = 1 for the system in linear response. Nevertheless, this can be overcome by tuning the initial weights K nm of the edges (n, m) ∈ E(G) in the network isolator. A possible choice is to tune the weights such that K new mn = K mn cos(ϑ * m − ϑ * n ) −1 . Note that there are several possibilities of tuning the weights to achieve rank(ã) = 1. However, any modification of the coupling matrix K will slightly alter the initial fixed point ϑ * → ϑ mod and thus slightly change the condition needed to fulfil rank(ã) = 1. Assume that we can find a tuned weight matrix K mod with slightly altered angles ϑ mod (see section appendix B for details). Then, we can apply theorem 1 to the modified set of equations wherẽ