Increasing evidence has shown that human activities have pushed the Earth’s system outside its resilient and sustainable state, which in turn poses an existential risk to humans [1]. Therefore, quantifying the impact of humans on the resilience of complex systems to environmental change has become an urgent yet challenging mission, known as Safe Operating Spaces [2], [3], as shown in Fig. 1(a). It quantifies safe boundaries for local stressors, such as exploitation rates, synchronized with climatic variations. Although conceptually powerful, a safe operating space is designed to model the properties of a highly abstract macroscopic system without considering the interactions between its dynamic components. However, real-world ecological systems, such as mutualistic networks (Fig. 1(b) [4]), contain interactions between many dynamic components governed by nonlinear dynamics [5], [6]. Consequently, the boundaries can be described as the border between the safe and collapse regimes [7], [8].

Studies on the boundaries between different operating spaces could help understand system resilience when considering human interventions such as exploitation and restoration [9], [10]. The intrinsic characteristics observed in numerous complex systems underscore their capacity to adapt their operations in order to preserve their fundamental functionality in the presence of errors and failures. Resilience loss is rarely predictable and is often irreversible, described by a sudden downgrade among alternative stable states [11], also known as regime shift [12], [13]. Regime shifts of ecosystems, such as mutualistic and food networks, are frequently triggered by disruptions such as species removal [14], [15], weight reduction [16], and individual trait variation [17]. Such critical transitions, in which a system abruptly shifts from one stable state to another, could fundamentally alter the structure and function of ecosystems and have profound implications. Consequently, investigating regime shifts and identifying tipping points is essential for understanding ecosystem stability and resilience. By identifying and understanding the drivers, patterns, and mechanisms of these regime shifts, ecosystem changes can be better predicted and managed to protect and maintain ecosystem functionality and sustainability. Such research also aids in formulating effective ecological policies and management measures to address the impacts of environmental change and human activities on ecosystems. Thus, studying regime shifts is crucial to achieving sustainable ecosystem management and biodiversity conservation.

Therefore, numerous efforts have been made to monitor system states, forecast extinction chains, and identify the most vulnerable species [18], [19] by leveraging empirical data and network theory. Gao et al. [16] developed a dimension-reduction method that derives an effective one-dimensional dynamic by making full use of the topologies and dynamics of large-scale networks. Subsequently, a reduced two-dimensional (2D) model of bipartite mutualistic networks was proposed to predict the tipping points [20]. Laurence et al. [21] established a spectral dimension-reduction framework with justifications for the necessary effective dimensions for studying critical breakdowns. Tu et al. [22] considered systems with heterogeneous dynamic mechanisms. Based on the interactions between network topology and dynamics, Morone et al. [23] discovered a k-core predictor of structural collapse in mutualistic ecosystems. Zhang et al. [24] introduced a scaling factor for recovery rates that places various systems with different dynamic parameters, network structures, and state perturbations on the same scale, enabling us to compare the distances to tipping points across multiple ecosystems.

Despite recent efforts, understanding the resilience of a networked system when a part of it is undergoing human overexploitation or other interventions remains an open question. This includes the intensive exploitation of certain fish species, high-value timber, and comprehensive development, all of which can result in significant environmental damage. This is because the nonlinear equation format is not identical and depends on the number of nearest neighbors when partial species are under exploitation. These limitations are rooted in a theoretical gap: Existing approaches usually consider a system that satisfies a mean-field approximation when all species are under a similar threat and will survive or become extinct simultaneously. However, depending on the exploitation rates and targeted ranges, the use of the mean-field approach to predict extinct species remains challenging.

In this study, we developed an analytical tool to quantify the response of an ecosystem to human exploitation based on a dimension-reduction method. This method allows us to quantitatively specify and predict species-level tipping points and extinction chains (Fig. 1(c)). Such analytical tools contribute to paving the way for studying the resilience of abundant natural, social, and engineering networks, predicting tipping points, and providing early warnings, thus guiding the conservation and restoration [25], [26] of real large-scale systems. By applying the tool to real-world ecological networks, we revealed a previously unrecognized regime shift encompassing five distinct states: safe, partial extinction, bistable, tristable, and collapse, which is significantly different from the qualitative or downscaled boundaries (Fig. 1(a)) in Refs. [3], [5]. Based on a phase diagram incorporating these five regimes, we developed an active reintroduction strategy to restore and revive degraded ecosystems. The ecological engineering approach, aimed at promoting human and ecological well-being, is further corroborated by the observed restoration data from eelgrass-fish mutualistic networks in the midwestern Atlantic coastal lagoons.

