Earlier collapse of Anthropocene ecosystems driven by multiple faster and noisier drivers

A major concern for the world’s ecosystems is the possibility of collapse, where landscapes and the societies they support change abruptly. Accelerating stress levels, increasing frequencies of extreme events and strengthening intersystem connections suggest that conventional modelling approaches based on incremental changes in a single stress may provide poor estimates of the impact of climate and human activities on ecosystems. We conduct experiments on four models that simulate abrupt changes in the Chilika lagoon fishery, the Easter Island community, forest dieback and lake water quality—representing ecosystems with a range of anthropogenic interactions. Collapses occur sooner under increasing levels of primary stress but additional stresses and/or the inclusion of noise in all four models bring the collapses substantially closer to today by ~38–81%. We discuss the implications for further research and the need for humanity to be vigilant for signs that ecosystems are degrading even more rapidly than previously thought. Current models, based on incremental changes in a single stress, have limited ability to anticipate abrupt ecosystem changes due to climate and human activities. Experiments on four models simulating ecosystems with a range of anthropogenic interactions show how much earlier abrupt change can happen.

For many observers, UK Chief Scientist John Beddington's argument that the world faced a 'perfect storm' of global events by 2030 1 has now become a prescient warning.Recent mention of 'ghastly futures' 2 , 'widespread ecosystem collapse' 3 and 'domino effects on sustainability goals' 4 tap into a growing consensus within some scientific communities that the Earth is rapidly destabilizing through 'cascades of collapse' 5 .Some 6 even speculate on 'end-of-world' scenarios involving transgressing planetary boundaries (climate, freshwater and ocean acidification), accelerating reinforcing (positive) feedback mechanisms and multiplicative stresses.Prudent risk management clearly requires consideration of the factors that may lead to these bad-to-worst-case scenarios 7 .Put simply, the choices we make about ecosystems and landscape management can accelerate change unexpectedly.
The potential for rapid destabilization of Earth's ecosystems is, in part, supported by observational evidence for increasing rates of change in key drivers and interactions between systems at the global scale (Supplementary Introduction).For example, despite decreases in global birth rates and increases in renewable energy generation, the general trends of population, greenhouse gas concentrations and economic drivers (such as gross domestic product) are upwards 8,9 often with acceleration through the twentieth and twenty-first centuries.Similar non-stationary trends for ecosystem degradation 10 imply that unstable subsystems are common.Furthermore, there is strong evidence globally for the increased frequency and magnitude of erratic events, such as heatwaves and precipitation extremes 11 .Examples include the sequence of European summer droughts since 2015 12 , fire-promoting phases of the tropical Pacific and Indian ocean variability 13 and regional flooding 11 , already implicated in reduced crop yields 14 and increased fatalities and normalized financial costs 9 .
The increased frequency and magnitude of erratic events is expected to continue throughout the twenty-first century.The Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report concludes that 'multiple climate hazards will occur simultaneously, and multiple climatic and non-climatic risks will interact, resulting in compounding overall risk and risks cascading across sectors and regions' 11 .Overall, global warming will increase the frequency of unprecedented extreme events 11 , raise the probability of compound events 15 and ultimately could combine to make multiple system failures Article https://doi.org/10.1038/s41893-023-01157-xas reinforcing feedbacks accelerate connections or human activities increase stress levels.However, extreme events could also counteract each other (for example, extreme droughts and extreme rainfall events) and interconnections could also have weakening effects (for example, where increased plant growth driven by increased CO 2 is counterbalanced by increased temperatures and droughts.To date, there is limited observational evidence showing that ecosystems have a record of tipping between alternate stable states 21 . Others 19 offer a mathematical tripartite classification of critical transitions that includes slow driver bifurcations, rate-induced (fast/ cumulative driver) and noise-induced (extreme event) tipping points.However, previous studies tend to focus on each of these categories individually.For example, there is a well-established body of physics and mathematical theory on 'mean exit times' 22 , with studies investigating the timing of tipping points in rate-induced [18][19][20] or noisy 19,23,24 systems.However, despite calls for more experimental evidence of the impacts of climate variability and extremes on ecosystems 25,26 , the relative importance or combined effect of fast drivers, multiple drivers and noisy system drivers on the collapse of real-world ecosystems is not known.Critical transitions driven by current pollution forcings such as greenhouse gas emissions 27 and nutrient loadings 28 are likely to be new, well beyond the envelope of natural variability.Hence, we avoid the use of the terms critical transition and tipping points, used formally in dynamical systems theory to represent shifts to alternative attractors and focus on abrupt threshold-dependent changes (ATDCs) that would be perceived by society as the quantitative (for example, fish and stock integrity) and/or qualitative (for example, ecosystem functions) collapse of a desirable system state 29,30 .more likely 16 .For example, there is a risk that many tipping points can be triggered within the Paris Agreement range of 1.5 to 2 °C warming, including collapse of the Greenland and West Antarctic ice sheets, die-off of low-latitude coral reefs and widespread abrupt permafrost thaw 17 .These tipping points are contentious and with low likelihood in absolute terms but with potentially large impacts should they occur.In evaluating models of real-world systems, we therefore need to be careful that we capture complex feedback networks and the effects of multiple drivers of change that may act either antagonistically or synergistically [18][19][20] .Prompted by these ideas and findings, we use computer simulation models based on four real-world ecosystems to explore how the impacts of multiple growing stresses from human activities, global warming and more interactions between systems could shorten the time left before some of the world's ecosystems may collapse.
Intuitively, stronger interactions between systems may be expected to increase the numbers of drivers of any one system, change driver behaviour and generate more system noise.As a result, we would anticipate that higher levels of stress, more drivers and noise may bring forward threshold-dependent changes more quickly.For any particular system (for example, the Amazon forest) it is possible to envisage a time sequence that starts with one main driver (for example, deforestation), then multiple drivers (for example, deforestation plus global warming), more noise through extreme events (for example, more droughts and wildfires), with additional feedback mechanisms that enhance the drivers (for example, diminished internal water cycle and more severe droughts).A vortex could therefore emerge, with drivers generating noisier systems as climate variability and the incidence of extreme events increases.Under worst-case scenarios, the circle becomes faster a, The four systems models simulated in this study (see section on Overview of systems models).b, Schematic representation of a system dynamics model (Lake Phosphorus model) with its external slow (blue and green) and noisy (red/orange) drivers depicted in colour (see Generation of future scenarios).c, Depiction of the four experiment types (section on Generation of future scenarios), ranging from changes in the primary baseline driver only (experiment 1), changes in all slow drivers and noise inputs simultaneously (experiment 4, where 'a' and 'b' represent noise profiles that are uncoupled or coupled to the primary driver trajectory, respectively): darker colours schematically represent steeper trajectories and/or higher noise levels.d, The two linear techniques used to check whether outcomes shift into a functionally different state (section on Time-series breakpoint detection)-the top panel is applied to Lake Chilika, Easter Island and TRIFFID, where the systems collapse from high quantitative outcome states to low quantitative outcome states and the bottom panel is applied to Lake Phosphorus (where lake phosphorus concentrations shift from low to high).e, Depiction of the time-series breakpoint date recognition (section on Time-series breakpoint detection).The Easter Island icon in a is made by Roundicons and the remaining three icons are made by Freekpik, as sourced from www.flaticon.com.

