Cascading Transitions in Coupled Complex Ecosystems

The Earth system is arguably one of the most complex systems in the known universe. Over 4.5 billion years it has self-organised into a state which features complex differentiated life that covers its surface and penetrates its crust. This widespread biosphere requires stability in the sense of maintenance of temperature and pressures on the surface that allow liquid water. Within this range there is significant scope for change, some of it dramatic. Glaciation to inter-glacial cycles are classic examples of planetary-scale transitions. This paper examines the role of biological feedback on planetaryscale transitions. We present a conceptual model in which stability emerges as a consequence of interactions between environment and life. These mechanisms not only lead to stable ecosystem configurations, but can also produce critical transitions which we characterise as cascading transitions. Such failures can interact and produce system-wide transitions. Our results would be of interest to those studying realworld ecosystems, tipping elements in the Earth system as well as theoretical studies on complex artificial life systems.


Introduction
Tipping points, where critical levels of stress poise a system on the precipice of a transition between its current state and another, are a common feature of many complex systems with at least bi-stability (Lenton et al., 2008).Examples are numerous across all scales and areas of science, from the earliest studies of phase transitions in physics to multi-agent social, ecological and other life-environment systems (Stanley, 1987;Holme and Newman, 2006;Chavez et al., 2013;Williams and Lenton, 2007).This paper is centred on relating local Earth system transition on the ecosystem scale, to planetary scale transitions.
It has even been proposed that, on the planetary-scale, the state of planet Earth may approach a tipping point where a critical point in a component of the Earth-system corresponds to a global tipping point such as by a domino effect, or 'tipping cascade', as the result of forcing effects from human activity (Barnosky et al., 2012).In such a case, the global ecology undergoes a catastrophic shift to some new, unpredictable composition.Indeed, Earth's history is punctuated by such transitions; the Snowball Earth hypothesis posits a massive glaciation event 650 million years ago, where it is thought that critical levels of sea-ice overcame the previously stable global conditions, causing a transition to an almost completely frozen state (Hoffman and Schrag, 2002).The hypothesis is aptly named as a description both of the state of the planet, and the runaway positive feedback mechanism thought to be responsible for such a catastrophic transition.
The principle of tipping points is less contentious when applied to individual Earth-system components and local ecosystems although it remains an ever more pressing issue (Lenton et al., 2008); pressure on the Earth system from humans, such as anthropogenic emissions, climbs ever higher, and the response of the components of the Earth system appear unpredictable.To what extent do we expect complex systems to maintain robustness in the face of increasing perturbations, how can we characterise the behaviour of nearcritical systems, and how do Earth-system components reorganise when pushed beyond their regulatory limits?In part, the difficulty in answering these key questions stems from the absence of any broad real-world mechanisms which govern the stationary and critical behaviour of Earth-systems.
It is clear that there is some degree of coupling between the biotic and abiotic components of Earth.This is particularly evident in large systems, such as Earth's atmosphere.While roughly stable, they exist far from chemical equilibrium; the levels attained from purely abiotic processes (Kleidon, 2011).The living elements of these systems carry out a crucial role in establishing this stationary behaviour.However, non-equilibrium systems require balance to establish stability, and disruptions to this either through external perturbation, or large internal fluctuations, may result in catastrophes; periods of rapid change for the composition of the biota, and the state of its environment (Lenton et al., 2008).This can happen when environmental conditions change abruptly, such that the biota exerts an overall effect in its environment, and this has a positive influence on the activity of the agents of change, the previously stationary behaviour is lost, and the configuration both of the biota and its environment approach a new attractor.In the case of lake eutrophication, the system can be shifted into an Iain S. Weaver (2015)

