Landscape quantifies the intermediate state and transition dynamics in ecological networks

Understanding the ecological mechanisms associated with the collapse and restoration is especially critical in promoting harmonious coexistence between humans and nature. So far, it remains challenging to elucidate the mechanisms of stochastic dynamical transitions for ecological systems. Using an example of plant-pollinator network, we quantified the energy landscape of ecological system. The landscape displays multiple attractors characterizing the high, low and intermediate abundance stable states. Interestingly, we detected the intermediate states under pollinator decline, and demonstrated the indispensable role of the intermediate state in state transitions. From the landscape, we define the barrier height (BH) as a global quantity to evaluate the transition feasibility. We propose that the BH can serve as a new early-warning signal (EWS) for upcoming catastrophic breakdown, which provides an earlier and more accurate warning signal than traditional metrics based on time series. Our results promote developing better management strategies to achieve environmental sustainability.

simulate forward using MATLAB's ode45 solver, which is a high-order algorithm for solving ODEs.The initial conditions involve randomly and uniformly distributing points within a feasible high-dimensional range (all species abundances are non-negative and bounded).We considered the system to have reached stability when the Euclidean distance between species abundance vectors in adjacent steps is less than 1e-7.As for the landscape, the weights for each attractor are obtained based on the proportion of points falling on each after uniform sampling.

C. Generalized dimension reduction approach for landscape
Here we show in detail the mathematical derivation of our dimension reduction approach and its advantages.Our approach is an improvement over previous work that proposed a general framework to perform a dimension reduction for landscape [3].We have modelled that the multidimensional stationary probability density function is where the symbol meanings are detailed in the main text.We initially applied this approach and selected the first two eigenvector directions for projections, to maximize the variance of population probability distribution after projection.However, we found that it is not very suitable for this bipartite system.First, the biological meanings of two coordinate variables are vague.The first principal component (PC1) and the second principal component (PC2) both are composed of plants and pollinators with equally important and indispensable weights (Figs Da-b and Ec-d).They further result that the reduced landscape can not adequately present the information we need.For example, the transition paths between stable states can provide the order information of species extinction and recovery, and we expect that the projected paths in the reduced landscape can roughly retain the information, at least about plant and pollinator population.But we cannot get this by applying this method.We also noted that PC2 is perhaps redundant, where its value of the intermediate state may be larger (Fig Dc) or smaller (Fig Dd) than that of the low state.This is related to the negative weights of the components that make up PC2 (Figs Da-b).Moreover, when we perturbed the value of parameter, we found that two coordinate variables (PC1 and PC2) exhibit a drastic change (Figs Da-b), which directly leads to a significant difference between their corresponding landscapes (Figs Dc-d).Here we changed the value of κ and fixed other parameters under default values.We showed the results of weights as well as their landscapes at κ = 1.09 (Figs Da and Dc) and κ = 1.1 (Figs Db and Dd).We concluded that the previous approach may undergo significant changes in variable weights and landscapes after small perturbations of parameters, i.e., it is not considered robust enough.
To address the above shortcomings and to take into account the plant pollinator-networks which are composed of two sets of nodes, we developed a novel dimension reduction approach for landscape.Inspired by hierarchical principal component analysis [4], the variables in this bipartite Network Can be naturally divided into two segments (plants and pollinators).We took the derivation of the plant population (denoted by P = [P 1 , . . ., P S P ] T ) as an example, and the conclusions of pollinator population (denoted by A = [A 1 , . . ., A S A ] T ) can be simply reached by replacing P with A. We first obtained its stationary probability density function December 20, 2023 2/34 through integrating A in Eq.S1: Then the expectation of component P i is derived as follows: ϕ j µ j (P i ) .
So the expectation of vector P is ϕ j µ j (P ), where µ j (P ) is the j th attractor's Gaussian mean of the plant segment.We are more concerned with the covariance matrix between variables, whose element σ pq (p th row, q th column) can be calculated by Cov (P p , P q ) = E [(P p − EP p ) (P q − EP q )] = E (P p P q ) − EP p • EP q December 20, 2023 where Σ j (P p , P q ) is the (p, q) element of the j th attractor's Gaussian covariance matrix of the plant segment (N P × N P submatrix in the upper left corner).The last equation comes from the definition of covariance: Σ j (P p , P q ) = E j (P p P q ) − µ j (P p ) µ j (P q ) Therefore, the covariance matrix is obtained from: where ⊗ is the Kronecker product.We can directly substitute A into B obtain the equations satisfied by pollinators.Then we separately selected the PC1 of plant covariance matrix and pollinator covariance matrix as new coordinates, denoted by plant PC1 and pollinator PC1 respectively.They are computed using singular value decomposition (SVD), and essentially the eigenvector corresponding to the largest eigenvalue of each two matrices.
Although current deviation form is similar to that in the previous approach [3], we will point out that they are fundamentally different.Our approach is more applicable to this bipartite system.First, for high-dimensional systems where both N P and N A are very large, it is more computable to solve the eigenvector corresponding to the largest eigenvalue of two lower dimensional matrices (N P × N P and N A × N A ), compared to solving the eigenvectors corresponding to the first two eigenvalues of a high-dimensional matrix ((N P + N A ) × (N P + N A )). Next, our novel coordinates are composed of linear combinations of plants and pollinators respectively, which ensures interpretability (Figs De, Df, Ea and Eb).We can roughly infer the overall plant behavior from plant PC1 and pollinators as well.Meanwhile, the almost non-negative weights indicate that, in the reduced landscape, large values of plant PC1 and pollinator PC1 correspond to the high state, while small values correspond to the low state.The same perturbation demonstrates the stability of our approach, where neither the weights in two coordinate variables nor the reduced landscapes change significantly (Figs De-h).
Furthermore, for wider range of parameters (κ ∼ U [1.01, 1.12]), we depicted changes in weights of each species that make up plant PC1 (Fig Ea  c).Meanwhile, they are more interpretable, which can describe overall evolution of plant and pollinator population but the individual PC1 can not.It is intuitive to choose two coordinates to visualize the landscape, which can provide as much detail as possible.We also noticed that PC2 varies dramatically, and thus not easy to choose the reference coordinates (Fig Ed).Indeed, by demonstrating the absolute dominance of plant PC1 in Σ(P ), pollinator PC1 in Σ(A), PC1 in Σ(P, A), combined with the little contribution rate of PC2 in Σ(P, A) (Fig Ee ), we reemphasized that PC2 does not provide valid information and also showed the superiority of our approach.

