A potential role for rare species in ecosystem dynamics

The ecological importance of common species for many ecosystem processes and functions is unquestionably due to their high abundance. Yet, the importance of rare species is much less understood. Here we take a theoretical approach, exposing dynamical models of ecological networks to small perturbations, to explore the dynamical importance of rare and common species. We find that both species types contribute to the recovery of communities following generic perturbations (i.e. perturbations affecting all species). Yet, when perturbations are selective (i.e. affects only one species), perturbations to rare species have the most pronounced effect on community stability. We show that this is due to the strong indirect effects induced by perturbations to rare species. Because indirect effects typically set in at longer timescales, our results indicate that the importance of rare species may be easily overlooked and thus underrated. Hence, our study provides a potential ecological motive for the management and protection of rare species.

The ecological importance of common species is relatively well documented [1][2][3] . Yet, the ecological and functional role of rare species, which make up the overwhelming majority of species in the world's ecosystems 4 , is less well known 5 . There are, however, studies clearly demonstrating the keystone role played by some rare species [6][7][8] and their significant contribution to the functional diversity of ecosystems [9][10][11][12] . Yet, the relative importance of rare versus common species for the stability of ecosystems in general is largely an open question.
Theoretically, the relative ecological role of rare and common species has received little attention. So far, however, in dynamical models of ecological networks, i.e. networks describing who interacts with whom in an ecosystem, enforced extinctions of common species has been shown to lead to larger community wide effects than extinctions of rare species 13,14 . It has also been found that functional extinctions -where an increased mortality rate and subsequently reduced abundance of a given species first leads to an extinction of another species in the community rather than the species itself -are more likely to occur if the perturbation affects a 'high biomass' (i.e. common) species compared to a 'low biomass' (i.e. rare) species 15 . Similarly, harvesting a species at maximum sustainable yield has been shown to lead to larger community wide effects if harvesting is applied to common compared to rare species 16 . All these studies indicate that when large and sustained selective perturbations affect common species it will lead to larger community wide effects than if such perturbations affect rare species. However, the methods used in these studies are not directly suitable for evaluating the relative dynamical importance of rare and common species, since these type of perturbations may be considered large from a community point of view when applied to some species, and small when applied to others. For example, selective removal of a high biomass species may be considered a much larger perturbation to a community than removal of a low biomass species. Thus, thoroughly comparing the relative dynamical importance of rare and common species in an unbiased way requires a more sophisticated approach. To this end, an alternative is to investigate how small perturbations, assumed to be of equal magnitude across species, affects the stability properties of a community. Analytical methods for addressing the response to such perturbations have been derived [17][18][19][20] , yet the relative effect of perturbing the density of rare and common species has not been thoroughly investigated [but see 20 ].
Here we use a theoretical approach, combining analytical and numerical methods, to investigate the relative importance of rare and common species for the stability properties of ecological network models. First, we investigate the relative role of rare and common species in governing the rate of return to equilibrium following generic www.nature.com/scientificreports www.nature.com/scientificreports/ pulse perturbations; that is, temporary perturbations simultaneously affecting all species in a community. Second, we investigate the recovery of ecological networks exposed to selective pulse perturbations (i.e. temporary perturbations that affect one species at a time) to rare and common species, respectively. Third and finally, we investigate community wide responses of small selective press perturbations, that is, permanent perturbations, to rare and common species. Overall, our analytical and numerical modeling results suggest that rare species may be of particular importance for the stability of ecological communities as communities recover more slowly following pulse perturbations to rare than common species, and because stronger indirect effects are induced if rare rather than common species are permanently perturbed.

