# Minimal fatal shocks in multistable complex networks

Jul 16, 2020

### Minimal fatal shock

The first step in identifying the MiFaS for a given system is to define a desired state (mathbf {X_0}). We then assume that, prior to perturbations, the system resides on (mathbf {X_0}) and that a shock—applied at (t=0)—kicks the system’s state instantaneously to (mathbf {X}(0)). A shock—now defined as (mathbf {x}(0) = mathbf {X}(0)-mathbf {X_0})—is said to be fatal if (mathbf {X}(0)) is located outside the basin of (mathbf {X_0}) and non-fatal if (mathbf {X}(0)) is located within the basin of (mathbf {X_0}). Accordingly, the MiFaS is a vector which displays the shortest distance between the desired state and its basin boundary and the corresponding direction in state space (Fig. 1a).

The second essential step, is defining a norm for the perturbation size. It is important to note that the use of a certain norm is not only a technical but also an interpretative decision. Throughout this work, we use the Euclidean distance to the desired state (mathbf {X_0}) to quantify the magnitude d of a perturbation

begin{aligned} d ; = , ||mathbf {x}(0) || , = , ||mathbf {X}(0) – mathbf {X_0}||. end{aligned}

(1)

To determine the MiFaS, we develop a search algorithm which is based on the minimal seed approach41 and which can be divided into two stages, the global random initialization (stage I) and the local non-random optimization (stage II).

In stage I, we randomly draw initial conditions from a shrinking subspace in state space to find a fatal shock with a preferably small magnitude d (see “Methods” and Supplementary Fig. S1). Stage II starts with the smallest fatal shock received from stage I (Supplementary Fig. S1). From this point on, we take two seemingly opposing steps. First, we adapt the direction of (mathbf {x}(0)) in order to move (mathbf {X}(0)) away from the basin of (mathbf {X_0}) while keeping d fixed. Second, we move (mathbf {X}(0)) towards the basin by reducing d by a step size (Delta d). By repeating these two steps iteratively, we attain smaller and smaller fatal shocks which finally converge towards a local MiFaS (see Fig. 1b and Supplementary Fig. S1). It is important to note that the outcome of the search—and thus the achieved local MiFaS—is dependent on the initialization in stage I. Accordingly, to attain the global MiFaS, we need to run the search algorithm multiple times and select the minimum of the local MiFaS as the global one.

The centerpiece of the outlined algorithm is the adaptation of the direction of (mathbf {x}(0)) during stage II, which aims at maximizing the distance between (mathbf {X}(0)) and the basin boundary of (mathbf {X_0}). However, since this distance is not easily accessible, it is approximated by an objective function which can be maximized within a constraint optimization. For the two applications we present here, the objective function can be thought of as the amplification of the shock over a preselected time T (see “Methods” for specific definition). The mechanism behind this is that trajectories close to the basin boundary stay close to it for long times as they move along the stable manifold of a saddle-type state while trajectories far off the boundary approach an alternative attractor faster and thus lead to earlier and stronger amplifications.

In summary, as a result of the optimization procedure we obtain the magnitude of the smallest distance to the basin boundary which can be utilized as a quantitative measure of global stability and the direction of the perturbation in the high-dimensional phase space.

### Plant–pollinator networks

In our first example, we consider a simple model of mutualism which captures the crucial aspects of a system of plants and their corresponding pollinators43,45. The mutualistic system is described as a bipartite network, with one set of nodes representing a number of (N_P) plant species and one set representing a number of (N_A) animal species whose dynamics are given by

