Ecological theory of mutualism: Robust patterns of stability and thresholds in two‐species population models

Abstract Mutualisms are ubiquitous in nature, provide important ecosystem services, and involve many species of interest for conservation. Theoretical progress on the population dynamics of mutualistic interactions, however, comparatively lagged behind that of trophic and competitive interactions, leading to the impression that ecologists still lack a generalized framework to investigate the population dynamics of mutualisms. Yet, over the last 90 years, abundant theoretical work has accumulated, ranging from abstract to detailed. Here, we review and synthesize historical models of two‐species mutualisms. We find that population dynamics of mutualisms are qualitatively robust across derivations, including levels of detail, types of benefit, and inspiring systems. Specifically, mutualisms tend to exhibit stable coexistence at high density and destabilizing thresholds at low density. These dynamics emerge when benefits of mutualism saturate, whether due to intrinsic or extrinsic density dependence in intraspecific processes, interspecific processes, or both. We distinguish between thresholds resulting from Allee effects, low partner density, and high partner density, and their mathematical and conceptual causes. Our synthesis suggests that there exists a robust population dynamic theory of mutualism that can make general predictions.

1985; Bronstein, 2015b;Raerinne, 2020), part of which we briefly describe below. We submit that ecology will benefit from integrating this coherent and robust body of theoretical work. Here, we contribute a first step toward such integration by presenting the ecological theory of mutualism available to the broader ecological community.
Specifically, we review its historical literature and synthesize generalities, both mathematical and conceptual, that can lay a foundation for a deeper understanding and integration of mutualism in ecology.
Foundational theory in ecology was initially developed using Lotka-Volterra models. In this framework, constant coefficients describe the positive or negative effects between two interacting species as a linear function of the other species' density. The Lotka-Volterra model for predation and competition predict stable cycles (oscillations, Lotka, 1925;Volterra, 1926) and competitive exclusion (Gause, 1934;Volterra, 1926), respectively, which stimulated fruitful empirical and theoretical work. Indeed, from the groundwork of Lotka-Volterra theory of predation came more general consumerresource theory, with useful and surprising results such as the paradox of enrichment (Rosenzweig, 1971) and a mathematical representation of seasonal cycling in lake food webs .
In contrast, Lotka-Volterra models for mutualism have been a less useful simplification than for predation and competition (Holland, 2015). The original model (Gause & Witt, 1935) and other formulations in which species benefit as a linear function of each other's density (Addicott, 1981) can predict unbounded population growth of both species. Additionally, the diversity of mechanisms by which species may benefit each other and the non-reciprocity of many of them, has cast suspicion on representing any "mutualistic" interaction as a simple exchange of positive effects (Bronstein, 2001a(Bronstein, , 2001b. Mutualisms are more likely to exhibit shifting net effects than other interaction types (Chamberlain et al., 2014;Jones et al., 2015), with several exchanges dipping, for example, into parasitism.
Despite all these interesting mechanisms and patterns ripe for study, mutualisms have been subjected to less theoretical study than predation and competition. Many have speculated on historical reasons (Boucher, 1985;Bronstein, 2015b;Raerinne, 2020), but we highlight two here. First, the terms used to identify interactions as "mutualism" have changed over time. Previous theory treated mutualism as a subset of facilitation, in which one species alters the environment to benefit a neighboring species (Callaway, 2007), or symbiosis, in which species coexist in "prolonged physical intimacy" (Bronstein, 2015b), or used those terms interchangeably.
Additionally, the terms "mutualism," "cooperation," and "protocooperation" have been used idiosyncratically for beneficial interactions within species as well as between them (Bronstein, 2015b). Furthermore, some mutualisms are "indirect," such that benefits to one partner can only be realized in the presence of an external species or environmental condition (Holland & DeAngelis, 2010). In this review, we limit our scope to mutualism defined as reciprocally beneficial interactions between two species (Bronstein, 2015b). We largely focus on direct mutualism or models that approximate the effects of indirect mutualism through two-species models, though we touch on some other cases (e.g., Thompson et al., 2006). Second, the mechanisms by which species benefit each other in mutualisms are extremely diverse. These mechanisms include, but are not limited to, habitat provisioning, deterrence of predators or competitors, increased growth, faster maturation, facilitated reproduction, improved digestion, parasite grooming, and resource consumption. Conceptual frameworks have attempted to organize this rich diversity, for example, by the types of benefits exchanged (nutrition, protection, or transportation), the mechanisms of exchange, or the obligacy of each partner (reviewed in Bronstein, 2015b;Douglas, 2015). This diversity of mechanisms makes the development of general but informative theory for mutualism more difficult than, for example, predator-prey theory, in which the interaction can be simply modeled as the consumption of individuals of one species by the individuals of the other species.
Here, we review ecological theory of mutualism, tracing authors' attempts to understand how mutualisms can persist stably overtime and synthesizing their results. We begin with an in-depth historical review of the theoretical study of mutualism, highlighting many now-obscure texts that have contributed to the field's current understanding. We focus exclusively on two-species populationdynamic models, leaving other aspects of historical mutualism research including game theory, biological market models, and ecoevolutionary dynamics to previous (excellent) sources (Bronstein, 2015b;Hoeksema & Bruna, 2000). We organize the development of the theoretical study of mutualism semi-chronologically, by its historical focus on the form of benefit either as linearly increasing with partner density or limited by intraspecific or interspecific density dependence, and its more recent incorporation into consumerresource theory (summarized in Table 1). After reviewing this rich and often overlooked body of work on the ecology of mutualism, we identify patterns in the predictions of these models that stand across systems and assumptions. In particular, we disentangle common terminology in order to clarify mechanisms that lead to predictable ecological dynamics. We find diverse, well-characterized ecological mechanisms that permit stable coexistence. We additionally find that mutualisms are characterized by thresholds in density that may cause system collapse, which can be explained by partner dependence and interaction strength. We argue that extant models make a robust set of qualitative predictions and that these predictions qualify as an ecological theory of mutualism.

