Foreseeing the future of mutualistic communities beyond collapse

Abstract Changing conditions may lead to sudden shifts in the state of ecosystems when critical thresholds are passed. Some well‐studied drivers of such transitions lead to predictable outcomes such as a turbid lake or a degraded landscape. Many ecosystems are, however, complex systems of many interacting species. While detecting upcoming transitions in such systems is challenging, predicting what comes after a critical transition is terra incognita altogether. The problem is that complex ecosystems may shift to many different, alternative states. Whether an impending transition has minor, positive or catastrophic effects is thus unclear. Some systems may, however, behave more predictably than others. The dynamics of mutualistic communities can be expected to be relatively simple, because delayed negative feedbacks leading to oscillatory or other complex dynamics are weak. Here, we address the question of whether this relative simplicity allows us to foresee a community's future state. As a case study, we use a model of a bipartite mutualistic network and show that a network's post‐transition state is indicated by the way in which a system recovers from minor disturbances. Similar results obtained with a unipartite model of facilitation suggest that our results are of relevance to a wide range of mutualistic systems.


S2 EXAMPLE: CRITICAL SLOWING DOWN IN A 4-SPECIES NETWORK
network's stability landscape. Attraction basins are shallow in between alternative stable states and the saddle 48 points on the thresholds that separate them. When approached by a threshold, the attraction basin of the initial 49 pristine state becomes increasingly shallow and the network increasingly slow when recovering from perturbations 50 in the direction of the saddle point on the approaching threshold. 51 For the here studied 4-species network ( Fig. 1) we found that the network's pristine state is initially accompa-52 Methods: To determine the rate at which pollinator abundances change as illustrated in Fig. 1 or ridges in the landscape. The stability landscape produced with this algorithm, is a useful tool to intuitively 79 illustrate the idea behind our method. As our system is non-gradient, it is not a way to determine the potential 80 energy of the system. To evaluate the performance of the first principal component, we determine the difference between the slope of 96 our indicator and the direction of the observed shift in abundance. We do this by determining the angle, θ, between 97 the direction of the indicator and the observed shift as follows: in which I is the indicator of a network's future state and ∆N (A) the observed shift in pollinator abundances.

99
I · ∆N (A) indicates that we take the dot product between these two vectors. To determine ∆N (A) , we take 100 the mean abundances over 200 time steps at 500 steps before the tipping point and subtract it from the mean 101 abundances 500 steps after the tipping point was found. Because we want to evaluate the accuracy of the first 102 principal component, and not whether points are also skewed in the right direction, we take −I as the input for 103 the formula above when we find an angle > π/2 (i.e. > 90 degrees). Both I and ∆N (A) are vectors of which To determine the aforementioned probability, we use the following probability density function: in which S (A) is the number of dimensions and h(θ) the probability density for a certain angle θ (ref. Cai et al. 115 (2013)). Our method may be interpreted as a test whether the null hypothesis that I and N are two random vectors 116 is true. This hypothesis is rejected when angle is found to be significantly smaller than the expected angle between 117 two random vectors, when the one-sided p-value is smaller than 0.01 (i.e. similarity > 0.99).

118
To evaluate the tendency of time points to be skewed in the direction of a network's future state, we determine 119 the skewness of the time points projected on the first principal component. When points are skewed in the direction 120 of the network's future state, we report a positive skewness. When points are skewed in the opposite direction, we 121 report a negative skewness. We consider a positive skewness as accurate and a negative skewness as inaccurate.

122
A strong positive or negative skewness is considered more accurate or inaccurate than a weak positive or negative 123 skewness. has changed. We, therefore, correct previously found indicator values such that there is no change larger than 90 138 degrees between two consecutive points at which the indicator's direction was determined. We assume the last 139 direction in which time points were found to be skewed to be the accurate one.

140
To determine whether there is a significant increase in the indicator's magnitude, we determine the Kendall rank 141 correlation coefficient, τ , for the last ten points at which the indicator's magnitude was computed. We consider 142 the increase significant when this coefficient was positive and its p-value < 0.05. Once a significant increase was 143 found, we tested whether the increase remained significant by determining Kendall's correlation for the last eleven 144 points the next time the indicator's magnitude is determined, for twelve points the time after that, and so on until 145 the tipping point is reached. We would again look at the last ten points when the increase was found to not be 146 significant anymore. By doing this, we could determine the range in conditions in which the indicator's magnitude 147 increased significantly.