begin{aligned} frac{mathrm {d} P_i}{mathrm {d} t} ,&= , alpha P_i , – , sum _{k=1}^{N_P} beta _{ik} P_i P_k , + , frac{sum _{j=1}^{N_A} gamma _{ij} A_j P_i}{1 + h sum _{j=1}^{N_A} gamma _{ij} A_j},nonumber \ frac{mathrm {d} A_j}{mathrm {d} t} ,&= , alpha A_j , – , sum _{l=1}^{N_A} {tilde{beta }}_{jl} A_j A_l , + , frac{sum _{i=1}^{N_P} {tilde{gamma }}_{ji} P_i A_j}{1 + h sum _{i=1}^{N_P} {tilde{gamma }}_{ji} P_i}, end{aligned}

(2)

where (P_i) denotes the abundance of plant species i ((i=1, ldots , N_P)) and (A_j) the abundance of animal species j ((j=1, ldots , N_A)). In Eq. (2), the parameter (alpha) gives the intrinsic growth rate, (beta _{ik}) (({tilde{beta }}_{jl})) the competitive pressure of plant (animal) species k (l) on plant (animal) species i (j), (gamma _{ij}) (({tilde{gamma }}_{ji})) the benefit plant (animal) species i (j) obtains from animal (plant) species j (i) and h the handling time for pollination. As a general principle, we assume the benefit a species gains from pollination to be obligatory for its own growth, an assumption which is necessary to obtain multistability in this model57. Therefore, we choose the net growth rate (alpha le 0).

In order to keep the parametrization as simple as possible, we set (alpha), (beta _{ii}) (({tilde{beta }}_{jj})) and h to be equal for all species. To reduce the complexity of the overall interaction pattern, we assume all-to-all coupling for the interspecific competition between species within one set, whereby (beta _{ik}=beta _0/(N_{P}-1)) for (i ne k) (({tilde{beta }}_{jl}=beta _0/(N_{A}-1)) for (j ne l)). By contrast, a mutualistic interaction between an animal and a plant species can either be absent, in which case (gamma _{ij}=0) (({tilde{gamma }}_{ji}=0)), or present, in which case (gamma _{ij}=gamma _0/kappa _i) (({tilde{gamma }}_{ji}=gamma _0/{tilde{kappa }}_j)), where (kappa _i) (({tilde{kappa }}_j)) denotes the degree or the number of mutualistic partners of plant (animal) species i (j). This formulation corresponds to a full trade-off between the benefit a species attains from one partner and the number of partners this species has45. An important aspect of the chosen parametrization is that species solely differ on account of their position in the mutualistic network. In the following, we determine the MiFaS for realistic plant–pollinator networks from the Web of Life Database58 representing networks from different geographic locations across various climate zones (see Supplementary Fig. S5 and Supplementary Table S2). With (alpha = -0.3), (beta _{ii}=1.0), (beta _0 = 1.0), (gamma _0 = 4.5) and (h=0.1), we choose the model parameters in a way that ensures that each of the studied systems possesses a state in which all species coexist. This ’desired’ state (mathbf {X_0}) is opposed to multiple ’undesired’ states in which one or more species are gone extinct (the MiFaS is actually fatal).

To interpret the results, it is useful to state some general considerations first. Due to the mutualism, the growth of a species depends on the abundance of its mutualistic partners. As the growth of these partners can also depend on further other partners, these further partners indirectly support the growth of the first species. We could continue building this chain of dependencies but essential is that species being close to each other within the network and especially those sharing partners benefit from each other. On the other hand, due to competition high abundances of one species directly impede the growth of all species within the same group (animals or plants). Hence, the net effect which an increase or decrease of a species’ abundance has on another species depends on the interplay between the two processes. The indirect benefits can either balance or enhance the negative effects due to competition depending on whether species are close (balance) or far apart (enhance).