Neuhauser and
Fargione (2004) Facultative only Graves et al. (2006) 10 Table S1 Thompson et al. (2006) Hermit crabs (N 1 )-Anemones (N 2 ) Closed system when I i = 0, i = 1 Obligate when Table S1 Holland and DeAngelis Fishman and Hadany (2010) Obligate only Kang et al. (2011) 16 Fungal garden (N 1 )-Leaf cutter ant (N 2 ) Obligate only Martignoni et al. (2020) 17 Below, we provide an in-depth description of this theoretical development (summarized in Table 1). We focus on the bulk of theory that conforms to the typical assumptions of population dynamic approaches (Gotelli, 2008). That is, we focus on models without immigration or emigration (i.e., closed systems), without age, stage, or genetic structure, and under the approximation that individuals encounter each other randomly with no spatial structure (mean field assumption). These models have tended toward increasing analytical F I G U R E 1 Characteristic dynamics for linear benefit models. In early models of mutualism, benefits were represented by a constant coefficient (interactions strength) multiplying a linear function of partner density. Benefits were modeled as affecting per-capita growth rate (low-density effect, Equation 4), equilibrium density (high-density effect, Equation 2), or both (Equation 1, see Table 2). When benefits have exclusively low-density effects, nullclines (curves of zero growth) are simply vertical (N 1 ) and horizontal (N 2 ) lines, always resulting in stable coexistence (qualitatively similar dynamics to those in a). Otherwise, the nullclines are linear, increasing curves, with different potential dynamics (a-d). When both partners are facultative mutualists (N i = K i > 0 when N j = 0), they display stable coexistence when benefits are weak (a) or grow without bound (unstable coexistence) when benefits are strong (c). When both mutualists are obligate upon their partner (N i = K i ≤ 0 when N j = 0) and benefits are weak, the system exhibits a threshold in density above which species exhibit unbounded growth and below which extinctions occur (b), whereas if benefits are strong, only extinctions occur (d). When mutualists are a facultative-obligate pair, any of the previous results can occur depending on relative interaction strength and obligacy. Benefit strength (weak or strong) is relative to intraspecific limitation. Arrows are vectors showing the "flow" of the system: arrow angle shows the direction of changes in density of N 1 (x-direction) and N 2 (y-direction) and arrow color shows the magnitudes of change in that direction (lighter colors are stronger changes). Nullclines are curves of zero change of density for one partner. Equilibria (colored or hollow dots) occur when both partners have zero change in density. Equilibria are locally stable (black dots) or unstable (red dots) if the system is attracted or repelled, respectively, the equilibrium after a small perturbation. Equilibria are half-stable "saddles" (hollow dots) if the system is attracted in some dimensions by repelled in others. Panels were generated using the model in Case 1.
Note: A full list of models cited in the main text is included in the supplementary information (Table S1). Equations largely follow the notation from the original citations. All parameters are positive (>0) unless otherwise specified. Models with unique mathematical forms are given unique equation numbers. We encourage the readers to refer to the original references for the model derivations and interpretation of parameters. Notes include inspiring system and obligacy, if specified by authors.
complexity as authors included more ecological mechanisms and system-specific realism (Table 2), leveraging numerical equation solvers. Accordingly, we use phase plane diagrams (Figures 1-4) to visualize the different qualitative dynamics of these models, as determined by species' curves of zero growth ("nullclines") and fixed points ("equilibria") of the system (summarized in Table 3). Gause and Witt (1935) proposed a model for "mutual aid" between a host and symbiont, inspired by Konstitzin (1934;Wolin, 1985).