Non-linear
Figure 1: The model biota consists of a large number, K, elements in the vector α.They are linearly coupled to their N environmental variables in vector E by the random, zero-mean coupling coefficients Ω.In turn, the environment determines their individual abundance by a simple niche idea; biotic elements are maximally abundant in the vicinity of their niche, µ.
anoxic state by increasing concentrations of otherwise limiting nutrients; while many plants and algae may flourish, the deoxygenation which results from their decomposition strongly favours the disruptive algae whose proliferation further alters the conditions to their relative favour.Lenton and Williams (2013) discuss conditions under which this mechanism may lead to to global transitions, dubbed 'planetary-scale tipping points'.These phenomena can occur when very large, and generally well-mixed systems undergo such a transition.Changes in the composition of the atmosphere, for example, is expected to have a rapid and global effect on the Earth system.Another mechanism by which the global state may abruptly shift is in response to some uniform, global pressure which causes a wide range of Earth-systems to transition simultaneously.Lenton and Williams (2013) point out the improbability of this explanation which would require a universal ecological threshold in spite of the inhomogeneity in the global ecology.The more likely alternative is that transitions in Earth systems on the local scale influence connected, codependent systems.By this process, the global system may self-organise such that local transitions can cause cascades, where failures of the regulatory mechanism in connected systems coincide not by chance but due to a causal or 'domino' effect.A sudden change in the environmental state of one system undergoing a transition will impact connected systems.If this is sufficient to destabilise the connected system, this can lead to an avalanche of transitioning systems, resulting in continentaland planetary-scale transitions.
Can transitions propagate through connected ecosystems in this way?The aim of this paper is to study a metasystem produced by the coupling of a large number of simple model ecosystems, whose coupling enables them to effect local transitions in one another.These models may be seen as a generalisation of Watson and Lovelock's (1983) Daisyworld; the two daisy types are replaced by a large and diverse biota while the simple, one-dimensional environment has been replaced by one of arbitrary dimensionality (Weaver and Dyke, 2013).We begin by describing the model in four components; i) abiotic model elements and their influence on the biota, ii) biotic elements and their effect on the environment, iii) external perturbing factors which affect the environment, iv) local coupling between adjacent systems defines the meta-system.
Following this, we examine the effect of local transitions on a network of complex ecosystem models with a range of connection strengths.The paper concludes with a discussion of the main results, along with a comparison to other complex systems which display a diversity of fixed points and internal transitioning behaviour.