We use a generic bipartite mathematical model to evaluate the dynamic evolution of species abundance driven by intrinsic growth, intraspecific and interspecific competition, mutualistic interactions between species [20], [27], [28] (e.g., plants and pollinators), and external exploitation stresses:where the $N$-dimensional system composes $N^{\mathrm{P}}$ plants and $N^{\text{A}}$ pollinators with abundance $P_i\left(i=1,...,N^{\text{P}}\right)$ and $A_i\left(i=1,...,N^{\text{A}}\right)$, respectively. $N^{\mathrm{P}}$ is the total number of the plants and $N^{\text{A}}$ is the total number of the pollinators. $f_i^{\mathrm{P}}$ and $f_i^{\mathrm{A}}$ are the exploitation rate on certain species i; ${\sigma }_i^{\mathrm{P}}$ and ${\sigma }_i^{\mathrm{A}}$ are the binary variables to indicate the exploited species. For simplicity without loss of generality, we apply the same exploitation rate $f$ to every exploited species. ${\alpha }_i^{\mathrm{P}}\mathrm{a}\mathrm{n}\mathrm{d}{\alpha }_i^{\mathrm{A}}$ are the intrinsic growth rates in the absence of any competition or external stresses; ${\beta }_{\mathit{ii}}\mathrm{a}\mathrm{n}\mathrm{d}{\beta }_{\mathit{ij}}$ represent the intraspecific and interspecific competition strength, respectively, depending on environmental conditions (typically, intraspecific competition weights are much higher than interspecific competition weights, ${\beta }_{\mathit{ii}}\gg {\beta }_{\mathit{ij}}$); ${\gamma }_{\mathit{ij}}={\epsilon }_{\mathit{ij}}{\gamma }_0/K_i^{\delta }$ is the benefit of mutualism regulated by a half-saturated constant $h$, denoting the essence of mutualistic effects. γ₀ is the level of mutualistic strength. ${\epsilon }_{\mathit{ij}}=1$ suggests that species $i$ directly interacts with species $j(j\ne i)$; otherwise, ${\epsilon }_{\mathit{ij}}=0$. $K_i$ is the number of species that benefit from the interaction with species $i$; $\delta$ determines the impact of the tradeoff between interaction strength and the number of interactions, that is, the network topology’s influence on mutualistic interactions, which is normalized in the range [0, 1]. Throughout the study, $t$ denotes time, the initial conditions of species abundance are randomly chosen, and the system is considered to have reached equilibrium when the difference between iterations is less than 0.00001. In this study, we apply the nonlinear dynamics in Eq. (1) to 26 benchmark mutualistic networks from the Web of Life database, composed of four categories: plant–pollinator, seed–disperser, anemone–fish, and plant–ant networks (Table S1 in Appendix A).

The concept of safe operating space for a mutualistic network is illustrated in Fig. 1(c). Within a safe regime (Fig. 1(c-i)), species abundance could recover after removing local stressors such as harvest or environmental change combinations. However, once local stressors or environmental conditions exceed the functional tipping point $f_{c1}$, the most endangered species become extinct, followed by an extinction chain as the local stressor continues to escalate (partial extinction regime). The system stays in the partial extinction regime because of the high local stressor and proper environment (Fig. 1(c-ii)), or the low local stressor and hostile environment (Fig. 1(c-iii)). When the local stressor, such as the exploitation rate, exceeds the ecological tipping point $f_{c2}$—the “point of no return”—the entire mutualistic network could be severely degraded or even collapse (Fig. 1(c-iv)) and cannot spontaneously recover, even after unloading the exploitation stress. Here, we defined the number of layers in an extinction chain as its depth and the largest number of species in one layer as its width.

We present a systematic tool for predicting tipping points and extinction in an ecosystem that considers both network structure and dynamics. The various fundamental processes highly depend on the exploitation range and rate; however, the beginning and endpoint are more discernible. Therefore, the general concept guiding our predictions is to retrace the extinction chains. The preliminary identification of species facing extinction is accomplished through network topology analysis, with particular emphasis placed on those detached from the largest cluster component following the removal of the exploited species. Subsequently, by employing dimension-reduction techniques, we estimate the system states at different stages along the extinction chains, which in this case refers to species abundances.