D. Calculation of transition action and MFPT and their consistency with barrier height (BH)
means that two stable points are directly connected with a straight line in high-dimensional space.On the other hand, larger L can lead to computational difficulties and potentially unstable calculations.Based on computational considerations, we opt for an appropriate value of L and keep it fixed for subsequent analysis.We demonstrate the impact of varying L on the transition action and computational time (Fig B ), and we choose L = 20 to balance the principle of least action and computation time.
Considering that the saddle point calculated in the high-dimensional system may not correspond well to that calculated in the reduced landscape, we chose RBH as a substitution indicator for BH in the main text.Moreover, the point with the local minimum potential energy in the high-dimensional system is nearly always projected onto the attractor (also with local minimum potential energy) in the reduced landscape.This correspondence implies that RBH can adequately describe the evolution of the system.Here we confirmed the strong relationship between ∆S and RBH (Fig Fa ).We can observe their positive correlation, in fact, they show very similar trends with κ varying (Figs 4A-F).The anomaly occurs at smaller κ, during which the weight of high state is much greater than that of the intermediate and low states.This complete dominance of the high state will make the calculation of RBH biased.
To calculate the MFPT in our multidimensional system, it is challenging to solve f (x) • ∇τ + d∇τ • I • ∇τ = −1.A feasible alternative solution is based on the trajectory simulation from equivalent high-dimensional Langevin equations dx(t) dt = f (x(t)) + ζ(t).We treated a supra-sphere as a proxy of an attractor, where its center is quantified from the stable-state mean in our TME method, and radius is the hyperparameter r we specified.For each experiment, we randomly picked an initial point in the supra-sphere and then simulated a trajectory through the Euler-Maruyama method.The timer of first passage time (FPT) is started when it initially leaves x i , and stopped when it first arrives in x j .We repeated 100 experiments and then averaged their FPT to approximate the MFPT.We simulated for long enough time, trying to capture the transition process.We recorded MFPT of three feasible transitions under different parameter κ (κ = 1.02 ∼ 1.07: intermediate to high, κ = 1.02 ∼ 1.06: low to high, κ = 1.08 ∼ 1.16: low to intermediate, all at 0.01 interval).Under Euclidean distance, we tuned the hyperparameter r, and showed the change of MFPT with varying κ (Fig Fb ).We chose r = 0.5 in our main text, but similar tendency indicates that its selection does not affect our general conclusions.As κ rises, the time required to transit to the state with higher abundance is longer.We clarified the approximate linear relationships between RBH and the logarithm of MFPT (log(MFPT)), as well as transition action and MFPT, choosing r = 0.3 (Fig Fc ) and r = 0.4 (Fig Fd).Similar linear correlations were found, again validating the robustness to the selection of hyperparameter r and the consistency of these three metrics.

