Modelling nonlinear dynamics of interacting tipping elements on complex networks: the PyCascades package

Tipping elements occur in various systems such as in socio-economics, ecology and the climate system. In many cases, the individual tipping elements are not independent from each other, but they interact across scales in time and space. To model systems of interacting tipping elements, we here introduce the PyCascades open source software package for studying interacting tipping elements (doi: 10.5281/zenodo.4153102). PyCascades is an object-oriented and easily extendable package written in the programming language Python. It allows for investigating under which conditions potentially dangerous cascades can emerge between interacting dynamical systems, with a focus on tipping elements. With PyCascades it is possible to use different types of tipping elements such as double-fold and Hopf types and interactions between them. PyCascades can be applied to arbitrary complex network structures and has recently been extended to stochastic dynamical systems. This paper provides an overview of the functionality of PyCascades by introducing the basic concepts and the methodology behind it. In the end, three examples are discussed, showing three different applications of the software package. First, the moisture recycling network of the Amazon rainforest is investigated. Second, a model of interacting Earth system tipping elements is discussed. And third, the PyCascades modelling framework is applied to a global trade network.

Abstract.Tipping elements occur in various systems such as in socioeconomics, ecology and the climate system.In many cases, the individual tipping elements are not independent from each other, but they interact across scales in time and space.To model systems of interacting tipping elements, we here introduce the PyCascades open source software package for studying interacting tipping elements (doi: 10.5281/zenodo.4153102).PyCascades is an object-oriented and easily extendable package written in the programming language Python.It allows for investigating under which conditions potentially dangerous cascades can emerge between interacting dynamical systems, with a focus on tipping elements.With PyCascades it is possible to use different types of tipping elements such as double-fold and Hopf types and interactions between them.PyCascades can be applied to arbitrary complex network structures and has recently been extended to stochastic dynamical systems.This paper provides an overview of the functionality of PyCascades by introducing the basic concepts and the methodology behind it.In the end, three examples are discussed, showing three different applications of the software package.First, the moisture recycling network of the Amazon rainforest is investigated.Second, a model of interacting Earth system tipping elements is discussed.And third, the PyCascades modelling framework is applied to a global trade network.

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 Langerhans Isles 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 Chapt.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 (Chapt.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 Chapt.4, we shortly summarise the functionalities of PyCascades.

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 (doi: 10.5281/zenodo.4153102).Furthermore, the code for each of the following fundamental features and the three applications is provided there.

Installation and package structure
PyCascades can be installed via the command line using the pip-command pip install pycascades==1.0.11 .
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 Tab. 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.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.
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.

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 where x i represents the state of the respective tipping element.τ i stands for a typical timescale of tipping.The direct interaction term g i pxq is assumed to be separable into the summands g i pxq " 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 preimplemented 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 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 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 crit, low .Eq. 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 crit, high " ´ccrit, low " a p4b 3 q{p27aq " a 4{p27 ¨4q « 0.19.The bifurcation diagram of this equation is shown in Fig. 2a.
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 with the parameters a and b.Here, the Hopf-bifurcation is given in polar coordinates with the radius r and the angle φ.Importantly, µ is the critical parameter and a bifurcation from a stable fixed point to a stable limit and an unstable fixed point occurs when µ 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 & Cane (1987) [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 If there is a connection between tipping element i and j, then A ij ‰ 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.

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 Chapt. 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 k " 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 " pk in i `kout i q{p ř mn a mn q, where k in i is the in-degree and k 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 k depends on the stochastic network.Therefore, it can happen that the actual average degree k of the network exceeds or falls below the desired average degree k des .While k ă k des , another link is added between two randomly selected nodes i and j.While k ą k 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.(2020) [30].Examples for a realisation of these three network types and an exemplary tipping cascade in those can be found in Fig. 3.

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 Here, σ denotes the noise level and dW {dt 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 doublewell potential is impacted by single strong extreme events from those α-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].
where σ is the standard deviation and the mean value µ " 0.

Applications
In this section, we show three examples of how PyCascades can be applied to realworld 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 5 000 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).

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ˆ0.5 ˝(appr.50 km) grid cells and assume that each is an interacting tipping element that can be described by the 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 pxq can be denoted as Here, R i is the rainfall in cell i, R is the average rainfall over the whole Amazon basin and R crit is the critical rainfall.Further, c 0 " a 4{27 if a " b " 1 and x 0 " 0. Lastly, δ 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 crit to be 1700 mm/yr 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. (2018) [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.

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), largecirculation 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 Here, ∆GMT is the increase of the global mean temperature, T 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 τ 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 [e.g.78].The typical tipping time scales τ i are chosen to be 4900, 2400, 300 and 50 years at 4 ˝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. ( 2020) [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 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 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.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 ∆GMT.For 2 ˝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.

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 [86,87].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) The state x i of a node i should represent its productivity that is between 0 (no production) to 1 (full production).
2) The element should react almost linearly to small perturbations as in the Acclimate model.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 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 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 log " 0.2 (see the blue line in Fig. 8a).
In order 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Ñi of a link from sector s in country c to node i is given as With this, we can write the coupling term for the differential equation as 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.Fig. 9a shows the number of tipped nodes, and Fig. 9b the average node state xxy "

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 Hopfbifurcation 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.4).c) Two cusp differential equations with the parameters a1 " a2 " 4, b1 " b2 " 1, x0, 1 " x0, 2 " 0.5 and c1 " 0.2, c2 " 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 panel c), and one tipping element that possibly could undergo a Hopf-bifurcation (a " 1, µ " ´1).In the panels c) and d), there is no interaction between the tipping elements, i.e., Aij " 0 @i, j. e) and f ) Same as in panel c) and d), but with A21 " 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 two fold-bifurcations or, respectively, a tipping cascade of one fold-and one Hopf-bifurcation can be observed.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.

Fig. 1 .Fig. 2 .
Fig. 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.

Fig. 3 .Fig. 4 .
Fig. 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, x0 " 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 ci " 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.

Fig. 5 .Fig. 6 .
Fig. 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ˆ0.5 ˝grid, where each grid cell represents one rainforest patch.The colourbar 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 gipxq " 0 @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.

Fig. 9 .
Fig. 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 xxy 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.

Table 1 .
File structure of PyCascades.Only important files are described below.Separate plot-files and extended readme-files for the Amazon rainforest and the climate tipping cascades are also supplied in the corresponding directory to facilitate the usage of PyCascades.
The distributions for Gaussian, Lévy and Cauchy noise in PyCascades are taken from python's SciPy libraries