There are four steps for tipping points and extinction prediction:

Step 1: Estimate the extinct species set according to network topology.

The extinction set $V=V_{\mathrm{exploited}}\cup V_{\mathrm{co}-\mathrm{extinct}}$ is obtained by excluding the exploited species set ${\mathit{V}}_{\text{exploited}}$, and identifying species detached from the largest cluster component as co-extinct species set ${\mathit{V}}_{\text{co}-\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{t}}$, including $\mathit{Z}$ species.

Step 2: Estimate the final system state when all species in extinct.

As all species $m$ $(m=1,...,Z)$ in $\mathit{V}$ extinct, that is ${{\sigma }_m^{\mathrm{P}}f}_m^{\mathrm{P}}{P}_m=0$, ${{\sigma }_m^{\mathrm{A}}f}_m^{\mathrm{A}}{A}_m=0$ in Eq. (1), we estimate the abundances of other species using the dimension-reduction method. Detailed information is provided in Appendix A, and we develop the reduced four-dimensional (4D) model as

Let $P_{\mathrm{e}\mathrm{f}\mathrm{f}}$ and $A_{\text{eff}}$$A_{\text{eff}}$ represent the effective abundances of plants and pollinators, respectively; $P_{\text{eco}}$ and $A_{\text{eco}}$ describe how different types of species influence each other and the strength of such interlayer effect. ${\alpha }^{\text{eff-P}}$, ${\alpha }^{\text{eff-A}}$, ${\alpha }^{\text{eco-P}}$, and ${\alpha }^{\text{eco-A}}$ are the intrinsic growth terms. The dynamical interactions of species within plants/pollinators are also governed by their nearest neighbors through the interaction terms; here, the intraspecific competition strength of species is represented by ${\beta }_{\mathrm{s}}^{\text{eff-P}}$, ${\beta }_{\mathrm{s}}^{\text{eff-A}}$, ${\beta }_{\mathrm{s}}^{\text{eco-P}}$, and ${\beta }_{\mathrm{s}}^{\text{eco-A}}$, which is analogous to ${\beta }_{\mathit{ii}}$ in the $N$-dimensional model. Effective interspecific competition strengths are ${\beta }^{\text{eff-P}}$, ${\beta }^{\text{eff-A}}$, ${\beta }^{\text{eco-A}}$, and ${\beta }^{\text{eco-A}}$, while interacting mutualistic effects among different systems denoted by ${\gamma }^{\text{eff-P}}$, ${\gamma }^{\text{eff-A}}$, ${\gamma }^{\text{eco-P}}$, and ${\gamma }^{\text{eco-A}}$. With the effective states, the system states $P_i$ and $A_i$—representing species abundance—could be estimated by decoupling the $N$-dimensional system. Eq. (1) could be converted as

We evaluate the equilibrium state when $\mathrm{d}{P}_i/\mathrm{d}t$ = 0 and $\mathrm{d}{A}_i/\mathrm{d}t$ = 0, ensuring the right-hand side of Eq. (3) is zero with $P_i\ge 0,A_i\ge 0$. Hence, the system states could be solved directly. This approach allows us to dynamically update the extinction set whenever a species’ abundance falls below the threshold $\Theta$. (For more detailed influences of parameters such as $\Theta$, refer to Fig. S1 in Appendix A) The solutions align well with those obtained from the original $N$-dimensional model, demonstrating a mean error less than 0.1, as shown in Appendix A Fig. S2. Additionally, a comparison of the system-state estimation with existing methods is presented in Appendix A Figs. S3 and S4.

Step 3: Derive the exploitation rate when the species in extinct.

We sort the species in $\mathit{V}$ based on their degree of endangerment, as indicated by their original abundances arranged in descending order. This sorting strategy underscores that species with higher abundances and positioned further from the extinction threshold have a better chance of survival. When the extinction procedure approaches an end, the exploitation rate is estimated with the following procedure:

(1) Final state approximation: Gradually remove species from $\mathit{V}$ until any effective state ($P_{\mathrm{e}\mathrm{f}\mathrm{f}}$ or $A_{\text{eff}}$) rises above 0, indicating an almost-final state of the system. This step helps in determining the critical juncture before system collapse.