E. The plant-pollinator network topologies can generate complex multistability
pollinators go extinct and the surviving species are typically generalists that have more connections with plants.Then we represent two additional scenarios for sub-networks, still obtained by removing pollinators from the orginal network (Hickling, Norfolk, UK, M PL 006) [5].
Network A (17 plants and 25 pollinators) is derived from removal of 60% pollinators (Table C), which also generates tristable behaviour (Fig Ga).A proportion of plants (P11, P13, P14, P15 and P17) are never pollinated by pollinators.The difference here is that more pollinators (40%) survive in the intermediate state.We also notice the abundance of some plants in the intermediate state is close to that in the high state, since higher pollinators richness and abundance promote plant diversity.This fact suggests that it may not be appropriate to project the landscape onto certain species.Not surprisingly, different patterns may emerge, including slightly higher abundance of some plants in the low state and the extinction of a few pollinators in the high state.These are due to the effect of intraspecific competition.
Network C (17 plants and 13 pollinators) is obtained from removal of 80% pollinators (Table E

F. RBH has better predictive efficiency than traditional metrics based on stochastic time series
We still focused on the variation of parameter κ, and first considered that it increases or decreases linearly with time.We were interested in the region where the phase transition occurs and simulated the trajectory of abundance change for each species during κ = 0.96 ∼ 1.2.We chose dt = 0.01, steps=2.4 × 10 5 , d = 5 × 10 −4 , κ = 0.96 + 10 −6 × i th step for collapse and κ = 1.2 − 10 −6 × i th step for recovery.Based on the parameter-varying Langevin equation, we simulated the abundance fluctuation of thirty species.Then it was projected onto the previously obtained coordinates (plant PC1 and pollinator PC1) and the occurrence of two phase transitions can be clearly seen (Figs 5A and B).
The early warning is concerned with the collapse process, also known as the transition to the low state.After we applied Gaussian filtering to each species trajectory (Figs Ia-b), we obtained the residuals' autocorrelation at lag 1 (AR(1)) and variance (Var), which are calculated in a time window (5% datapoints) to represent the value of a single point.The coefficient of variance (CV) and fano factor (fano) are both December 20, 2023 6/34 quantified from original species abundance, since we considered that it leads to a smaller denominator after filtration.These indicators show an upward trend before the system undergoes the phase change [6,7].Our procedures are shown and illustrated in terms of AR (1).For each species corresponding trajectory, the goal is to find the smallest κ 0 that satisfies AR(1)(κ 0 ) < AR(1)(κ) for all κ ∈ (κ 0 , κ * ), where κ * represents the tipping point of collapse.It is equivalent to consider slope AR(1) (κ) > 0 for all κ ∈ (κ 0 , κ * ) [7].Based on least-squares linear regression over 20% datapoints of AR(1), we estimated their slopes and found the minimum κ 0 according to the above criterion.We calculated the minimum κ 0 corresponding to all species.Then the procedures are repeated for other metrics (Var, CV, fano) (Fig 5D).
Based on these metrics, the early warning is possible for some species, especially for pollinators surviving in the intermediate state (pollinator 2).For three species (plant 3, plant 15 and pollinator 2), we showed the variation of AR( 1 Ie).However, we found that some species (like plant 6, pollinator 3 and pollinator 8) cannot provide accurate early-warning signals (EWSs) (Figs Id and If), i.e., no so-called κ 0 exists (grey points in Figs Id and If).A possible explanation is that most pollinators become extinct before the intermediate state.Their degrees are generally lower, but that's not always the case, like Network A and Network C. Their extinction does not mean the collapse of the system, and they do not exhibit obvious characteristics in terms of system collapse.The traditional perspective is that species with lower degrees are feasible to anticipate critical transitions because they go extinct earlier [7].We did find this pattern, but here the real breakdown is marked by the extinction of pollinators surviving in the intermediate state (generally with larger degrees), and their monitoring is more effective.
Derived from our reduced landscape, the RBH can be seen as a sequence with varying κ.We desired to find the critical point when the whole sequence up to this point exhibits nonlinearity.The Brock-Dechert-Scheinkman (BDS) test is applied to detect it [8].We eliminated linear trends using the first difference, and then performed the test on its residuals.They first need to accept the hypothesis that they are white noise, by the Ljung-Box Q-test (LBQ test).The collapse transition induced by strong nonlinear responses is expected to reject the BDS's null hypothesis that residuals are independent and identically distributed [6,8].After comparisons with traditional metrics obtained before, the RBH's critical point is much earlier (smaller κ).We also performed above analysis for Network A (Figs Gb-e Besides, when κ exhibits nonlinear changes, the predictive effectiveness of RBH remains significantly better than traditional indicators (Figs J and K).We also explored the situation when κ exhibits weak autocorrelated fluctuations (κ(t + 1) = ω(t + 1)+ 1. 1, where ω(t + 1) = 0.3 × ω(t) + ξ(t) and ω(0) = 0, ξ(t) ∼ N (0, 1/40)), the BDS test does not reject the null hypothesis, indicating that the EWS can to some extent mitigate the issue of false positives (Fig L ).

G. Global sensitivity analysis identifies the indispensable role of species surviving in the intermediate state
Here we further illustrate the importance of intermediate state species for mitigating collapse, by analyzing the relative change of transition action based on global sensitivity analysis.As we did in our main text, we selected a set of parameters in the domain of multistability, and made suitable perturbations, trying to ensure the presence of intermediate states.Then we calculated the relative change in the transition action, to quantitatively portray the evolution of the system.Despite the network discussed in the main text, the only surviving pollinator in the intermediate state do have the largest degree.But it is not always the case.For example, for Network A, ten pollinators survive in the intermediate state, in which the degree of I9 ( 1) is less than the degrees of extinct pollinators I2, I3 and I6 (2).For Network C, three pollinators survive in the intermediate state, in which the degree of L6 ( 3) is less than the degree of extinct pollinator L9 (4).
We perturbed parameters of these three types (chosen by degree, generalist and surviving in the intermediate state), and marked the superscript 'degree', 'generalist' and 'our' respectively.For comparison we again did the case of perturbing all parameters, where no superscripts are labelled.We applied suitable (small) perturbations (Network A: 1%, Network C: 0.5%), trying to avoid the phase change.We calculated the relative changes of transition actions of intermediate to Surviving species are associated with the complex network topology, and they cannot be directly inferred from the magnitude of the degree.The most generalist pollinator plays a certain but relatively limited role.However, for full recovery (intermediate to high transition, blue bar in Fig M ), the effectiveness of perturbation on above three types becomes smaller, especially for the individual generalist pollinator.As we discussed in our main text, recovery to the high state is a joint result of interspecific collaboration, and the role played by a few pollinators is restricted.December 20, 2023 8/34 For Network C, although four stable states can appear, we noticed that the intermediate II appears in a small region during our perturbation experiments.Also considering the consistency of studies in Network A, we focused on the transition actions from intermediate I to low state and intermediate I to high state, and tentatively referred to intermediate I as the intermediate state.We found that, in the case of increasing and decreasing 0.5% of γ 0 , the intermediate state disappears, which confirms the influence of the reciprocity parameter on the system once again (Fig N ).We noted that the intermediate state is fully converted to the low state at decreasing 0.5% of γ 0 , while to the high state at increasing 0.5% of γ 0 .For the former case, the transition action from intermediate to low state is 0 and its relative change is -1; the action from intermediate to high state is infinity and its relative change is infinity.For the latter case, the relative change is infinity for the intermediate to low transition, and is

H. Linear stability analysis demonstrates that intermediate states only appear in multidimensional frameworks
The dimension reduction approach is much challengeable since we always have to strike a better balance between simplified modelling and detailed descriptions.Prior research has noted how to predict tipping points in plant-pollinator mutualistic networks through reduced two-dimensional (2D) equations [9].Nevertheless, the main weakness of the paradigm is that we cannot identify stable intermediate states via all suitable parameters.To the best of our knowledge, this is the first time to analyze the linear stability of the 2D reduced models in detail.
We again showed the 2D equations with respect to P ′ and A ′ [9]: where definitions for most parameters (α, β, h, µ, κ) is the same as in multidimensional model.Additionally, ⟨γ P ⟩ = γ 0 , in which k i and d j are nodal degrees of the i th plant and j th pollinator respectively.To solve the stationary points, the right-hand sides of Eq.S2 are set equal to zero.Then we obtained: where the sign ± can be either positive or negative.Thus we had four possible combinations (++), (+−), (−+), (−−).We presented the simplified forms with µ ≈ 0: December 20, 2023 9/34 . However, this does not mean four pairs of points.(P * , A * ) (P * > 0, A * > 0) actually corresponds to two pairs of points, as A * is the solution of a quadratic equation (Eq. 3 in the main text).Therefore, we can analyze the linear stability of above four pairs by their corresponding Jacobian matrices.Since the intrinsic growth rate remains positive (α > 0), at least one eigenvalue of J (P1,A1) = α 0 0 α − κ has a positive real part.By Lyapunov's first method, (P 1 , A 1 ) is always unstable.As for (P 2 , A 2 ), the eigenvalues of and κ − α.They cannot be both negative, which results in the unstable (P 2 , A 2 ).However, (P 3 , A 3 ) can become stable on condition that κ > . Indeed, it corresponds to the low state we considered in our multidimensional model.Here the reduced model is in good agreement with the original multidimensional model (Fig 1D).All pollinators went extinct (A 3 = 0) and all plants showed similar patterns of low abundance which can be approximated by P 3 = α β .This also implies that increasing κ will inevitably make the low state stable, both for the reduced and multidimensional models.
Although we have analyzed (P * , A * ) in the main text, next we show more details.First, its Jacobian matrix is 1+h⟨r P ⟩A ′ is always negative, (P * , A * ) is stable if and only if the Jacobian determinant is positive.Then we derived the condition (Eq.5).By simulation of the effects of increasing κ from zero (in the main text), we have proved that (P 5 , A 5 ) corresponding to larger root of Eq. 3 is always stable, and A 5 has a lower bound (Eq.6): Here our perspective is that the critical point (lower bound) of the high state (P 5 , A 5 ) is an intrinsic property, which depends on the internal parameters and network, instead of external environmental parameter κ.We solved the critical value of κ which corresponds to the lower bound of A 5 : It is derived from ∆ = q 2 2 − 4q 1 q 3 = 0, in which the larger root December 20, 2023 10/34 is discarded because its corresponding abundance is negative, without any ecological meaning.It can also explain the appearance of bistability and hysteresis loops.When κ approaches its critical point κ 1 , pollinator spcies abundance is close to its lower bound.Continuous increase in κ leads to the disappearance of the high state, and before that, the low state has met the conditions for reaching stability.It should be pointed out that the exact abundance in the high stable state is influenced by κ.We currently denote the critical point of κ by internal parameters, with the aim of better explanation and simpler forms.Other representations may result in lengthy calculations, but do not affect the results.
To facilitate analysis, we investigated the effect of controlling individual parameter.We have already discussed the case of κ in the main text.As for other parameters, we first took full advantage of κ 1 's form, which can be regarded as a function with respect to α, β, h and γ 0 .The emphasis here is the sign of the partial derivative has not changed in the domain we cared about.Here ∂ ∂α κ 1 and ∂ ∂γ0 κ 1 are positive, but ∂ ∂β κ 1 and ∂ ∂h κ 1 are negative.We added the weak assumption (β < ⟨γ A ⟩ < ⟨γ P ⟩) for convenience in the proof, which are always satisfied in our cases.The conclusions can be reached through arithmetic-geometric mean inequality, and neglecting small terms O(h 2 ) since the default value of h is 0.2.We calculated that This implies to us that changing parameters in one direction causes monotonic effects in κ 1 .Worse conditions (smaller α, γ 0 or larger β, h) make the κ 1 smaller, and the high stable state disappers when κ 1 is smaller than current κ.
The monotonic impacts of increasing parameters on other variables (A 5 lower bound, the axis of symmetry, A 4 * A 5 and A 5 ) can be found similarly (Table B), and we omit the lengthy details.Here o indicates that increasing value of the variable has no effect on the first-row variable, + represents increased parameters have positive effects on either increasing value or moving right, − represents negative effects and / denotes the effect is not verified.An unexpected finding is the opposite effect of α and γ 0 on the lower bound, where it decreases when α increases, but increases when γ 0 increases.December 20, 2023 11/34 These continuous and monotonic changes ensure that A 5 is the only stable point over a large parameter range.Numerical simulations also show the expected and consistent results (Fig O).We considered cases of varying individual parameter and fixed other parameters that do not affect the shape of the function.Varying multiple parameters is also a matter of changing individual parameters in sequence.Overall, worse conditions cause A 5 to approach its lower bound, also the dividing point between stability and instability, but A 4 never goes beyond it.So far, we have reported a basic understanding of the reduced 2D system.It has at most two stable states, (P 3 , A 3 ) and (P 5 , A 5 ), over an appropriate and widespread range of parameters.

I. A new network to replicate our findings
We reproduce our conclusions in a new network in the Web-of-Life database [10]

p
ss (P 1 , . . ., P NP , A 1 , . . ., A N A ) dA 1 • • • dA N A j ss (P 1 , . . ., P N P , A 1 , . . ., A N A ) dA 1 • • • dA N A 1 , . . ., P N P , A 1 , . . ., A N A ) dA 1 • • • dA N A = M j=1ϕ j p j ss (P 1 , . . ., P N P ) ), pollinator PC1 (Fig Eb), PC1 (Fig Ec) and PC2 (Fig Ed).During the variation period, we also showed the distribution of variance contribution rate of plant PC1 in Σ(P ), pollinator PC1 in Σ(A), as well as PC1 and PC2 in Σ(P, A) (Fig Ee).Our two coordinates (plant PC1 and pollinator PC1) are much more robust, and their combination contains the information of PC1 (Figs Ea- ) (Fig Ic) and variance (Fig Ie) with varying κ.For each species and metric, we also marked the location of the minimum κ 0 (red points in Figs Ic and ) and Network C (Figs Hb-e), the hysteresis loops always appear (Figs Gb-c and Hb-c).In general, RBH serves as the earliest warning signal of collapse, except for choosing p < 0.001 in Network C (red line in Fig He).By the way, through collapse and recovery trajectory in Network A, we found that the collapse process goes through the intermediate state (Fig Gb), while recovery does not (Fig Gc).
low transition and intermediate to high transition (Network A: Fig M, Network C: Fig N).For Network A, we found that, in the case of increasing 1% of γ 0 and decreasing 1% of κ, the intermediate state disappears and is fully converted to the high state (Fig M).The transition action from intermediate to high state is 0 and its relative change is -1; the action from intermediate to low state is infinity and its relative change is infinity.Remarkably, we again emphasized that, for intermediate to low transition (gold bar in Fig M), perturbing species surviving in the intermediate state ('our' in Fig M) has almost the same effect as perturbing all species (no superscripts in Fig M).Good ecological management of these species can essentially mitigate the collapse.It is also more effective than perturbing the same number of species based on degree ('degree' in Fig M), as well as singly perturbing the most generalist species ('generalist' in Fig M).
-1 for the intermediate to high transition.The important role of species surviving in the intermediate state for collapse mitigation is discovered again (Fig N).The conclusion that controlling a few pollinators has a limited role in the transition to the high state is also consistent (Fig N).
(Fig Pa).Similarly, by removing a subset of pollinators, we identified a subnetwork that exhibits an intermediate state (Network B: Fig Pb, Table D), enabling us to describe the abundance of each species in stable states (Fig Pc).Next, we plotted the energy landscape (Fig Qa) in reduced coordinates (Fig Qb) and calculated the transition paths between stable states (Fig Qc), once again revealing a tendency to transition from the high state to the low state through the intermediate state.We have also demonstrated the landscape change with parameter variation (Fig R), the consistency between RBH and transition action (Fig S), and RBH as an EWS (Fig T).Sensitivity analysis once again highlights the role of pollinators surviving in the intermediate state in preventing system collapse (Fig U).Here, the pollinators surviving in the intermediate state happen to be the two largest-degree pollinators.