Results
the relative role of rare and common species in governing the recovery of ecological networks following generic pulse perturbations. Consider the dynamics of a system of n interacting species: where dN i /dt is the rate of change in density of species i, N i is the density of species i, r i is the intrinsic growth rate of species i, α ii is the per capita interaction strengths within species and α ij is the per capita interaction strength between species j and i (i.e. the effect of species j on species i). Following a linearization around an equilibrium point of the non-linear system (Eq. 1) the dynamics of small pulse perturbations (i.e. temporary perturbations) can be studied (Methods). If such perturbations eventually decay, the equilibrium point is called locally stable. Here we investigate the recovery process towards such locally stable equilibria (See Supplementary Fig. 1 for a conceptual illustration of this approach), using both analytical and numerical methods, and investigate the relative dynamical importance of rare and common species during the recovery process.
If it is assumed that a generic pulse perturbation, on average, affects each species in a community to an equal extent, the median initial recovery rate (i.e. just after the perturbation has been imposed [at t = 0 + ]) is given by (See Methods & Supplementary Fig. 1): here R Init is the median recovery rate of the system at time t = 0 + , α ii is intra-specific competition and N i is the equilibrium density of species i. If the per capita intra-specific competition among species in a community is assumed to be approximately equal, this equation illustrates that the median initial return rate is primarily governed by common species. Yet, if there is a strong negative correlation between the density of species and their intra-specific competition coefficients this result does not hold.
Furthermore, the initial return rate to equilibrium following a pulse perturbation is only an initial snap-shot (i.e. at time t = 0 + ) of the whole recovery process (See Supplementary Fig. 1). Other species may therefore be dynamically important during the rest of the recovery process. As a pulse perturbation decays it will eventually (asymptotically) decay at a rate given by the dominating eigenvalue, λ, of the Jacobian matrix evaluated at an equilibrium of the system in question (See Supplementary Fig. 1). This return rate -one measure of a system's resilience -has frequently been used as a measure of ecological stability [20][21][22] . For systems close to instability, that is, systems that are particularly vulnerable to perturbations, λ will be close to zero and resilience, Λ = −Re(λ), can be approximated by (See Supplementary Information Section 1 for derivation): Here D is the determinant of the Lotka-Volterra interaction matrix A with elements α ij (Eq. 1), and D -i is the determinant of the reduced matrix A −i with the ith row and column removed from A. Equation (3) shows that resilience is a weighted harmonic mean of species equilibrium densities, and therefore is primarily governed by the densities of rare species. Now, if species 1 is very rare compared to the other species, resilience can be further approximated by Λ ≈ − − DN D / 1 1 (See Supplementary Information Section 1 for details). This shows that resilience is most sensitive to changes in the equilibrium density of the rarest species in a community. Yet, if the dominating eigenvalue is complex, resilience is instead most sensitive to changes in the equilibrium density of the second rarest species in the community (See Supplementary Information Section 1). It is further worth noting that the dominating eigenvalue is not only important for species dynamics at longer timescales. In relation to other eigenvalues it also has the largest overall influence on species densities throughout the recovery process, and its relative importance increases over time (See Supplementary Information Section 2).
We numerically investigated our analytical results, relaxing the assumption of closeness to an instability, using models of seven real food webs (with antagonistic trophic links [Methods & Supplementary Information Section 3]) and a broad range of synthetic bipartite networks with a mixture of mutualistic, antagonistic and competitive interactions (Methods & Supplementary Information Section 4). Across replicates of the food webs we found the median initial return rate (i.e. just after a perturbation has been imposed) often to be related to the equilibrium biomass of the most common species (R init in Fig. 1). However, median initial return rates also depend on the intra-specific competition being experienced by species in a community (Eq. 2). In the bipartite networks, where intra-specific competition coefficients are different for each species in the network (Methods), there is no clear association between species density and median initial return rate (R Init in Supplementary Fig. 2). In fact, across www.nature.com/scientificreports www.nature.com/scientificreports/ replicates of these networks the abundance of the most common species is not significantly correlated to median initial return rate (rightmost bar R Init in Supplementary Fig. 2). Thus, our numerical analyses illustrate that species density in a community cannot be used as a general predictor of initial return rate to equilibrium. Yet, at longer timescales there is a somewhat less ambiguous relationship between species density and return rate across the investigated systems. For recovery rates at intermediate timescales, i.e. where half of the perturbation has recovered (See Supplementary Fig. 1 & Methods), rare species densities are more strongly associated with return rates than common species (R Interm in Fig. 1 & Supplementary Fig. 2), and at even longer timescales the significance of rare species is even more pronounced (R Asympt in Fig. 1 & Supplementary Fig. 2). show the Spearman rank correlation (Rho) for the association between return rates (R Init -median return rate immediately after the pulse perturbation is being imposed; R Interm -median return rate for the case when half of the perturbation has recovered; R Asympt -asymptotic return rate [Resilience], i.e. return rate as t -> ∞) and equilibrium biomass of a species with a given biomass rank, across 100 food web replicates. Bars are from rarest (left) to most common (right) species, with rank being based on equilibrium biomass within a food web replicate. Red dashed horizontal lines show the two tailed threshold level of the correlation coefficient at a significance level of 5%. The correlation between the equilibrium biomass of a species, with a given biomass rank within replicates, and return rates are thus significantly positively or negatively correlated if bars cross any of the two threshold levels (i.e. the red horizontal lines), meaning that these species affect return rates. See Supplementary Fig. 1 for a graphical illustration of the different return rates being used.
www.nature.com/scientificreports www.nature.com/scientificreports/ Recovery dynamics following selective pulse perturbations. In order to disentangle the relative effect of perturbations to common and rare species, respectively, we conducted selective pulse perturbations. Initially, i.e. just after a selective perturbations has been imposed (at time t = 0 + ), the rate of recovery is given by (See Supplementary Information Section 5 for derivation): t ii i 0 where α ii is the intra-specific competition coefficient for species i and N i its equilibrium density. Equation (4) shows that communities in general recover more slowly initially, if selective perturbations affect rare rather than common species. If the initial rate of recovery to equilibrium is considered a relevant measure of ecological stability (a system is more stable if it is approaching the equilibrium faster), communities should be more sensitive to perturbations of rare than common species. Yet, again this result is contingent on the intra-specific competition of species in the community. If there is a strong negative correlation between intra-specific coefficients and the equilibrium densities of species this result does not hold.
As already mentioned, the initial recovery rate following a pulse perturbation is only an initial snap-shot and may not be representative of the whole recovery process. We therefore numerically investigated the time it takes until a selective perturbation decays to half of its initial magnitude (i.e. a perturbations half-life; See Supplementary Fig. 1 & Methods). Across all systems studied there is a strong and negative association between species equilibrium density and perturbation half-life ( Fig. 2 & Supplementary Fig. 3). Indeed, within replicates of food webs as well as bipartite networks there is always a negative relationship between a species density rank and perturbation half-life (small barplots in Fig. 2 & Supplementary Fig. 3).
Selective press perturbations. To further investigate the relative dynamical importance of rare and common species we conducted small selective press perturbations; that is, small and permanent perturbations (See Methods & Supplementary Fig. 1). These analyses show that asymptotic return rate (i.e. resilience) is more sensitive to changes in the density of rare species, than to changes in the density of common species ( Fig. 3 and Supplementary Fig. 4; See also Supplementary Figs 6, 7). Indeed, within network replicates there is a strong and significant negative correlation between the density rank of a species (the rarest species having rank 1) and the effect of a change in its equilibrium density on resilience (small barplots in Fig. 3 and Supplementary Fig. 4; See also Supplementary Figs 6, 7). This result is probably due to the larger changes in species densities that arise from selective press perturbation to rare than common species (Fig. 4 & Supplementary Fig. 5). Indeed, within replicates of the bipartite networks and the seven food webs there is almost always a negative association between species density and the overall direct and indirect effects caused by press perturbation to a given species (small barplots in Fig. 4 & Supplementary Fig. 5). Analytical results further show that this numerical result is general for the class of models studied here, i.e. generalized Lotka-Volterra models (See Supplementary Information Section 6 for derivation).