(2) States estimation for surviving species: At the end of the extinction chain, we estimate the abundances of surviving species by solving Eqs. (2), (3). Here, the abundances of all species $m\in \mathit{V}$ are set to 0, focusing solely on those species that persist.

(3) States approximation for extinct species: For species $m\in \mathit{V}$, we approximate their abundances while holding the abundances of surviving species constants.

(4) Exploitation rate estimation: Calculate exploitation rate $f$ with Eq. (1), utilizing the approximated final effective states and system states. Considering potential inaccuracies arising from extinction order and system-state estimation, we define $f$ as the maximum of all calculated rates for both plants and animals, $f=\mathrm{m}\mathrm{a}\mathrm{x}[f_1^{\mathrm{P}},...,f_{N^{\text{P}}}^{\mathrm{P}},f_1^{\mathrm{A}},...,f_{N^{\text{A}}}^{\mathrm{A}}]$, and then round it up as $f=\lceil f\rceil$. In single-species exploitation scenarios, $f=f_{c2}$.

Step 4: Peel back to get the extinction chain from the final degraded state.

Gradually decrease the exploitation rate from, and estimate the abundances for species in set. The surviving species’ states are set as constants, thus releasing the computational complexity, now that only species’ abundances need to be tracked. The revived species, if any, comprise the latest layer of the extinction chain with the corresponding exploitation rates as their tipping points and then obtain the new extinction set.

Repeat until all species in are assigned in the extinction chain from the bottom up with their tipping points.

In addition, we predict extinction chains based on network structure and dynamics such as species abundance, average vertex strength, the mutualistic species ranking method MusRank [29], and nodal contribution to nested architecture [19]. Species abundance is obtained by the four-dimensional model free from exploitation. The average vertex strength $(\overline{s_i})$ is defined as the average of the species’ link strengths, $\overline{s_i}={\sum }_{j\ne i}^N{\epsilon }_{\mathit{ij}}(K_i+K_j)/{\sum }_{j\ne i}^N{\epsilon }_{\mathit{ij}}$. The accuracy of predicting extinction chains using such criteria when capturing one or two species is shown in Appendix A Figs. S5 and S6. The results of the proposed framework, as shown in Appendix A Figs. S5–S9 and Tables S2–S4, certainly outperform those methods.

We investigated the tipping points and extinction chains of the ecological network denoted as M_PL_017. The network comprises 25 plants, 79 pollinators, and 299 interactions obtained from Bristol, United Kingdom. The network topology is illustrated in Fig. 2(a). To investigate how an ecosystem reacts to overexploitation, we introduce the exploitation term to a certain fraction of the species until the system reaches a new stable state and the abundance of extinct species falls below the threshold $\Theta$. Subsequently, we withdraw the exploitation rate and observe whether the abundance recovers.

Mathematically, we start from the initial condition, $P_i\left(0\right)\gg 0$ and $A_i\left(0\right)\gg 0$ (healthy high state), numerically solve the ordinary differential equations (ODE) equations in Eq. (1) with given $f$, and then we set all ${\sigma }_{\mathit{i}}$ to be 0 (no exploitation) and obtain its final steady state. Figs. 2(a)–(d) show the final steady states as a function of the exploitation rate $f$ if we only exploit a single species. When $f$ is small, the system always recovers its original state. However, co-extinction occurs as $f$ increases. For example, when capturing species 1, ten species spontaneously go extinct if $f\ge f_{c1}=3.7$, while the rest of the species can still recover to a high state close to their original state. As we continue to increase $f$, species 1 becomes extinct. The widths and depths of the extinction chains are 10 and 2, respectively. Similarly, we show the extinction chains when we exploit species 4 and 15 in Figs. 2(c) and (d), respectively. The tipping points and shapes of the extinction chains vary with the specific exploited species. More importantly, the generalists, such as species 1, with more interactions and higher abundance, often cause more co-extinction than others, such as species 4 and 15, not to say specialists.

Our numerical findings suggest that an extinction chain reveals interdependencies for species survival. Identifying these interdependencies could provide insights into why certain species face extinction, even in the absence of direct exploitation. Consequently, detecting extinction chains could guide us in considering the most endangered species rather than only the exploited ones. Note that this process generalizes resilience studies focusing on node/link removal or weight reduction [16], as the node removal process is analogous to the special case when $f>f_{c2}$. Such continual and gradual exploitation mimics the more widespread exploitation of the ecosphere and validates that the extinction chain can be triggered by both local extermination of species or a significant decline in abundance, corresponding to functional extinction [18].