This model was a modification of the Lotka-Volterra competition
equations with positive (instead of negative) interaction coefficients (Equation 1; see Table 2 for numbered equations). Benefits increased linearly with increasing partner density, while the strength of negative (intraspecific) density dependence arising from processes external to the mutualism also increased linearly with the density of the species receiving the benefit (i.e., the recipient species; Figure 1). In this formulation, mutualism has two effects: it increases the lowdensity growth rate of the recipient and the highest density at which the recipient can persist (typically, the equilibrium density). The second effect has been called an increase in "carrying capacity," but we reserve that term for density in the absence of the mutualistic F I G U R E 2 Characteristic dynamics for saturating benefit models. Density-dependent benefit functions stabilize linear benefit models ( Figure 1). Benefits may saturate (decrease in strength) with increasing recipient density ("intraspecific density-dependence," Case 2.1), increasing partner density ("interspecific density-dependence," Case 1.2), or both (Case 2.2), resulting in stable coexistence (see Table 3). Specifically, when paired with a partner with linear (a, b) or saturating (e, f) benefits, feasible systems exhibit the same qualitative dynamics: stable coexistence at densities higher than either partner could achieve alone (off-axes black point), and potential or guaranteed threshold effects when one or both partners are obligate mutualists. Under a certain threshold (red dashed line), one population is at too low density to support its partner, collapsing the system (b, f). This threshold causes extinction of obligate partners, even if initially highly abundant (e.g., follow lighter colored trajectories in panel f). These dynamics of coexistence and threshold effects are robust across models of mutualism with saturating benefits, regardless of the mechanism by which benefit saturates (Cases 1.2, 2.1, 2.2). Benefits may also increase in strength with increasing recipient density (also called "intraspecific density-dependence," Case 3.2), causing unbounded growth in the absence of other limitations. Specifically, feasible systems between two facultative partners of this form exhibit unstable coexistence (c, d) and a potential threshold under which the system exhibits stable coexistence at low density or explodes with unbounded population growth at high density (d). Panels were generated using models in Case 1.1.1 (N 1 only, a, b), Case 1.2 for (N 2 only a-b, both e-f), or Case 3.2 (c, d) of Table 3 (a) (b) Characteristic dynamics for shifting net effects and consumer-resource models. Models that investigated shifts in net effects as a balance of costs and benefits ("context-dependency") led to a synthesis of mutualism into a consumer-resource framework. Models with saturating benefit functions and linear costs (a-b) tend to display stable coexistence (a) and threshold effects (b) like earlier models ( Figure  2). Stable coexistence is "mutualistic" if the nullclines intersect such that both species achieve higher density than they would alone, or if increasing the density of one species from equilibrium permit growth of its partner. Otherwise, the interaction is "parasitic." Linear costs can make the coexistence equilibrium a stable spiral, with damped oscillations toward equilibrium (b, d, f, g). Models with unimodal benefit response that allow negative effects (net costs) at high density (c, d) or that include both separately saturating costs and benefits (e, f) display more complex dynamics. Depending on its parameterization, the mutualism-competition model by Zhang (2003) displays mutualistic stable coexistence (not shown), competitive exclusion (c), or competitive dominance (d), with dominant species dependent on initial densities (i.e., system initialized to the left or right of the separatrix). The consumer-resource model by Holland and DeAngelis (2010) also displays a range of dynamics depending on parameterization (e, f), including multiple stable coexistence equilibria (f). Mutualistic coexistence occurs when the ratio of consumers to their resources is not above a certain threshold (i.e., to the left of the left separatrix, or below the bottom separatrix). Otherwise, consumers overexploit their resources (causing more costs than provided benefits), leading to system collapse. Recent works use a consumer-resource approach with system-specific mechanisms (g, h), but often exhibit the simpler qualitative dynamics of saturating benefit models ( Figure 2) with the potential for oscillations (g). Panels show the following models: (a-b) Neuhauser (Vandermeer & Boucher, 1978), those that can persist at positive density ("carrying capacity", K) in the absence of their partner (K > 0). Gause and Witt also commented that increasing the strength of mutualism ( ij , Equation 1) increases both species' equilibrium biomass until they pass to infinity, but that infinite populations are obviously unreasonable and microcosm studies suggest that interaction strength should decrease as species grow.
These two studies (i.eGause & Witt, 1935;Kostitzin, 1934) initiated theoretical research on what we now call mutualism around the same time as theoretical research on predation and competition, but then paused for nearly 40 years.
Beginning in the 1970s, mutualism received attention as a destabilizing force in ecological networks represented as random community matrices (May, 1972(May, , 1973, with the unbounded growth in the Lotka-Volterra models of mutualism being called a "silly solution" (May, 1976). Using Lotka-Volterra models, authors better characterized the conditions that lead to unbounded growth found by Gause and Witt's original model of mutualism (Albrecht et al., 1974;Goh, 1979;Travis & Post, 1979;Vandermeer & Boucher, 1978).
Other forms of linear benefits were investigated such as those that increase per-capita growth rate, equilibrium density, or both ( Figure 1). Whittaker (1975) introduced a model in which mutualism increases the equilibrium density of one partner and both the equilibrium density and per-capita growth rate of the other partner.
This model accommodates "obligate" mutualists like symbionts living on a host that cannot persist in the absence of that host, that is, have zero carrying capacity (K = 0) in the absence of their partners. The mutualistic symbiont-host interaction linearly increases the carrying capacity for the symbiont (Equation 2) while benefiting the host population by increasing its low-density growth rate and its equilibrium density (Equation 1). Later, Addicott (1981) introduced a model in which mutualism only increases the per-capita growth rate (Equation 4), inspired by the ant-aphid mutualism described in Addicott (1979). Addicott emphasized that these different linear benefit models could be used in a mix-and-match style to accommodate different types of benefits exchanges. Vandermeer and Boucher (1978) proposed the groundbreaking idea that mutualistic partners may exist along continuums of obligacy and interaction strength. The authors defined facultative mutualists as those with positive carrying capacity in absence of their partner. Obligate mutualists were defined more abstractly with zero or negative carrying capacity in absence of their partner (K ≤ 0), which represents the demographic drawdown that mutualism must exceed to allow persistence of the population. Negative carrying capacity arises mathematically when a population has a negative "intrinsic" growth rate, as is the case when its per-capita death rate exceeds its per-capita birth rate (K i = r i ∕a ii < 0, where a ii > 0 is a self-limitation coefficient, Table 2). This choice is useful both mathematically and ecologically because it allows the strong demographic pulldown when death rates exceed birth rates to be represented, without introducing numerical issues due to zero carrying capacity. Vandermeer and Boucher's analysis of Gause and Witt (1935)'s model found that obligate-obligate partnerships would either F I G U R E 4 Distinguishing characteristic dynamics. N 1 (x-axis) is obligate mutualist and N 2 (y-axis) is facultative in all panels. (a) Threshold effects: N 1 goes extinct when the density of N 2 is below a threshold (separatrix). The system achieves stable coexistence when N 2 is above the threshold, and both species achieve higher densities than either would attain alone. (b) Overexploitation dynamics: the system collapses above a threshold in the ratio of consumer (N 2 ) to resource (N 1 ) species density. At low density, both partners will grow due to benefits from mutualism until they reach stable coexistence at higher density than either species could achieve alone. Above a threshold of N 2 density (separatrix), both populations will grow but N 2 will increase to such an extent that it exerts more costs than benefits it provides (exploitation). N 1 will begin to decline at low density while N 2 continues to grow, eventually leading to both going extinct. At even higher initial densities of N 2 , N 2 will immediately overexploit N 1 and both species will go extinct, without even acquiring enough benefits to allow its own population to grow. (c) Allee effects: N 1 will go extinct if its density is under a threshold of its own density (left side of N 1 non-trivial nullcline) because it becomes too rare to receive benefits from the mutualistic interaction. The system tends toward stable coexistence at higher density than either partner could achieve alone when N 1 is above a threshold of its own density (separatrix). Note that threshold effects induced by partner decline (a) cause Allee effects in both species because at low density they cannot support a sufficient partner population density to allow their own population growth. Overexploitation (b) by the high-density consumer (N 2 ) also induces an Allee effect in the resource species (N 1 ) where lower resource density causes lower benefits from the interaction.  S35: Accelerating negative densitydependence; "K-selected," sedentary, & stage-structured organisms, e.g., flowering plants (Moore et al., 2018) 1.2 …as a function that saturates with increasing partner density  Rio, 1981;Wright, 1989) S34: limited by rewards availability (Type I, on saturating plant rewards Revilla, 2015) S27: Mortality declines due to protection or deterrence by partners (Thompson et al., 2006) Case 2: Intraspecific density-dependence in mutualism only: benefits saturate with increasing recipient density Benefits accrue directly to per-capita growth rate… 2.1 ……with increasing recipient density Plant reproduction is a function of pollinator visitation… S7: Type II, on plants (Soberón & Martinez del Rio, 1981) S33: Type I, on saturating plant rewards (Revilla, 2015) Also see S3 ) S10: Plant reproduction is a function of pollinator visitation (Type II), limited by ovule availability (Wells, 1983) S11: Pollinators forage on plants (Type II), limited by search time (Wells, 1983) Also see S4, S31 Case 3: Benefits of mutualism reduce intraspecific density-dependence in population dynamics Benefits reduce negative density-dependence… 3.1 …via increasing carrying capacity as a linear function of partner density  (May, 1976;Whittaker, 1975) S28: Partners supply substrate or habitat, e.g., domatia for aphids (Thompson et al., 2006) Also see S12 …via decreasing self-limitation 3.2 ……as a linear function of partner density  (Wolin & Lawlor, 1984) 3.3 ……as a function that saturates with increasing recipient & partner density Note: Description of nullcline geometry, qualitative dynamics, and empirical assumptions under which seven generic models of mutualism may arise. In all models, benefits of mutualism are a function of partner density (N j ). All models also include a form of intraspecific density dependence, that is per-capita growth rate is dependent upon recipient density (N i ). To better interpret the historical literature, we categorize models into three cases of intraspecific density dependence (see text). Only Case 2 yields feasible dynamics in the absence of self-limitation (i.e., when s i = 0). Intrinsic (per-capita) growth rate determines obligacy in all models (r i is facultative), with one exception. Case 3.1 uses the (deprecated) historical convention in which carrying capacity directly determines i is facultative). All other parameters are assumed to be positive. Nullcline geometry is restricted to the ecologically relevant region (N 1 ≥0, N 2 ≥0). Only feasible dynamics are listed: "SC" is stable coexistence, "UC" is unstable coexistence," "UC/E threshold" is a threshold dividing the plane into unstable coexistence at higher density or extinction at lower density, "HD" is high density, etc. Alternative qualitative dynamics (listed on separate lines) are possible based on parameterization of the models.