At first, we compute the minimal fatal shock (MiFaS) for an exemplary network from Morant Point in Jamaica (Fig. 2a). The topology of this system is characterized by an asymmetric division into a small tree-like part and a large core, i.e. a large mostly well connected component. This topological division is mirrored in the direction of the MiFaS which is visualized by the color-coding. A small negatively perturbed part consisting of the tree-like periphery (nodes within the yellow shaded region in Fig. 2a) plus its single non-peripheral neighbor is opposed to the rest of the network which is positively perturbed. This division exemplifies how the mutualistic and competitive interactions between species shape the system’s response to perturbations. In the tree-shaped part of the network, all species are close to each other but far away from most other species. Furthermore, due to the sole connection between the two characteristic structural parts of the network, the share of partners between the two is minimal. As a result, the interdependency of species within the tree-shaped part is extremely high. Accordingly, the loss of abundance of any species in the tree-like structure—as it is the case in the MiFaS (Fig. 2)—significantly affects all other species in this tree-like periphery. On the contrary, the competitive stress due to species within the large component is high as it is not balanced by the indirect benefits. It is actually even enhanced as the increase of abundance of one species boosts the growth of its partners which again enhances the competive stress on the peripheral tree-like structure.

After the system has been hit by the MiFaS, all ten species within the tree-like periphery are lost in the long run (Fig. 2c and yellow shaded region in Fig. 2a). The remaining species—except for the single neighbor of the periphery—tend to higher abundances as the competitive pressure on them is relaxed. Accordingly, the new asymptotic state (Fig. 2c) again shows that the net impact of the peripheral species on most other species has been negative. Apart from the new asymptotic state, the transient leading there (Fig. 2b,c) is of interest as well. In fact, the transient behavior is typical for an initial state close to the basin boundary which is made up by the stable manifold of a saddle point. The transient at first moves towards the saddle fast (Fig. 2b), stays in its vicinity for some time as the repulsion is weak and finally settles on an attractor which, in this case, is the undesired state of partial extinction (Fig. 2c).

Overall, we examine the MiFaS for a total of 59 plant–pollinator systems, each being based on one of the real-world network topologies. For comparison, we order the networks from sensitive to robust according to the magnitude of their respective MiFaS and depict the direction of the MiFaS for five further exemplary systems (Fig. 3).

Some characteristics found for the MiFaS of the exemplary network (Fig. 2) prove to be generally valid. For each system, the division of the MiFaS into a small negatively perturbed part and a larger but weaker positively perturbed part displays how mutualistic interdependency and competition shape the system’s response to perturbations. In this context, the negatively perturbed part marks the weakest point of the network at whose outer edge the extinction occurs. Speaking in ecological terms, we find these weak points always being associated with specialization and the distribution of negative perturbations depends on the nature of the caused interdependency: in the exemplary system (network 1 in Fig. 3), where the specialization among all species within the tree-like structure is rather mutual, all involved species are significantly perturbed (the same for network 13 and partly for network 4, Fig. 3). However, the more asymmetric the specialization gets—meaning that many specialists are connected to a single generalist—the stronger the negative perturbation focuses on this generalist (networks 4 (rightarrow) 26 (rightarrow) 27 (rightarrow) 49, Fig. 3). This perturbation structure proofs to be efficient as the dependency of the generalist on each single specialist is low but its cumulated dependency on all specialized partners is high. A perturbation at the generalist therefore induces a negative feedback whose strength also depends on the number of connections the generalist has to other-non-specialized species. Accordingly, network 49 is much more robust than network 26 as the decisive generalist is highly connected to the core.

The positive contribution to the overall MiFaS marks the impact of competitive forces which depends on the global interdependency among species. In the case of a single well-connected core and a periphery which only consists of specialists being directly connected to this core, indirect positive effects between species balance competive effects as all species are close and well connected. Accordingly, we do not find any significant contribution of positive perturbations to the overall MiFaS (networks 37, 49, Fig. 3). The contrary is the case if the core is not well build, meaning that only a few connections between important hub nodes exist (networks 4, 26) or if—due to strong reciprocal specialization—a larger peripheral structure exists (networks 1, 13). In such cases, positive perturbations at rather central core-species contribute significantly to the overall MiFaS and thus to the extinction of peripheral species. In summary, a strong global interdependency among all species favors a system’s robustness whereas a strong local interdependency paired with a weak global interdependency depicts the worst case scenario.