Fig A .
Fig A. In the case of fixed parameters, noise-induced transitions between stable states.(a-b) Taking κ = 1.12, with initial points in the low state (a) and high state (b), noise induces the system to transit to the intermediate state, presented by plant 1 and pollinator 1. (c-d) Taking κ = 1.02, with initial points in the low state (c) and intermediate state (d), noise induces the system to transit to the high state, still presented by plant 1 and pollinator 1. (e) With initial points in the low state without external force, the plants recover earlier than the pollinators, consistent with our calculated transition paths.(f) Under conditions of external intervention (keeping pollinator 2 at high abundance), the plants still recover earlier.
Fig E. The superiority of our approach is further illustrated by the variation of coordinate weights and variance contribution rates.(a -d) We vary different parameters and observe changes in the weights of each species that make up the coordinates: (a) plant PC1 and (b) pollinator PC1 in our approach; (c) PC1 and (d) PC2 in previous approach [3].We use the same perturbations of parameters to ensure that two approaches are comparable.(e) The distribution of variance contribution rate of plant PC1, pollinator PC1, PC1 and PC2 With parameters varying.

FigF.
Fig F. The consistency among the barrier height, transition action as well as MFPT, which does not depend on the selection of the hyperparameter r (left: intermediate and high states, middle: low and high states, right: low and intermediate states).(a) ∆S and RBH exhibit the same trend.(b) The line graphs of MFPT from one attractor to another with varying κ, and different radii r are selected as 0.3, 0.4 and 0.5, where the error bars correspond to stand error.(c-d) Approximate linear relationships between RBH and log(MFPT) (blue line), as well as between ∆S and log(MFPT) (green line) under different r (c: r = 0.3, d: r = 0.4).