Discussion
The dynamical importance of common species, which make up the core of ecological communities, is unquestionably due to their high abundances 1-3 . Yet, the relative dynamical importance of common vs. rare species is less well understood 2 . Theoretically, we here investigated this fundamental ecological question, and found that the importance of rare species may easily be overlooked and thus underrated. Small selective perturbations (assumed to be of equal magnitude across species) to rare species generally led to larger dynamical effects than similar perturbations to common species. This contrasts with previous studies showing that large perturbations, such as species extinctions, led to larger dynamical effects when common species, rather than rare species, were eradicated [13][14][15][16] . These apparently conflicting results illustrate that the specific type of perturbation being imposed will affect which specific species that are identified as major contributors to ecological stability 23 . Moreover, the identification of keystones in a community also depends on the specific stability metric being studied 24,25 . Here we used the rate of return to locally stable fixed points as a measure of ecological stability. This definition of stability assumes a static equilibrium point, while empirically based stability metrics are often based on some form of population or community variability 24,25 . Thus, in order to relate our theoretical results to empirical results based on the latter view, it needs to be explored if and how our theoretical stability metrics can be related to population variability. Although this objective warrants future research focus, it is worth noting that there is a direct relationship between the dominating eigenvalue of interaction matrices and population variability in stochastic linear discrete models 26 . This points to the possibility that our results, based on deterministic models, may also hold information for other empirically based stability measures. We have further assumed that species trophic interactions can be described by linear functional responses. The extent to which our results are generalizable to systems with non-linear functional responses therefore remains to be explored. However, it is first worth noting that several theoretical studies have compared model results based on systems with linear and non-linear functional responses, and found that many aspects (such as the probability of observing pyramidal patterns of abundance 27 , the risk of secondary extinctions 15,28,29 , and the probability of local stability of ecological networks 30 ) were qualitatively similar for the two scenarios. Secondly, whether linear or non-linear functional forms is a more realistic model of species interactions, in the range of species densities observed in multispecies systems, is an open question [31][32][33] . Still, to validate if our results holds also under non-linear functional responses, linearization using numerical methods is straight forward 34,35 , and the methods used here could thus easily be adopted to other non-linear systems in future work.
Another important aspect that should be acknowledged is the fact that equilibrium densities are functions of the interaction strengths in a model. One might thus argue that the findings presented here, that is, results www.nature.com/scientificreports www.nature.com/scientificreports/ illustrating the relative importance of rare vs. common species for the stability of ecological communities, may be confounded by a structured pattern of species interaction strengths (something that may emerge from body size scaling of parameters). Still, although we observe some weak relationships between species interaction strength www.nature.com/scientificreports www.nature.com/scientificreports/ and equilibrium biomass in the model realizations of the empirical networks ( Supplementary Figs 8-10), these relationships are of variable signs for the different networks (positive relationships for some systems and negative for others), indicating that it is not differences in species interaction strengths that are driving the patterns we www.nature.com/scientificreports www.nature.com/scientificreports/ see in the model food webs. Moreover, for both the bipartite networks, where there are no relationships between equilibrium densities and interaction strengths ( Supplementary Fig. 11), and the analytical results, which are to a large extent independent of the specific parametrization used, we found similar results with respect to the relative www.nature.com/scientificreports www.nature.com/scientificreports/ importance of rare vs. common species. Altogether, we thus think that the results presented here are fairly robust to the specific model parameterization being applied.
Our contribution partly builds on the elegant and informative tools developed by Arnoldi et al. 20 to study short-to long-term recovery of communities following pulse perturbations. Although the major contribution by Arnoldi et al. 20 was to develop tools for studying recovery dynamics of communities following generic pulse perturbations, they also began to explore the relative role of abundant and less abundant species for ecosystem recovery. To this end, they derived an analytical result showing that the rare species in a two species community (where a negative unidirectional effect of a common species on a rare species is assumed) governs asymptotic recovery rate, and found numerically that the rarest species, in a system with competitive interactions only, governs asymptotic recovery rate. In contrast, our study disentangles the relative role of rare and common species for ecological stability at short-to long-time scales, in arbitrarily complex systems (i.e. specious systems with any interaction types), following generic pulse perturbations, selective pulse perturbations as well as selective press perturbations. Using a broad set of analytical as well as numerical methods we have developed the theoretical understanding of the potential relative role of rare and common species for ecological stability by significantly extending the findings of Arnoldi et al. 20 , specifically showing that: (i) the long term recovery rate of arbitrarily complex generalized Lotka-Volterra model systems (i.e. with any combination of interaction types) that are close to an instability is governed by rare species and that the contribution of rare species to the recovery process can be significant also at shorter time scales; (ii) press perturbations to rare species lead to an overall larger effect on ecological stability than what perturbations to common species do; (iii) communities are more sensitive to selective pulse perturbations to rare than common species since communities recover more slowly when rare rather than common species are perturbed.
One might argue that our result, showing that communities recover slowly following selective pulse perturbations to rare species, only depends on the fact that rare species approaches equilibrium slowly. By extension, this might be interpreted as if perturbations to rare species only have negligible impact on community dynamics, and only at very long timescales as the rare species recover. Yet, in general we found this not to be true. In fact, (i) rare species have a strong effect on return rates already at intermediate time scales and (ii) stronger indirect effects (i.e. effects on other species mediated via chains of interactions) appear if rare species are perturbed compared to if common species are perturbed (See also Supplementary Information Section 7). Since indirect effects typically set in at intermediate to long timescales 36,37 , our results indicate that the effect of perturbations to rare species may be hidden in ecological dynamics (i.e. not immediately evident and manifested only when indirect effects unfold as time progresses), but the response of rare species perturbations may nevertheless be strong. Therefore, our study indicates that rare species may be just as ecologically important as common species, but the dynamical effects caused by perturbations to rare species may be hard to detect as well as difficult to predict (due to the indeterminacy of ecological interactions involving indirect effects 19 ). Thus, our results provides a potential ecological motive for the protection and management of rare species.