### Great Britain power grid

As a second example we consider a coarse-grained model of a power grid which exhibits synchronization dynamics. In this framework, a power grid is described as a network of Kuramoto-like13 second order phase oscillators whose dynamics are given by

begin{aligned} frac{mathrm {d} phi _i}{mathrm {d} t}&= omega _i nonumber \ frac{mathrm {d} omega _i}{mathrm {d} t}&= P_i – alpha omega _i + sum ^N_{j=1} K_{ji} , sin (phi _j-phi _i), end{aligned}

(3)

where (phi _i) and (omega _i) denote the phase and frequency deviation of oscillator i from a grid’s rated frequency (which will hereinafter be referred to as phase and frequency). The parameters (alpha) and (P_i) are the grid’s damping constant and the net power input/output of oscillator i, respectively. The capacities of the transmission lines and therefore also the topology of the grid are contained in the matrix K, with (K_{ji}=K_{ij}>0) if oscillators i and j are connected and (K_{ij}=0) otherwise.

As an example, we consider the Great Britain power grid which consists of 120 nodes and 165 transmission lines59. For reasons of simplification, we assume one half of the oscillators to be generators ((P_i=+P_0)) and one half to be consumers ((P_i=-P_0)) whose distribution within the grid we draw randomly (see Fig. 5). Furthermore, we choose the same maximum capacity for all transmission lines, either (K_{ij}=K_0) or (K_{ij}=0). In a realistic parameter setting of this model, one ’desired’ synchronized state ((phi _i=const) and (omega _i=0) for all i) representing stable operation competes with several ’undesired’ non-synchronized states. With (alpha =0.1), (P_0=1.0) and (K_0=5.0), we choose the model parameters accordingly. In this setting, the MiFaS represents the smallest perturbation to the synchronous state which induces a shift to one of the non-synchronous states interpreted as a power outage.

The combination of frequencies and phases is actually problematic when determining the MiFaS since they differ in units. We therefore only take into account perturbations in the frequencies (omega). In this context, choosing the frequencies (omega) instead of the phases (phi) seems reasonable as disturbances usually occur due to fluctuations in the power generation or consumption60. Such parametric disturbances would first affect the frequencies via (mathrm {d}omega /mathrm {d}t) (Eq. 3). Furthermore, considering only frequencies allows a clearer depiction of the MiFaS, since the corresponding vector contains exactly one entry per node of the power grid.

Examining one random realization of the power grid (Fig. 4a), we find that, like in the exemplary plant–pollinator network, the MiFaS is associated with a tree-like structure including the most peripheral nodes of the network (according to the resistance centrality proposed by61, see Supplementary Fig. S7). In fact, the same structure is highlighted by some of the eigenmodes of the graph Laplacian (see Supplementary Fig. S8). However, apart from the observation that the MiFaS is orthogonal to a neutral perturbation affecting all oscillators in the same way which is equivalent to its first eigenmode, we find no simple connection to the graph Laplacian (see Supplementary Information).

In order to understand the effectiveness of the MiFaS, it is instructive to have a closer look at how the desynchronization occurs after the system has been hit by the MiFaS (Fig. 4c,d). The desynchronization is triggered by an overload on the transmission line which connects the seven northermost oscillators to the rest of the grid (Fig. 4b). Due to the accumulation of consumers within this tree-like structure (5 consumers towards 2 generators), already in the unperturbed state, the load—(K sin (phi _j – phi _i)) for the line connecting nodes j and i—on the ’trigger transmission line’ is comparatively close to its maximum capacity K (see Fig. 4c). Intuitively, a strong deceleration of oscillators inside plus an acceleration of oscillators outside the tree-like structure seems to be an efficient way to induce an overload. Indeed, we find the strongest negative perturbations at the seven oscillators within (Fig. 4b) as well as positive perturbations at several oscillators outside the tree-like structure. However, in the northern part of the grid, the overall MiFaS roughly follows a broad gradient distribution with negative perturbations on both sides of the trigger transmission line and the strongest positive perturbations at rather distant nodes in the northwest of Great Britain. This distribution is efficient as the perturbations in frequencies first have to be transferred into phase deviations to induce an overload. A relatively smooth gradient ensures that the arising phase deviations are balanced slowly and thus a large transmission load can build up.