We also present a common scenario in which a portion of the ecosystem is under exploitation. Fig. 2 shows examples of the exploiting 20% and 80% of all species. We discovered that apart from the exploited species, other species without direct exploitation may also become extinct [31] owing to the loss of mutualistic interactions, forming extinction chains with various patterns. The co-extinct species are highlighted in the network (Figs. 2(a), (e), and (i)). When only a few species are exploited, extinction chains may not appear, and the ecosystem could maintain a safe or partial extinction state, even under high perturbation rates. By contrast, as the exploitation range increases, for example, to 80%, as shown in Figs. 2(i)-(l), the system may collapse even under a small exploitation rate. We present additional findings involving diverse networks characterized by varying growth rates and exploitation levels in Appendix A Figs. S10-S14 and Table S2.

By examining the exploitation process of 26 benchmark mutualistic networks, we encountered universal co-extinction phenomena and extinction chains in mutualistic ecosystems of different types, scales, and geographical locations. However, within the same network, exploiting various species results in distinct sizes and shapes of extinction chains. Moreover, their depths and widths varied significantly among the various networks (Figs. S15-S17 in Appendix A). Although ecosystems show resilience and revival under minor exploitation rates and ranges, they cross the safe operating space boundary and fall into undesirable regimes under augmented stress. Enlarging the exploitation rate triggers co-extinction, whereas an increased exploitation range scales up ecosystem degradation, decreases tipping points, and expands the width of extinction chains. Indeed, the shape of the extinction chain is vital for early warning of co-extinction and corresponding conservation strategies. For example, a flat but wide extinction chain requires more urgent safeguards than a deep but narrow chain because the collapse of the entire system may be abrupt and unexpected.

To effectively predict diverse exploitation patterns in complex systems, we extend the dimension-reduction approach [16] to bipartite networks that model the topology and dynamics of mutualistic networks. By estimating the influence of the exploited species on the overall species abundance, we derive an analytical method to accurately predict the co-extinction species and their tipping points in extinction chains (Figs. S1, S5, S10-S14, S18, and S19 in Appendix A). For simple cases, such as single species (Figs. 2(a)-(d)) or two-species exploitation scenarios (Figs. S6 and S20 in Appendix A), the order of the species in the extinction chain and its shape can be accurately predicted. In the multiple-species cases (Figs. 2(e)-(l) and Fig. S7 in Appendix A), extinction chain prediction becomes more challenging as the exploitation range or rate increases, and nontrivial phenomena appear, such as bistable and tristable regime shifts. Despite the extreme complexity of multiple-species exploitation, we precisely forecast the most endangered species (Fig. S7(e)), which can be used as an early warning signal of cascading extinction in mutualistic networks. All extinction chains presented in Fig. 2 were correctly predicted. The error estimation and analysis are shown in Appendix A.

To reveal the regimes of exploited ecosystems and define a safe operating space, we methodically expand the exploitation rate and range of mutualistic systems. Regime shifts occur as ecosystems are pushed from the “basin of attraction” of states, leading to “folded” multiple alternative stable states separated by unstable equilibrium [32]. When a part of the system is exploited, the final system state is determined not only by the exploitation rate and range but also by the specific exploited species. To reveal the connections among the exploitation rate, exploitation range, and system state, we performed extensive simulations exploiting randomly selected species, given specific exploitation rates and ranges (Fig.3; Figs. S5, S16, and S17). Different simulations led to distinct system states owing to the variations in the exploited species. From all the simulation results, we calculate the probabilities of the system being in a state of safety, partial extinction, or collapse.

Under mild perturbations with low exploitation rates and ranges, the system remains in the safe or partial extinction regime even when all species are exploited (Fig. 3). In particular, given a certain exploitation range and rate, the system may fall into a multistable area (Fig. 3(b)). Fig. 3(e) shows the corresponding degradation process: After decreasing from a single stable state, such as the safe or partial extinction (one basin of attraction), the system could enter into either a partial or a majority extinction state (bistable regime)—on some occasions, even collapsing with an increasing exploitation range (tristable regime)—and then gradually degrading to either the majority extinction or collapse state (bistable regime).