148
As a measure of a 'regime shift' we determined whether there was a change in abundance of more than 1.5 149 over a period of 1% of the entire time series (200 time steps). We did this by taking the mean abundances over a 150 period of 200 time steps before this period and 200 time steps after this period and determining Euclidean distance 151 between these two mean abundances. To make sure that this large shift in abundances was not a temporal large 152 deviation from the species' mean abundances, we added as a second criterion that the abundance of at least one 153 species should be near extinction, i.e. below 0.1. 154 We did not apply any preprocessing to handle trends in the time series. We expect the indicator to be relatively 155 robust against such trends, because trends only alter the direction of the first principal component when their effect 156 on this direction is stronger than the effect of critical slowing down. Not applying any preprocessing is a good way 157 to test this robustness. When using the indicator as part of a different study it may, however, be worth considering 158 to apply a preprocessing method (see ref. Dakos et al. (2012)). It may improve the performance of the indicator, 159 especially when trends are strong.

S5 ADDITIONAL INFORMATION BIPARTITE MUTUALISTIC NETWORKS
Nontrivial equilibrium abundances,N , competitive interaction strengths, c, mortality rates, d, and saturation terms, 162 h, are randomly sampled from predefined probability distributions, and the total amount of resources received by 163 species i at the system's nontrivial equilibrium, R i (N (P ) ), are assigned such that the rate at which abundances 164 change at the system's nontrivial equilibrium, dN (P ) /dt, is zero: .
The total amount of resources provided at the system's nontrivial equilibrium, R i (N (P ) ), is thus approximately 166 the same for highly specialized and more generalist species, provided that their losses due to competition, c, and 167 mortality rates, d, and their nontrivial equilibrium abundances,N , are similar.

168
The extent to which species are saturated is determined by the total amount of resources provided, R i (N (P ) ), 169 and the rate at which species become saturated as determined by saturation term h i . In our simulations, we assume 170 nontrivial equilibrium abundances,N , and inter-and intraspecific competition, c ij and c ii , to be similar for all 171 species. Highly saturated species are, therefore, the ones with a high h i . Species are saturated relatively quickly, 172 and, according to equation S4, the total amount of resources provided at the system's nontrivial equilibrium is high 173 when species have a high h i .

174
Parameters are assigned such that there are substantial differences in the extend in which species are saturated 175 by drawing saturation terms, h i , from a scaled beta distribution with range ∼ (0.05, 0.35) and shape parameters 176 α = 1 and β = 5. Due to this distribution, there are few highly saturated species, i.e. h i close to 0.35, and 177 many non-saturated species, i.e. h i close to 0.05. Strong mutualistic interactions between non-saturated species 178 lead to strong positive feedbacks. Non-saturated species thus need to obtain a relatively large share of resources 179 from a few, highly saturated species for the network to be stable. Relative mutualistic benefits at initial conditions, 180 θ 0,ik , are therefore ordered such that larger benefits are obtained from the more saturated species. To make sure 181 that the sum of all relative benefits is one, we take relative mutualistic benefits, θ 0,ik , from a symmetric Dirichlet 182 distribution. The distribution's concentration parameter, α, determines the extent in which species are specialized 183 and is, for each species, taken from a uniform distribution between zero and one.

189
Delayed negative feedbacks become stronger as the strength and variability of interspecific competition increases.

200
As conditions change, either a single eigenvalue or a pair of complex conjugate eigenvalues goes to zero. In the 201 first case we are dealing with a saddle-point approaching the network's initial state, caused by a positive feedback.

202
In the second case, we are dealing with a Hopf bifurcation caused by a delayed negative feedback.

203
Data sets consist of 100 initial networks. For each network, 10 final distributions of relative mutualistic ben-204 efits, θ f inal,ik , were drawn, allowing us to determine the extent in which a community's future state depends on 205 the specific way in which relative mutualistic benefits are changed. Parameters were assigned such that this depen-206 dency is high. Networks were discarded from a data set when they were unstable at initial conditions, M = 0. We 207 determined the frequency at which this occurred as a measure of how difficult it is to find a stable solution for the 208 initial networks of a given data set. The final distribution of relative mutualistic benefits was redrawn either when 209 the network would become unstable within the range of conditions M = (0, 0.5), or when a network would still 210 be stable at M = 1.

211
To test whether the indicator also works when equilibrium abundances change, we analyzed networks of 10 212 plant and 10 pollinator species of which the final equilibrium abundances are different. We do this by changing the 213 nontrivial equilibrium abundances of species as follows: in whichN 0,i is the initial,N f inal,i the final, andN * i the actual nontrivial equilibrium abundance of species i.