Methods
theory on recovery dynamics. The dynamics of a non-linear system, with time derivatives dN dt i , close to an equilibrium point, N , can be studied using a first order linear approximation: where C is the Jacobian matrix of the system evaluated at an equilibrium x is the displacement from the equilibrium, i.e. = −x N N i i i , and x(0) is an initial displacement from the equilibrium at time t = 0. The system of differential equations (Eq. 5) has the particular solution: and if the real part of the eigenvalues of C are all negative, the initial perturbation eventually decays; that is, e Ct → 0 as t → ∞. Such systems are called locally (or asymptotically) stable. Further, for a locally stable system, a small initial perturbation will eventually (asymptotically) decay at a rate determined by the dominating eigenvalue of C, i.e. the eigenvalue with largest real part. This rate of return is sometimes called the resilience of a system: Return rates following generic pulse perturbations. In a recent study, Arnoldi et al. 20 developed a set of different return rate metrics that are useful for studying short to intermediate term responses to pulse perturbations of locally stable communities. It is known that resilience is an intrinsic property of a locally stable equilibrium and thus independent of the specific pulse perturbation being imposed. However, at shorter timescales recovery dynamics of any given community depends on the direction of the perturbation. Therefore, Arnoldi et al. 20 developed an approach for investigating the median response across all perturbation directions, for different perturbation scenarios (i.e. assumptions on the relative strength of the initial perturbation on different species). One of these perturbation scenarios assumes Gaussian distributed perturbations that are on average of equal magnitude for all species in a community. We used this specific perturbation model as a null-model to investigate www.nature.com/scientificreports www.nature.com/scientificreports/ the relative importance of rare and common species for the different parts of the recovery process (from the initial to the asymptotic response) of locally stable ecological models exposed to random pulse perturbations.
Median initial return rate (R Init ). Initially, just after a pulse perturbation is being imposed to a stable equilibrium (i.e. at t = 0 + ), the median return rate across random perturbation directions R Init can be described by 20 : (where α ii is Lotka-Volterra intraspecific competition terms and N i the equilibrium density of species i); and secondly assuming that intraspecific competition, α ii , is approximately equal among species in a community; it can be seen that the initial return rate of a community is primarily governed by the abundance of common species. If the latter assumption is not true, other species may be important. Yet, a strong negative correlation between abundance and the strength of intraspecific competition is needed for common species not to have a major impact on the initial response to a perturbation.