Ultimately, an irrational exploitation range leads to collapse. A sufficiently high perturbation rate (Figs. 3(c) and (f)) leads to more severe damage and deteriorates the system state to easily shift from safe to partial extinction, bistable (either majority extinction or collapse), and eventually to a certain collapse. The detection of bistability and tristability could help us understand the resilience patterns and regime shifts of mutualistic networks, which have profound implications for the restoration of degraded ecosystems [33], [34], [35] and for building sustainable systems.

When ecosystems suffer from drastic changes or hostile environments, such as those resulting from eutrophication, excessive interspecific competition leads to species extinction, even in the absence of exploitation stress (Figs. S21–S23 in Appendix A). The conventional result would be that more species decline, which gradually leads to a loss in biodiversity under fierce competition (large ${\beta }_{\mathit{ij}}$) as the exploitation rate and range [19], [29], [36], [37] increase. Counterintuitively, under specific scenarios (Fig. 4(a)), the decline in species diversity is mitigated at certain exploitation rates. Moreover, there are specific counterintuitive regions in the phase diagram where exploiting more species could even help increase ecosystem biodiversity.

One plausible explanation for this phenomenon is rooted in ecological principles, according to which a reduction in the abundance of dominant species could mitigate interspecific competition within the ecosystem [38]. In many ecosystems, dominant species exert strong competitive pressure, thereby reducing overall biodiversity. However, strategically exploiting these dominant species could alleviate competitive pressure by decreasing their abundance. This creates opportunities for other species to thrive and boosts the overall biodiversity. This phenomenon has also been observed in various other ecological scenarios, including food webs and predator-prey relationships [39], [40]. Protecting top predators is crucial because they regulate lower trophic levels, similar to how humans conduct harvesting activities.

To elucidate the mechanisms underlying these counterintuitive phenomena, we conducted an in-depth analysis of the evolution of species abundance with an increasing exploitation range. In this analysis, specialists, characterized as leaf nodes within the system, are considered indicators because of their lower degree of connectivity and heightened sensitivity to critical transitions [41]. As shown in Fig. 4, the average species abundance of leaf nodes is positively correlated with cluster size, whereas the average species abundance of non-leaf nodes is negatively correlated with the number of co-extinct species. In counterintuitive cases, the average species abundance of leaf nodes shows a similar initial decrease followed by an increase in cluster size as the exploitation range expands with losses in average species abundance, and we determined that the anomalous twist occurs when the average species abundance with exploitation surpasses that under exploitation-free conditions. We predict the counterintuitive area by focusing on system responses when all species are being exploited. More specifically, the counterintuitive area is predicted when the average species abundance of the leaf nodes surpasses 105% of the original scenario before exploitation. When most of an ecosystem is exploited in a hostile environment, a regime shift between degradation and disastrous collapse becomes more evident (Figs. S21-S23 in Appendix A).

We mapped multiple simulations and acquired a complete phase diagram of the M_PL_017 network (Fig. 5). The ecosystem exhibited first-order phase transitions on the partially extinction-bistable, bistable-tristable, and tristable-collapse boundaries, indicating tipping points from safe to degraded. Moreover, regardless of the initial conditions, the final status of the ecosystem in the exploitation rate-competition intensity space can be predicted solely according to the combinations of exploitation rate $f$ and competition intensity $\beta$. When ecosystems are already degraded or have collapsed owing to overexploitation, we seek practical strategies to restore them to prosperity using the revealed five-phase diagram. Using the derived exploitation rate-cluster size phase diagram, we provided systematic schemes for ecosystem restoration, which include both the reintroduction of extinct species and reduction in the exploitation range or rate. Furthermore, we simulated active restoration strategies on network M_PL_017, demonstrating the ongoing reintroduction of 5% of the original biomass for all extinct species. Threats were effectively contained by decreasing the exploitation range (Fig. 5(b)), exploitation rate (Fig. 5(f)), or both (Figs. 5(c)–(e)). These strategies successfully steered the damaged system back to safety through diverse trajectories, demonstrating trends similar to those of the degradation process caused by exploitation.