This transfer can be observed in the first stage of the transient following the MiFaS (Fig. 4c,d). In this stage, the system evolves rather smoothly towards a point where the frequency deviations of all oscillators are close to zero but where, at the same time, the transmission load on the trigger line (red line in Fig. 4) has passed its maximum capacity. The system subsequently enters a stage in which both transmission loads as well as frequencies oscillate erratically until the oscillations suddenly collapse and the system settles on an undesired attractor. It is remarkable that the final overload (green line in Fig. 4) is not located on the line which triggered the desynchronization but on a line deeper in the tree-like structure (Fig. 4c). The final overload is similar to a cutoff of two consumers from the rest of the grid, as the frequencies in the two departed components evolve more or less independently. It is however important to note that this particular undesired state represents only one of several possible outcomes. Indeed, already the slightest variation (smaller than the finite precision of the search algorithm) of the initial perturbation can lead to a different non-synchronous asymptotic state, although the trigger transmission line is always the same. Such high sensitivity is often an indicator for complexly intervowen basins of attraction, characteristic to many highly multistable systems62.

In order to gain more insights into how certain topological features harm a power grid’s stability against shocks, we examine some of the local MiFaS inducing power outages (Fig. 5). These local minima correspond to different outcomes of the applied optimization scheme for the same network topology and parametrization and thus represent further close but less crucial distances between the desired state and its basin boundary. As we are interested in distinct topological weak points of the grid, we take into account only those local minima which differ in the involved trigger transmission line (highlighted edges in Fig. 5).

The local MiFaS, and in particular the examination of the associated trigger transmission lines, reveal two mutually reinforcing sources for the emergence of weak points. Firstly, desynchronization events are triggered on transmission bottlenecks which result from the loose connection between a peripheral subgraph and the rest of the grid. Four out of five of the shown local MiFaS (Fig. 4 and Fig. 5a–c) are actually related to the most pronounced case of such a bottleneck which is a bridge, i.e. a single edge connecting two subnetworks. Secondly, the accumulation of oscillators of the same type within a subgraph induces a local mismatch between power generation and consumption (Fig. 4 and Fig. 5a–d). We find each of the shown local MiFaS to be related to such a local mismatch. Already in the unperturbed state, this mismatch has to be balanced by a high initial load on the connecting transmission line(s) which in turn results in a low threshold for an overload (Fig. 5d). This overload is then triggered by the MiFaS by reinforcing the generation/consumption imbalance between the two subgraphs. Accordingly, all fatal shocks involve strong frequency perturbations with a sign according to the already established power mismatch in the peripheral subgraph and frequency perturbations in the opposite direction in adjacent areas of the grid. However, as in the global MiFaS, the boundary between positive and negative perturbations is not sharp but more (Fig. 5a,c,d) or less (Fig. 5b) follows a kind of gradient.

Of particular interest is the local MiFaS shown in Fig. 5c as its underlying topological motif is quite common in the network: a node with degree 1, also termed ’dead end’32. Apart from the two dead ends within trees (Fig. 5a,b), the portrayed dead end is the one being most sensitive to perturbations despite or seemingly because it is connected to a rather central node of degree 6 (see also Supplementary Fig. S7). For none of the surroundings of the other dead ends, which are all adjacent to lower degree nodes, we find a local MiFaS of similar low magnitude. Accordingly, we conclude that a rather central position of the node from which the peripheral subgraph branches off might actually harm its robustness against particular perturbations.