Median intermediate return rate (R Interm ).
Intermediate return rates were found from numerical simulations of food webs and bipartite networks (See below). This was done by recording the time it takes until the magnitude of an initial pulse perturbation recovers to a proportion c of its initial size (See Supplementary Fig. 1): where T(c) is the time it takes until the magnitude of an initial perturbation has recovered and remains below a certain proportion c of its initial size. For 100 initial random perturbations given by: and ′ x i is a standard normal variable with a mean of zero and standard deviation of unity (i.e. ′ ∼ x N(0, 1) i ), numerical simulations were performed (Eq. 5) and the perturbation half-life (i.e. T(c) for c = 0.5) was recorded for each of these perturbations. The median across these was thereafter used to calculate an average return rate, given by 20 (See Supplementary Fig. 1): Interm 0 where x 0 is the pulse perturbation and x(t) are species' displacement from equilibrium at time t.
Asymptotic return rate (R Asympt ). Asymptotic return rate, i.e. resilience, was calculated using the dominating eigenvalue of the Jacobian matrix evaluated at stable equilibria (i.e. using Equation (7)).
Numerical analyses of the association between median return rate and species density. Numerical analyses were performed in order to disentangle the relative role of rare and common species for the different parts of the recovery process. First; for each of the 100 replicates of a given system (See below) median initial, intermediate and asymptotic return rates were derived as described above. Second; Spearman rank correlations were calculated for the association, across replicates of a given network, between equilibrium density of a species with a given density rank within a replicate and median return rate. Fig. 1; Eq. 9) was derived by simulating the recovery dynamics following a selective perturbations of magnitude ||x 0 ||=1, where the ith element of x 0 is equal to one, that is x 0 (i) = 1, and all other elements are zero.

