Dynamics of Tipping Cascades on Complex Networks

TL;DR

This paper explores how different network topologies influence the spread of tipping points in complex systems, using simulations on various network models and real-world data, revealing that clustering and spatial organization increase vulnerability.

Contribution

It introduces a comprehensive analysis of how network structure affects tipping cascades, combining model simulations with real-world data from the Amazon rainforest.

Findings

Clustering and spatial organization increase network vulnerability.

Real-world networks like the Amazon are more susceptible to cascades.

Network topology can inform system design to enhance robustness.

Abstract

Tipping points occur in diverse systems in various disciplines such as ecology, climate science, economy or engineering. Tipping points are critical thresholds in system parameters or state variables at which a tiny perturbation can lead to a qualitative change of the system. Many systems with tipping points can be modeled as networks of coupled multistable subsystems, e.g. coupled patches of vegetation, connected lakes, interacting climate tipping elements or multiscale infrastructure systems. In such networks, tipping events in one subsystem are able to induce tipping cascades via domino effects. Here, we investigate the effects of network topology on the occurrence of such cascades. Numerical cascade simulations with a conceptual dynamical model for tipping points are conducted on Erd\H{o}s-R\'enyi, Watts-Strogatz and Barab\'asi-Albert networks. Additionally, we generate more realistic networks using data from moisture-recycling simulations of the Amazon rainforest and compare the results to those obtained for the model networks. We furthermore use a directed configuration model and a stochastic block model which preserve certain topological properties of the Amazon network to understand which of these properties are responsible for its increased vulnerability. We find that clustering and spatial organization increase the vulnerability of networks and can lead to tipping of the whole network. These results could be useful to evaluate which systems are vulnerable or robust due to their network topology and might help to design or manage systems accordingly.

Figures

Equations

$$\frac{dx}{dt}=-a(x-x_0{)}^3+b(x-x_0)+r,$$

$$\frac{d{x}_i}{dt}=-a(x_i-x_0{)}^3+b(x_i-x_0)+\underset{{\tilde{r}}_i(x_1,x_2,...,x_N)}{\underbrace{r_i+d\sum _{j=1,j\ne i}^N{a}_{ij}{x}_j}},$$