Model formulation
The model can be viewed as a simplification of Watson and Lovelock's (1983) Daisyworld model which foregoes the need to precisely prescribe the behaviour of the small biota, instead favouring a large ensemble of randomly parametrised elements, similarly to Dyke (2010).The model components are summarised by Fig. 1, and explained in detail in this section.The biota is represented by a number K of independent and randomly parametrised biotic elements.We avoid making any assumption about the form the biota takes; be they individual organisms, species, ecotypes or populations.Furthermore, the biota is considered well-mixed and possesses no organisation or structure in and of itself; a cogent feature of some other artificial ecosystems (Chavez et al., 2013).Their individual abundance, metabolic activity or overall influence is denoted in the elements of the vector α.Importantly, the biota lacks any form of selfinteraction through processes such as competition or predation.These elements interact collectively with a simple environment, comprised of N variables denoted by the vector E. As before we do not prescribe a physical interpretation of the environmental dimensions beyond asserting they are variables which are both influenced by the biota, and influence the composition of the biota.The effect of individual biotic elements is unidirectional and monotonic; Figure 2: The activity of an element of the biota is a function of its environment E and niche position µ are shown in Fig. 2a for an N = 2-dimensional environment.The precise choice of function has been shown to be arbitrary with only the characteristic width, σ E being important (Weaver and Dyke, 2013).As the essential range is populated, the total activity shown in Fig. 2b approaches uniformity.
each species has only an increasing or decreasing effect on each aspect of the environment which is directly proportional to their biotic activity.These minimal ingredients, encapsulated in Fig. 1, have been shown to produce a range of interesting behaviours (Dyke and Weaver, 2013).Most significantly, homeostatic fixed points in the environmental variables emerge over a wide range of assumptions merely from the effects of a randomly assembled biota, foregoing the need for parametrisation or tuning.Feedback mechanisms exist within the environment in real systems; chemical weathering of silicate rocks is a negative environmental feedback where increases in global temperature, for example by increases in atmospheric CO 2 levels, lead to increased silicate rock weathering, a process which captures CO 2 .However, this model lacks such environment-environment or lifelife feedbacks; any coupling between biotic elements is mediated through the environment and vice-versa.Such interactions are excluded to avoid prescribing which limiting factors may be relevant.Competition is an example of feedback within the biota but is only one of many potential factors; on the biospheric level, the majority of biotic elements may not directly interact at all and so relaxing the assumption of strong interactions between biotic elements broadens our results and analysis.
i) Environment affects life The first component of Fig. 1 we discuss is the influence of the state of the environment E on the composition of the biota α.We implement a very simple niche idea whereby each element of the biota has an optimal environment, a point in the space of environmental variables where its activity is maximised.The time evolution of each element of the biota is, in the simplest case, a linear relaxation towards some steady-state activity α * with characteristic time scale τ α .
where α * is the fixed point in the time evolution of α, dictated by the state of the environment and a niche parameter unique to the each biotic element µ, which determines the most advantageous environmental conditions for each biotic element.In the vicinity of µ, α * is maximised and decays by some function as the environment departs the niche.Niches are randomly distributed uniformly in a finite volume of environmental variables called the essential range, which imposes limits on the ability of the biota to prosper in extreme environmental conditions.Again, in the simplest case this range is equal in all environmental variables and denoted σ µ .For the purpose of this paper, the precise choice of α * can be shown to be arbitrary (Weaver and Dyke, 2012;Dyke and Weaver, 2013) and we choose a N -dimensional Gaussian for visualisation purposes where appropriate although a parabola, step-function or some multi-modal function would equally suffice providing it has a well defined variance.
where σ E is the niche width, the characteristic width of the function, illustrated in Fig. 2a.As the essential range be- Figure 3: Two biotic elements with an increasing (ω 1 > 0) and decreasing (ω 2 < 0) effect on one environment variable are illustrated by Fig. 3a along with their niche positions.As the number of such elements, K, becomes large, the essential range is populated by a diverse biota and the net effect maintains zero mean, but significant variance, illustrated by Fig. 3b.
comes populated by a diverse biota, the mean activity inwith K, while the deviations increase with √ K such that in the limit of a very diverse biota of many unique elements, the total steady-state activity is approximately uniform over the essential range, shown by Fig. 2b.
ii) Life affects environment To be alive necessitates having some influence on the local environment through consumption or excretion of resources to maintain a metabolism although organisms may influence other environmental factors such as temperature, surface albedo and the abundance of other chemicals through their influence on other abiotic processes like weathering and erosion.Examples of such processes on Earth are abundant over nearly all spatial scales; from local effects such as changes in soil composition to modifications of atmospheric composition globally, by oxygenic photosynthesis (Wilkinson et al., 2009;Sanders and van Veen, 2011;Goldblatt et al., 2006).These impacts are encapsulated by concepts of niche construction and ecosystems engineering (Jones et al., 1994;Odling-Smee et al., 2003).
In the simplest case, the effects are linear and simply proportional to the activity of the biota, implemented by assigning each biotic element a unique influence on each aspect of its environment independently, stored in the matrix Ω such that the summed effect on each environmental variable, F may be found simply from the matrix product (3) The effects in Ω are chosen randomly with zero mean such that the model has no propensity for positive or negative feedback.The steady-state effects of two opposing biotic elements are shown by Fig. 3a, while the net effect of a large population is shown by Fig. 3b.Unlike the steady-state activity of a large population, the net effect of a population has zero mean and with appropriate scaling (discussed shortly) the variance in the surface is significant.The time evolution of the environmental variables is, in the simplest case, linearly driven by contributions from the biota F and external influence from adjacent systems, N , and global, abiotic influences P , whose form is discussed shortly.
The linear time evolution of α ensures this system has a 2Ndimensional phase space (Weaver and Dyke, 2013) though this can be reduced to N -dimensions if the time scales of processes associated with changes in the biota, α and the environment, E may be assumed separated, such that changes to the ecology occur on very much shorter time scales than those to the environment, or τ α τ E .In this instance, the activity of the biota may be replaced by its steady-state value α * .
It is useful to ensure that important model characteristics are invariant with the arbitrary model parameters, those which are neither random nor fundamental model components.This enables direct comparison of model behaviour over a range of parameters which might otherwise be obfuscated.We rescale the effect of the biota by introducing the normalisation constant where we have used the constant C to normalize the total biotic effect, F (E) such that its variance is independent of our choice of variance of the random coupling parameter, σ Ω , the fundamental niche width, σ E and the number of biotic components, K.
iii) Ecosystem connections Our model consists of a number of these systems connected into a network of interacting sub-systems, each with their own set of homeostatic states.
In the absence of any coupling, where N = 0, such a model behaves as a large number of independent systems whose large-scale statistics have already been established in depth (Dyke and Weaver, 2013).In this model, we couple adjacent systems by allowing the environmental variables to diffuse across edges with the following equation where D E parametrises the coupling strength and nn.denotes the set of n neighbouring vertices.This coupling can be seen as a site being perturbed by its neighbours in direct proportion to the difference between the state of its environment and that of the neighbouring site and the factor 1 n ensures this forcing is not more or less influential on vertices with higher or lower than average degree.This term is similar to the 'leakage' or heat conductance described by Harvey (2004), a modification to the original Daisyworld model where populations of biotic elements have distinct but coupled environments.Indeed, the imposition of a network structure coupling distinct artificial ecosystems is not new (Punithan and McKay, 2013) although a key difference in this work is our focus not on large-scale pattern formation, but the propagation of local transitions through a connected meta-system.
To establish invariance between the behaviour of our coupled ecosystem models and our parametrisation, it is useful to define the coupling strength D E in terms of the width of the effect function, σ E .We have previously ensured the variance of the effect function to be unity, and so we define D E in terms of the expected difference between homeostatic fixed points at a certain perturbation level By this definition, a system undergoing a transition with D E = 1 E increases the magnitude of its diffusion to neighbouring sites by approximately unity.
It is important to note that we have not explicitly included the diffusion of our biota.The reason for this is that if we choose such a diffusion process to be linear with the biotic activity then the inclusion of another diffusion process can be shown to do nothing beyond modify the existing diffusion process.To show this, we start from Eq. ( 1), taking the product with the coupling matrix, Ω as shown in Eq. ( 3), modifying the equation to include a diffusion term exactly as Eq. ( 7).In the steady state, we have where F * (E) is the steady state value of F (t) is the absence of diffusion, where D F = 0. Solving Eq. ( 9) for F (t) gives where n is the number of neighbours to the system, and we have eliminated the term including the summation which evaluates to zero at the steady state.Substituting this into the steady state condition for Eq. ( 4) gives which is identical to Eq. ( 4) with the substitution of F (E) for its steady state value and D E for the effective diffusion rate (1 + D F + nD F )D E .In summary, if a linear relaxation is chosen for the change in the composition of the biota towards its steady statue value, and the time scale of this process is much shorter than changes to the environmental state, any diffusion term in this process may be absorbed into the environmental diffusion process giving an effective diffusion rate, D E , of where D F and D E are diffusion rates of the biota and environment respectively.