We also verified the proposed model (Eq. (1)) and an active restoration strategy using actual revival data [30] for the midwestern Atlantic coastal lagoons (Fig. 6) over the last 20 years. The dataset includes information on eelgrass cover, reintroduced amounts of eelgrass, and fish biomass. We focus on reintroduced species of eelgrass (Zostera marina) and ten fish species (Bairdiella chrysoura, Anchoa mitchilli, Lagodon rhomboides, Leiostomus xanthurus, Sygnathus fuscus, Orthopristis chrysoptera, Menidia menidia, Centropristis striata, Paralichthys dentatus, and Chilomycterus schoepfii) most recorded over the years. Reintroducing rates are obtained using annual seeded area and eelgrass cover, as in Fig. 6(b), while parameters such as intrinsic growth rate $\alpha$ and intraspecific and interspecific competition strengths $\beta$ are selected using a genetic algorithm, fitted to the time series of seagrass cover and fish biomass in the initial five years. The parameters are processed to minimize the mean square error (MSE) compared with real data. Initial abundance of fish as well as the reintroduced species of eelgrass are randomly chosen from (0, $\Theta$), where the threshold for species extinction is $\Theta =0.001$$\Theta =0.001$. The system cannot automatically recover without reintroduction at a given initial abundance. As we conduct the reintroduction strategy, the normalized numerical simulation results match the real data quite well in terms of eelgrass cover (Fig. 6(c), correlation coefficient: 0.9848; MSE: 0.0042) and fish biomass (Fig. 6(d), correlation coefficient: 0.9222; MSE: 0.0166). The simulation results for seagrass cover and fish biomass recovery aligned well with the data collected from the restoration project, thus, the feasibility of the proposed model is verified. In addition, we found that, in two of the 26 networks, adding a certain leaf node [42] partly enhanced the resilience (Fig. S24 in Appendix A). This suggests the potential application of active reintroduction or alien invasion in preventing massive extinctions (Table S5 in Appendix A).

In summary, we constructed a dynamic mutualistic model that incorporates internal driving mechanisms, interspecies interactions, human activities, and external environmental disturbances. We propose a dimension-reduction method designed to accurately predict tipping points and extinction chains under human exploitation, effectively bridging the gap between complex systems theory and practical predictions of natural ecosystem outcomes.

Although ecosystems exhibit resilience and revival under minor exploitation rates and ranges, they cross the safe operating space boundary and shift to undesirable regimes under augmented stress. This study introduces a comprehensive phase diagram that categorizes the exploitation rate range space into five distinct regions: safe, partial extinction, bistable, tristable, and collapse. The present study sheds light on the analysis of extinction evolution, the identification of extinction factors, and understanding of the mechanisms behind the counterintuitive functioning of exploited ecosystems.

Furthermore, our research extends beyond a theoretical exploration to offer practical applications in ecological engineering. From a human intervention perspective, the revealed five-phase diagram guides systematic restoration schemes for ecosystems already degraded or on the edge of collapse due to overexploitation. Particularly, we focus on transitioning ecosystems from zones at risk of extinction to safe areas through strategic control measures. This includes the reintroduction of extinct species and a targeted reduction in exploitation rates and ranges to recover vulnerable ecosystems.

This study mapped simulations to thoroughly understand phase transitions in the networked M_PL_017 ecosystem. To validate the approach, the proposed model and active restoration strategy were quantitatively verified using data from midwestern Atlantic coastal lagoons over the last 20 years. This not only underscores the practical relevance of the study but also reinforces its real-world applicability. The mutualistic system, involving eelgrass and ten fish species, showed simulation results that aligned well with data collected from the restoration project, particularly in terms of seagrass cover and fish biomass recovery. This robust, application-oriented approach contributes significantly to ecosystem management and restoration.

With the continuous expansion of open-source datasets for diverse ecosystems, the proposed prediction and rescue methods can be refined and expanded. This includes broader coverage of larger and more varied mutualistic ecosystems, as well as potential applications to different ecosystem types characterized by commensalism, parasitism, competition, and beyond. Such endeavors promise to deepen our understanding of ecological conservation and restoration, paving the way for more effective ecosystem management strategies in the future.

This work was supported by the National Key Research & Development Program of China (2022ZD0119601), the National Natural Science Foundation of China (62225306 and U2141235), and the Guangdong Basic and Applied Research Foundation (2022B1515120069). Jianxi Gao acknowledges the support of the US National Science Foundation (2047488) and the Rensselaer-IBM Artificial Intelligence Research Collaboration. Guangwei Wang acknowledges the support of China Scholarship Council (202206160043).

Guangwei Wang, Xueming Liu, Ying Xiao, Ye Yuan, Linqiang Pan, Xiaohong Guan, Jianxi Gao, and Hai-Tao Zhang declare that they have no conflict of interest or financial conflicts to disclose.