Fig G.
Fig G.The analysis is repeated for Network A, reproducing the superiority of RBH as an EWS.(a) The abundance of 17 plant species and 25 pollinator species (indicated by each column, see Table C for details) in low, intermediate and high stable states using the default parameters.The default values are also set as α (P ) i = α (A) j = 0.3, β (P ) ii = β (A) jj = 1, β (P ) ik = β (A) jm = 0.01(i ̸ = k, j ̸ = m), κ = 1.07, γ 0 = 1, δ = 0.5, h = 0.2, µ (P ) = µ (A) = 0.001.(b-c) The process of system collapse (b) and recovery (c) is simulated from the Langevin equations, by modelling linear change in κ.The system goes through the intermediate state for the collapse, while directly transits to the high state for the recovery.(d) Calculated RBH between low and high states with increasing κ.Three marked points are with the same meaning as the main text.(e) For complete collapse, the RBH again serves as the earliest warning signal compared with AR(1), variance (Var), coefficient of variation (CV) and fano factor (fano).

Fig H .
Fig H.The analysis is repeated for Network C, reproducing the superiority of RBH as an EWS.(a) The abundance of 17 plant species and 13 pollinator species (indicated by each column, see Table E for details) in low, intermediate and high stable states using the same default parameters.(b-c) The process of system collapse (b) and recovery (c) is simulated from the Langevin equations.(d) Calculated RBH between low and high states with increasing κ.Three marked points are with the same meaning as the main text and Fig Gd.(e) For complete collapse, if we choose p < 0.05 (blue line) or p < 0.01 (gold line), the RBH can serve as the earliest warning signal compared with AR(1), variance (Var), coefficient of variation (CV) and fano factor (fano).
Fig I. Interpretation of the time-series metrics.(a) Species abundance fluctuations and shifts are shown under the parameter κ increasing.Gaussian filtering algorithms are used for preprocessing, presented by plant 1 and pollinator 1.(b) The residual is obtained from the difference between the original abundance and the Gaussian filtering results, presented by plant 1. Traditional time-series indicators rely on the residuals to determine the occurrence of critical points, which is identified for the continued increase of the indicators.(c-f) We examine whether all species can achieve early warning.Some species can, like plant 3, plant 15 and pollinator 2 (red points in c and e); while some species cannot, like plant 6, pollinator 3 and pollinator 8 (grey points in d and f).