TA B L E 3 (Continued)
collapse to extinction when benefits are weak or exhibit a threshold population size under which they go extinct and above which they grow unboundedly when benefits are strong (Figure 1b,d). They also found that facultative partners are likely to coexist stably when benefits are weak or exhibit unbounded growth when benefits are strong (Figure 1a,c, also see Wolin, 1985). Remarkably, Vandermeer and Boucher (1978; also see Christiansen & Fenchel, 1977) anticipated the qualitative dynamics generated by extending these models to saturating benefit responses. However, the authors emphasized that unbounded growth was still an ecologically relevant result because it indicates persistence of the two-species mutualistic system. Indeed, they argue that persistence (whether species persist or go extinct) is a more biologically useful metric than neighborhood stability (whether the system returns to equilibrium after a small perturbation). Subsequent authors also emphasized other properties of stability of mutualism such as return time to equilibrium (Addicott, 1981;Wolin, 1985), the domain of attraction to equilibrium (Benadi et al., 2013a), species persistence (Valdovinos et al., 2013(Valdovinos et al., , 2016(Valdovinos et al., , 2018, maintenance of diversity (Benadi et al., 2013b), and biomass variability (Hale et al., 2020).

| Saturating benefit models
The earliest models that incorporated saturating benefits within mutualism invoked unspecified (intraspecific) environmental constraints that limit population growth in the presence of a mutualist (Dean, 1983;May, 1976;Whittaker, 1975;Wolin & Lawlor, 1984).
For example, Whittaker (1975) assumed extrinsic, intraspecific limiting factors to the benefits a host could receive from its symbiont (Equation 3, Figure 2a). This is the first of many models that exhibit thresholds (sensu Vandermeer & Boucher, 1978), where the low density of one partner pushes the whole system to collapse (sometimes called "Allee thresholds," e.g., Johnson & Amarasekare, 2013).
This focus on extrinsic limits to benefit was epitomized by Wolin and Lawlor (1984). They derived models for five different ways in which mutualism could affect per-capita birth or death rates as functions of recipient density. For example, they compared models in which mutualism reduces intraspecific density-dependent limiting factors only in per-capita birth rates (Equation 6, Figure 2c Figure 1a). These models were classified as describing mutualisms with effects primarily at high versus low self-density.
Wolin and Lawlor concluded that low-density effects (i.e., primary effects on per-capita growth rate as opposed to equilibrium density) are stabilizing in terms of faster return times and the existence of a feasible, stable equilibrium. These models of "intraspecific densitydependence" (so-called by later authors, Holland, 2015) lacked biological mechanisms or reference to clear ecological examples, which perhaps pivoted the field away from this otherwise fruitful approach. In contrast, Soberón and Martinez del Rio (1981) proposed a detailed pollination model in which plant benefits are a function of pollinators' visitation rate, modeled as a saturating Type II functional response. Thus, benefits to plants saturate as a function of their own density (intraspecific density dependence), but due to factors intrinsic to the mutualism (that is, time constraints for pollinators handling flowers during foraging visits). Such an approach has seen a resurgence in recent literature (see Consumer-resource approach, below) but was largely abandoned at the time.
Starting in the late 1980s, authors began to focus on "interspecific density dependence," which has been considered more consistent with other theories of interspecific interactions (Holland, 2015). Wright (1989) proposed a model in which per-capita benefits saturate in terms of partner density analogously to consumers foraging on resources due to limitations of consumer handling of resources or uptake rate (Figure 2e,f). In the mutualistic case, benefits are assumed to saturate with increasing partner density, often as an additive, first-order term to per-capita growth rate following a Holling Type II functional response (Bazykin, 1997;Hale et al., 2021;Holland & DeAngelis, 2010;Thompson et al., 2006;Wright, 1989;Wu et al., 2019). On the other hand, Thompson et al. (2006) proposed a theoretical framework that organized both terrestrial and aquatic mutualisms into those that affect birth rate, death rate, habitat acquisition, or a combination of these benefits for each partner. Other authors have used different mathematical forms for analytical tractability (García-Algarra et al., 2014;Pierce & Young, 1986). Regardless, these assumptions result in both an increase in low-density growth rate and an increase in equilibrium density in the presence of mutualists.
These studies using the interspecific density dependence approach included more ecological justification for mechanisms that limited benefit accrual. However, phenomenological accounts of environmental conditions limiting population growth were still present with most models via an undiscussed intraspecific limitation term (see discussion by Johnson & Amarasekare, 2013). That is, authors assumed that at least one partner was limited by negative density dependence to ensure curved nullclines and stable coexistence in the mutualism (see Intraspecific density-dependence, below). Recently, Moore et al. (2018) introduced nonlinearities in intraspecific limitation while maintaining linear benefits (Table 3, Case

1.1.2-3). Mutualism is stable when density dependence accelerates
with increasing recipient density. Ecologically, this means that the growth rate of the population receiving the benefit decreases faster and faster at higher density, which has been observed empirically (Moore et al., 2018). This result highlights the importance of investigating the effect of more realism in intraspecific limitation on the dynamics of mutualism, which has been largely underexplored.
Other authors derived models with benefits limited by both inter-and intraspecific density dependence (Fishman & Hadany, 2010;Johnson & Amarasekare, 2013;May, 1976May, , 1978Wells, 1983; Table 3). This added complexity was usually justified by systemspecific considerations (e.g., May, 1976;Wells, 1983), but it also emerges from individual-level mechanisms in plant-pollinator systems (Fishman & Hadany, 2010) or intraspecific competition for food or services (Johnson & Amarasekare, 2013). In general, these limitations emerge when systems are limited both by availability of service providers (e.g., pollinators) and by the substrates that receive benefit (e.g., flowers to be pollinated, seeds to germinate, or individuals to protect from predators; Hale et al., 2021).

| Cost-benefit models and shifting net effects
Empirical work bloomed in the 1980s, revealing that mutualisms are not only more (omni)present than previously expected but also context-dependent (Bronstein, 1994;Chamberlain et al., 2014;Herre et al., 1999;Thompson, 1988). That is, the net effect of these interactions often shifts between mutualism and parasitism or competition due to the relative balance of costs and benefits of participating in the interaction (also called "context-dependency").
Moreover, costs and benefits themselves may be strongly varying across space, time, and other abiotic conditions. Early investigations of this topic used models that could accommodate different types of interactions through smooth transitions in parameter values (Gilpin et al., 1982;Pierce & Young, 1986;Whittaker, 1975). For example, Pierce and Young (1986) did not provide a specific mathematical form but used a geometric argument to investigate the dynamics of an ant-lycaenid butterfly interaction in which lycaenids may be mutualistic, commensalistic, or parasitic to tending ants. Neuhauser and Fargione (2004) explored the mutualismparasitism continuum using the classical predator-prey (or hostparasite) Lotka-Volterra model with the additional possibility of the parasite benefiting the host (Figure 3a,b). The model includes both benefits and costs, and it was applied to study plant-mycorrhizae interactions across gradients of soil fertility. The authors assumed that mycorrhizal fungi not only increase host-plant equilibrium density (benefits) but also linearly increase plant death rate due to exploitation (costs). This and other cost-benefit models can exhibit coexistence equilibria that are stable spirals, meaning that the population densities will oscillate toward a fixed point (see Patterns from Theory). Zhang (2003) also modified a Lotka-Volterra model to accommodate mutualism but chose the competition instead of the predator-prey version of the model (Figure 3c,d). The modified model assumed that the interaction between species was competitive at high density and mutualistic at low density, modeled phenomenologically as parabolic nullclines. This model can predict competitive exclusion, competitive coexistence where one partner dominates depending on initial density, thresholds in which low density of one partner drives the system to collapse, or "mutualism" according to the criterion that species coexistence stably at higher density than either could have achieved alone. Unfortunately, it is difficult to understand which of the diverse dynamics this model can exhibit are most ecologically relevant because interpretation is not provided for its parameters. A mechanistic derivation that achieves similar dynamics could be useful future work (but also see Gross, 2008 for a similar approach on an explicit resource).
Other models also described different outcomes depending upon relative species' density (Hernandez, 1998;Holland et al., 2002;Tonkyn, 1986;Wang, 2019). In an important advance, Holland et al. (2002) proposed a suite of models in which different net effects result from the difference between increasing benefit functions and linear, saturating, or decreasing cost functions (see Figure 1 of Holland et al., 2002). Their approach balances out different mechanisms that cause net effects of the interaction to shift as the relative densities of the populations change over time.
In seeking to represent the phenomena or mechanisms of shifting interaction outcomes, cost-benefit models revealed a much more complex set of potential dynamics for mutualism than had been previously reported. Saturating costs bend species' nullcline toward the partner's axis at high partner density, curving it back around toward the origin into a lobe shape (Figure 3c-f). This is because high partner density exerts high saturating costs on the recipient, which may exceed the benefits that can be acquired.
Up to five nontrivial equilibria occur when coexistence is feasible.
Moreover, separatrices running through saddle points define basins of attraction that lead to extinction or potential single-species persistence for facultative species. This ensures instability when one population is of substantially higher density than the other due to overexploitation of the rare partner ( Figure 4b). These dynamics contrast with the threshold effects ( Figure 4a) wherein the low-density partner benefits from mutualism but cannot provide sufficient reciprocal services. When the low-density partner becomes even rarer, it experiences an Allee effect, leading to its extinction ( Figure 4b). The high-density partner will also go extinct if it is obligate upon the low-density partner.
This much more complex set of potential dynamics that emerges from cost-benefit models exemplifies the criticism of mutualism theory as either too system-specific or too abstract to provide general insight into patterns and processes in mutualism (Bronstein, 2001a;Holland, 2015). Additionally, the field had not clearly connected the costs and benefits observed for individuals participating in a mutualism to potential population-level effects. The time was ripe for a conceptual synthesis.

| Consumer-resource approach to mutualistic interactions
In a landmark work, Holland and DeAngelis (2010) formalized a consumer-resource approach to mutualism, providing a bridge between mutualism and the ecology of other interspecific interactions. In their framework, mutualisms may be "unidirectional" or "bidirectional" consumer-resource interactions, in which one or both partners benefit from consuming costly resources provided by the other (Figure 4b, Figure 3e,f, respectively). Such framework accommodated the shifting net effects of previous models (Holland & DeAngelis, 2009, previous section), and formalized the concept of ecological costs and benefits as modifications to demographic rates due to resource provisioning and nutrient or service consumption. Notably, this framework allowed mutualisms to be modeled as a dynamic continuum along a spectrum of other interspecific interactions, such as predator-prey and competitive interactions (Holland, 2015;Holland & DeAngelis, 2009). This was possible by clarifying the "currency" of the effects of mutualism as energy or biomass exchanges that manifest in changes to percapita growth rate (or its components: birth, death, immigration, etc.). This framework stimulated recent development of theory for more specific systems (e.g., Kang et al., 2011;Martignoni et al., 2020).
Holland and DeAngelis (2010) modeled specific study cases similarly to previous studies (see Saturating benefits, above), but with costs defined separately from benefits via saturating interspecific functions, accrued through provisioning resources. In contrast, service-provisioning by consumers is assumed to incur only fixed costs that can be accounted for in parameter values, like increased handling time when foraging for resources. The nonlinear costs cause lobe-shaped nullclines allowing up to five coexistence equilibria. Like the earlier Zhang (2003) model, many dynamics are possible including mutualistic stable coexistence and oscillations. However, instead of the competitive exclusion and competitive coexistence outcomes of Zhang's model, "parasitism" by one partner is due to exploitation by a high-density partner that outweighs the benefits it provides to the lower density partner. In most dynamics of the Holland and DeAngelis model, parasitism collapses the system to extinction instead of allowing a stable but exploitative interaction like in Zhang's model.

Following Holland and DeAngelis' publication, authors began to
investigate accounting for resource dynamics in consumer-resource mutualisms more mechanistically. Resource dynamics were also considered in some earlier literature investigating mutualistic exchange of resources and between guild members sharing resources (bidirectional consumer resource), largely in the context of investigating coexistence mechanisms (e.g., Gross, 2008;McGill, 2005;Meyer et al., 1975). However, Benadi et al. (2012) and Valdovinos et al. (2013) proposed consumer-resource models for pollination networks (unidirectional consumer resource) in which consumption was on nectar "rewards" rather than individuals of the resource populations directly (but also see Scheuring, 1992 for a similar stage-structured model). These models separated the dynamics of the plants' vegetative biomass from the dynamics of the plants' floral rewards either implicitly (Benadi et al., 2012(Benadi et al., , 2013a or explicitly (Valdovinos et al., 2013). Explicitly separating vegetative and rewards dynamics introduces complexity but allows (1) Revilla, 2015;Wang, 2019) and community (Benadi et al., 2013b;Hale et al., 2020;Valdovinos et al., 2016) scales. For example, Revilla (2015) assumed rewards achieve steady state compared to changes in population density and derived models in which the linear consumption rate on rewards mediates benefits to the resource species. Hale et al. (2020) considered that pollinator visits can be approximated by consumption of floral rewards, and assumed that benefit to both plant and pollinator species is proportional to consumption rates on floral rewards. Hale et al. (2021) further specified whether benefits should be proportional to per-capita consumption rate (as may be the case for animal-dispersed plants) or to total consumption rate (as may be the case for animal-pollinated plants, which require obligate outcrossing). The latter leads to emergent Allee effects (Courchamp et al., 2018) for obligately animal-pollinated plants, explained by the plants' inability to attract pollinators at low density.

| PAT TERN S FROM THEORY
Historically, theory in mutualism has been focused on understanding how mutualisms can stably persist. Here, we broaden our scope to ask, what dynamics does the theory predict mutualisms will exhibit, and are they dependent upon ecological system or model formulation? We found that predictions for the population dynamics of mutualisms are qualitatively robust across the models reviewed, despite differences in level of detail, types of benefit, and inspiring systems.
We synthesize these general findings below.
We found that theoretical investigation of pairwise mutualism has repeatedly and robustly shown that mutualisms are stable.
The pattern of stable coexistence of mutualists at high density is robust across mechanisms that limit benefit (Figures 2 and 3, Table 3). Both inter-and intraspecific density dependence in saturating benefit functions lead to the same qualitative dynamics when they are present in at least one partner (also see Thresholds, below).
However, intraspecific density dependence and its effect on stability have been a source of confusion in the mutualism literature for decades.

| Intraspecific density dependence
We found that authors described their models as exhibiting intraspecific density dependence in three (not necessarily distinct) cases. In the first case, authors are referring to the negative density dependence term in a simple population dynamic model (Case 1 of Table 3). This term causes the decline in per-capita growth rate with increasing population density, and historically was modeled through a carrying capacity function (− N i ∕K i in Table 3). It is now typically modeled through a "self-limitation" term (− s i N i in Table 3), though it may represent any form of negative density dependence such as the Janzen-Connell effect, not just intraspecific competition for limited resources. To display a nullcline in the relevant ecological quadrant, it is necessary for mutualism models to include nonzero negative density dependence unless they include some other source of dependence on recipient density (e.g., s i can be zero in Case 2 of Table 3 because benefit saturates in terms of recipient density). Moore et al. (2018) found that one species having an accelerating negative density dependence term is also sufficient to allow stable coexistence if per-capita benefits accrue linearly (Case 1.1.2). However, the form of negative density dependence (accelerating, decelerating, or constant) does not typically affect nullcline geometry if per-capita benefits saturate (e.g., does not affect the qualitative dynamics of Cases 1.2, 2.1, 2.2 of Table 3).
In the second case, authors refer to intraspecific density dependence in their models when benefits from mutualism increase percapita growth rate directly (i.e., affect density-independent rates such as increased per-capita birth rate or decreased per-capita death rate), but benefits saturate with increasing recipient density (Case 2 of Table 3). This emerges when benefits are a function of the partner's visitation rate on the recipient or consumption rate on rewards provided by the recipient or when the recipient has limited substrate with which to convert interactions into benefits. This may generally be the case when mutualists provide reproductive or protective services (e.g., Soberón & Martinez del Rio, 1981, Thompson et al., 2006, Johnson & Amarasekare, 2013, but also see nutritional exchanges in Parker, 2001, Martignoni et al., 2020. In the third case, authors refer to intraspecific density dependence when benefits from mutualism reduce negative density dependence (Case 3 of Table 3) so that the effect of mutualism is most prominent at high recipient density (Wolin & Lawlor, 1984). Authors have chosen this approach when mutualists provision habitat (e.g., Thompson et al., 2006), reduce density-dependent mortality such as seed predation via the Janzen-Connell effect (e.g., Hale et al., 2021), or in the case of symbionts, which live within host populations (e.g., Whittaker, 1975). Here, benefits may be mediated through carrying capacity (Case 3.1) or through a self-limitation term (Cases 3.2, 3.3), with different resulting nullcline geometries.
Linear increases in carrying capacity or decreases in self-limitation rate can yield unbounded population growth (Table 3). More generally, even models with saturating benefits can exhibit unstable behavior when benefits accrue directly to a term that represents intraspecific density dependence, which decreases per-capita growth rate at high density (not shown). If mutualism decreases negative density dependence to such an extent that it induces positive density dependence at high partner density, the recipient population will begin accruing increasing benefit with its own increasing density (Case 3.2). Then, the system can display unbounded growth ( Figure 2c,d) unless benefits are additionally limited by extrinsic or intrinsic factors such as the number of seeds that can germinate after seed dispersal or the number of ovules that can be pollinated by pollinators (Case 3.3, Figure 3h).
Though all three of the above cases have been called "intraspecific density-dependence" in the mutualism literature, they refer to different ecological phenomena and have different implications for the dynamics of mutualism. All models must include some form of per-capita dependence on recipient density for feasible nullclines, but this may be manifest through a self-limitation term or through per-capita benefit functions that decrease with increasing recipient density. Models in which benefits reduce negative density dependence in a recipient population tend to allow unbounded population growth unless there are additional limits to benefits accrued. In contrast, models in which per-capita benefits saturate with increasing recipient density are stable, and exhibit the robust dynamics of high density stable coexistence and a low-density threshold observed in models with benefits that saturate with increasing partner density (i.e., interspecific density dependence).

| Mutualisms exhibit thresholds when at least one partner is obligate
Nearly all models that predict stable coexistence at high density also predict destabilizing thresholds at low density when one or more partners are obligate upon the mutualism (Figure 2a,b,e,f, Figure 3a,b,g,h). Specifically, if either species dips below a critical threshold in population density, the obligate partner(s) will go extinct, even if initially at high density (Figure 4a). This collapse occurs because, under the threshold, the low-density species cannot provide sufficient benefits to its higher density partner. Threshold effects occur in systems with interaction strengths high enough to allow feasible coexistence, but with per-capita growth rates small enough (very negative for obligate partners, near-zero for facultative partners) that a partner can potentially achieve densities low enough for long enough that its obligate partner will go extinct.
Understanding threshold dynamics provides rich insight into interaction strength, obligacy, and positive feedbacks in mutualistic interaction. By definition, obligate mutualists have negative per-capita growth rate in the absence of their partner. Thus, obligate mutualists can only be saved from population decline by benefits from mutualism that exceed their own negative intrinsic growth rate, that is, via strong mutualistic interactions. If both partners are initially at high enough density, obligate mutualists can achieve positive population growth, resulting in stable coexistence. However, if an obligate mutualist is at high density but its partner is at low density, the obligate mutualist will decline quickly due to both its negative intrinsic growth rate and strong intraspecific limitation at high density. The low-density partner may be growing due to mutualistic benefits, positive intrinsic growth, or release from intraspecific limitation. However, under the threshold, its population cannot recover fast enough to provide sufficient benefit to cancel out the negative intrinsic growth rate of the obligate partner and save it from decline. On the other hand, facultative partners can rely upon their own positive intrinsic growth rate to recover from low density, even after declines due to strong intraspecific competition or insufficient benefits provided by its partner. Thus, destabilizing threshold effects do not occur when both partners are facultative. However, highly nonlinear models can exhibit similar thresholds in facultative partnerships where coexistence occurs below the threshold at low, rather than high densities ("bistable coexistence," Hale et al., 2021;Parker, 2001).
Threshold dynamics emerge from the unique nature of mutualism and are potentially characteristic of this interaction. In predator-prey interactions, a low-density predator may benefit from a higher density prey population that is declining, but negative feedback in the system also limits the growth of the predator population at high density and subsequently allows the recovery of the prey population from low density. In competition interactions, the higher density partner exerts stronger and stronger negative effects on the rare population, causing the rarer population to go extinct if interspecific competition exceeds intraspecific competition for at least one of the competitors. In contrast, the positive feedback in the mutualistic system requires that both partners can provide sufficient benefits to the other to maintain the interaction.
Notably, threshold effects also occur in models that take very different approaches than those reviewed here. For example, Ingvarsson and Lundberg (1995) observed threshold effects dependent upon the ability for pollinators to find flowers in a modified disease model for mutualism, while Wang (2019) showed that the thresholds observed in Revilla's (2015) model more precisely occur between pollinator and rewards density rather than pollinator and plant density directly. This further emphasizes the potential generality of thresholds in mutualisms.

| Allee effects
Allee effects are a form of threshold where the population exhibits negative per-capita growth rate when rare. Here, we use "Allee effects" to refer specifically to strong, demographic Allee effects (Kramer et al., 2009) that emerge from the mutualism (i.e., are not hard coded into the population dynamics, Courchamp et al., 2018). Allee effects can emerge from many mechanisms, but we distinguish between a few proximal causes that suggest differing management recommendations for driving a collapsing system to high-density stable coexistence. The most obvious case is also the least common form of threshold observed in mutualism models: Allee effects driven by the inability of a population to support itself.
This type of Allee effect has also been observed in food chains that include protection mutualism (Morales et al., 2008) and in models of sequential colonization of patches by plants and mobile mutualists (Amarasekare, 2004). As mentioned above, Hale et al. (2021) find Allee effects in obligate plants when they become too rare to attract sufficient visitation from pollinators (Figure 4c). From a management perspective, it would be necessary to supplement the population experiencing the Allee effect (the declining, low-density partner) to prevent its extinction (Figure 4c). The partner-induced threshold described above also leads to Allee effects, wherein species decline when their partner is too low in density to support positive growth. In this case, it would also be necessary to supplement the low-density species, though it may already appear to be recovering due to positive population growth and high partner density.
Indeed, from a management perspective, this would achieve the counter-intuitive goal not of saving the low-density population, but rather its high-density partner from extinction ( Figure 4a). Finally, Holland and DeAngelis (2010) find Allee effects in animal populations induced by overexploitation from another consumer mutualist.
In this case, the management recommendation would be to equalize partners' population densities to avoid overexploitation (Figure 4b).

| Strong interactions are needed for obligate mutualists to persist
Research on mutualistic interactions has yet to firmly define interaction strength (Valdovinos, 2019). In Lotka-Volterra models, interaction strength is simply defined by the benefit coefficient (α ij in Equations 1, 2, 4). However, as authors have gained deeper mechanistic understanding of mutualism, it has become clear that interaction strength is a more complex topic related to the "effectiveness" of mutualistic partners (Schupp et al., 2017;Vázquez et al., 2015). Schupp et al. defined the effectiveness of a population for providing mutualistic benefits to its partner as the product of the "quantity" and "quality" of benefits provided. The term "quality" accounts for the species-specific and interaction-specific traits, as well as the environmental context that determine how much benefit a partner can receive from a unit of benefit "quantity". Examples of such benefit quality are the nutrition acquired from a foraging visit or the probability of a seed recruiting after being removed by a disperser.
The parameters that determine the quality of the mutualistic interaction are useful for understanding the criteria for stable coexistence and thresholds. Weak interactions between facultative partners in Lotka-Volterra models are considered stabilizing because they ensure stable coexistence instead of permitting unbounded growth. Specifically, mutual benefits must be weaker than species' intraspecific limitation (Gause & Witt, 1935;Travis & Post, 1979). However, stable coexistence always occurs between facultative mutualists in models with saturating nullclines regardless of interaction strength. Conversely, in saturating systems with at least one obligate partner, interactions must be sufficiently strong to overcome the negative intrinsic growth rate of the obligate partner for coexistence to be feasible (Bazykin, 1997).
In this case, destabilizing threshold effects can occur not because of interaction strength, but due to the low intrinsic growth rate of the partner. Overall, stronger interactions stabilize systems with threshold effects by decreasing the threshold in population density that causes the system to collapse, which allows positive growth from lower densities.

| Effects of mutualism varies between low and high population density
Empirical work has shown that the effects of mutualism vary with both recipient (Wolin & Lawlor, 1984) and partner density (Holland, 2015), and models show that this can lead to different ecological dynamics. When benefits are strongest at low recipient density, we can expect the robust dynamics of stable coexistence and threshold effects described previously (Figure 2). When benefits are strongest at high recipient density, models predict unbounded growth unless Early syntheses reported that mutualism with the strongest effects at high recipient density is less likely to be stable than that with the strongest effects at low recipient density (Addicott, 1981;Wolin, 1985). At that time, authors represented high-density effects of mutualism as direct modifications to species' carrying capacity (Equations 2, S9, S16; Wolin & Lawlor, 1984). Authors now represent the effects of mutualism exclusively through changes in demographic rates (Holland, 2015) unless explicitly representing habitat provisioning, for example, corals or plants with domatia and their animal partners (Thompson et al., 2006). Mutualism may still have the strongest effects at high density (e.g., if benefits reduce negative density dependence due to intraspecific competition or the Janzen-Connell effect), but this would be represented by modifying intraspecific limitation due to mutualism.
Categorizing mutualisms by their relative magnitude of costs and benefits at low versus high density of recipients versus partners is still a profitable approach that could lead to a next-generation theoretical framework that organizes mutualism by their population dynamics. Additionally, separating out the specific demographic rates affected by mutualistic interactions (as in Thompson et al., 2006 andHale et al., 2021) will likely clarify the differences and similarities between mutualisms. Even if the population dynamics of most models of mutualisms are qualitatively robust, the details of the low-density dynamics and the criteria for collapse can provide insight for system-specific mechanisms and patterns among them Wu et al., 2019).

| Costs of mutualism can cause damped and undamped oscillations
Models that incorporate costs to the mutualistic interaction can exhibit the same qualitative dynamics described above. That is, they are stable when incorporating limiting factors to benefits and self-limitation, exhibit thresholds when at least one partner is obligate, and need strong interactions for obligate partners to persist.
Additionally, these models can produce oscillations. Linear costs can result in damped oscillations when the equilibrium is a stable spiral (Figure 3b,g; Kang et al., 2011;Neuhauser & Fargione, 2004). Undamped oscillations occur when overexploitation by the consumer causes an Allee effect in the resource, which does not necessarily lead to extinction (Figure 3f). After depleting their resource population, the consumer population also declines, eventually allowing the resource to receive sufficient benefit compared to losses due to consumption. The system thus recovers, and coexistence is maintained in this region via a limit cycle (i.e., oscillations) around a stable center (left-most stable equilibrium, Figure 3f). This outcome is not seen in simpler models without cost terms, which predict stable coexistence at a nonoscillatory node (Figure 2), or with linear cost terms, which can predict damped oscillations in a stable spiral (Figure 3b,g).
Note that oscillation has been considered an important dynamic for mutualism models to reiterate, as justified by observations that mutualist populations can vary in space and time (Holland, 2015).
However, such variability need not necessarily be driven by the underlying population dynamics. Far simpler models of mutualism can produce oscillations when accounting for discrete time dynamics (e.g., Gilpin et al., 1982). Additionally, population oscillations observed in nature may be caused by external factors, such as environmental variation. This emphasizes that introducing explicit cost terms into mutualism should be adequately justified at the population level. Regardless, the models in question suggest that oscillations can be induced predictably, for example, by decreasing the density-dependent mortality of an obligate symbiont (Neuhauser & Fargione, 2004, e in Equation. 9), which could potentially be tested empirically by using different fungal strains in a plant-mycorrhizal system (Martignoni et al., 2021).
Ecological theory of mutualism has been criticized as sparse, largely consisting of models that are either too abstract to be useful or too case-specific to reveal general patterns (Bronstein, 2015a). This is an accurate description of many of the models we reviewed; however, remarkably, nearly all these models conformed to the same dynamics. We found that many historical models make similar qualitative predictions despite their different derivations, mechanisms, and inspiring systems. When feasible, coexistence is stable, and populations grow with bound. Mutualisms with at least one obligate partner exhibit thresholds, under which the low density of one partner destabilizes the system. If a species sustains nonlinear, population-level costs from mutualism, it may be overexploited to extinction by its partner. These patterns suggest that there exists a robust population dynamic theory of mutualism that can make general predictions. With this groundwork of theory laid, authors can now focus on how relaxing the assumptions of current models affects their predictions. For example, spatial and transmission models reiterate the threshold predictions of models that conform to the mean-field assumption (Ingvarsson & Lundberg, 1995;Mohammed et al., 2018) as do models with explicit rewards dynamics compared to those that approximate rewards dynamics as at steady-state (Revilla, 2015;Wang, 2019).

| Avenues for future research
Future work should understand how predictions from pairwise models scale to the network level. Threshold effects only occur when at least one partner is an obligate mutualist. Most species have multiple potential partners and thus are not truly "obligate" in the sense that only a specific pairwise interaction can allow positive population growth. Instead, most mutualists are likely to be facultative, engaging in diffuse interactions with many potential partners. However, it is likely that mortality exceeds reproduction in the absence of mutualistic interactions for many species. In this sense, species may be obligate mutualists even though they have multiple partners. Additionally, species are likely to have critical (cumulative) thresholds to allow population growth. For example, Valdovinos and Marsland (2021) identify the quality of visits needed from pollinators for plants to persist. Below such threshold, the plant species and the animals depending on those plants go extinct. Understanding how destabilizing thresholds may emerge or be ameliorated due to obligate mutualists in a network setting is an important goal for future work. Moreover, emphasis on consumer-resource approaches with a common "currency" of energy or biomass flows (Holland, 2015) makes mutualisms amenable to integration into interspecific network models such as food webs (e.g., Hale et al., 2020). Such integration can illuminate how context mediates interaction outcomes between potential mutualists, for example by shifting interactions into overexploitation or competition regimes. Indeed, understanding the structure and dynamics of these "multiplex" ecological networks that include multiple types of interactions has been identified as a primary goal in ecology (Kéfi et al., 2012).
Future work should also interrogate the assumptions and predictions of these models empirically. A main assumption is that mutualisms have population-level impacts. However, most empirical studies quantify the benefits and costs of mutualisms at the individual level in terms of fitness or even by using a single proxy for fitness (Bronstein, 2001a;Ford et al., 2015). Those effects do not necessarily imply population-and community-level impacts of mutualism (Flatt & Weisser, 2000;Ford et al., 2015;Palmer et al., 2010;Williamson, 1972). Therefore, empirical work is of foremost importance to evaluate whether mutualisms affect the population dynamics of mutualistic partners. Among the predictions of these models (stable coexistence, threshold effects, overexploitation), threshold effects have received the most attention (Latty & Dakos, 2019), but more empirical work is still needed. Wotton and Kelly (2011) and Kang et al. (2011) observed threshold effects directly in frugivory systems and in ant-fungal gardens, respectively, although the authors did not identify their results as such. Hale et al. (2021) showed that threshold effects in obligate plants may be swamped out by Allee effects (e.g., Forsyth, 2003), which suggests that targeted experiments to explore population trajectories should consider the criteria for observing different dynamics (Figure 4).
One difficulty of empirical applications is that an out-of-thebox consumer-resource approach following Holland and DeAngelis' (2010) framework can be logistically overwhelming. Nonlinear cost and benefit functions generate so many dynamics that they are nearly intractable analytically (but see numerical toolkit by Wu et al., 2019). Moreover, with up to four separate functional responses to parameterize, this framework requires an extremely high number of parameters to estimate empirically. This level of detail may be necessary to describe some two-species mutualism but is likely not general. Simplifications like approximating costs and benefits as proportional to consumers' foraging rate Revilla, 2015;Soberón & Martinez del Rio, 1981) can facilitate integration between theoretical and empirical approaches. Additionally, costs that scale with rewards construction can be approximated as fixed reductions to benefit, and thus accounted for in the measured parameters Revilla, 2015;Figure 3h). Systems with these complementary saturating benefits and fixed costs are likely to display much more limited dynamics than those shown in Figure 3c-f. For example, Kang et al. (2011) and Martignoni et al. (2020), Martignoni et al. (2021) adapted Holland and DeAngelis' approach to specific empirical systems, leading to models that predict the threshold and stable coexistence dynamics of simpler saturating benefit models (Figure 3g).
Reviewers for an earlier version of this manuscript commented that our results cement the idea that pairwise models of mutualism have been "pushed…as far as they will go," that "this literature has limited usefulness for motivating the theory of the future," and that it may be "the nature of mutualism" that its dynamics are "not very interesting…for a broad audience in ecology and evolution." Though we cannot speak to whether mutualism is of interest to specific individuals, we do believe that this attitude may have contributed to the long-term stagnation and repeated loss and rediscovery of theory in mutualism. A clear summary of the population dynamics of pairwise mutualisms (as we presented here) is an important groundwork for directing research into modules and networks including mutualistic interactions, the evolutionary origins of mutualism, and, pressingly, directing conservation efforts across systems (Figure 4). Both within the discipline and more broadly, there is an impression that theory is lacking. But it is simply not the case that ecological theory of mutualism is incoherent or underdeveloped: we find here that it is remarkably self-consistent despite the diversity of inspiring systems and modeling frameworks. It is not a mystery how pairwise mutualisms can persist stably, at least theoretically. Mutualisms are highly stable at high density, and the network setting may diffuse the risk of low density-thresholds leading to population collapse. A similar set of empirical literature to support or dispute the models' results has yet to accumulate, but we hope that by clearly outlining dynamical expectations of mutualistic theory, such work will be more accessible to empiricists.