$$a_{ij}=\{\begin{array}{ll}1,& \text{if there is a directed link from element i to j}\\ 0,& \text{otherwise}\end{array}.$$

Full text

††preprint: APS/123-QED

I Introduction

In the last decades the study of tipping elements has become a major topic of interest in climate science. Tipping elements are subsystems of the Earth system that may pass a critical threshold (tipping point) at which a tiny perturbation can qualitatively alter the state or development of the subsystem Lenton et al. (2008). However, tipping points also occur in various complex systems such as systemic market crashes in financial markets May et al. (2008), technological innovations Herbig (1991) or shallow lakes Scheffer and van Nes (2007) and other ecosystems Scheffer et al. (2001). Understanding their dynamics is thus crucial not only for climate science but also for other disciplines that use complex systems approaches.

Many tipping elements are not independent from each other Brummitt et al. (2015). In such cases, if one tipping element passes its tipping point, the probability of tipping of a second tipping element is often increased Kriegler et al. (2009), yielding the potential of tipping cascades Steffen et al. (2018) via domino effects with significant potential impacts on human societies in the case of climate tipping elements Cai et al. (2016). In this study, we investigate the dynamics of complex networks of interacting tipping elements. A tipping element is described by a differential equation based on the normal form of the cusp catastrophe which exhibits fold-bifurcations and hysteresis properties. The interactions are accounted for by linear coupling terms. Many environmental tipping points can be described as fold bifurcations Lenton (2013) and prototypical conceptual models that exhibit fold bifurcations have been developed for the Thermohaline Circulation Wright and Stocker (1991), the Greenland Ice Sheet Levermann and Winkelmann (2016), or tropical rainforests Staal et al. (2018a) among others. Coupled cusp catastrophes have already been studied in detail for two or three subsystems Abraham et al. (1991); Brummitt et al. (2015); Klose et al. or in combination with Hopf bifurcations Dekker et al. (2018). On the other hand, threshold models for global cascades on large random networks have been investigated Watts (2002).

Here, we study cascades in complex systems with continuous state space that are moderate in size, yet large enough for statistical properties of the complex interaction networks to become relevant. Cascades in complex systems with continuous state space have been investigated for example for power grids Yang et al. (2017); Schäfer et al. (2018). We use a paradigmatic coupled hysteresis model based on the normal form of the cusp catastrophe. Employing different network topologies such as Erdős-Rényi-, Watts-Strogatz- and Barabási-Albert-networks as well as networks generated from moisture-flow data of the Amazon rainforest, we investigate the effect of topological properties of the network. We find that networks with a large average clustering coefficient are more vulnerable to cascading tipping and discuss how this is connected to the occurence of small-scale motifs such as direct feedback and feed-forward loops. We consistently observe networks with spatial organization like small-world or the Amazon networks are more vulnerable than strongly disordered networks.

II The Model

II.1 System

In our conceptual model, a tipping element is represented by a (real) time-dependent quantity $x(t)$ that evolves according to the autonomous ordinary differential equation

[TABLE]

where $r$ is the control parameter and $a,b>0$. The parameters $a$ and $b$ control the strength of these effects respectively and $x_{0}$ controls the positon of the system on the x-axis. The equation has thus one stable equilibrium for $|r|>r_{\mathrm{crit}}$ and a bistable region for $-r_{\mathrm{crit}}<r<r_{\mathrm{crit}}$ (see the bifurcation diagram depicted in the box in Fig. 1).

We describe the characteristic behaviour of Eq. 1: If the system state is initially in the lower stable equilibrium ($x\approx 0$) and $r$ is slowly increased, eventually at $r=r_{crit}$ a tipping point is reached and a critical transition to the upper stable equilibrium ($x\approx 1$) occurs. If $r$ is afterwards decreased, the system state stays on the upper branch and only at $r=-r_{crit}$ tips down to the lower branch again. Equation 1 is a minimal model for ecosystems with alternative stable states and hysteresis Scheffer et al. (2001) but can as well be used to conceptualize other systems with similar properties such as the Thermohaline Circulation and ice sheets Stommel (1961); Levermann and Winkelmann (2016).

Next, we consider a directed network of $N$ interacting tipping elements as a linearly coupled system of ordinary differential equations

[TABLE]

where $d>0$ is the coupling strength and

[TABLE]

For simplicity, we use the same parameters $a$ and $b$ for all tipping elements in the network. An illustration of such a system with several tipping elements is depicted in Fig. 1. Similar systems have already been studied with diffusive coupling focusing on hysteresis effects Eom (2018).

We briefly review the behaviour of two tipping elements with unidirectional coupling ($X_{1}\rightarrow X_{2}$) Brummitt et al. (2015). The elements of the adjacency matrix are $a_{21}=1$ and $a_{12}=0$ which means that element 1 has an effect on element 2 but there is no effect in the other direction. As $r_{1}$ is slowly increased, it approaches its tipping point at $r_{\mathrm{crit}}$ and eventually tips from $x_{-}$ to $x_{+}$. The effective control parameter $\tilde{r}_{2}$ is thus increased by $\Delta\tilde{r}=d(x_{+}-x_{-})$. For $r_{2}=0$, a tipping event in the second element is induced if $\Delta\tilde{r}>r_{\mathrm{crit}}$ and therefore if the coupling strength exceeds a critical threshold of $d_{\mathrm{c}}=\frac{r_{\mathrm{crit}}}{x_{+}-x_{-}}$.

II.2 Network Models

To investigate the effect of the network topology on tipping cascades we use different network models: We use three well-known models, the Erdős-Rényi model (ER) Bollobas (2001), the Watts-Strogatz model (WS) Watts and Strogatz (1998) and the Barabási-Albert model (BA) Barabási and Albert (1999). We slightly extend the two last models such that we are able to generate and compare directed networks with controllable average degree ${\langle k\rangle=\langle k_{in}+k_{out}\rangle}$. Furthermore, we use models to control the reciprocity and average clustering coefficient as well as a directed configuration model and a stochastic block model. All network models are shortly discussed in the following paragraphs:

(i) The ER model is a simple random network model, where a directed link between two elements $i$ and $j$ is added with probability $p$. The resulting average degree is $\langle k\rangle\approx p(N-1)$.

(ii) The WS model is usually used to generate networks with large clustering coefficients, but small average path lengths to resemble the small-world phenomenon Milgram (1967). We implement a directed WS model as follows: Initially, a regular network is generated where each node $i$ is connected in both directions to its $m$ nearest neighbors, e.g., nodes $i+1,i-1,...,i+\frac{m}{2},i-\frac{m}{2}$. Therefore, $m$ has to be an even integer and the average degree of the resulting regular network is equal to $m$. In order to generate networks with arbitrary average degree, $m$ is chosen such that the average degree of the resulting regular network is larger than the desired average degree. Then, until the average degree of the network matches the desired average degree, links are randomly deleted. Finally, each of the remaining links is rewired with probability $\beta$, similar to the usual WS model Watts and Strogatz (1998). With increasing rewiring probability $\beta$ the generated network becomes more and more random and its properties approach the properties of ER networks for $\beta\rightarrow 1$.

(iii) The BA model is used to generate scale-free networks, i.e. networks with a power-law degree distribution. We implement a directed BA model as follows: We start with two bidirectionally coupled nodes. Every additional node is in both directions connected to an already existing node $i$ with probability $p=\frac{k^{\mathrm{in}}_{i}+k^{\mathrm{out}}_{i}}{\sum_{m,n}a_{mn}}$. When the specified network size $N$ is reached, the average degree $\langle k\rangle\approx\frac{\sum_{m,n}a_{mn}}{N}$ is compared to the desired average degree. If the average degree is smaller than the desired average degree, links between randomly selected nodes $i$ and $j$ are added with probability $p=\frac{k^{\mathrm{in}}_{i}+k^{\mathrm{out}}_{i}+k^{\mathrm{in}}_{j}+k^{\mathrm{out}}_{j}}{2\sum_{m,n}a_{mn}}$ until the average degree matches the desired average degree. Else, if the average degree is greater than the desired average degree, links are randomly deleted as in the WS model.

(iv) To generate networks with arbitrary reciprocity $R$, we initially generate an ER network where all links are reciprocal ($R=1$). Afterwards, links are randomly chosen and rewired until the desired reciprocity is achieved.

(v) The procedure to generate networks with arbitrary average clustering coefficent $\mathcal{C}$ is similar. Initially a network with only reciprocal triangles between three randomly chosen nodes is generated. Afterwards links are randomly chosen and rewired again until the desired average clustering coefficient is achieved. That way, we are able to generate networks with an average clustering coefficient between $\mathcal{C}=0.05$ and $\mathcal{C}=0.35$. Note that the reciprocity is also large for networks with a large average clustering coefficient.

(vi) A directed configuration model can be used to generate networks with arbitrary average in- and out-degree. Links are randomly assigned to node pairs where the corresponding in- and out-degree has not been reached before Newman et al. (2001).

(vii) Finally, stochastic block models (SBM) are used to generate networks with community structures. For each (directed) combination of communities there is a seperate link probability which is usually high within the community and small between two different communities Holland et al. (1983).

II.3 Simulation Procedure

We use the system given in Eq. 2 and conduct cascade simulations on different network topologies. The parameters of the equation are chosen such that $r_{\mathrm{crit}}=0.183$ and for $r=0$ the two stable equilibria are $x_{-}=0$ and $x_{+}=1$ for all elements. The resulting parameters are $a=4$, $b=1$ and $x_{0}=0.5$. Consider a network with $N$ tipping elements and a topology, that is described by the adjacency matrix $A=(a_{ij})$. Initially, $r_{i}=0$ and $x_{i}=0$ for all $i=1,...,N$. The algorithm of a cascade simulation is the following:

Choose a random starting node $m$ of the network. 2. 2.

Slowly increase $r_{m}$ ($r_{m}\rightarrow r_{m}+\Delta r$). 3. 3.

Let the system equilibrate, e.g., integrate the ODE system until $\dot{x}_{i}<\varepsilon$ for all $i=1,...,N$. 4. 4.

Check if at least one element tipped. If not jump back to step 2. Otherwise, count the number of tipped elements.

We normalize the number of tipped elements to the number of nodes that can be reached on a directed path from the starting node (the size of the out-component), where we do not take the starting node into account. We call the resulting number cascade size $L$. The ODE system was integrated with the function scipy.integrate.odeint from the scipy python package Oliphant (2007). In all simulations, $\Delta r=0.01$ and $\varepsilon=0.005$ was used. Examples of tipping cascades with size $L=1$ are shown in Fig. 2 for ER networks with different size $N$.

III Results and Discussion

III.1 Cascades on Generic Network Topologies

We start with cascade simulations on networks generated with the ER model. For any parameter combination we generate $100$ different networks and simulate one cascade on each network. We use the average cascade size from these simulations as a measure of the vulnerability of the corresponding network structure ranging from robust ($\langle L\rangle=0$) to highly vulnerable ($\langle L\rangle=1$) networks. The dependence of the average cascade size with respect to the coupling strength is shown in the upper panel of Fig. 3 for random networks with a fixed average degree $\langle k\rangle\approx 10$. For low coupling strengths ($d\lesssim 0.1$) the network is not affected by the externally induced tipping of one element and the average cascade size remains zero. With increasing coupling strength, a transition robust to vulnerable networks is observed. From the analysis of the unidirectional system, a sharp transition at $d\approx r_{\mathrm{crit}}$ would be expected for all networks. However, only for $N\rightarrow\infty$ the transition becomes more and more steep and approaches $r_{\mathrm{crit}}$. For networks of finite size, the onset of the transition is shifted to lower coupling strengths with decreasing network size. We hypothesize that the reason for this is two different effects: The first effect is the destabilization of the system by feedback loops ($X_{1}\leftrightarrows X_{2}$) which can lead to a decrease of the tipping point $r_{\mathrm{crit}}$ of certain nodes. The second effect is due to the gradual change of the state of a tipping element $X_{3}$ that is coupled to another element ($X_{1}\rightarrow X_{3}$). When the element $X_{1}$ tips, the state of the element $X_{3}$ will be slightly altered even if it does not tip. If it is coupled to another element $X_{2}$ however ($X_{2}\rightarrow X_{3}$), the effective control parameter of element $X_{3}$ will be slightly increased by an increment of the order $\Delta\tilde{r}\sim d^{2}$. Therefore an additional indirect coupling with one intermediate node, called feed-forward loop, will decrease the critical coupling strength $d_{\mathrm{c}}$ of the target node. With this we can explain the size dependency of the transition which is shown in the lower panel of Fig. 3. With increasing network size while fixing the average degree, the relative density of these motifs decreases and for $N\rightarrow\infty$, the destabilizing effect of the motifs vanishes.

To test this hypothesis, cascade simulations on networks with different reciprocities and average clustering coefficients are conducted. Simulation results for different reciprocities $R$ can be seen in the left panel of Fig. 4. As expected, for networks with high reciprocity, the transition region is shifted to lower coupling strengths. As can be seen, however, the dependence on the reciprocity is rather weak. Simulation results for networks with different average clustering coefficient $\mathcal{C}$ are shown in the right panel of Fig. 4. It can be clearly seen that the vulnerability to tipping cascades is significantly increased for high average clustering coefficients. There are eight motifs that contribute to the average clustering coefficient in a directed network, two (indirect) feedback loops and six feed-forward loops Milo et al. (2002). We suspect that the effect of indirect feedback loops is smaller than the effect of direct feedback loops for $d<1$. Therefore, we conclude that feed-forward loops are mainly responsible for the increased vulnerability of networks with large average clustering (see Fig. 4).

We also observe a transition of the average cascade size when the coupling strength is held constant at $d=0.15$ and the average degree is varied (Fig. 5). In that case the transition is shifted to higher average degrees when the network size increases, because a higher average degree is necessary to yield the same relative density of destabilizing motifs.

Cascade distributions for $\langle k\rangle\approx 5$ and selected coupling strengths at the onset, in the center and at the end of the respective transition region are shown in Fig. 6. We find a bimodal distribution of very small cascades ($L\approx 0$) and very large cascades ($L\approx 1$). For networks with small-world and scale-free topology generated with the WS model with $\beta=0.1$ and the BA model, respectively, we observe similar transitions of the average cascade size. For the scale-free topology, the large cascades are distributed around an average size $\langle L\rangle<1$. This can be explained by the preferential attachment mechanism. Through this mechanism a large number of weakly connected elements develop which can only be tipped when the coupling strength is very large ($d\gtrsim r_{\mathrm{crit}}$).

Now we focus on the effect of the network topology. For all network models, the transition from robust to vulnerable networks is shifted to lower coupling strengths, when the average degree is increased (Fig. 7). The topology of the network has a significant effect on this shift of the transition region for sparse networks ($\langle k\rangle\approx 5$). For networks with small-world and scale-free topology, the transition is shifted to lower coupling strengths compared to the simple random topology generated with the ER model. For the scale-free topology the transition width is also significantly increased for $\langle k\rangle\approx 5$. For denser networks ($\langle k\rangle\gtrsim 19$), the differences between the network topologies are less pronounced.

We further investigate in which way the rewiring in the WS model decreases the vulnerability of the network. In Fig. 8 the shift of the transition region to higher coupling strengths with respect to the rewiring probability $\beta$ can be clearly seen. The increase of the critical coupling strength mainly occurs between $\beta=0.1$ and $\beta=1$. The lower panel of the figure again demonstrates how this corresponds to the decay of the average clustering coefficient $\mathcal{C}$. Thus, we again conclude that tipping networks with an increased average clustering coefficient such as small-world networks (but also spatially structured networks Wiedermann et al. (2016), see III.2) are especially vulnerable to cascades and that the average clustering coefficient is a good indicator for the vulnerability of a network topology.

III.2 Cascades on Spatial Network Topologies from Moisture-Flow Data

To investigate the effects of spatial organization of the network on vulnerability with respect to tipping cascades, we apply our model to network topologies generated from data of atmospheric moisture flows between different forest cells in the Amazon. On a local-scale, the Amazon may exhibit alternative stable states between rainforest and savanna, with tipping points between them depending on rainfall levels Hirota et al. (2011); Staver et al. (2011); Xu et al. (2016); Ciemer et al. (2019). Models that capture the basic mechanisms also reveal a bifurcation structure with hysteresis and saddle-node bifurcations with rainfall level as control parameter, comparable to our conceptual model H. van Nes et al. (2014). On a regional scale, the forest enhances rainfall through the ”transpiration” of groundwater to the atmosphere; local-scale tipping may thus increase the vulnerability of remote forest patches by allowing less local precipitation to be passed on to other patches because the transpiration capacity of savanna is lower than that of forest. Therefore, the Amazon can be thought of as a spatial network of local-scale tipping elements. Note that the Amazon as a whole is often viewed as a tipping element Cox et al. (2004). In our framework, vulnerable regimes where tipping of single cells induces large cascades correspond to such threshold behaviour of the large-scale Amazon system. Complex-network approaches such as a cascade model inspired by the Watts-model Watts (2002) have been applied to observation-based data of Amazon forest patches Zemp et al. (2017). Here we analyze the effect of the network structure of transpired-moisture flows for the Amazon that were calculated by Staal et al. Staal et al. (2018b), aggregated to a single year (2014) on 1 degree spatial resolution.

As our analysis will be focused on the effect of the network topology, we neglect the actual moisture-flow values and use a homogeneous coupling strength analogous to the above simulations. This makes the simulation results less realistic and applicable, however, we do not aim to draw conclusions about the Amazon system. Rather, we want to compare the network topology to common random networks and identify topological effects on the vulnerability of tipping networks with respect to tipping cascades.

To generate and compare networks with arbitrary average degree, similar to the random network topologies above, we calculate a moisture-flow threshold from a specified average degree. Only when the moisture flow between two cells exceeds the threshold, these cells are connected with a link in the corresponding direction. If a large average degree is specified, the threshold becomes small and the resulting network will be dense. That way we are able to generate networks with arbitrary average degree from the data. An example network with $\langle k\rangle=5$ is depicted in Fig. 9.

The average cascade size is calculated by conducting one cascade simulation with each node of the generated network as starting node and averaging over the cascade size. We generate networks from data with $1\times 1\degree$-grid ($N=567$) and with $2\times 2\degree$-grid ($N=160$) and $\langle k\rangle=5$. The average cascade size of ER networks with the same size is shown for comparison (upper panel of Fig. 10). For the Amazon network, the onset of the transition from robust to vulnerable networks is shifted to a lower coupling strength of $d\approx 0.08$ compared to the ER network. In contrast to the ER networks there is no strong size dependency. However, a small shift to lower coupling strengths is observed.

Additionally to the Amazon moisture-flow network obtained by thresholding, we generate networks with a directed configuration model Newman et al. (2001) and a stochastic block model (SBM) Holland et al. (1983) to isolate the effects of the degree sequence and the community structure of the network, respectively. For the directed configuration model, we specify the joint degree sequence of the Amazon network. For the stochastic block model, we apply a Girvan-Newman algorithm to the original Amazon network Girvan and Newman (2002). The algorithm progressively removes edges with the highest edge betweenness, i.e., those rare links that connect seperate communities. When the network breaks into two components, we calculate the elements of the probability matrix (fraction of links over possible links for the corresponding combination of components). With the probability matrix and the component sizes, we then generate a random network with the stochastic block model.

In the lower panel of Fig. 10, the transition of the configuration model and the SBM is compared to the original Amazon network and the ER network with $N=160$. Although the vulnerability of the network is increased in both cases compared to the ER model, none of the topological properties alone, degree sequence or community structure, sufficiently explains the early onset of the transition in the original Amazon network.

Cascade distributions for the coarse resolution ($2\times 2\degree$-grid) are depicted in Fig. 11. They show that already for values of $d\approx 0.1$ cascades with two typical cascade sizes occur for the original Amazon network. With increasing coupling strength the frequency of these cascades increases and the cascade size is shifted to higher values. Comparing this observation to the network in Fig. 9 suggests that these cascades correspond to the two subclusters in the north and south-west regions of the Amazon rainforest. These subregions form clusters that are much more strongly connected than the rest of the network and are thus much more vulnerable to tipping cascades. Interestingly, seperate tipping of subclusters is not observed for the networks generated with the SBM implying that some relevant topological property of the spatially structured Amazon network, for example the anisotropy of the link direction due to atmospheric wind patterns, might still be missing.

The robust and vulnerable regimes of the networks are shown in Fig. 12. Consistent with the above results, we observe a shift of the transition to lower coupling strengths with increasing average degree $\langle k\rangle$ where the transition is smooth for the Amazon network and steep for the configuration model and the SBM. Similar to the random network topologies, the differences are only relevant for the sparse regime below $\langle k\rangle\lesssim 19$.

IV Conclusion

The aim of our study was to assess the effect of the network topology on the vulnerability of tipping networks to cascades. This is not only important for understanding the effect that the tipping of potential tipping elements in the climate system might have on the complete Earth system, but also of high relevance for other fields that use complex system approaches. We found that networks with large average clustering coefficients and spatially structured networks are more vulnerable to tipping cascades than more disordered network topologies. This implies that the risk of a cascade to be triggered could be surprisingly high for real-world networks where large clustering is common. Furthermore, we found that the effect of the network topology is relevant only for relatively sparsely connected networks. The analysis of the Amazon network suggests that the structure of the forest-climate system in the Amazon might yield subregions that are especially vulnerable to tipping cascades. A detailed study using actual moisture-flows could investigate the question if the Amazon rainforest consists of separate sub-regional-scale tipping elements. Generally, heterogeneity in the parameters, for example the temporal and spatial scales or the coupling strengths of the ODE system stated in Eq. 2, could have a further influence on the results Rocha et al. (2018).

Acknowledgements

The authors thank M. Wiedermann and J. Heitzig for fruitful discussions. N.W. acknowledges support from the the IRTG 1740/TRP 2015/50122-0 funded by DFG and FAPESP. N.W. is grateful for a scholarship from the Studienstiftung des deutschen Volkes. R.W. and J.F.D. are thankful for support by the Leibniz Association (project DominoES). A.S. and J.F.D. acknowledges support from the European Research Council project Earth Resilience in the Anthropocene (743080 ERA). A.S. and O.A.T. thank support from the Bolin Centre for Climate Research. J.F.D. thanks the Stordalen Foundation (via the Planetary Boundaries Research Network PB.net) and the Earth League’s EarthDoc program for financial support. The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Land Brandenburg for supporting this project by providing resources on the high performance computer system at the Potsdam Institute for Climate Impact Research.