1 Introduction

In the recent years, complex systems research has increasingly focused on the matter of tipping points [1,2,3] since they occur in many different systems including ecosystems, over economics, the Earth’s climate system and social systems [4,5,6,7,8,9,10]. Tipping points are the critical thresholds of tipping elements, where a small perturbation can be sufficient to invoke a qualitative change of the whole system. Whether such qualitative changes can be seen as something desirable or undesirable depends a lot on the context: for instance, a potential transition of climate tipping elements towards a potential “hothouse” state might be dangerous for humanity [11, 12], while a rapid transition towards a sustainable future lies well within the scope of desired tipping events [13]. However, oftentimes tipping elements do not exist in isolation, but interact across scales in time and space [14, 15] such as connected lakes in ecology [16, 17], in the adoption of new technologies in the economy [18] or the climate tipping elements in the Earth system [19]. Since several decades, networks are an established tool for the description of complex systems [e.g., 20, 21]. Complex networks are structures that represent certain entities as their nodes and their interaction as their edges. They have been used, for example, to model oscillators in power grids [22], food webs [23], interactions of climate system components [24] and the collaboration network of scientists [25]. Critical behaviour has also been revealed on the network level. For instance, it has been shown that the likelihood of developing diabetes depends of the criticality of excitable tissue in the Islets of Langerhans of the pancreas [26].

Since there is increasing interest in modelling interacting tipping elements within the context of complex systems [27,28,29], we bring these two strands of research together since tipping elements on networks can not only tip themselves but also imply tipping of neighbouring systems or even the network as a whole. Building upon recent developments in studying interacting nonlinear dynamics on complex networks [30,31,32,33,34], we here introduce the unified Python package PyCascades.

In Sect. 2, we describe how PyCascades can be installed and what the package contains (Sect. 2.1). Further, we describe the general structure of our package (Sect. 2.2), the building blocks of nonlinear dynamical systems, namely the tipping elements and their interaction structure (Sect. 2.3) as well as the network types natively included in the package (Sect. 2.4) and lastly, the extension to several types of stochastic tipping elements (Sect. 2.5). Thereafter, we apply our modelling framework to three different examples (Sect. 3). First, we use our model to simulate tipping cascades in the Amazon rainforest, which is connected by a network of atmospheric moisture flows (Sect. 3.1). Second, we show how PyCascades can be extended to large scale Monte Carlo ensemble studies such that many uncertainties can be propagated (Sect. 3.2). Third, we exchange the fundamental differential equation that has been used in the two earlier examples to model tipping cascades in an economic example of a global trade network (Sect. 3.3). Lastly in Sect. 4, we shortly summarise the functionalities of PyCascades.

2 Methods