215
The total amount of resources provided at the system's nontrivial equilibrium, and the strengths of mutualistic 216 interactions are determined by equations 4 and S4. We tested three scenarios. One in which the nontrivial equi-217 librium abundances of species tend to increase,N f inal,i ∼ U (2, 3), one in which they stay the same on average 218N f inal,i ∼ U (1.5, 2.5), and one in which they tend to decreaseN f inal,i ∼ U (1, 2). Competitive interaction 219 strengths were taken from the following distributions: c ii ∼ U (0.9, 1.1) and c ij ∼ U (0.02, 0.08). Changing 220 abundances affect all relationships as described by the Jacobian matrix. The main effect of a decline in abun-221 dance is, however, a reduction of the direct negative effects of species on themselves which undermines resilience.

222
Increasing abundances tend to promote resilience.

223
To test whether the indicator may accurately indicate the future state of larger networks, we analyzed networks species. We assigned competitive interaction strengths such that the rate at which species lose in abundance 226 due to competition, , is approximately the same for different numbers of species, as well 227 as the relative difference between intra-and interpecific competition, c ij /c ii . When a species group consisted 228 of 10 species we assumed c ii ∼ U (0.9, 1.1) and c ij ∼ U (0.02, 0.08). When a group consisted of 20 species one. To make sure that the rate at which abundances change at the nontrivial equilibrium, dN i /dt, is zero, we 242 assign the intrinsic growth rates, r, as follows: The contribution of species to the overall resilience of a network is determined by critical abundance A i .

244
Species with a high critical abundance, A i , collapse more easily and the overall resilience of the community is 245 highest when such species are facilitated by species with a low critical abundance. A change from such a distri-246 bution to a more random distribution of facilitative interaction strengths will undermine resilience. To generate 247 time series in which the resilience of the here described facilitative communities is undermined, we assume that 248 conditions, M , affect facilitative interactions as follows: in which γ 0,ik is the initial, γ f inal,ik the final, and γ * ij the actual facilitative interaction strength. Conditions, 250 M , change from zero to one over time. We assume that the total amount of facilitation received, S j=1 γ ij N j, 251 remains equal as conditions change. We therefore determine the final facilitative interaction strength as follows: in which θ ij is the fraction of the total facilitation received by species i from species j. 253 We assign parameters such that there are substantial differences in the critical abundances of species by draw-254 ing critical abundances, A i , from a scaled beta distribution with α = 5 and β = 1 and range ∼ (0, 1.5). Due to 255 the beta distribution, there are few highly vigorous species (i.e. A i close to 0) and many non-vigorous species (i.e.

256
A i close to 1.5). The initial facilitative interaction strengths are taken from the following uniform distribution: Initial facilitative interaction strengths are ordered such that species receive most facilitation, 258 i.e. highest γ 0,ij , from species with the lowest A i . We assume that as conditions change, the strength of some 259 facilitative interactions increases strongly while others approach zero. Final relative facilitative benefits, θ f inal,ik , 260 are therefore selected with a probability of 0.75 and set to zero. To the remaining interactions, relative benefits are 261 assigned by taking them from a uniform Dirichlet distribution (α = 1). As with the model of mutualistically inter-262 acting species, we chose for this distribution of critical abundances, A i , and facilitative interaction strengths γ ij , 263 because it leads to a high variety in potential future states to which a network may shift. Other parameters and equi-264 librium abundances are taken from the following uniform distributions: Simulations were made with networks of 10, 20 and 40 species. As for the bipartite model of mutualisti-267 cally interacting species, we assign parameters such that the rate at which abundance is lost due to competition, 268 S j=1 c ij N j /K i , remains approximately the same for different species numbers, as well as the relative difference 269 between intra-and interpecific competition (see main text). Carrying capacities, K i , were therefore taken from 270 respectively K i ∼ U (5, 6), K i ∼ U (7.63, 9.15), and K i ∼ U (12.89, 15.47), depending on the number of species.

271
The amount of noise, determined by standard deviation δ, is assumed to be equal for all species. For the results

272
shown in this document we assume standard deviation δ = 0.05. As with the model of mutualistically interacting 273 species time series had a length, T , of 20.000 time steps.

S7 SUPPLEMENTARY RESULTS
Independent of the parameter ranges chosen, we found that regime shifts were preceded by a substantial period 276 in which the indicator's magnitude increases significantly, i.e. the 'critical range'. Our indicator would, provided 277 that the future state is indicated accurately, point towards a network's future state during a substantial part of this 278 period (Fig. S7). In Fig. S8 we provide information about the critical ranges as observed in a single data set interspecific competitive interaction strengths, c ij , were taken from ∼ U (0.02, 0.08), we found that such likely 287 cascading, full network collapses took up a bit more than 12% of the data set. For specific parameter ranges not 288 tested by us, this percentage may be higher.