Selective pulse perturbations. Perturbation half-life. The half-life of a perturbation (Supplementary
Selective press perturbations. Sensitivity of resilience. Sensitivity analyses were applied on both food webs and computer-generated bipartite networks to study the effect of a small change in the equilibrium density of a species (assumed to be caused by changes in the intrinsic growth rates of species) on network resilience. Sensitivity of resilience (S Λ ) to a small change in the equilibrium density (N k ) of species k is given by ΛŜ N Supplementary Information Section 8 for details). Here α kj is the per capita effect of species j on the per capita growth rate of species k, v and w are the left and right eigenvectors corresponding to the dominant eigenvalue of the Jacobian matrix of the network, v k and w j are elements k and j of the respective vectors, v k is the complex conjugate of v k , and <w, v> is the scalar product of w and v. Elasticity of resilience (E Λ ) to a change in the equilibrium density (N k ) of species k is given by ΛÊ www.nature.com/scientificreports www.nature.com/scientificreports/ The inverse Jacobian matrix. We used the inverse Jacobian matrix, a metric that has been extensively used elsewhere 18,19,34 , to investigate the direct and indirect effects of a small continuous density addition to one species on the equilibrium densities of all species in a community (See Supplementary Information Section 9 for derivation).
We investigated the total effect of a small continuous density addition to species j on all species equilibrium densities, given by:  Table 1 for descriptive statistics) were included in the analysis. These food webs were assembled from different ecosystems, and we used information on average body size of trophic groups to parameterize them. Food web networks, i.e. the binary network describing who feeds on whom, and body size information were retrieved from different sources (Food web networks: Baltic Sea 17 Here N i is the biomass density of species i in a community with n species, r i is the intrinsic growth rate of species i and α′ s ij are the strengths of interaction between and within species. For systems with linear functional response there is only one internal equilibrium point with all species present. For each natural food web we generated 100 feasible (all species have positive equilibrium densities) and locally stable replicates. Those replicates were generated by varying some of the model parameters (See below and reference 17 ). The replicates have a decent variability in equilibrium biomasses ( Supplementary Fig. 12) and the difference in equilibrium biomass between the most common and the rarest species in each replicate, across all food webs, is approximately three orders of magnitude  , m j is the average body mass of species j and P ij is the preference of species j for species i 17 . In the Baltic Sea and Lake Vättern the preferences are based on relative gut content data 17 and for the remaining food webs consumers are assumed to show equal preferences for all their prey species. The per unit biomass effect of prey i on consumer j is positive and related to the effects of consumers on resources by α α = −e ji ij where e is the conversion efficiency drawn from the uniform distribution [0.05, 0.2]. In order to generate feasible models of the natural food webs intraspecific competition (self-limitation) in species above the lowest trophic level, α ii , was set individually for each web (Baltic Sea, α ii = −0.0005; Broadstone Stream, α ii = −2.1; Lake Vättern, α ii = −0.005; Montane Forest, α ii = −0.01; Skipwith Pond, α ii = −1.5; Trelease Woods, α ii = −0.03; Tropic Sea, α ii = −0.005). For species at the lowest trophic level strength of intraspecific competition was set to −1. Species at the lowest trophic level were also assumed to directly compete with each other and the strength of this interspecific competition, α ij , was weak compared to that of the intraspecific competition and drawn from the uniform distribution [−10 −12 , 0]. Species at higher trophic levels competed with each other indirectly through shared resources.

Bipartite model networks.
Bipartite model networks were generated as in Mougi and Kondoh 30 . The networks consist of a mixture of competitive, mutualistic and antagonistic links, where relative fractions of interaction types are specified before adding dynamics 30 . Network connectance, C (defined by C = L/L max , where L is the number of links and L max is the maximum number of possible links in a bipartite network) is either low or high (C = 0.3 and C = 0.7, respectively); the proportion of mutualistic links, i.e. the proportion of the realized links that are mutualistic rather than antagonistic in the networks, is either low (0.3) or high (0.7); and the degree of intraspecific competition among species is either low (ω = 0.5), intermediate (ω = 1) or high (ω = 2). The dynamics is represented by continuous time models with type 1 functional response (See Supplementary Information Section 4 for a detailed model description). These bipartite model networks have a decent variability in equilibrium biomasses (Supplementary Fig. 13) and the difference in equilibrium density between the most common and the rarest species in each replicate, across all networks, is approximately two orders of magnitude www.nature.com/scientificreports www.nature.com/scientificreports/ code Availability All calculations were performed using Matlab 2018a, and the recovery dynamics was numerically simulated using the stiff ode solver ode23s. Figures were produced using the R software (version 3.4.3). The code generated during the current study are available from the corresponding author on reasonable request. www.nature.com/scientificreports www.nature.com/scientificreports/