Analysis
The life-environment systems detailed in the previous section have proven themselves to possess interesting behaviour in their own right; they generate a number of homeostatic fixed points in the environmental state where the configuration of the biota opposes changes in external perturbations to maintain the environmental state.This remains true even for a complex, high-dimensional environment (Dyke and Weaver, 2013).In this section, we thoroughly analyse the behaviour of coupled pairs of systems with respect to the value of the coupling parameter, D E , which correlates the state of their environments.By doing so, we detail metrics to shed light on the coincidence of transitions in connected systems.We perturb one system from its initial stable state to an adjacent one (whose basin of attraction borders that of its initial state, if any).The behaviour of the model is trivial for the limiting cases of D E = 0 and D E → ∞.In the  Figure 4: A system perturbed from the basin of attraction of a stable point undergoes a transition to the next stable fixed point.Fig. 4a shows how step sizes between fixed points are distributed in systems coupled to a stable environmental state.Fig. 4b shows, to first order, the feedback between connected systems; when a transition in one system does not coincide with a neighbour's transition, it produces only a small change in the environmental state inversely proportional to the gradient of F (E). Deviations from this covariate distribution are caused by cascading transitions.
first case, systems are uncoupled, while in the second case, the state of the environment is uniform across coupled systems, and the behaviour of the model is identical to a single system whose biotic effects consists of the sum of all systems' biota.In this case, correlations in changes between neighbouring systems is complete, as transitions cannot be confined to a single system.
To begin, we examine the probability distribution of transition sizes, |∆E| in a system coupled to another with an environmental state equal to its initial state.Previous work determined the expected number density of stable fixed points which can be used to find the mean transition size in an uncoupled system, given by Eq. ( 8).The simplest way to determine the form of this distribution however, is through Monte Carlo simulation whose results are given by Fig. 4a, which shows increasing coupling strength more strongly correlates environmental states in coupled systems, reducing the mean transition size slightly.
Next, we determine the response in the neighbouring system to the abrupt change in environmental diffusion which, to first-order (neglecting feedback), is inversely proportional to the gradient of its biotic effect function, F (E), in the direction of the change in environmental diffusion, ∆E.As the sum of the independent random contributions from the biota, the function is normally distributed with the variance fixed at unity by Eq. ( 6).The derivative ∇ ∆E F (E) (abbreviated to F (E)) is trivially shown to be independent of F (E) and similarly distributed with variance 1 2σ 2