289
In Fig. S10 we provide examples of two cascading collapses and one immediate network collapse. Species that 290 collapsed a bit later, were also the ones for which the indicated loss in abundance was smallest, suggesting that 291 the indicator indicates the initial regime shift accurately. The amount of time in between two consecutive partial 292 network collapses can be extremely small. Also when cascades are not clearly visible, we suspect therefore that 293 the inaccurate prediction of a full network collapse is caused by the occurrence of a cascading collapse.

294
In Fig. S13 we provide an example of a network for which the future state is hard to predict because it may 295 shift to several alternative future states. When making five simulations in which relative mutualistic benefits, θ ik , the several future states to which this system may shift is the fact that this system is approaching a Hopf bifurcation, 301 leading to oscillating (Fig. S14), chaotic or other complex dynamics (Fig. S15). Such dynamics may explain a 302 high sensitivity to perturbations in more than one direction.

303
In Fig. S18-S20, we show examples of time series in which not only the relative benefits, θ ij , change over time.

304
The nontrivial equilibrium abundances,N i , and thus the total gain from mutualistic interactions, R i (N i ), changes 305 as well. We found that a change in abundance over time does not have a strong effect on the performance of the 306 indicator (Fig. S17). In comparison to data sets in which abundances stay (on average) the same, full network 307 collapses are much less frequent when abundances increase and much more frequent when abundances decrease.

308
Quite a large fraction of full network collapses is indicated accurately when abundances decrease. Cascading 309 collapses may occur less frequently because all species experience a similar loss in resilience as a consequence of 310 a decline in abundance. Another difference is that the length of the critical range tends to be a bit shorter when 311 abundances in-or decrease.

312
In Fig. S21 and S22, we show that the indicator performs well, also when we apply our method to net- , increases substantially as the number of species increases. Systems with 317 many species may, therefore, be comparable with smaller networks in which interspecific competition is relatively 318 strong. In those networks we also observed that full network collapses were less frequent. Increasing numbers of 319 species did not have clear effect on the length of the critical range, nor on the fraction of the critical range in which 320 the future state was indicated accurately by the slope of the indicator (Fig. S22). We did, however, found some    . Far from the tipping point, in window I and II, deviations from the species' mean abundances are relatively small. Close to the tipping point, in window III, the distribution of points in the network's phase space is highly asymmetrical. Deviations from the mean abundances in time window III usually involve a simultaneous increase in the abundance of species A1 and a relatively larger decrease in the abundance of species A2, suggesting that this will also be the direction in which the network will shift once a threshold is passed. Figure S4 The measures of asymmetry together forming our indicator as they were determined for window III in Fig. S3. In this example, we found a large negative score (-0.79) indicating a relatively large decline in abundance for the pollinator on the x-axis and a relatively smaller positive score (0.32) indicating a relatively smaller increase in abundance for the pollinator on the y-axis. The length of the indicator corresponds to the amount of variance explained by the first principal component.

Figure S9
The number of pollinator species collapsing to extinction as observed in data sets of 1000 regime shifts. Each panel shows results when sampling competitive interaction strengths from a different parameter range (see ranges indicated).
In the extreme case where there was no competition (top left panel), we found almost exclusively full network collapses (i.e. all ten pollinator species collapsed to extinction). As the strength of competition increases, full network collapses become less frequent. Partial network collapses tend to be small independent of the strength of competition, i.e. the most common partial collapse led to the extinction of only one single pollinator species. The fraction of regime shifts for which the change in abundance was not well indicated is shown in red. The fraction accurately indicated by the first principal component, i.e. the slope of the indicator is accurate, but not by the direction in which time points are skewed is shown in light blue. Fully accurate predictions are indicated in dark blue.

Figure S10
Two cascading collapses and one immediate collapse. (a) Example of a cascading collapse that eventually leads to the collapse of four pollinator species. Three species (blue, green and purple) collapse to extinction rapidly. A fourth (black) species collapses as well, but remains for a short while at a lower abundance before collapsing to extinction (red arrow, a.I).
Out of the four species that collapse to extinction, the black species is also the one for which the indicated loss in abundance is smallest (red circle, a.II). (b) Example of a cascading collapse that eventually leads to a full collapse of the network (i.e. the most common outcome of a cascading collapse). Two species (black and yellow) collapse to extinction rapidly. The other species collapse as well, but remain for a short while at a lower abundance before collapsing to extinction (red arrow, b.I).
The indicated loss in abundance of the rapidly collapsing species is much bigger than the loss indicated for the species that collapse a bit later (red circles, b.II). (c) Example of a full network collapse that was accurately indicated. All species collapse at approximately the same time (c.I). All species were indicated to lose in abundance (c.II).