This chapter describes the basic features that are supplied by PyCascades from the installation and the structure of the package to the fundamental features that have been developed. Here, a tutorial can be found that guides the interested reader through the most important first steps to simulate tipping cascades on interacting tipping elements (https://doi.org/10.5281/zenodo.4153102). Furthermore, the code for each of the following fundamental features and the three applications is provided there.

2.1 Installation and package structure

PyCascades can be installed via the command line using the pip-command

pip install pycascadesFootnote 1.

Alternatively, the package can directly be downloaded via the website following the zenodo-doi: 10.-5281/zenodo.4153102. The layout of the file structure of PyCascades can be found in Table 1. Important files, which led to the outcomes of this work, are listed and described there. A dedicated tutorial has been developed, guiding the interested reader through some important first steps and features of the software package. For the Amazon rainforest application and the climate tipping elements application, further readme-files have been added in the respective directory. There, it is explained how the respective simulations can be started and evaluated. Additionally, the plot scripts for these two applications are deposited.

Table 1 File structure of PyCascades
Fig. 1
figure 1

UML class diagram of the core of PyCascades that depicts structure and dependencies of PyCascades’ functionalities separated in the different python classes. The class tipping_network is derived from the DiGraph class of the NetworkX package [36]. It aggregates instances of the general classes tipping_element and coupling. The evolve class is associated with one instance of the tipping_network class and simulates the evolution of the complex dynamical system which is implemented by the concrete tipping_element and coupling objects with their specific parameters. For simplicity, only classes and class members important to the understanding of the PyCascades core are shown.

2.2 Structure of the core of PyCascades

PyCascades provides a convenient framework to solve differential equations on complex networks, i.e., it describes the dynamics of states of nodes in such a network as well as their interactions. The basic assumption is that the dynamics of tipping elements can be separated into one part for the isolated dynamics of the tipping element and another part representing the interaction terms (see Sect. 2.3 for more details). For that, it builds on SciPy differential equation solvers [35] for the dynamics and on NetworkX [36] to generate the underlying network.

Fig. 2
figure 2

The bifurcation diagrams of the pre-implemented tipping elements are shown, which possess a Cusp- or a Hopf-bifurcation, respectively (a and b). Further, two examples of tipping events are shown, once for the case where no tipping cascade emerges (c and d), and once where a tipping cascade emerges (e and f). a Bifurcation diagram of a fold-bifurcation (see Eq. 3). b Bifurcation diagram of a (supercritical) Hopf-bifurcation (see Eq. 4). c Two cusp differential equations with the parameters \(a_1 = a_2 = 4\), \(b_1 = b_2 = 1\), \(x_{0,\ 1} = x_{0,\ 2} = 0.5\) and \(c_1 = 0.2\), \(c_2 = 0\). Thus, the first tipping element is slightly pushed over its upper critical value and a state transition occurs. d One cusp, initialised as in c, and one tipping element that possibly could undergo a Hopf-bifurcation (\(a = 1\), \(\mu = -1\)). In c and d, there is no interaction between the tipping elements, i.e., \(A_{ij} = 0\) \(\forall i,j\). e and f Same as in c and d, but with \(A_{21} = 0.5\) such that the second tipping element (Cusp-2, Hopf) is coupled to the state of the first element (Cusp-1, Cusp). Therefore, in the lower panels a tipping cascade of twofold bifurcations or, respectively, a tipping cascade of one-fold- and one Hopf-bifurcation can be observed

The core of PyCascades is structured as follows (see Fig. 1). It provides the two classes tipping_element and coupling that implement the two described types of dynamics. From these classes that can be viewed as references, concrete classes for tipping elements and interactions can be derived. Currently, PyCascades provides the classes cusp and hopf derived from tipping_element and linear_coupling derived from coupling. Other types of tipping elements or couplings can be implemented in an analogous way. The class tipping_network which is derived from the DiGraph class of NetworkX is used to combine different tipping_element and coupling objects into a network and identify each tipping_element object with a node and each coupling object with a link. Finally, an evolve class is provided with methods to integrate the resulting ODE system or to trigger tipping events.

2.3 Different types of tipping elements and interactions

Through the tipping_element class in PyCascades different types of tipping elements can be defined and coupled together. Each tipping element can be described by its individual dynamics \(f_i\) and the interaction term \(g_i\), i.e., the coupling to other tipping elements. This yields

$$\begin{aligned} \begin{aligned} \tau _i \frac{\text {d}x_i}{\text {d}t}&= f_i(x_i) + g_i(x), \end{aligned} \end{aligned}$$
(1)

where \(x_i\) represents the state of the respective tipping element. \(\tau _i\) stands for a typical timescale of tipping. The direct interaction term \(g_i(x)\) is assumed to be separable into the summands

$$\begin{aligned} g_i(x) = \sum _j g_{ij}(x_i,x_j), \end{aligned}$$
(2)

linking the tipping elements i and j.

In principle, any kind of tipping element can be supplied in the tipping_element class of PyCascades, but as of now, there are two kinds of tipping elements predefined that are ready to be used and implemented. These two tipping elements are elements that possess a Cusp-bifurcation or a Hopf-bifurcation [37]. The first pre-implemented tipping element is the Cusp-differential equation, which has been used in many contexts before to model nonlinear transitions between two alternative stable states [15, 38]. The normal form of its differential can be written as

$$\begin{aligned} f_\text {Cusp} (x) = \frac{\text {d}x}{\text {d}t} = -a\left( x-x_0\right) ^3 + b\left( x-x_0\right) + c. \end{aligned}$$
(3)

Here, \(a, b > 0\) and \(x_0\) represents a shift on the x-axis. The parameter c is the critical parameter, which invokes a shift from a lower stable state to an upper stable state as soon as the critical value \(c_\text {crit, high}\) is surpassed. The other way round, when c is diminished, a state transition from the upper to the lower stable state occurs at \(c_\text {crit, low}\). Equation 3 has the normal form of a fold-bifurcation and has, as a paradigmatic model, been applied in many different areas such as systems in ecology, climate science and economics [15, 29, 31, 39, 40]. For the special case that \(a=4, b=1\) and \(x_0=0.5\), the two stable states are located at \(x_{1}= 0\) and \(x_{2} = 1\) for \(c=0\). The critical parameter lies at \(c_\text {crit, high} = - c_\text {crit, low} = \sqrt{(4b^3)/(27a)} = \sqrt{4/(27\cdot 4)} \approx 0.19\). The bifurcation diagram of this equation is shown in Fig. 2a.

Fig. 3
figure 3

The three paradigmatic and pre-implemented network types are shown together with an example of a tipping cascade in this network. Exemplary network structure with 15 nodes and an average degree of 3 for a an Erdős–Rényi network, b a Watts–Strogatz network (rewiring probability: 0.15) and c a Barabási–Albert network. Exemplary tipping cascades for a network of 15 nodes, an average degree of 3, where each node is represented by a Cusp-tipping element (see Eq. 3) with \(a=4\), \(b=-1\), \(x_0=0.5\) for all nodes. The couplings between the nodes are alternately set to 0.2 and −0.4 for interaction 1, 2, 3, etc. Then at \(t=0\) one randomly chosen node i is set to the upper state by choosing \(c_i=0.2\) such that tipping cascades can emerge. The number of tipped elements are shown for d an Erdős–Rényi network, e a Watts–Strogatz network (rewiring probability: 0.15) and f a Barabási–Albert network

The second tipping element that is provided by PyCascades is a Hopf-bifurcation. The normal form in polar coordinates of this bifurcation can be written as

$$\begin{aligned} \begin{aligned} f_\text {Hopf, r}(r) = \frac{\text {d}r}{\text {d}t}&= \left( \mu - r^2\right) \cdot r a \\ f_{\text {Hopf, } \phi } (\phi ) = \frac{\text {d} \phi }{\text {d}t}&= b \end{aligned} \end{aligned}$$
(4)

with the parameters a and b. Here, the Hopf-bifurcation is given in polar coordinates with the radius r and the angle \(\phi \). Importantly, \(\mu \) is the critical parameter and a bifurcation from a stable fixed point to a stable limit and an unstable fixed point occurs when \(\mu \) crosses zero from below. The bifurcation diagram is shown in Fig. 2b. Applications of Hopf-bifurcations have been found, for instance, in predator–prey cycles in Lotka–Volterra systems or in the Hodgkin–Huxley model of neurons [41, 42]. In the climate system, there exist conceptual models that represent the El-Niño Southern Oscillation as a Hopf bifurcation [43, 44] based on a model by Zebiak and Cane [45].

Next, for the interactions, any type of coupling can in principle be used and implemented in PyCascades. However, for the moment, only linear interactions are considered

$$\begin{aligned} g_i(x) = \sum _{j=1}^N A_{ij} x_j. \end{aligned}$$
(5)

If there is a connection between tipping element i and j, then \(A_{ij} \ne 0\), otherwise \(A_{ij} = 0\). In Fig. 2c–f, we show an example how tipping cascades can emerge from the coupling between two tipping elements for the case of two cusp-differential and for the case of one cusp coupled to the normal form of a Hopf-bifurcation.

Fig. 4
figure 4

Noise-induced tipping events with respect to different types of noise (Gauss, Lévy, Cauchy) that are available in PyCascades. Simulation of one Cusp-tipping (\(a=4\), \(b=1\), \(c=x_0=0\)) element with stochastic noise (see Eq. 6) of the following type: a Gaussian noise (\(\sigma =0.25\), see Eq. 7), b Lévy noise (\(\sigma = 1.0\), see Eq. 7) and c Cauchy noise (\(\sigma = 1.0\), see Eq. 7)

2.4 Paradigmatic network types of interacting tipping elements

With PyCascades it is possible to investigate the dynamics of tipping and tipping cascades in larger directed networks. These types of networks can either be explicitly spatially embedded (see Sect. 3) or well-known predefined network models such as the Erdős–Rényi model, the Barabási–Albert model or the Watts–Strogatz model [46,47,48]. Originally, the network models that are inbuilt in python’s network package NetworkX are undirected for Watts–Strogatz networks and Barabási–Albert networks, while we require directed networks. Additionally, it might also be helpful or necessary to be able to determine a certain average degree. Therefore, a generalisation of these two networks types has been developed. (i) Watts–Strogatz network: a regular network is created where each node i is connected to its m closest neighbours in both directions. m must be an even integer and the average degree \(\left\langle k \right\rangle = m\). m is chosen in such a way that the average degree of the resulting network is larger than the desired average degree. Then links are randomly deleted until the desired average degree is reached. Lastly, each of the remaining links is rewired with the desired rewiring probability as in the usual Watts–Strogatz model. (ii) Barabási–Albert model: first, two nodes are bidirectionally coupled. Each further node is, again, bidirectionally coupled to one already existing node i with the probability \(p = (k^{\text {in}}_{i} + k^{\text {out}}_i)/(\sum _{mn} a_{mn})\), where \(k^{\text {in}}_{i}\) is the in-degree and \(k^{\text {out}}_i\) is the out-degree of node i. With \(a_{mn}\), the sum of all edges in the network is denoted. In the end, the average degree \(\left\langle k \right\rangle \) depends on the stochastic network. Therefore, it can happen that the actual average degree \(\left\langle k \right\rangle \) of the network exceeds or falls below the desired average degree \(k_\text {des}\). While \(\left\langle k \right\rangle < k_\text {des}\), another link is added between two randomly selected nodes i and j. While \(\left\langle k \right\rangle > k_\text {des}\), a link between two randomly selected nodes i and j is deleted. For comparison of the construction of these network models, see also Krönke et al. [30]. Examples for a realisation of these three network types and an exemplary tipping cascade in those can be found in Fig. 3.

2.5 Stochasticity in tipping elements

In the real world, systems often underlie fluctuations, which under certain circumstances can cause critical transitions, called noise-induced tipping. Numerous prominent examples can be found in dynamical systems such as in electronics, optics or neurons, but also in ecology and in the Earth’s climate system [49,50,51,52,53]. Therefore, we decided to create a class for a stochastic version of the cusp tipping element (Eq. 3) for additive noise

$$\begin{aligned}&f_\text {Cusp, stoch.} \nonumber \\&\quad = {\text {d}}x = \left[ -a\left( x-x_0\right) ^3 + b\left( x-x_0\right) + c \right] {\text {d}}t + \sigma {\text {d}}W. \nonumber \\ \end{aligned}$$
(6)

Here, \(\sigma \) denotes the noise level and \({\text {d}}W/{\text {d}}t\) describes the Wiener process or Brownian motion. In the case of random white noise (Gaussian white noise) as used here, W is sampled from a Gaussian distribution. To implement stochastic differential equations, python’s SciPy function odeint has been replaced by sdeint [54]. sdeint has several algorithms implemented, which are able to solve stochastic differential equations. Here, and in the provided version of PyCascades, an order 1.0 strong stochastic Runge–Kutta algorithm is employed [55].

Furthermore, Gaussian noise distributions are not necessarily able to describe all types of fluctuations in real-world systems since in reality noise might be correlated or not be standard normally distributed. Besides Gaussian noise, PyCascades allows to compute systems with Lévy and Cauchy noise (see Fig. 4). These types of noise (Lévy, Cauchy) may be more suitable for describing extreme events than Gaussian noise; however, in the implemented form, they still remain uncorrelated. It has been found that the probability of jumping between the two stable states in a double-well potential is impacted by single strong extreme events from those \(\alpha \)-stable noise distributions [56]. For instance, it has been proposed that this might have been of relevance for climate system states on a millennial time scale during the last glacial period as was observed in ice-cores [57]. Also, transitions triggered by extreme events emerging from Lévy-distributions in other nonlinear climate system components such as the Amazon rainforest or the thermohaline circulation have been investigated, as well as transitions in gene expression processes in molecular biology [58,59,60].

Fig. 5
figure 5

Tipping cascades in a conceptual model of the Amazon rainforest connected via an atmospheric moisture recycling network. a Sketch of the network of interacting rainforest patches in the Amazon rainforest. Precipitation rains down over some parts of the rainforest and parts of it are re-evapotranspirated and transported further by the wind (atmospheric moisture transport). b–d Exemplary tipping experiment on a 0.5\(\times 0.5^\circ \) grid, where each grid cell represents one rainforest patch. The colour bar represents the likelihood of tipping averaged over the years 2003–2014. We show a comparison between b coupling switched on (see Eq. 8), c coupling switched off (see Eq. 8 with \(g_i(x) \equiv 0\) \(\forall i\)) and d the difference between the panels b and c. For each year in the study period (2003–2014), we performed one such tipping experiment, and the results shown are an average over this period

The distributions for Gaussian, Lévy and Cauchy noise in PyCascades are taken from python’s SciPy libraries

$$\begin{aligned} \begin{aligned} p_\text {Gauss}(x)&= \frac{1}{\sqrt{2 \pi \sigma ^2}} \cdot \exp \left( - \frac{x^2}{2 \sigma ^2} \right) \\ p_{\text {L}\acute{\mathrm{e}}\text {vy}}(x)&= 0.5\cdot p_{\text {L}\acute{\mathrm{e}}\text {vy},\ \text {pos.}} + 0.5 \cdot p_{\text {L}\acute{\mathrm{e}}\text {vy},\ \text {neg.}} = \\&= 0.5 \cdot \sqrt{\frac{\sigma }{2 \pi x^3}} \cdot \exp \left( -\frac{\sigma }{2x} \right) \\&\quad + 0.5 \cdot \sqrt{\frac{\sigma }{2 \pi \left| x \right| ^3}} \cdot \exp \left( -\frac{\sigma }{2\left| x \right| } \right) \\ p_\text {Cauchy}(x)&= \frac{1}{\sigma \pi } \cdot \frac{\sigma ^2}{\sigma ^2 + x^2}, \end{aligned} \end{aligned}$$
(7)

where \(\sigma \) is the standard deviation and the mean value \(\mu \equiv 0\).

3 Applications

In this section, we show three examples of how PyCascades can be applied to real-world systems. The first application is the moisture-recycling network of the Amazon rainforest, where we introduce PyCascades on a spatially embedded network. In the second application in a subset of interacting climate tipping elements, we combine the PyCascades modelling framework with a large-scale setup of Monte Carlo simulation to show how numerous uncertainties in parameters can be propagated systematically. The third application, a global trade network of more than 5000 nodes and 400,000 links, complements our analysis by simulating tipping cascades with a modernised, economically motivated differential equation (see Eq. 10) replacing the Cusp-differential equation (see Eq. 3).

3.1 The Amazon rainforest

It is suspected that the Amazon rainforest is a tipping element in the Earth’s climate system [4], which might approach a tipping point due to various anthropogenic pressures including climate change, fires and land-use change [61,62,63]. The Amazon rainforest might exhibit multistability at certain rainfall levels, as suggested by conceptual models and observational data [64,65,66,67,68]. This implies that rainforest patches may transition to a savannah-like state when the rainfall drops below a certain critical level. These rainforest patches depend on each other, as rain is re-evaporated by the trees and thus preserved in the system through atmospheric moisture recycling [69, 70] (see Fig. 5a). This means that the Amazon rainforest is an excellent example of how tipping cascades can travel through a system, which can be modelled with PyCascades. We divide the Amazon into 0.5\(\times 0.5^\circ \) (approximately 50 km) grid cells and assume that each is an interacting tipping element that can be described by Eqs. 3 and 5. For simplicity, we chose \(a_i=b_i=1\) and \(x_{0,i}=0\) for all tipping elements and further assume that the critical parameter is only dependent on the rainfall a rainforest cell receives, which tips in case the received rainfall is less than the critical rainfall. Then the critical parameter \(c_i\) and the coupling \(g_i(x)\) can be denoted as

$$\begin{aligned} \begin{aligned} c_i&= c_0 \cdot \frac{R_i - \left\langle R \right\rangle }{R_\text {crit} - \left\langle R \right\rangle } \\ g_i(x)&= \frac{1}{2} \frac{c_0}{R_\text {crit} - \left\langle R \right\rangle } \sum ^N_{j=1} \delta ^\text {Rain}_{ij}. \end{aligned} \end{aligned}$$
(8)

Here, \(R_i\) is the rainfall in cell i, \(\left\langle R \right\rangle \) is the average rainfall over the whole Amazon basin and \(R_\text {crit}\) is the critical rainfall. Further, \(c_0 = \sqrt{4/27}\) if \(a=b=1\) and \(x_0=0\). Lastly, \(\delta ^\text {Rain}_{ij}\) is the moisture transport in mm/yr from cell j to cell i. Since the distance between the two stable states is 2, a prefactor of 1/2 is required to re-normalise the coupling. We choose the critical rainfall \(R_\text {crit}\) to be 1700 mm/year for all cells, which is approximately the value below which the alternative savannah state becomes more resilient than the rainforest state [71]. The atmospheric moisture recycling simulations used in this work were performed by Staal et al. [72] for the years 2003–2014 and assembled into a network by Krönke et al. (2020) [30]. In this simplified example, we assume that if a forested grid cell tips, moisture recycling via that cell stops. We performed a tipping experiment for each year between 2003 and 2014 and averaged the results over this period. We find tipping events in several parts of the Amazon basin which cascade to other forest patches (Fig. 5b–d). This analysis illustrates how PyCascades can be applied to simulate tipping events and cascades in a real-world network of interacting tipping elements.

Fig. 6
figure 6

Construction of the simulation ensemble to cover a large part of the phase space. a Uncertainties in the critical temperatures [73] and b interaction strengths [19] are put into a Latin hypercube sampling algorithm [79] to construct suitable initial conditions c that cover a larger part of the state space than normal random sampling would

3.2 Climate tipping elements

Apart from the Amazon rainforest, there exists a range of processes and systems in the Earth’s climate system that exhibit threshold behaviour [4]. These tipping elements contain biosphere components (e.g. Amazon rainforest, coral reefs), large-circulation patterns (e.g. Atlantic Meridional Overturning Circulation, monsoon systems) or cryosphere components (e.g. Greenland Ice Sheet, West Antarctic Ice Sheet). Under ongoing global warming, many of them are at risk of transitioning into an alternative, tipped state at lower levels of global warming than previously though [12, 73]. Such transitions would have dangerous consequences for humanity and biosphere integrity in the Earth system [11, 12]. There is an additional risk that tipping elements are strengthened by reinforcing, positive feedbacks within the climate system such that cascades might be triggered, potentially up to a planetary-scale tipping cascade that could push the Earth towards a “hothouse” state [11]. Moreover, the tipping elements in the climate system are interacting and there is a subset of five tipping elements where the interaction structure has been made explicit by an expert elicitation [19]. This network and their interactions have been used by some studies to investigate the risk of tipping cascades in the climate system, but also to quantify economic damages exerted by interacting tipping elements [27, 33, 74].

Here, we show how PyCascades can be used to simulate tipping events in four of the five aforementioned tipping elements: the Greenland Ice Sheet (GIS), the West Antarctic Ice Sheet (WAIS), the Atlantic Meridional Overturning Circulation (AMOC) and the Amazon Rainforest (AR). For these four tipping elements, there exist conceptual models of their nonlinear behaviour with respect to a forcing parameter [64, 66, 75,76,77], which can be traced back to increases in levels of global warming above pre-industrial [73]. Therefore, we can arguably model these four elements by

$$\begin{aligned}&\frac{{\text {d}}x_i}{{\text {d}}t} = \left[ \underbrace{- x_i^3 + x_i + \frac{\sqrt{4/27}}{T_\text {crit, i}} \cdot \varDelta \text {GMT}}_{\text {Individual dynamics term}} + \underbrace{\sum _{\begin{array}{c} j\\ j \ne i \end{array}} \frac{s_{ij}}{4} \left( x_j + 1 \right) }_{\text {Interaction term}} \right] \frac{1}{\tau _i} \nonumber \\&\text {with } i=\left\{ \text {GIS, WAIS, AMOC, AR} \right\} .\nonumber \\ \end{aligned}$$
(9)

Here, \(\varDelta \)GMT is the increase of the global mean temperature, \(T_\text {crit, i}\) the critical temperature at which a certain tipping element transgresses its baseline state, \(s_{ij}\) the interaction strength between the tipping elements and \(\tau _i\) the time a certain tipping event needs. Each \(s_{ij}\) has a certain physical meaning, for instance, the freshwater entry from the GIS weakens the AMOC, while a weaker AMOC cools the northern hemisphere at the same time [78, e.g.]. The typical tipping time scales \(\tau _i\) are chosen to be 4900, 2400, 300 and 50 years at 4 \(^\circ \)C above pre-industrial levels of global warming for GIS, WAIS, AMOC and AR, respectively. For more details, see Fig. 6 and please be referred to Wunderling et al. [31]. The parameter uncertainties and a potential interaction structure are shown in Figs. 6 and 7. In Eq. 9, there are many parameters with uncertainties, for instance at which temperature \(T_\text {crit, i}\) a critical transition occurs or how strong the interactions \(s_{ij}\) are. While upper and lower limits are given in the literature [19, 73], their uncertainty need to be propagated thoroughly. For this purpose, we use a large-scale Monte-Carlo ensemble based on the latin hypercube sampling (LHS) method pyDOE [79]. The LHS is a sampling method that generates initial conditions that can be used in a Monte Carlo ensemble. They cover the state space of all uncertain parameters to a higher degree than random sample generation and are, therefore, better suited to create Monte Carlo ensembles in higher dimensional systems. In Figs. 6 and 7, we demonstrate that constructing a large-scale Monte Carlo ensemble can be combined with simulating tipping cascades with PyCascades. In the critical temperatures \(T_\text {crit, i}\) and the interaction strengths \(s_{ij}\) are 11 parameters with uncertainties (see Fig. 6). Upon that, we construct an ensemble of 1000 initial conditions.

Fig. 7
figure 7

Tipping cascades in a conceptual model of climate tipping elements. a Interaction structure of four tipping elements (Greenland Ice Sheet: GIS, West Antarctic Ice Sheet: WAIS, Atlantic Meridional Overturning Circulation: AMOC, AR: Amazon rainforest). Destabilising interactions are depicted in red, stabilising interactions are depicted in light blue and uncertain interactions are shown in black. b Likelihood for the respective tipping element to transgress its stable branch computed by running Eq. 9 into equilibrium. The error bar shows the standard deviation arising from the nine different possibilities of constructing the network. There are two uncertain links since their direction of interaction is unclear, meaning they could be stabilising, destabilising or zero, i.e., \([-,0,+]\). Permuting these three options, gives nine different network structures and for each of them, we simulate 1000 ensemble members. We chose a scenario, where \(\varDelta \text {GMT} = 2\) \(^\circ \)C above pre-industrial levels

In Eq. 9, we assume that the interaction term is 25% as important as the individual dynamics term. Thus, the interaction strength \(s_{ij}\) is divided by 4 in Eq. 9. While this poses a hypothetical scenario, it allows us to estimate the likelihood of tipping of certain element at a certain increase of the global mean temperature \(\varDelta \text {GMT}\). For 2 \(^\circ \)C, we find that the likelihood of tipping is around 50% for the GIS and WAIS, while it is significantly lower for the AMOC (around 25%) and the AR (less than 5%). There is a relatively high likelihood that GIS and WAIS tip since their critical temperature is lowest and there is a strong interaction link from GIS to WAIS. Therefore, the likelihood of tipping is lower for the AMOC, but the uncertainty is higher due to the strong negative feedback loop with GIS. Lastly, the AR has a very low likelihood of tipping since it is only connected to the other tipping elements via one uncertain link from AMOC.

3.3 International trade network

In the third example, we apply the PyCascades framework of interacting tipping elements to the International Trade Network. We construct the network from the EORA multi-regional input–output (MRIO) database [80] as also done in other studies [81, 82]. The database, which has also been subject to static analyses [83], consists of 188 countries with 27 economic sectors each, and includes the annual monetary flows between these sectors and regions. We interpret the individual sectors in each country as nodes of a network, and the flow \(f_{ij}\) in the MRIO table as the weight for each directed link from node j to i. In our analysis, we use the data for the year 2012. Following previous analyses [81, 82, 84], we also use a threshold of 10\(^6\) US-$ such that we exclude unrealistically small flows.

Propagation of economic losses on the trade network has previously been studied, for instance, with the Acclimate model [84]. This model interprets the economic sectors in each country as firms producing a commodity specific for the respective sector. Each firm does so using other commodities as inputs with specific, fixed proportions according to a Leontief production function [85] as also used in simpler input–output models. These fixed proportions are taken from the multi-regional input–output (MRIO) table underlying the construction of the trade network, which constitute the baseline state (untipped state) of the model. If, for instance, the transportation sector in a country receives an input of ten billion US-$ from the oil sector and 90 billion US-$ from the machinery sector, it might have an output of 110 billion US-$ according to the MRIO table such that the created surplus would be 10 billion US-$. However, the according sector always produces according to these proportions using nine times as much “machinery commodity” as “oil commodity”, and produces ten percent more “transportation commodity” than the sum of its inputs. If a firm receives only a certain fraction of the baseline state of a commodity due to some perturbation, the firm’s output is limited to the same fraction of the baseline output. However, in the Acclimate model firms have idle capacities, i.e. the ability to temporarily produce more than their baseline output, if they have the necessary inputs and demand is high. The dynamics of this anomaly model are focused on perturbations around the baseline state with each agent aiming for maximum profit and minimum costs under local circumstances. After a shock or perturbation ceases the model returns back to the baseline state, which constitutes an equilibrium of the model’s dynamics.

Tipping is not at the centre of the Acclimate model whose scope, as an anomaly model, vanishes for very large perturbations such as bankruptcies. We here, thus, define a simple dynamic for tipping on the trade network while keeping the linear Leontief production assumption for small perturbations, whereas nonlinear dynamics are assumed for larger perturbations. Nonlinear behaviour and tipping is common in economic networks, for instance in the banking sector [5, 86]. The nodes in the trade network are only to be perceived as representative firms, i.e., aggregates of national sectors, which consist of a variety of connected actors—so each node represents a network itself and might show nonlinear as well as tipping behaviour. The standard form of a tipping element defined by the Cusp-differential equation (see Eq. 3) is not well suited for this purpose. Instead, we are here looking for a new differential equation with the following properties:

  1. 1.

    The state \(x_i\) of a node i should represent its productivity that is between 0 (no production) to 1 (full production).

  2. 2.

    The element should react almost linearly to small perturbations as in the Acclimate model.

  3. 3.

    For large perturbations there should be a collapse of the productivity, including tipping and hysteresis.

To meet these criteria, we propose the differential equation

$$\begin{aligned} \begin{aligned} \frac{{\text {d}}x_i}{{\text {d}}t}&= r_i - x_i - a\sqrt{x} \cdot \exp \left( -bx_i\right) \\&\quad + w_\text {log} x_i \cdot \left( 1-x_i \right) , \end{aligned} \end{aligned}$$
(10)

where a and b are parameters and \(r_i\) is the limiting relative input as in the Leontief production function. The bifurcation diagram is given in Fig. 8a. The first two terms in equation 10 represent a linear response to perturbations, similar to the Acclimate model (see the dotted line in Fig. 8a). The third term is responsible for the nonlinear behaviour, causing tipping and hysteresis (see the dashed line in Fig. 8a). However, an economic tipping element defined by these three terms would be inherently unstable. Even small perturbations would finally lead to a collapse of the node since perturbations are almost always growing due to the structure of the network. However, we know that the trade network is not that fragile. Therefore, we add a logistic growth term to the differential equation to stabilise the network on the individual node level with the weight \(w_\text {log}\). Here, we argue that a certain flexibility in substituting inputs exists. Within limits, it is therefore possible to return to the original production due to the logistic growth term. To illustrate an example, we choose \(a = 4, b = 10\) for the parameters in Eq. 10. Therefore, the two bifurcation points lie at \(r_1=0.4\) and \(r_2=0.6\). The strength of the logistic growth term is chosen as \(w_\text {log} = 0.2\) (see the blue line in Fig. 8a).

Fig. 8
figure 8

Tipping cascades in a conceptual model of interacting countries in the International Trade Network. a Bifurcation diagram for the proposed economic tipping element defined by Eq. 10. b Tipping cascade between three different countries starting at the red start node at \(t=0\). The colour of the nodes represents the time at which a certain sector in that country tips and the colour of the arrows indicate the targeted country

To calculate the input term represented by \(r_i\), every flow is normalised to the sum of flows from nodes of that sector to the target node. So the new weight \(w_{c,s\rightarrow i}\) of a link from sector s in country c to node i is given as

$$\begin{aligned} w_{c,s\rightarrow i} = \frac{f_{c,s\rightarrow i}}{\sum _{k} f_{k,s\rightarrow i}}. \end{aligned}$$
(11)

With this, we can write the coupling term for the differential equation as

$$\begin{aligned} r_i = \min _{s \text { in sectors}}\left\{ \sum _{c\text { in countries}} w_{c,s\rightarrow i} \cdot x_{c,s}\right\} . \end{aligned}$$
(12)

To simulate cascades, we start with all nodes in the untipped state, here \(x_i=1\) for all nodes i. We select a random starting node and tip it by setting its productivity to zero and then evolve the system with PyCascades. We exemplify this for a cascade between three countries, where one node has been tipped artificially (see Fig. 8b). The graph illustrates how the cascade propagates within and across the different countries forming densely connected network communities. Once the cascade reaches a country, most of that country’s nodes tip almost at the same time. However, this gradual and sequential propagation of a tipping event is only one pattern of cascading behaviour observed. In Fig. 9, we show cascades for 30 different start nodes, chosen such that a wide range of different tipping cascades can be observed. Figure 9a shows the number of tipped nodes, and Fig. 9b the average node state \(\langle x\rangle = \frac{1}{N}\sum _{i}^N x_i\).

Fig. 9
figure 9

Tipping cascades with the economic tipping element defined by equation 10 on the normalised trade network. The cascade depicted in Fig. 8 is plotted in red. a Number of tipped nodes over the course of the simulation time. Cascades either lead to tipping of almost all nodes (global tipping) or the dynamics of the tipping cascade stops growing after very few nodes are tipped. b Average node state \(\langle x\rangle \) over the course of the simulation time. At the onset of the tipping cascade, the average stage shows a sharp decline. Since a decline in the average state does not automatically mean that nodes were tipped, some average state time series stabilise while others show a global cascade

4 Conclusion

In this work, we have outlined the software package PyCascades, which is designed for simulating nonlinear dynamics, in particular tipping behaviour of interacting systems. For that purpose, two different types of tipping elements (Cusp and Hopf-bifurcation type models) are provided in PyCascades as well as different paradigmatic complex network types (Erdős–Rényi, Barabési–Albert, Watts–Strogatz networks) and a stochastic version of the tipping elements supplying Gaussian, Lévy and Cauchy noise. PyCascades is written in the programming language Python and is written with an object-oriented architecture such that it remains flexible and can easily be adapted or extended to further applications or theoretical problems.

However, a distinct limitation is, as of now, that only systems can be investigated, where the individual dynamics part can be separated from the interaction part. We also suspect that there is considerable potential for improvement in some technical details. For instance, more interaction types or multiplicative noise could be implemented. Another distinct constraint of PyCascades is that only paradigmatic dynamics of tipping elements are implemented. In particular, it would be highly desirable to develop process-based tipping elements depending on the respective application.

All in all, due to modular setup, PyCascades has the potential to contribute to relevant questions about the emergence of tipping cascades in various contexts, ranging from economics, ecology, climate science and beyond.