The balance of interaction types determines the assembly and stability of ecological communities

What determines the assembly and stability of complex communities is a central question in ecology. Past work has suggested that mutualistic interactions are inherently destabilizing. However, this conclusion relies on the assumption that benefits from mutualisms never stop increasing. Furthermore, almost all theoretical work focuses on the internal (asymptotic) stability of communities assembled all at once. Here, we present a model with saturating benefits from mutualisms and sequentially assembled communities. We show that such communities are internally stable for any level of diversity and any combination of species interaction types. External stability, or resistance to invasion, is thus an important but overlooked measure of stability. We demonstrate that the balance of different interaction types governs community dynamics. A higher fraction of mutualistic interactions can increase the external stability and diversity of communities as well as species persistence, if mutualistic interactions tend to provide unique benefits. Ecological selection increases the prevalence of mutualisms, and limits on biodiversity emerge from species interactions. Our results help resolve long-standing debates on the stability, saturation and diversity of communities. A model of community assembly that is sequential and has saturating benefits of mutualism produces communities with internal stability for any level of diversity, and for which external stability depends on the balance of interaction types.

stone of ecology 2 . However, the field was given a jolt by the seminal theoretical work of May, who showed that complex ecological communities in fact tend to be less stable, in the sense of returning to equilibria from small perturbations to existing species 3,4 . This finding set off the great complexity-stability debate, as nature abounds with complex communities containing large numbers of species. A large theoretical and empirical literature developed studying different notions of complexity, diversity and stability 2,[5][6][7][8][9] . Yet despite this enormous attention, theoretical work has left important aspects of the question understudied.
The first understudied aspect is accounting for the diversity of interaction types between species (e.g., mutualism, competition and exploitation), which is expected to affect the stability and assembly of ecological communities. May's original analyses considered random interactions among species, while subsequent studies initially considered species networks connected only by a single interaction type, focusing on purely competitive, exploitative (predator-prey) or mutualistic communities. Only more recently has theoretical work started to incorporate multiple interaction types and ask how the balance of interaction types within a community affects its stability [10][11][12][13][14][15][16][17][18] . Most of these studies find that mutualistic interactions have destabilizing effects on communities, which are stable only when mutualistic interaction strengths are weak or asymmetric [13][14][15][16] . These results lead to the conclusion that mutualisms must play relatively small roles in complex communities.
However, this conclusion rests largely on the neglect of a second important aspect of the complexity-stability question: the role of nonlinear functional responses. It is easy to see that the destabilizing effect of mutualisms stems from the fact that unbounded positive feedbacks can arise from Lotka-Volterra models with a linear functional response, where the per capita effect of one species on another is independent of their population sizes. But linear functional responses are biologically unrealistic when studying mutualisms [19][20][21][22][23][24] . The benefits that a focal species receives per capita from a mutualist partner cannot increase without limit as the partner species' population size increases; they must saturate eventually 20 , as the focal species become limited by resources or services other than the one the mutualistic partner is providing. Nonlinear functional responses with saturating effects for mutualistic interactions are therefore more realistic models to consider the effects of mutualisms on community stability [20][21][22][23][24] . Yet they remain critically underused, having been used to model mutualistic communities in only a handful of studies 10,11,22,[25][26][27][28] .
Exploitative interactions are also expected to exhibit nonlinear saturating functional responses 24,[29][30][31] , since predators can become satiated with high levels of prey (either literally-that is, eating to their physiological limit-or time budget-wise, owing to the handling time for each prey). Basic theory shows that the effects of saturating functional responses in pairwise mutualistic and exploitative interactions should be opposite, stabilizing in the former case while destabilizing in the latter. Thus, the overall effect of incorporating saturating responses on community stability probably depends on the balance of different kinds of interactions.
Finally, most of the theory on the relationship between diversity and stability is concerned with the internal stability of communities-that is, whether a collection of species that coexist at an equilibrium will return to this equilibrium if perturbed (also called asymptotic stability). Theoretical models typically approach this question with random communities assembled all at once, where an interaction matrix of a certain size is populated from some distribution, and its stability determined. However, communities in nature assemble sequentially, as new species enter and existing ones go extinct. This means that community matrices drawn all at once may not be representative of natural communities, and that the external stability of communities (that is, whether new species can invade and existing species go extinct) plays an important role in determining their structures. With a few exceptions 32,33 , the relationship between external stability and diversity has been severely understudied relative to internal stability.
Here, we present an approach that jointly addresses these three neglected aspects of the complexity-stability debate. We construct a model of community assembly where communities emerge from the sequential invasion of randomly selected species, rather than studying community matrices assembled all at once. We vary the proportions of competitive, exploitative and mutualistic interactions in the invader pool to ask how the balance of interaction types affects the structure and stability of assembled communities. We use saturating functional responses for mutualistic and exploitative interactions. We show that in this setting, the problem of internal stability effectively vanishes, while the problem of external stability-Elton's original question-comes into focus.