Article
https://doi.org/10.1038/s41893-023-01157-x We have selected a range of system dynamic models that have been previously used to demonstrate generalizable findings (for example, with regard to safely overshooting ATDCs 27 ) and can be externally manipulated to simulate internal emergent ATDCs at local and regional scales-as if they were impacted through stronger connections to other systems.Reflecting modern ecosystems, these models show varied anthropogenic interactions, ranging from social-ecological systems with strongly coupled human-nature feedbacks to ecological systems with predominantly one-way interactions where ecosystems are influenced by the external impacts of people.The ability of these models to capture feedback loops, delays and interactions between components is well established 31,32 and has motivated their use in various recent studies of sustainability and resilience 21,[33][34][35] .Therefore, guided by the ref.19 typology of tipping points, we aim to generalize the dynamics of increasing the numbers of drivers, their rates and variability (as proxies for stronger interactions between systems and noise) on the speed at which ATDCs are reached in four ecosystem dynamics models (Fig. 1): Lake Chilika lagoon fishery 21,33 , Easter Island 36 , Lake Phosphorus 28,37 and a modified version of The Hadley Centre Dynamic Global Vegetation Model (TRIFFID) of forest dieback 27,38 .

Results
As described in the Methods, the four models each have a primary (baseline) slow driver (Fig. 2, grey boxplots), where linear changes in their trajectories over time can initiate ATDCs in their respective outcome variable (Lake Chilika, fish population; Easter Island, human population; TRIFFID, tree coverage; Lake Phosphorus, lake phosphorus concentration).When the strength of the primary slow driver in each model is increased, the modelled systems collapse sooner-as defined by a statistical breakpoint in their temporal trends (section on Time-series breakpoint detection).Increasing the strength of multiple drivers with additional secondary and tertiary drivers further reduces the breakpoint date (Fig. 2), with variation around these median responses determined by the relative strength of the additional drivers-with addition of a weak secondary driver bringing forward the start of system collapse substantially less than the addition of a strong secondary driver (Supplementary Fig. 2-1).
In addition to earlier breakpoint dates, extra drivers can also cause ATDCs at levels where it would be resilient to the primary slow driver in isolation (Supplementary Section 2).For example, across the 1,000 timesteps of the Lake Phosphorus model, the system is stable at normalized baseline driver rates up to 0.348 (that is, Lake Phosphorus

Article
https://doi.org/10.1038/s41893-023-01157-xconcentration does not go through a breakpoint; Supplementary Fig. 2-4d).However, the addition of a single secondary driver of normalized strength 0.3 can lead to breakpoints occurring at normalized primary driver strengths 0.312 (reduction from baseline: 0.036 (10.3%);Supplementary Fig. 2-4d) and the addition of an extra tertiary driver with normalized strength 0.3 can lead to breakpoints at normalized primary strengths 0.270 (reduction from baseline: 0.078 (22.4%);Supplementary Fig. 2-4d).With all additional drivers, 12.3% of breakpoints observed in the Lake Phosphorus model occurred at primary driver strengths below the minimum threshold required to result in a breakpoint when the primary driver is acting in isolation (Lake Chilika, 1.2%; Easter Island, 14.8%; TRIFFID, 7.7%; Supplementary Table 2-1).
Next, for each of the four models, the trajectories of the primary slow drivers were randomly perturbed by the addition of noise (section on Generation of future scenarios).Noise was generated within the system dynamics software used to run the models (STELLA 39 ) by randomly sampling per timestep from a normal distribution with a mean value of 0 and standard deviation of σ, meaning that random perturbations on the system could work in both positive (σ > 0) and negative directions (σ < 0).The value of σ was randomly sampled once per simulation to explore the effects of different noise scales on the time to reach the breakpoint date (section on Generation of future scenarios).The addition of high noise (normalized σ > 0.666) shows that increasing the variability of the primary slow driver (in isolation) across all four models can bring forward the date of system collapse (Fig. 3).
The effects outlined above are synergistic-combining multiple drivers with noise further reduces the breakpoint date beyond the effects of either multiple drivers or noise acting alone (Fig. 4).For example, at a normalized slow baseline driver strength of 0.3 in the Easter Island model (Fig. 4b), the addition of low uncoupled noise (normalized σ ≤ 0.333) with all possible additional drivers switched on with normalized strengths of over 0.666 ('high' secondary and tertiary trajectories) brings the median breakpoint forward from timestep 1,179 to timestep 426 (63.8% reduction), whereas high noise levels (defined as normalized σ > 0.666) brings the breakpoint forward from timestep 1,179 to timestep 225 (80.9% reduction).The finding that the breakpoint date is most advanced by the combination of high noise and high secondary trajectories is consistent across the other three models, with the median breakpoint date at a normalized slow baseline driver strength of 0.3 changing from year 2047 to year 2035 (37.5% reduction) for Lake Chilika, timestep 238 to timestep 92 (61.3% reduction) for TRIFFID and timestep 848 to timestep 388 (54.2% reduction) for Lake Phosphorus.Across all combinations of noise and multiple drivers, 1.7%, 7.5%, 6.6% and 8.9% of modelled breakpoints occurred at primary driver strengths below the minimum threshold required to result in a breakpoint when Model timestep units and boxplot dimensions are the same as in Fig. 2; see Supplementary Table 3-1 for the number of model simulations underpinning each boxplot.

Article
https://doi.org/10.1038/s41893-023-01157-xacting in isolation for Lake Chilika, Easter Island, TRIFFID and Lake Phosphorus, respectively (Supplementary Table 2-4).All results presented above are robust to different modelling and monitoring decisions.For example, these results are consistent regardless of whether the noise is coupled to (allowed to grow with) the magnitude of the primary slow driver or uncoupled and sampled from a constant distribution (Supplementary Figs.2-2 and 2-3 and Supplementary Tables 2-3 to 2-5) and irrespective of whether linear, nonlinear or threshold-type boundaries 40 are used to define the breakpoints (Supplementary Section 4 and Supplementary Figs.4-1 to 4-6).

Discussion
Previous findings have supported the idea that Earth's subsystems may interact to the extent that an abrupt shift in one raises the probability that a shift may occur in another [41][42][43] .In this paper, we have explored, through four ecosystem models, how these interactions may alter the timing of ATDCs through the effects of strengthened drivers, multiple drivers and higher internal variability or noise.The potential effects are substantial with combinations of a strengthened main driver, an additional driver and noise giving at least 38-81% reductions in the future date of a predicted ATDC compared to estimates for a non-interacting system with a constant single driver and no noise.Importantly, the effect per unit time on bringing forward an ATDC is greatest at low driver trajectories, which further strengthens the suggestion that abrupt Earth system changes may occur sooner than we think (Supplementary Introduction).Our findings also show that 1.2-14.8% of ATDCs can be triggered by additional drivers and/or noise below the threshold of driver strengths required to collapse the system if only a single driver were in effect.
Overall, we find that, as the strength of a main driver increases, the systems collapse sooner.Adding multiple drivers brings collapses further forward, as does adding noise, and the two effects can be synergistic.However, the relative importance of these changes varies across systems.For the Chilika fishery, the most influential driver is captured as the primary driver and so additional drivers have limited influence, with the addition of noise in the primary driver bringing the breakpoint date much closer to the present.For Easter Island, TRIFFID and Lake Phosphorus, the opposite is true-the addition of high levels Note, the Lake Phosphorus model (d) did not produce any outcomes between the 0.65-0.75primary driver range within the 'high trajectory, high noise' scenario; however, the median breakpoint date of the adjacent range (0.55-0.65) is 346.Model timestep units and boxplot dimensions are the same as in Fig. 2; see Supplementary Table 3-1 for the number of model simulations underpinning each boxplot.

Article
https://doi.org/10.1038/s41893-023-01157-x of noise in the primary driver advances the date of system collapse far less than additional drivers.Thus, while the earliest collapses in all the systems are found when both additional drivers and noise are applied, an important implication for real-world governance is that the precise importance of individual driver trajectories and noise is system-dependent.

Earlier occurrence of abrupt threshold-dependent changes
Our results show that systems do not collapse at a constant level of cumulative stress (that is, total stress built up over time) irrespective of the rate of stress change (Supplementary Section 5) but rather underline the importance of rate over accumulated stress [18][19][20] .Simulations where the primary, secondary or tertiary drivers change more rapidly tend to shift earlier and are less able to absorb cumulative stress (Supplementary Fig. 5-1).Thus, the same ecosystem can collapse as a result of sustained/cumulative pressure of a slower driver but will probably collapse faster if the rate of change is increased [18][19][20] .Increasingly fast driver rates will eventually overwhelm the ability of balancing feedback loops to compensate for increased stress on the system; thus, signifying a loss of resilience.In the absence of strong balancing loops, a fast driver allows reinforcing feedback loops to grow (Supplementary Section 6).
The driver may also re-energize dormant reinforcing feedback loops or allow new coupled, reinforcing feedback mechanisms to emerge (compare ref. 44).For the Easter Island, TRIFFID and Lake Phosphorus models, as the balance of feedback loops shifts towards reinforcing loops, the probability that the system will be driven out of its attractor into an ATDC increases (Supplementary Section 6).Additional drivers limit further the balancing ability of balancing feedback loops and increase the probability of collapse.For Lake Chilika, the pre-ATDC phase is dominated by reinforcing feedback loops driving fisher population growth towards dangerous levels, with collapse coinciding with the growth of balancing feedbacks in the form of reduced fish populations.These rebalance the system by limiting the effectiveness of the fisher population's fishing efforts (Supplementary Fig. 6-1).
In our analysis, the rise in driver stress is continuous over time.Where the stress is applied in discrete events (for example, wildfire events), the same response can be expected where elapsed time between events is insufficient for balancing feedback loops to rebalance the system or where large stress events motivate additional amplifying loops.This is similar to the impact of extreme events (that is, noise; Figs. 3 and 4), which has the ability to push a system out of its attractor temporarily or permanently; an effect that strengthens as the system becomes increasingly sensitive to perturbations close to a potential ATDCs 19,23 .However, sequences of extreme events from multiple drivers, such as extreme drought followed by extreme rainfall, may only act antagonistically where sufficient time allows for the system to repair the extreme impacts.Our study only looks at driver noise; there could coincidentally or equally be natural 'state' change/ noise (vertical axis on phase-plot figures)-for example, natural tree mortality, natural lake infilling, fluctuating populations in ecosystems, or ageing population, behavioural/psychological changes in the social domain-all of which could alter the probability of ATDCs even in the absence of, or changes in, the external drivers 19,23 .

Moving forward
These results have research implications for further developing and applying models of ecosystems to study ATDCs.Whilst our findings derive from models based on real-world systems, the greater complexity of reality may limit the transferability of our results.The Lake Chilika model is the most complex of the four models, with upwards of 100 model variables capturing hydroclimatic, ecohydrological, fishery and socio-economic dynamics interacting to create four balancing loops and seven reinforcing loops-and is validated against historical data 33 .Of all the models, it shows the least dramatic reductions in the date of any ATDC (Supplementary Introduction).Therefore, it is plausible that more complex systems will have stronger regulating mechanisms that stabilize the system through sets of balancing feedback loops 44 , constraining the more extreme of our findings.
Mechanistically, in simpler models, such as the Lake Phosphorus model, regime shifts may be triggered by a single feedback loop.In more complex models (and probably ecosystems), our analysis of feedbacks strengths shows evidence for an instability cascade through the system via multiple feedback loops.For example, the collapse in the Easter Island human population reflects the cumulative effects of several feedback loops triggered by overharvesting the tree population.Growing instability weakens the balancing feedbacks for the tree population, rat population and agricultural carrying capacity (Supplementary Fig. 6-2), allowing the reinforcing loop for the decline in human population to strengthen.In general, increasing driver strengths can trigger these mechanisms earlier, whereas additional drivers have the ability to shift the nature of the cascade (for example, including/excluding different feedbacks; Supplementary Figs.6-5 to 6-8).However, in spatial terms, multiple interacting feedback mechanisms may lead to spatial re-organization which slows the rate of collapse 45,46 , with stochasticity promoting temporal stability -particularly in local regions with small populations 24 .There is the possibility, too, that interconnections could have weakening effects and, where the impacts are slower than the system response, extreme events could counteract each other.Thus, our quantitative findings could be viewed as representing worst-case scenarios for the different ecosystems 7 .
Nevertheless, the finding that additional stress produces qualitatively similar emergent phenomena in a range of simulation models should not be dismissed lightly 47,48 .The consistency across models representing varying processes, interactions and contexts may indicate that equifinality makes the accurate representation of internal system dynamics less important than the external drivers/stresses in simulating complex realities 49 .Clearly, model development is required to better capture the diversity of system elements, interactions and feedbacks for more complex systems and, in particular, more realistic coupling of human decision-making and ecological/environmental dynamics.With the exception of Lake Chilika 33 , each model in this study was originally created to study the impact of a primary driver influenced by predominantly external anthropogenic processes, presumably the driver perceived as the most impactful.Our results show that this assumption may not be the case (for example, Easter Island) and models should include a range of plausible drivers and scenario combinations if they are to avoid underestimating the risk of ATDCs.Moreover, new ecosystem models should allow for the growth of feedback loops and long-term simulations to observe the mechanisms that underpin ATDCs 48,50 .For example, more realistic social-ecological coupling may lead to shifts in the human decisions capable of either shifting an ATDC much closer to the present or avoiding it completely.Monitoring of real-world systems should therefore capture multiple plausible drivers, their variability and their feedbacks to social systems.More ATDCs will occur unexpectedly if the focus on perceived main drivers ignores other drivers that increase cumulative stress and gradually reduce the resilience of systems, as exemplified in the lake water regime shift at Erhai, western China 28 .There, abrupt lake eutrophication was initially perceived to have been driven by transgression of a threshold in nutrient enrichment driven by agricultural runoff but historical analysis has shown that the shift was also affected by lake water-level management, seasonal climate and fish farming 44 .
Substantial research has focused on identifying early warning metrics linked to critical slowing down theory which applies primarily to 'equilibrium' system states with single, slow drivers 51 .If, as we indicate, real-world tipping elements are more likely to be driven by multiple, fast drivers and extreme events, it is less likely that early warning signals in the frequency domain will be observed 20,51 for noise-induced thresholds.Certainly, excluding noise from model systems, whilst a potentially useful simplification for theoretical understanding, risks Article https://doi.org/10.1038/s41893-023-01157-xcreating a false sense of security by overestimating the distance remaining before critical thresholds are breached in the real world where multiple drivers and noise are abundant 27,52 .Therefore, alternative approaches to identifying resilience loss in real systems before ATDCs through structural metrics [53][54][55] and early warning signals generated by agent-based models 50 should be considered more widely.
Previous studies of interactions between tipping elements have focused on large-scale systems and suggest substantial social and economic costs from the second half of the twenty-first century onwards 42,56 .Our findings suggest the potential for these costs to occur sooner.For example, it is not clear whether the IPCC estimate for a tipping point in the Amazon forest before 2100 (ref.11) includes the possibility for interacting drivers and/or noise; if not, our findings suggest that a breakdown may occur several decades earlier (Supplementary Introduction).This would occur where local-scale failures in elements (such as species populations, fish stocks, crop yields and water resources) combine with more extreme events (such as wildfires and droughts) to precondition the large-scale system, already vulnerable to the influence of other large-scale tipping elements, to collapse earlier-a meeting of top-down and bottom-up forces (Supplementary Introduction).This vertical integration of forces is reinforced by the rising trend in global warming that already represents a spatial integrator which may be expected to strengthen before it subsides.Clearly, climate economics need to incorporate these synergistic and cumulative effects that are occurring at local and regional scales into larger scale models where they are currently lacking 57,58 .The dominance of accelerating trends in global time series of economic consumption (for example, refs.9,59) makes our finding that ramping up the main driver is the easiest way to bring forward an ATDC particularly worrying.Similarly, the implication for regions experiencing more extreme events is that an ATDC may occur even before the main driver has ramped up.
The commonality of findings across four well-studied ecosystems has potentially profound implications for our perception of future risks associated with the climate and ecological crises.While it is not currently possible to predict how climate-induced ATDCs and the effects of local human actions on ecosystems connect across temporal and spatial scales, our findings show the potential for each to reinforce the other.The ability of present policy and practice to prevent an ever-deepening vortex of degradation in local and regional ecosystems requires urgent investigation 7 .

Overview of systems models
Here, we briefly describe the four previously published models used to investigate the effects of multiple drivers and noise upon the timing of ATDCs.Each model was replicated and simulated within the system dynamics software STELLA Architect v.1.6.1 (ref.39), with outputs exported into CSV files as time series and analysed in the statistical software R v.4.1.0(ref.60).The models, example data and code used in the analyses are available via https://doi.org/10.5281/zenodo.7946433.
The Lake Chilika fishery model 21,33 is a social-ecological model designed to simulate the future fish population and catch trajectories of the Chilika lagoon, Odisha, India.The model is able to explore the impacts of many slower drivers (fisher population growth and increased rainfall and temperatures under climate change) and many faster drivers (abrupt changes in fish prices and fishing gear) on the sustainability and resilience of the fish population until 2100.As described in detail in ref. 33, the model includes coupling between many social and ecological components of the system.First, the efficiency of fish catch efforts is proportional to the fish population density within the lagoon (as fish density declines, catch per unit effort also decreases).Second, as a form of environmental carrying capacity, the fisher population growth is proportional to the total number of livelihoods supportable by the overall fishery value, which is derived from the total fish catch in any given month.Third, fishers may invest their fishing revenues into more intensive fishing gear (motorboats), which also has implications for fish catch and fish stock health over time.The model is also able to simulate many natural resource governance approaches (for example, fishing quotas and alternative livelihoods), although the model runs conducted here are all under the baseline governance scenario 33 (the tidal outlet between the lagoon and the Bay of Bengal is re-opened every 10 years to rejuvenate fish migration and lagoon salinity).The model has been previously validated against empirical data through standard behaviour-matching techniques and Monte Carlo sensitivity analysis 33 .The Lake Chilika model is run for a total of 1,536 timesteps (months), with each time series aggregated to the annual scale (about 1973-2100).Future trajectories, detailed in the section below on (Generation of future scenarios), activate from timestep 504 ( January 2015) after the completion of the historical parameterization and validation periods 33 .
The Easter Island model aims to explore alternative hypotheses behind the collapse of the Easter Island civilization 36 .The initial parameterization of the model here is the same as the 'ecocide' configuration detailed in ref. 36.The main internal social-ecological feedback driving model dynamics is the balancing feedback between human population growth, tree coverage and land clearance, whereby the overharvesting of the primary resource (palm forest) can lead to overshoot dynamics and the eventual demise of the human population (ecocide).As noted in ref. 36: 'While it is obvious that the islanders were not directly living from palm trees, the forest provided several valuable and difficult to substitute ecological services, including food from fruits and palm nuts, timber to construct houses and sea-going canoes for fishing'.In addition to this main internal social-ecological feedback, many external variables can be modified to change the speed of human population growth, including the tree clearance rate per capita, the maximum carrying capacity of the agricultural system (to help support human population growth) and the mortality rate of trees (representative of potential disease outbreaks).The model is run for 1,500 timesteps (years), with future scenarios active from the first timestep (Generation of future scenarios).
The TRIFFID model is a modified version of The Hadley Centre Dynamic Global Vegetation Model, originally developed by ref. 38 to explore the effects of atmospheric CO 2 concentrations on the rate of Amazon dieback.Here, we simulate the modified version developed by ref. 27, which is based around a central Lotka-Volterra equation describing the change in vegetation coverage as the primary external driver (local atmospheric temperatures) increases.On any given timestep, the change in vegetation coverage (dv/dt) is driven by a temperature-dependent growth term and an externally set disturbance rate: Where v is the vegetation coverage (ranging from 0-100 %), T f is the temperature forcing parameter (see below), g is the vegetation growth rate (the increase in vegetation coverage [v]), g 0 is the maximum growth rate (set at 2 % per year), y is the disturbance rate (see below), T l is the local temperature, T opt is the optimal temperature (28 °C), β is the half-width of the growth versus temperature curve (10 °C) and α is the difference in temperature between surface bare soil and forest (5 °C).Therefore, the growth term is assumed to be parabolic with the local temperature (equation (1b)), meaning that once the local temperature increases beyond the optimal temperature, negative tree growth ensues (that is, additional tree mortality 27 ), which in turn leads to an Article https://doi.org/10.1038/s41893-023-01157-xincrease in temperature (equation (1c)), which may eventually produce the runaway loss in tree coverage.Although the meaning of the disturbance rate is not specified by ref. 27, it may proxy human-induced ecosystem stresses such as deforestation for agricultural land and disease-driven forest dieback.The model is run for 500 timesteps, with future trajectories active from the first timestep (Generation of future scenarios).The Lake Phosphorus model is a simplified version of the original 'lake response to P input and recycling' model 37 , as modified by ref. 28.The model is designed as a simple ecosystem model, with lake water phosphorus concentration driven by a generic external phosphorus input (which may, for example, proxy external inputs from agricultural runoff, sewage and industrial discharges from factories, construction sites and urban areas) 61 .In turn, lake water phosphorus is recycled back into the system as an ecological reinforcing feedback loop, proportional to the lake phosphorus concentration on any given timestep.Phosphorus is also removed from lake waters via sedimentation, where the volume removed in sediment is proportional to the phosphorus concentration of the lake.Therefore, on any given timestep, the change in lake phosphorus concentration (dP/dt) equals: Where P is phosphorus concentration, α is phosphorus input rate, r is the maximum recycling rate, s is the phosphorus loss rate, n is the strength of the recycling response to phosphorus concentrations (n = 8) and t is time (see below).The model is run for 1,000 timesteps (unitless), with future scenarios active from the first timestep (see below).Given the simplicity of this model, an area for future research could be expanding the original model to explore how adaptive management mechanisms may help to avoid ecosystem thresholds, for example, by linking government fertilizer incentives to lake phosphorus levels as the ecosystem approaches a threshold.

Generation of future scenarios
Using the above models, we performed four in silico experiments (presented visually in Fig. 1): • Experiment 1: only the primary slow driver in each model changes over time and all other drivers remain constant (Fig. 2 baseline).
• Experiment 2: multiple slow rates, with up to two additional (secondary and tertiary) slow trajectories on top of the primary driver changing over time (Fig. 2 multiple drivers).• Experiment 3: the addition of noise to the primary trajectory (Fig. 3), with all other drivers held constant.The magnitude of noise may be either coupled or uncoupled from the trajectory of the primary driver.• Experiment 4: the addition of noise to the primary driver, with up to two additional slow drivers (Fig. 4).The magnitude of noise may be either coupled or uncoupled from the trajectory of the primary driver.
To survey a wide range of future trajectories and generate a sufficient number of simulations that collapsed (Time-series breakpoint detection), each of the models were run for the following number of iterations (including both coupled and uncoupled settings): • Lake Chilika fishery: 70,000 • Easter Island: 70,000 • TRIFFID: 70,000 • Lake Phosphorus: 120,000 In turn, to maximize computational efficiency both in STELLA and in R, the following logic was applied to the inbuilt Monte Carlo function in STELLA to automatically generate the four different experiment types described above (the baseline primary driver always remains 'on/active'): • If µ 1 > 0.4 then secondary driver active else secondary driver remains at default value.• If µ 2 > 0.4 then tertiary driver active else tertiary driver remains at default value.• If µ 3 > 0.4 then noise active else noise level remains at zero.
Where µ 1 , µ 2 and µ 3 represent 'on switches', with values randomly sampled from uniform distributions between 0 and 1 per simulation.The number of simulations per model experiment which showed ATDCs are detailed in Supplementary Table 3-1.Whilst some insights could be obtained deterministically 62 , this is not possible for all models (for example, Lake Chilika) nor for all experiments (those involving additional noise).Thus, undertaking these model runs and analyses of the outputs (below) is the most consistent, feasible approach suitable across all models and experiments, allowing for comparisons across experiments, as well as investigation of synergistic impacts-fulfilling our primary aim of investigating the impact of the interaction of fast drivers, multiple drivers and system noise on the timing of tipping points in ecosystems.
To investigate experiment 1, each of the four models has one primary baseline driver which changes from its default value in every simulation: • Lake Chilika fishery: fisher population growth rate (net difference between the birth rate per 1,000 population and the death rate per 1,000 population) • Easter Island: tree clearance rate (trees per person per year) • TRIFFID: local temperature (°C) • Lake Phosphorus: phosphorus input rate (unitless) Baseline outputs were generated with the primary driver active and the secondary and tertiary drivers remaining at its default value and the noise level remaining at zero (Supplementary Table 3-2).In turn, the Monte Carlo sensitivity analysis function in STELLA randomly samples a future change trajectory for the primary slow driver per simulation (as plotted on the horizontal axes of Figs.2-4).The primary trajectory is sampled between the lower and upper limits of uniform distribution bounds, meaning that there is a uniform likelihood of selecting any given trajectory between the bounds (Supplementary Table 3-2).A future change trajectory of '0' would cause no change from the default value; the maximum trajectory change limits for each of the models can be seen in Supplementary Table 3 The built-in STELLA 'TIME' function generates future scenario trajectories that change linearly over time (with a constant gradient over the model horizon).Therefore, the value of the primary driver at any given timestep equals:

Maximum trajectory value i
Total number of timesteps in model where i is the simulation number and t is the timestep (for example, t = 1, 2, 3… total number of timesteps in model).Using the Easter Island model as an example: if a maximum tree clearance value of 7 has been sampled for the given simulation, then its value after 500 timesteps would be equal to 500 × (7/1,500) = 2.333.The plausible trajectory funnels for each of the primary drivers are plotted in Supplementary Fig. 3-1.
To simulate experiment 2, secondary and tertiary driver trajectories are also activated using the following logic: • Secondary: primary driver active and secondary driver active and tertiary driver remains at default value and noise level remains at zero or Article https://doi.org/10.1038/s41893-023-01157-x • Tertiary: primary driver active and secondary driver remains at default value and tertiary driver active and noise level remains at zero or • All: primary driver active and secondary driver active and tertiary driver active and noise level remains at zero For each model, this specifically involved the following variables (Supplementary Table 3-2): • Lake Chilika fishery: (1) annual rainfall totals and mean near-surface air temperatures, as per IPCC (2013) climate change projections for the east coast of India; (2) price of fish per unit (Indian rupee per kg), leading to a more commercially oriented fishery, with an increasing number of fishers able to upgrade from traditional fishing boats to more intensive motorboats 33 .• Easter Island: (1) agricultural carrying capacity of the system, which enables a higher human population to be supported per unit of land cleared for agriculture; (2) the mortality rate of trees as a proxy for a disease-spread event.
• TRIFFID: (1) temperature-independent disturbance rate of vegetation coverage, that is, causes of forest clearance which are not directly linked to temperature changes (for example, deforestation).Note: due to the small size of the model, TRIFFID does not have a tertiary driver.• Lake Phosphorus: (1) rate of phosphorus recycling within the lake environment, (2) rate of phosphorus removal from the lake via sedimentation.
For the Lake Chilika and Easter Island models, these additional drivers are external forcings (similar to the primary driver).However, since the TRIFFID and Lake Phosphorus models are designed with only a single external forcing, additional drivers were also generated internally by altering parameters that operate on state variables.Whilst mathematically, internal and external forcings are fundamentally different things, both potentially impact the state of the system and, ecologically, changing internal model parameters can act as a proxy for an external process causing that change.For example, in the Lake Phosphorus model we have a system with a bifurcation in one dimension of slow external forcing (α) and we additionally vary internal parameters of the system (P recycling rate and P removal rate) as a proxy for, for example, anthropogenic disturbance impacting the species composition within the lake 63 .
Each of the additional driver trajectories are produced via the same approach as in equation ( 3): the Monte Carlo sensitivity analysis function in STELLA randomly samples a trajectory between the lower and upper bounds of a uniform distribution for each driver (Supplementary Table 3-2); in turn, the TIME function linearly increases the value of the driver from its default value to its sampled trajectory value by the final timestep of the model horizon.
To produce one secondary trajectory per simulation in the Lake Chilika model, the sampling of rainfall and temperature trajectories are connected along a linear gradient, ranging from no change to a combination of +30% rainfall change and +4.5 °C temperature change by 2081-2100 relative to 1986-2005 (as per representative concentration pathway 8.5 projections for the east coast of India 64 ).In STELLA, this is operationalized by the model variable 'climate change trend', with Monte Carlo sensitivity analysis in STELLA randomly sampling a value between 0 and 1 per simulation.As an illustration, if a value of 0.6 was to be sampled, then the change in rainfall by 2081-2100 (relative to 1986-2005) would equal 0.6 × 30(%) = 18%, whilst the change in temperature would equal 0.6 × 4.5(°C) = 2.7 °C.
To investigate experiments 3 and 4, the value of each primary slow driver is perturbed per timestep by randomly generated noise.We simulate a standard Wiener process to generate noise, equal to √dt × N(0, 1), where 'dt' equals change in time and 'N(0,1)' is a normal distribution with a mean of 0 and standard deviation of one.In turn, the product of the Wiener process is multiplied by the scaling factor σ, providing an overall level of noise to be added to the value of the primary driver on any given timestep.As per the future trajectories, the magnitude of σ is randomly sampled once per simulation from uniform distributions, with lower and upper limits shown in Supplementary Table 3 Therefore, building on equation (3) above, the value of a primary driver at any point in time in experiments 3 and 4 equals:

Maximum trajectory value i
Total number of timesteps in model Equation ( 4) as detailed above only refers to the 'uncoupled' noise simulations.Therefore, to explore the effects of 'coupled' noise, whereby the magnitude of noise increases with the growth in the primary driver, 20,000 simulations were run per model spread evenly between experiments 3 and 4, with the magnitude of noise coupled to the magnitude of the primary driver trajectory.Given the differences in the shape of the noise spectrums, we do not directly compare outcomes from the uncoupled and coupled noise simulations in this study.Instead, the purpose of modelling coupled noise is to ascertain whether worsening magnitudes of extreme events over time are associated with earlier ATDCs.In the coupled simulations, equation ( 4) is modified to: For experiment 3 (single slow driver plus noise), the runs were generated in STELLA 39 with the following logic: primary driver active and secondary driver remains at default value and tertiary driver remains at default value and noise active.For experiment 4 (noise plus multiple slow drivers), the logic used included: • Primary driver active and secondary driver active and tertiary driver remains at default value and noise active.• Primary driver active and secondary driver remains at default value and tertiary driver active and noise active.• Primary driver active and secondary driver active and tertiary driver active and noise active.

Time-series breakpoint detection
The identification of the timing of the ATDCs in the model runs was a two-step process.First, to ensure that we were only analysing model runs that had transitioned (collapsed) to quantitatively and qualitatively functionally different states (for example, fishery collapse, civilization collapse, forest dieback or lake eutrophication), we assessed whether model outcomes had crossed a predefined threshold at any point over the model horizon.For the three models which observe collapses in the outcome variable (Lake Chilika fishery, Easter Island and TRIFFID), model runs were considered to have reached a collapsed state if the outcome variable reached a value beneath 20% of its initial value at any point during the simulation.This demarcation is therefore representative of type-1 boundaries defined by ref. 40, with the numerical value of the boundary in line with the concept that fish stocks may be considered collapsed once their biomass falls beneath 20% of the biomass needed to maintain maximum sustainable yield 65,66 .In the case of the Lake Chilika fishery model, which has inbuilt social-ecological feedbacks that may trigger the recovery and later recollapse of the fishery 21,33 , we subset the time series to the period before the first timestep beneath 20% of the initial fish population.As we are only interested in the initial collapse, not Article https://doi.org/10.1038/s41893-023-01157-xsubsetting this time period would risk capturing subsequent collapses and recoveries in the analysis.
With lake eutrophication caused by an increase in phosphorus concentrations, a linear threshold beyond which we could be confident that the model had entered a qualitatively different state could not be identified.Therefore, as per the approach taken by ref. 67 for identifying abrupt events in global climate models, we adopted a ref.40 type-2 boundary to include only simulations which reached lake phosphorus concentrations greater than four times the standard deviation (s.d.) of the pre-ATDC time series.Therefore, runs of the Lake Phosphorus model which did not exceed this 4 × s.d.threshold were not considered to reach phosphorus concentrations sufficiently outside of the pretransition envelope of variability and were therefore excluded from our analysis.
The second stage of time-series breakpoint detection used the optimal breakpoint function of the R package strucchange v.1.5-2(ref.68) to identify ATDC dates in the time series that had successfully met the above qualifications (that is, reached an alternative state).As described in ref.21, the optimal breakpoint function finds the most substantial deviation from stability in classical regression models (Supplementary Fig. 3-2), whereby regressions coefficients shift from one regime to another.Therefore, the breakpoint date is taken as the most substantial deviation of the outcome variable en route to its qualitatively and quantitatively alternative state.

Boxplots and output graphs
For each of the experiments (Generation of future scenarios), boxplots were generated to visualize the distribution of breakpoint dates for each of the slow driver and noise level combinations (Figs.2-4).To standardize the comparisons between experiments, the normalized magnitude (0 → 1) of the primary trajectories (Supplementary Table 3-2) for each model was plotted on the horizontal axes.In turn, to visualize how the breakpoint dates change with the addition of secondary or noisy stresses over the range of the primary trajectories, model outcomes that tipped (Time-series breakpoint detection) were subset in the statistical software R between normalized primary trajectory values of 0.25-0.35,0.45-0.55 and 0.65-0.75.From here, boxplots for each of the driver combinations (for example, 'primary only' and 'primary and secondary') and primary driver subsets (for example, 0.25-0.35and 0.45-0.55)were graphed in R using the package ggplot (v.3.3.5;ref. 69).

Inclusion and ethics statement
This research is global in scope, using models that have been appropriately cited throughout.Roles and responsibilities were agreed amongst collaborators ahead of the research.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Fig. 1 |
Fig.1| Schematic overview of the framework developed to explore the influence of slow driver trajectories and/or noise on the timing of ATDCs.a, The four systems models simulated in this study (see section on Overview of systems models).b, Schematic representation of a system dynamics model (Lake Phosphorus model) with its external slow (blue and green) and noisy (red/orange) drivers depicted in colour (see Generation of future scenarios).c, Depiction of the four experiment types (section on Generation of future scenarios), ranging from changes in the primary baseline driver only (experiment 1), changes in all slow drivers and noise inputs simultaneously (experiment 4, where 'a' and 'b' represent noise profiles that are uncoupled or coupled to the primary driver

Fig. 3 |
Fig. 3 | The relationship between the breakpoint date and the primary slow driver (grey) for varying levels of uncoupled noise in the primary slow driver.Where normalised primary slow driver (σ) values ≤ 0.333 signify 'low noise' (yellow), normalised σ values > 0.333 and ≤ 0.666 signify 'mid noise' (orange), and normalised σ values > 0.666 signify 'high noise' (red; section on Generation of future scenarios).The normalized primary driver trajectories are apportioned into three discrete ranges: low-0.25-0.35,mid-0.45-0.55 and high-0.65-0.75.a-d, Subplots: Chilika model outputs, primary slow driver-fisher population growth (a); Easter Island model outputs, primary slow driver-tree clearance (b); TRIFFID model outputs, primary slow driver-temperature change (c); Lake Phosphorus model outputs, primary slow driver-phosphorus input (d).Model timestep units and boxplot dimensions are the same as in Fig.2; see Supplementary Table3-1 for the number of model simulations underpinning each boxplot.

Fig. 4 |
Fig.4| The relationship between the breakpoint date and the primary slow driver (grey) when weak and strong multiple driver trajectories are combined with weak and strong uncoupled noise (N).Where normalized driver trajectories (T) ≤ 0.333 are classified as weak and T > 0.666 as strong, similarly, normalized noise (σ) values ≤ 0.333 and > 0.666 are clasified as weak and strong, respectively.The normalized primary driver trajectories are apportioned into three discrete ranges: low-0.25-0.35,mid-0.45-0.55 and high-0.65-0.75.a-d, Subplots: the Chilika model, primary slow driver-fisher population growth, additional driver-climate change and fish price (a); the Easter Island model, primary slow driver-tree clearance, additional drivers-