4 𝜿
Fig K.The scenarios of exponential and logarithmic changes in κ.(a) Corresponding to 5 ○ in Fig Ja (the exponential changes in κ), (b-c) we showed the RBH variation (b) and identified its earliest warning (c), which implies that the rapid change contributes to achieving earlier warnings.(d) Corresponding to 6 ○ in Fig Ja (the logarithmic changes in κ), (e-f) we showed the RBH variation (e) and identified that RBH's warning occurs later than the exponential scenario, but it is still superior to the traditional indicators on average (f).
Fig L. To show the RBH's inability to predict transitions that have not occurred, we perturbed κ in an autocorrelated manner.(a) We simulate κ fluctuating over time with weak autocorrelation.(b) The corresponding RBH varies with changes in κ, and the color bar represents the time progression.(c) The RBH evolves over time, and there is no statistically significant point in the sequence to reject the BDS test.
Fig M. Global sensitivity analysis is repeated for Network A, to emphasize the indispensable role of pollinators surviving in the intermediate state.We perturbed 1 % of each parameter from the default value to guarantee the existence of the intermediate state.The left panel corresponds to the perturbations in the direction of collapse, and the right corresponds to the opposite perturbations.We reported the relative changes in transition actions of intermediate to low transition and intermediate to high transition.No superscript (α (A) , β (A) , γ 0 , κ) means perturbation on all species parameters.The superscript 'our' means only perturbation on the parameters of the pollinator species surviving in the intermediate state (I1, I4, I5, I8, I9, I11, I12, I14, I16, I17).The superscript 'degree' means only perturbation on the parameters of the same number of pollinator species with larger degrees (I1, I2, I3, I4, I5, I6, I8, I11, I12, I16).The superscript 'generalist' means only perturbation on the parameters of the most generalist pollinator species (I1).Increasing γ 0 or decreasing κ makes the relative change be -1 and infinity, and for visualization, we denoted them as -0.75 and 0.75 respectively.