results
We first give an overview of our model (see Methods for a detailed model description). We model sequentially assembled communities of species that can have competitive, mutualistic or exploitative interactions with each other. We use saturating (Holling type II) functional responses 20,29 with the half-saturation constant h for mutualistic and exploitative interactions. This setting raises a question that does not come up with linear functional responses: whether each mutualistic or exploitative interaction is unique, or whether they are interchangeable with each other. For mutualisms, if each partner provides a unique type of benefit (different nutrients or services) to a focal species, the interaction strength should saturate in that partner's density only. We call this scenario the unique interactions model (UIM). At the other extreme, if all mutualisms provide the same type of benefit, then the functional responses should saturate in the total density of all mutualist partners. We call this scenario the interchangeable interactions model (IIM). Similar considerations apply to exploitative interactions, where prey and predators may be unique or interchangeable with each other. Both types of models have been considered in previous studies 10, 28 , though their results have not been directly compared. We would expect interactions in nature to be at some intermediate between the completely unique and completely interchangeable models. We thus provide results for both the UIM and IIM, focusing on the UIM to introduce our main conclusions, and later comparing the results with the IIM. In Supplementary Section 1, we consider a model that is intermediate between the UIM and IIM, which shows that the results from the UIM are more robust. We assume that competitive interactions retain the oft-assumed linear (type I) functional response, and hence that the issue of interchangeability does not arise for competitive interactions.
Each community starts with a small number of species and grows by the repeated invasion of new species. Each time a new species arrives in the community, its interactions with existing species are generated randomly, with average connectivity (that is, the probability of any interaction) C and frequencies P m , P c and P e for mutualistic, competitive and exploitative interactions, respectively. When we introduce a new species, it can either have a positive growth rate when rare and be able to invade, in which case it is added to the community, or have a negative growth rate when rare, and thus is not added to the community. In the former case, we simulate the dynamics of the community until it settles to an ecological steady state, which might be an equilibrium or a limit cycle. We introduce a new species to this steady state, reflecting a separation-of-timescales assumption [34][35][36] . At each equilibrium, we generate new species until an invasion occurs and track the probability of invasion. Species can also go extinct. We repeat this process until the community assembly reaches a steady state in species richness, quantified using a unit root test. In the main text, we present the results setting C at 0.5, but all our findings are robust across different values of C and h (Supplementary Information).
In our model, new species are introduced into communities without directly controlling the species richness. In other words, the biodiversity of a community is constrained by emergent limits on community complexity. We find that the relative proportions of interaction types govern the emergent species diversity in the assembled communities ( Fig. 1). Notably, as long as communities have nonzero competition, they saturate as large numbers of invasions occur and the species richness converges to some steady state (Fig. 1a). The steady state is characterized by small, stochastic oscillations around some baseline value, so that a sequence of consecutive invasions is offset by a subsequent cascade of extinctions, even when thousands of invasions occur. This holds for both the UIM and IIM. In both models, communities with more mutualism and less competition show a higher diversity at this steady state (Figs. 1a,c and 6a), with the mutualism effect more pronounced in the UIM, and the competition effect more pronounced in the IIM. Regardless of C or h, the community diversity saturates (Supplementary Figs. 8 and 9 for UIM; Supplementary Figs. 25 and 26 for IIM). The lower the value of C, the higher the steady-state species richness, but the balance of interaction types yields similar qualitative patterns regardless of C and h.
Asymptotic stability provides information on whether communities return to an equilibrium state after small perturbations 3 . In our simulations, the population dynamics almost always converge to an equilibrium, in rare cases settling to stable oscillations ( Supplementary Fig. 10). Populations never grow boundless, owing to the saturating functional responses for mutualistic and exploitative interactions as well as intraspecific density dependence. The probability of oscillations is very low for all proportions of interaction types for both the IIM and UIM ( Fig. 1d and Supplementary Fig. 27), and this holds for almost all combinations of C and h (Supplementary Figs. 11 and 27). In particular, we find that mutualistic interactions do not result in asymptotically unstable communities even when they are strong and abundant, consistent with past results in purely mutualistic communities with saturating functional responses 22 . Instead, internal stability is a ubiquitous property of sequentially assembled communities with saturating functional responses in mutualistic and exploitative interactions.
In contrast, we find that the external stability, or the resistance of communities to invasion, is a critical determinant of community structure and function. Even though this was the original question motivating Elton 1 , and many authors noted that invasibility can be a measure of community instability 2,8,9,37,38 , resistance to invasion from the outside has largely been overlooked in the theoretical debate on the complexity-stability relationship.
The species richness of communities is mediated through the balance of invasion and extinction rates, which in turn depend on the balance of interaction types. Figure 2a,b shows that in the UIM, the probability that a new species can invade decreases with higher P c , as one might expect. However, it also shows that, contrary to what one might expect, higher values of P m also lead to lower probabilities of invasion. This is due to the fact that highly mutualistic communities allow resident species to achieve high population sizes, such that an invader with a competitive interaction with a resident will face strong competition. Because the benefits from mutualistic interactions saturate but the cost of competitive interactions does not, even a few competitive interactions can be enough to keep an invader out. Mutualisms with saturating functional responses thus increase the competitive burden to invaders and render communities externally stable.
More specifically, for communities with moderate P e , increasing P m while decreasing P c monotonically enhances external stability. At low P e , external stability is maximized at intermediate P m . Conversely, while holding P c constant, external stability increases with increasing P m (except when P c is very low). For low and moderate P m , increasing P c (and thereby decreasing P e ) augments stability, whereas when P m is high, decreasing P c increases stability.
These trends are observed both throughout the entire community history (Fig. 2a) and during the steady state after the diversity has converged (Fig. 2b). A community's stability does not change as it leaves its transient phase (before the species diversity converges) for the steady state (Fig. 2c). Any changes in the probability of invasion are small, and there is no clear pattern throughout the different communities. These trends are qualitatively robust for all C and h .
In contrast to the UIM, mutualisms do not markedly increase the external stability of communities in the IIM ( Fig. 6a and Supplementary Fig. 28). Under the IIM, invasion probabilities are generally low, show no strong pattern with P m and weakly decrease with increasing P e . Looking at the mean population sizes in the IIM confirms that this contrasting result is due to the fact that varying P m and P e in the IIM does not cause the population sizes to change substantially (Fig. 6d). The mechanism that stabilizes communities against invasion in the UIM is thus unavailable in the IIM. This shows that mutualisms can be externally stabilizing as long as they allow resident species to achieve high densities.
Our model also illustrates how species interactions can modulate the relationship between species richness and resistance to invasion in different ways. Under the UIM, we can observe both a positive and negative relationship, depending on which interaction proportions are varied. When holding P m constant and trading off P e against P c , the species richness and invasion probability are positively correlated, producing what has been termed the invasion paradox 39 (Figs. 1c and 2a,b, and Supplementary Figs. 9, 12 and 13). However, when keeping P c constant and trading off P m against P e , we get a negative relationship between species richness and invasion probability. Here, more mutualistic communities have more species but are less likely to be invaded. In contrast, under the IIM, the relationship between invasion probability and species richness is largely negative, and modulated by P c (Fig. 6b,c).
The other determinant of community richness is the probability that an existing species goes extinct. Our model shows that this probability also depends on the balance of interaction types. We define the probability of extinction simply as the number of species that went extinct divided by the total of number of species that existed during Each line corresponds to a community over 1,500 species invasions that occur at equilibria. The same trend is observed for all P c > 0. b, A schematic on how to read the ternary plots 64 . The three lines show the directions in which each parameter is constant. These are in the same directions as the corresponding axis tick marks. Each hexagon represents a single combination of the three interaction types, given by the three axes. c, A ternary plot showing the steady-state species richness of different communities with varying P c , P m and P e . The colours represent the species richness. Grey communities have no competition and are not considered since they are biologically unrealistic. d, The probability that an equilibrium in the community history shows oscillations in population dynamics. This probability is defined as the number of times oscillations occur divided by the total number of equilibria. The colours represent the probabilities.
the relevant part of the community history. In the UIM, for any constant level of P c , decreasing P e and increasing P m decrease the probability of extinction ( For constant moderate or high P e , extinction is minimized for low P c and high P m (for high h, this is true for all levels of P e ). With constant P m , increasing P e or decreasing P c decreases extinction. As expected, these trends are the opposite of those seen in Fig. 1c for steadystate species richness. In the IIM, the extinction probability again is mostly a function of P c ( Fig. 6c and Supplementary Fig. 31), without very strong dependence on P e or P m . In both the UIM and IIM, community saturation and the emergence of a limit on biodiversity occur due to changes in the extinction probability as opposed to invasion, which stays constant before and after the convergence of species diversity to the steady state ( Fig. 2c and Supplementary Figs. 14 and 30), whereas the probability of extinction increases after the steady state is reached ( Fig. 3c and Supplementary Figs. 17 and 33). The steady-state species richness is maintained when the invasion of new species precipitates the extinction of older species, rather than through a resistance to invasion by new species. As a result, the distribution of species persistence is skewed right for all communities, and most species do not survive very long ( Fig. 4a and Supplementary Fig. 18). To quantify how persistence varies with species interactions, we define relative persistence as the number of equilibria in which the species remains extant normalized to the number of equilibria in the respective simulation. In the UIM, high P m and low P e both favour slow species turn-over and longer species persistence ( Fig. 4 and Supplementary Figs. 19 and 20) as they both reduce the invasion and extinction probabilities. In particular, keeping P c constant and increasing P m lengthen the tail of the persistence distribution so that some species survive much longer than others (Fig. 4a). In contrast, species persistence in the IIM again depends mostly on P c (Supplementary Figs. 34 and 35), which modulates both invasion and extinction probabilities.
Finally, we consider the ecological selection (or filtering) of interaction types during community assembly. Our model throws at communities randomly generated species that have their own values of P m , P c and P e (we call this the extrinsic distribution). Although candidate invader species are generated randomly (see Methods), the ones that are able to invade and persist in a community may be biased towards one or the other interaction type. Indeed, we observe the signature of such selection acting on the distribution of interaction types in a community in both the UIM and IIM. We find that selection acts to decrease P c , increase P m and slightly increase P e , with a more dramatic increase in P e at lower P m ( Fig. 5 and Supplementary Fig. 36). This selection is caused by a combination of biases in the invasion and extinction rates. Successful invaders generally show the same bias for a higher P m and lower P c (Fig. 5 and Supplementary Fig. 37). This is not surprising since an invader that benefits from existing species will have an easier time invading than one that suffers competition from them. One exception is in the UIM, when P e is high and P c is low, where invaders tend have a higher P c . Extinction generally is biased in the same way as the overall community composition ( Fig. 5c and Supplementary Fig. 38), except in the high-P e , low-P c cases in the UIM. The differences in interaction types between species that invade and those that go extinct are generally small ( Fig. 5d and Supplementary Fig. 39), except again in the UIM for high-P e , low-P c communities, where extinction is biased towards lower P c . The finding that ecological selection increases P m is generally robust to variation in C or h (Supplementary Figs. 21-24 and 36-39), though the latter can have slightly differing effects on other interaction types.

Discussion
Our model incorporates two key biologically realistic assumptions: saturating functional responses for mutualistic and exploitative interactions, and sequential selection of candidate invader species during community assembly. We show that under these assumptions, the internal stability of communities is always guaranteed (in the form of either equilibria or stable limit cycles) with all combinations of interaction types, species diversity and community connectivity. Part of this result can be explained by the fact that, by construction (and for biological realism), we restrict our attention to feasible communities of species that actually have positive abundances at equilibria. An early but relatively overlooked result by Roberts 40 shows that feasible equilibria tend to also be internally stable. Consistent with this, Stone 41,42 showed that as the variability of interaction strengths increases, communities become unfeasible (some species go to negative abundance at equilibrium) sooner than they become asymptotically unstable, and that in a broad range of cases, feasibility implies asymptotic stability. Feasible communities thus tend to be stable. In addition, we find that mutualistic interactions with saturating functional responses promote internal stability, even if saturation happens only at high population densities, also consistent with previous results 28 . It is worth noting that we assume that all species experience some (relatively weak) self-regulation (that is, within-species negative density dependence). We share this assumption with the large body of theoretical work on community stability and diversity, although we differ from most of these studies in considering the sequential assembly of communities. Recently, Barabás et al. 43 showed that in empirically constructed food-web models, a high proportion of self-regulation is necessary and sufficient for asymptotic stability even with weak self-regulation, which is concordant with our finding in sequentially assembled communities. Our results contribute to the growing realization that internal stability is a ubiquitous property of communities regardless of their complexity and composition of species interactions, and therefore unlikely to be a differentiating feature of natural communities. However, external stability (that is, the ability of a community to resist invasion) emerges in our model as a crucial factor that varies with community structure and is pivotal in determining species richness. Yet, external stability has received much less attention than internal stability 1,8,9,32,33,37,38 . A notable exception is a model by Law and Morton 33 where they consider the sequential assembly of a purely exploitative community (that is, a food web) from a finite species pool. By considering the endogenous assembly of communities from a large source pool with different types of interactions, we are able to answer a major outstanding question of ecology: how does the balance of species interactions drive community assembly and structure? Our model reveals how the external stability of communities, their steady-state richness and the extinction and persistence of organisms are determined by the balance of interaction types.
Our results show that mutualisms can play a key role in determining the assembly and stability of communities. In the UIM, we find that when holding either P c or P e constant, mutualism augments internal stability, decreases extinction, increases species persistence and thereby decreases species turn-over. In the IIM, however, where mutualistic (and exploitative) interactions saturate in the total density of all interaction partners, these effects are absent or weaker. In particular, mutualisms no longer clearly increase external stability in the IIM. This is because adding more mutualistic interactions in the IIM does not actually lead to the generation of more gross mutualistic benefit in the community, as most functional responses tend to be in the saturated ranges already, and adding new mutualistic partners does not noticeably increase the mutualistic benefit that a focal species gets. More mutualistic interactions in the IIM thus do not lead to (much) higher population sizes for resident species, which takes away the mechanism by which mutualisms stabilize communities against external invaders in the UIM. Importantly, results from models that are intermediate between the UIM and IIM suggest that the patterns we document in the UIM obtain even when only a small fraction of interactions are unique (see Supplementary Section 1). Our model therefore predicts that to the extent that each mutualistic interaction provides some unique benefits, adding more mutualistic interactions will increase resident species' densities, and will render communities more externally stable. In nature, it is likely that at least some benefits from each mutualistic interaction are unique. Many mutualisms provide fundamentally different benefits (defence, nutrients, pollination and so on), which will not be interchangeable. Even within the same kind of interaction, adding more partner species can result in higher benefits to a focal species 44 . These kinds of unique or complementary benefits were shown experimentally in plant-mycorrhizal interactions 45,46 , antplant symbioses 47 and marine defensive mutualisms 48,49 . These studies suggest that mutualistic partners even from the same functional category can provide some unique benefits.
We find that mutualisms also raise the emergent limits on diversity, with a stronger effect in the UIM, allowing the community to saturate at a higher species richness. Accordingly, selection acts to increase mutualism during community assembly. Our demonstration that the ratio of P m to P c is an important determinant of community diversity, with a high ratio leading to high diversity, is consistent with a previous empirically derived model 17 . However, our findings contradict generalizations of May's original stability-complexity results 3 , which show not only that mutualism is destabilizing but also that it is even more destabilizing as the product of species richness and C increases 13,14,50 . Such a destabilizing effect of mutualism contradicts empirical observations of frequent mutualistic interactions in communities with high diversity 50,51 . In our model, the internal stability is not compromised with increasing species richness since sequentially assembled communities are selected for feasibility and the benefits from mutualisms saturate at high densities, both realistic features of communities in nature. The first feature means that communities can get richer only if they remain feasible, and therefore no species blows up in density, which would cause some other species to go extinct. The second means that the blowing up in density is unlikely in the first place, as the positive feedback between mutualist densities is cut short at high enough densities.
Previous work by Mougi and Kondoh 10-12 also showed that mutualisms may not be inherently destabilizing, but these authors focused exclusively on the internal stability of fixed-size communities. Furthermore, they assumed that the total interaction strength of each species for each interaction type is constrained. This means that as mutualisms become common, each mutualistic link gets weaker, in addition to interactions saturating in total species abundance (similar to our IIM) in their saturating functional response model. In contrast, Kawatsu and Kondoh 28 considered a case that is closer to our UIM, and showed that saturating functional responses can stabilize fixed-size communities as in our model. We add to these results by showing that saturating responses with ecological selection is enough to almost guarantee internal stability in both the UIM and IIM. Another recent line of work by Butler and O'Dwyer 52,53 explicitly modelled species limited by resources that can be produced by other species. This model allows saturating functional responses in mutualisms to emerge endogenously from resource limitation. Consistent with our results, they found that feasible mutualistic communities are stable even if mutualistic links are strong, provided that the exchange of resources is reciprocal between species. More generally, resource-mediated interactions tend to produce saturating functional responses as opposed to linear ones 54 . An important future direction is to model the external stability of communities with resource-mediated interactions.
One of our main findings is that all of the communities in our model show an emergent limit on species richness. This is true regardless of C or the distribution of interaction types (as long as there is non-zero competition, which is certainly true in natural communities) or whether interactions are unique or interchangeable. In other words, all communities become saturated when mutualism and exploitation have saturating effects and communities are sequentially assembled. This result contributes to an ongoing debate that has centred on the relative importances of local biotic interactions and regional processes in determining community richness and saturation [55][56][57][58][59][60] . We show that the composition of local species interactions can generate an intrinsic limit on community diversity. The exact effect of the balance of interactions on species richness depends somewhat on whether interactions are unique or interchangeable, but diversity generally is higher with lower P c and higher P m . This saturation comes about even though we effectively model an infinite species pool with interactions randomly drawn according to the desired proportions, rather than a finite regional or metacommunity pool. Regardless of how many species attempt to invade the community, the species richness fluctuates around some steady-state value that depends only on the distribution of interaction types. Our results thus suggest that species interactions are sufficient for community saturation on ecological timescales.
Our results on species richness and its equilibration add a new dimension to MacArthur and Wilson's classical theory of island biogeography 61 , which states that the species richness of an island is controlled by immigration and extinction patterns. We show that these patterns are in turn governed by the balance of species interactions. We therefore posit that the distribution of interaction types is a fundamental determinant of island biodiversity, and of community assembly in general. MacArthur and Wilson investigated how an island's area and distance to the mainland affect the rates of immigration and extinction, and they predicted an equilibrium of species diversity when these two rates are equal. We extend their analysis by considering how species interact with one another. We show that all communities can exhibit a dynamic equilibrium of species richness. Their equilibration model requires the extinction rate to increase with species diversity, and we provide insight into the underlying proportions of interactions that allow this to happen. In fact, we reveal that an increasing extinction rate (rather than a decreasing invasion rate) is indeed what drives the convergence of diversity. Our finding that species interactions can explain the distribution of island species is consistent with previous results 32 .
The patterns we uncover also shed new light on the relationship between the diversity and invasibility of a community, where evidence for both a negative and positive relationship exists 39 . Elton's hypothesis that diversity decreases invasibility implicitly assumes that competition is the primary driver of community composition and invasibility. This reasoning was confirmed by Case 32 , who modelled the invasion of new species and extinction of resident ones in purely competitive communities. We show that accounting for different interaction types may help explain the invasion paradox. Rather than a direct correlation between resistance to invasion and species richness, we find that species richness (Figs. 1 and 6a) and the probability of invasion (Figs. 2 and 6b) are independently modulated by the distribution of interaction types. Our results thus renew the focus on the balance of interaction types as a unifying theme for explaining the invasibility of ecological networks 62 .
The role of species interactions in community assembly has long been a central concern of ecology 55,56 . It has also recently become a focus of studies on community stability 10,13,14,18 . Our results suggest that the balance of species interactions is a fundamental driver of community dynamics. Numerous studies have examined the relationships between community stability, invasibility, diversity and saturation, but we argue that species interactions operate at a more elementary level and are the determinants of such other, higher-order properties of communities. In particular, mutualisms may play a more important role than suggested by previous studies that assume the linearity of its functional response. We show that considering sequentially assembled ecological communities with biologically realistic nonlinear functional responses allows a resolution of the complexity-stability debate and brings into focus other important properties of communities.

Methods
We construct a model of community assembly to study how communities with desired properties are formed over time. The species richness of a community changes over time and is denoted S. The community is initialized with S 0 species and grows in size as successful invasions by new species occur. One of the key features of our model is the ability to consider all combinations and proportions of all interaction types in the construction of the interaction matrix A of the community, similar to Coyte et al. 14 . We likewise limit our model to mutualism, competition and exploitation for simplicity, though our model can easily be generalized to include commensalism and amensalism. As model parameters, we impose a desired C of A and desired values of P m , P c and P e . The normalization constraint is simply P m + P c + P e = 1.
For the UIM, the population dynamics of species i are given by and for the IIM, the population dynamics of species i are given by where X i is the population of species i, r i is its intrinsic growth rate, a ij is the interaction coefficient between species i and j, s i is the density-dependent selfregulation term (and is negative) and K i is the carrying capacity. C i I is the set of interactions between i and its competitors. Likewise, M i I is the set of interactions between i and its mutualistic partners. E þ i I refers to the set of interactions between i and the species that it exploits, while E À i I refers to the set of interactions in which i is getting exploited. In equation (1), the separation of the summation over these sets is simply the standard implementation of the predator-prey type II functional response, where interaction effects for both predator and prey saturate with increasing populations of prey. In both the UIM and IIM, we can account for any variation in P c (both a ij and a ji are negative), P e (a ij and a ji are opposite signs) and P m (both a ij and a ji are positive) throughout the community.
The UIM and IIM are implemented identically except for the population dynamics given above. Each model begins with S 0 species whose population dynamics are governed by equation (1) or (2). Their initial populations are set equal to a constant x 0 multiplied by a random variable drawn from a uniform distribution on (0, 1). To focus on the effects of interaction types, we set the intrinsic growth rates, density-dependent self-regulation terms and carrying capacities of all species to be equal, so r i = r, s i = s < 0 and K i = K for all i. We start with an S 0 by S 0 interaction matrix A 0 , which stores the interaction coefficients of equation (1) or (2). The diagonal values of A 0 are set to s. For each pair of species (that is, each (i, j) combination where i ≠ j), the probability that an interaction exists is given by C; that is, with probability 1 − C, we set the interaction coefficients a ij = a ji = 0. For each pair i and j that does interact, the interaction is mutualistic with probability P m , in which case a ij and a ji are both drawn from a half-normal distribution N(0, σ), where σ is the scale parameter. Likewise, with probability P c , i and j compete with each other, in which case both a ij and a ji are drawn from a negative half-normal distribution −N(0, σ∕K), where the scale parameter is scaled by the carrying capacity. This is such that the average effect of competition on population dynamics is comparable to that of mutualism. Finally, with probability P e = 1 − P m − P c , i and j engage in an exploitative interaction. In this case, we assign one of i or j randomly as the exploited species, while the other becomes the exploiter. If i is the exploiter, a ij is drawn from the halfnormal distribution N(0, σ) and a ji from −N(0, σ), and vice versa if j is the exploiter.
The construction of A 0 follows the approach developed by Coyte et al. 14 . We then add community assembly to the model and use a different manner of simulating population dynamics. After A 0 is created, we simulate the population dynamics of all S 0 until the populations come to an equilibrium. The integration of the Lotka-Volterra dynamics occurs with adaptive step sizes using the Dopri5 integrator 63 . The equilibrium of populations in our model can be achieved in two ways. First, we impose a population change threshold δ, so that if in a single time step every single population fluctuates by less than δ, then the populations are determined to have equilibrated. Second, for simulation purposes we impose a time limit t l so that if the populations still have not equilibrated after t l time steps, we use the current population sizes as the equilibrium values. We choose a large enough t l so this does not occur unless populations have entered a stable limit cycle and are oscillating. As species populations change in size, some species may go extinct. We define an extinction threshold ϵ so that any species whose X i dips under ϵ is considered extinct and removed from the community at the next equilibrium.
After the populations have equilibrated, we add a new species z into the community, so that A increases in size from S 0 by S 0 to S 0 + 1 by S 0 + 1 (in the case that no species have gone extinct in the first iteration of the model). Entries in the new row and column that are added to the interaction matrix are drawn in the exact same manner as described above for the initialization of A 0 , so the desired properties C, P m , P c and P e are preserved as the community changes in size. Importantly, not every species introduced randomly in this manner will be able to invade the community. For a successful invasion to occur, the growth rate of z must be positive when a small but non-zero number of z individuals are suddenly added into the community. If the invasion condition is not met and a failed invasion occurs, we replace the failed invader with a new species whose interaction coefficients are drawn anew in the same aforementioned manner. We repeat this process until a new species is able to invade the community. We track the number of failed invasions. For computational tractability, we impose a failed invasion threshold β so that if over β new species fail to invade a community at any given equilibrium, the community is deemed externally stable and the simulation ends. This scenario is not reached in the communities we simulate.
If the invasion condition is met, the community grows by one species. The new species is initialized with a small population size equal to a constant x z multiplied by a random variable drawn from a uniform distribution on (0, 1). We then simulate the population dynamics until the next equilibrium, at which point we introduce another new species, and so on. In this manner, at each equilibrium a new species is added into the community. The species richness and probabilities of invasion and extinction are tracked. We impose a limit n on the species richness of a community for computational tractability so that if S exceeds n, the simulation ends (this does not occur in our simulations). As mentioned in the main text, communities exhibit a convergence of species richness. The presence of a steady state of species diversity is assessed using an augmented Dickey-Fuller test, a unit root test whose null hypothesis is that the species richness time series has a unit root (that is, it is non-stationary and has a time-dependent structure) and whose alternative hypothesis is that species richness does not have a unit root and is stationary over time. We analyse the species richness of the most recent w 1 equilibria of the community and reject the null hypothesis when the P value is below some threshold ρ. The rejection of the null hypothesis indicates that the species richness has entered a steady state. If a community enters the steady state, the simulation runs for w 2 more equilibria before it is terminated. Simulations with communities that never enter the steady state can terminate in one of three ways: by reaching n, by reaching the limit on the number of equilibria (κ) or by reaching t l , though in our simulations these limits are never reached.
For each combination of C and h, we simulated communities by varying P c , P m and P e from 0 to 1 in intervals of 0.1, subject to the condition that their sum is equal to 1; there are 66 such combinations of interaction types, and we simulate each combination with five replicates to get 330 communities. Since there are 9 combinations of C and h, there are a total of 2,970 communities. We discard communities with P c = 0 during our analysis, since communities at this edge case do not saturate in species richness. Since real communities are expected to have at least some competition, this does not affect our conclusions. Supplementary Table 1 summarizes the parameters we used in our simulations. We choose s and σ such that the magnitude of the intraspecific interaction is generally larger than the magnitude of interspecific interaction. This biologically realistic assumption easily arises from the fact that members of the same species compete more strongly for the same set of resources and space than do members of different species. In the main text, we use C = 0.5 and h = 100, while in the Supplementary Information we extend our analysis to other values of C and h. We perform the full analysis for both the UIM and IIM.
All data shown in the figures in the main text are averages over five replicates, except for Figs. 1a and 4a, which show data from representative communities, and Fig. 5. To study selection (Fig. 5), we simulate communities until the species diversity converges to a steady state. We take advantage of the fact that the communities appear to exhibit ergodicity (as suggested by Fig. 1 and Supplementary Fig. 8). We continue to randomly create new species and allow the community to come to 100 S s additional equilibria, where S s is the steady-state species diversity. At each equilibrium, species are introduced until a successful invasion occurs. We sample the community 100 times at an interval of every S s equilibrium and calculate the average distribution of interaction types (Fig. 5a). After the steady state is reached, we save the interaction types of every tenth successful invader and every tenth species that goes extinct. The average interaction proportions are calculated and plotted (Fig. 5b-d).
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Data sharing is not applicable to this article as no datasets were generated or analysed during the current study. All simulation results can be reproduced using the code provided.

code availability
The model of community assembly was implemented in Python. The simulation code is available at https://github.com/erolakcay/CommunityAssembly.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted Give P values as exact values whenever suitable.

For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings
For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated Our web collection on statistics for biologists contains articles on many of the points above.

Software and code
Policy information about availability of computer code Data collection All results presented here are from numerical simulations, written in python 3.7. The code is available on the corresponding author's github page (linked in the manuscript)

Data analysis
All results presented here are from numerical simulations, written in python 3.7. The code is available on the corresponding author's github page (linked in the manuscript) For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability All data are simulation data, and can be replicated using the code provided.

Field-specific reporting
Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

Life sciences
Behavioural & social sciences Ecological, evolutionary & environmental sciences