E
. In order to produce a fixed point, the function must cross the plane F (E) = 0, and therefore the likelihood of the function pro-ducing a fixed point in a small interval increases linearly with the magnitude of the gradient, and the gradient distribution becomes To find the feedback from a transition in a neighbouring system, we require the reciprocal distribution which may be found from the moment generating function of the product

Results
A comparison between the statistics of our non-transitioning model with first-order feedback and a simulated coupled system is provided in Table 1.Our analysis in the previous section has shown that, to first-order, the coupling between a transitioning system and its neighbour is related to the distribution of the effect function gradient.This gradient is responsible for the homeostatic properties of individual systems as it dramatically reduces the magnitude of changes to its state in response to changes in forcing from neighbouring systems or external sources; significant increases in these perturbations cause only marginal changes to the state of the model, providing the homeostatic fixed point behaviour prevails.This can be seen by examining the correlation coefficient for the non-transitioning system, corr |∆E 1,2 | , which measures the linear dependance of changes to the state of the coupled systems.In contrast, our simulated system produces much larger values, corr sim |∆E 1,2 | ; the increase in linear dependence is therefore a result of coincidence of large transitions in the environmental state in the coupled system.This difference is mirrored by the integrated residuals, R 2 , which quantifies the difference between the non-transitioning and simulated distribution of environmental changes.Dyke and Weaver's (2013) simple ecosystem model exhibits a diversity of stable environmental configurations which emerge from a randomly parametrised biota.By connecting pairs of such systems, this paper has presented a minimal mechanism by which ecosystem-level transitions may cause cascades, which on a more extensive lattice of systems leads into large, potentially spanning transitions.Fig. 5 shows a lattice of one-thousand systems into a much larger metasystem.Changes in the coupling coefficient introduce longer range correlations in the environmental state and, when perturbed, lead to cascading transitions across parts the network.
The degree of connectivity, and the strength of connections between real Earth systems is unclear although these questions are ever more pressing; the environment and its biotic inhabitants are subjected to ever increasing anthropogenic pressures both locally, such as land use change, and globally, such as emissions of greenhouse gasses (Steffen et al., 2015).Investigating the robustness and limitations of regulatory mechanisms on Earth is as important as understanding the influence of a transition both in terms of alternative stable states, and the influence of such a dramatic change on coupled Earth systems (Williams and Lenton, 2010).This work is an early step towards understanding the collective behaviour of coupled multi-stable complex systems and lays the foundation for a more thorough and general analysis of correlations and transitions within the network.

Increasing
Iain S.Weaver (2015)  Cascading Transitions in Coupled Complex Ecosystems.Proceedings of the European Conference on Artificial Life 2015, pp.612-619 (a) The activity of a biotic element with a niche in the center of the essestial range.(b)The total activity of many biotic elements witha diversity of niches.
Iain S.Weaver (2015)  Cascading Transitions in Coupled Complex Ecosystems.Proceedings of the European Conference on Artificial Life 2015, pp.612-619(a) The effect function generated by two biotic elements with opposing coupling values.(b) The total effect function of many randomly parametrised biotic elements.
(a) Distance between adjacent fixed points with varyied coupling strength.(b) Codistribution of changes in environmental state without cascades.

Figure 5 :
Figure5: A lattice of connected systems mirrors the studied behaviour of connected pairs of systems; increasing coupling strength increases the correlation between the environmental state of groups of systems.Along with correlations comes the likelihood of local transitions causing cascades, dramatically changing the environmental state of large sections of the network.
Cascading Transitions in Coupled Complex Ecosystems.Proceedings of the European Conference on Artificial Life 2015, pp.612-619

Table 1 :
and ρ sim are the probability distributions computed in Fig.4band generated by simulated systems respectively.Large values of R 2 indicate greater departure from Iain S.Weaver (2015)Cascading Transitions in Coupled Complex Ecosystems.Proceedings of the European Conference on Artificial Life 2015, pp.612-619We measure changes to the environmental state of two coupled systems when one undergoes a transition into a new stable state.We quantify the discrepancy between the codistribution where there are no coincident transitions, where only first order coupling is considered, and a simulation which allows such events to occur.thenon-transitioningdistribution and can be attributed to small higher-order feedback, and an increasing coincidence of transitions in connected systems.In addition to this, we can compare the correlation of |∆E 1 | and |∆E 2 | from the simulated systems to that predicted by Fig.4b.A higher correlation indicates the state changes in connected systems have a greater degree of linear dependence, corresponding to increases in the coincidence of transitions.