FigFigjj = 1 ,FigFig S .
Fig N. Global sensitivity analysis is repeated for Network C, to emphasize the indispensable role of pollinators surviving in the intermediate state.We perturbed 0.5 % of each parameter from the default value, and the presentation follows Fig M. The superscript 'our' means only perturbation on the parameters of the pollinator species surviving in the intermediate state (L1, L2, L6).The superscript 'degree' means only perturbation on the parameters of the same number of pollinator species with larger degrees (L1, L2, L9).The superscript 'generalist' means only perturbation on the parameters of the most generalist pollinator species (L1).Increasing or decreasing γ 0 makes the relative change be -1 and infinity, and for visualization, we denoted them as -0.5 and 0.5 respectively.
), taken as an example of the emergence of two intermediate states (intermediate I & II, I 1 &I 2 ) (Fig Ha).Both of them have characteristics mentioned previously, and the one closer to the high state (intermediate II, I 2 ) is featured by more abundant and more diverse pollinator communities (Fig Ha).In fact, only one species (L6) differs significantly, which survives in I 2 but goes extinct in I 1 .

Table B .
Monotonic by parameters changed in one direction κ 1 A 5 lower bound the axis of symmetry A 4 * A 5 A 5

Table C .
Network A: more pollinators surviving in the intermediate state