Unifying mutualism diversity for interpretation and prediction

Coarse-grained rules are widely used in chemistry, physics and engineering. In biology, however, such rules are less common and under-appreciated. This gap can be attributed to the difficulty in establishing general rules to encompass the immense diversity and complexity of biological systems. Even when a rule is established, it is often challenging to map it to mechanistic details and to quantify these details. We here address these challenges on a study of mutualism, an essential type of ecological interaction in nature. Using an appropriate level of abstraction, we deduced a general rule that predicts the outcomes of mutualistic systems, including coexistence and productivity. We further developed a standardized calibration procedure to apply the rule to mutualistic systems without the need to fully elucidate or characterize their mechanistic underpinnings. Our approach consistently provides explanatory and predictive power with various simulated and experimental mutualistic systems. Our strategy can pave the way for establishing and implementing other simple rules for biological systems.


3
How populations interact to generate various outcomes is a key question in biology. Mutualism, where two or more populations provide reciprocal benefit, is an essential type of ecological interaction (1).
In marine ecosystems, coral reefs are based on mutualistic interactions between coral and algae that provide ecosystem services for humans and habitats for diverse organisms (2). Plant-bacterial mutualism is estimated to generate 60% of the annual terrestrial nitrogen input (3). Cross-feeding in microbial communities is also a mutualistic interaction that influences community structures and is the cornerstone of various microbial metabolic tasks (4,5). Although mutualistic coexistence is beneficial in maintaining the biodiversity, function and stability of ecosystems, under some conditions mutualistic systems can collapse, where one or more mutualistic partners is lost. This perturbation can further trigger the extinction or invasion of other populations and alter ecosystem functions (6)(7)(8)(9). A framework to interpret and predict mutualistic outcomes is useful to prevent undesirable system behaviors and provide guidance for modulating and engineering synthetic mutualistic systems.
Quantitative rules have been developed to elevate our understanding and provide predictive power for various biological systems (10)(11)(12)(13). However, such a framework is not yet available for determining mutualism outcomes. Main barriers in developing such a framework are the diversity of mutualistic interaction mechanisms and the complexity of underlying dynamics. Indeed, even engineered mutualistic systems that are by-design capable of cooperation, may not coexist. For example, it is still difficult to predict a priori whether an engineered microbial auxotrophic pair can persist or not (14)(15)(16). Previously, theoretical criteria in the form of model parameter inequalities have been developed for specific mutualistic systems such as cross-feeding mutualisms (17), plant-pollinator mutualisms (18), seed-dispersal mutualisms (19), ant-plant mutualisms (20), and plant-mycorrhizal mutualisms (21). These criteria depend on the underlying mechanisms assumed in the models and are not applicable to other types of mutualistic systems. General criteria have been developed (22-25), such as the classic criterion which states that intraspecific competition must be greater than mutual benefit for a mutualistic system to be stable (26). However, these criteria usually describe transitions between stable coexistence and unbounded growth, and fail to address the transitions between coexistence and collapse and other mutualism dynamics (27, 28) (Fig. S1, Supplementary Text section I).

Abstraction reveals a general rule
To reveal any commonality of mutualistic systems, we first summarized the logic of mutualism ( Fig. 1a). Mutualism can be defined as the collective action of two or more populations, where each population produces benefit ( ) that reduces the other's stress ( ) at a cost ( ) to itself. , , and are universal features of mutualistic systems (Table 1). In addition to benefit and cost, which are conventionally considered as the driving forces of mutualistic outcomes (29-31), it has been shown that stress factors such as nutrient limitations (32), rising temperature (33, 34), rising CO2 levels (35), invasive species (36) are also determining factors of mutualistic outcomes. Although stress is not always explicitly acknowledged in previous models, it has been implemented as decrease in growth rate (37), increase in linear turnover rate (32,(38)(39)(40)(41), or increase in intraspecific competition (22,25). Here we use stress as an umbrella term to capture the effects of stress factors in reducing baseline fitness of individual populations to provide a more complete picture of mutualism (see section II.A of Supplementary Text for the detailed reasoning).
To reflect the diversity of natural mutualistic systems, we systematically generated a total of 52 ordinary differential equation models based on the basic logic of mutualism with various implementation details (see Materials and Methods and section II.B of Supplementary Text for model assumptions). These implementation details are designed to comprehensively cover the various common and plausible forms of kinetic models that have been adopted in previous studies (see section II.C of Supplementary Text for a summary). Specifically, the models all revolve around the logistic growth equation but differ in the locations of , , and and in the complexities the models capture, such as competition, partner-density- 9 We cocultured the two strains starting from the same initial density with different concentrations of IPTG and aTc, which are the two dimensions of . The outcomes of coexistence and collapse are evident in the bimodal distribution of optical density (OD) at 32 hours of culturing (Fig. S6c). can be quantified by treating monocultures with the same set of [IPTG] and [aTc]. We used for M2 since it has a wider dynamic range (Fig. S6d). Using these data (Fig. 3b), we obtained a calibrated ( ). The confidence of ( ) is evaluated by the consistency of the top 5 ( ) and relative standard deviation of each ( ) (Fig.   S6e). Consistent with the circuit logic, ( ) increases with increasing [aTc]. The calibration reveals that [IPTG] also modulates ( ), which indicates unintended system complexities, such as QS cross-talk and unequal fitness of the two populations (Fig. 3c). We used cross-validation to evaluate how well new observations can be predicted. We found that ( )⁄ provides a cross validation accuracy of 96.8% for coexistence versus collapse and it is also predictive of total final density (Fig. 3d).
We then applied our procedure to data on a pair of Saccharomyces cerevisiae auxotrophs that is previously published (32). In this system, one strain cannot produce tryptophan (Trp) and the other cannot produce leucine (Leu). The mutualistic interaction of this system is realized by the exchange of the two amino acids in cocultures (Fig. 3e). Because [Leu] is maintained as 8 times of [Trp], we use [Trp] as one of the dimensions of to represent overall concentration of supplemented amino acid. The ratio of initial densities of the two strains composes the other dimension of (Fig. 3f, Fig. S7a). All top 5 ( ) reveal that intermediate ratios of initial density and increasing amino acid concentrations elevate ( ) (Fig. 3g). However, at the highest level of supplemented amino acid ([Trp] = 16nM, [Leu] = 128nM), top-ranked ( ) have qualitatively different trends, indicating a low confidence of ( ) at high concentration (Fig.   S7b). Interestingly, this high variability coincides with the system transitioning into competitive exclusion (32). Nevertheless, ( )⁄ is still predictive of final densities (Fig. 3h). Furthermore, we explored using the concentration of supplemented amino acid as a single system variable. ( )⁄ in this case can also predict the probability of coexistence as the ratio of initial densities varied (Fig. S7c).
In the third example, we applied our framework to previously published measurements of 14 engineered auxotrophic E. coli strains that compose 91 pairwise mutualistic systems (48) (Fig. 3i). The genetic context of the two partners vary while the growth environment is kept the same. The classification of coexistence versus collapse is based on the bimodal distribution of total density (Fig. 8a). of each auxotroph is determined based on final cell densities of monocultures when supplemented with different concentrations of its corresponding amino acid (Fig. S8b). We sorted the auxotrophs by the number of partners they coexist with to convert categorical indices into an ordinal scale. Thus, is composed of ordinal rankings of the two strains and measurements of coexistence versus collapse and are both arranged accordingly (Fig. 3j). We used strain 1 as the reference strain for the calibration. The calibrated ( ) generated a cross validation accuracy of 91.8% and we verified that ( )/δ is predictive of final total density ( Fig. 3k-l, Fig. S8c). We noticed a relatively high level of variability of total density when ( )/δ > 1, which can be due to system-specific properties that are not fully accounted for by mutualistic interactions.

Generalization of the framework
In nature, mutualism can occur among three or more partners (49). Thus, we tested our framework with simulations and experimental measurements of N-mutualist systems. Here we show the calibration procedure with simulated data from a 5-mutualist system and found that the quality of the calibration results is well-maintained ( Fig. 4a; Fig. S9a; Supplementary Text section VI.A). The same 14 auxotrophs (48) presented above were also used to construct all possible three-member systems. Using the same procedure with a 3-dimensional , where each dimension represents the genetic context of one strain, we found the predictor ( )⁄ provides an 89.3% cross validation accuracy and remains predictive of the total density, which indicates the scalability of the framework in experimental settings (Fig. 4b, Supplementary Text section VI.B). We hypothesized that ( ) calibrated for pairwise interactions can be used to directly construct a metric for 3-member systems, since theoretical analysis shows that n-member ( ) can be approximated by pairwise ( ) (Table S2). We assume of a 3-member system is the average of for all 3 combinations of its underlying 2-population systems and the same is true for . The constructed ⁄ for 3-member systems explains 80.8% of system outcomes (Fig. S9b). This result suggests the possibility of directly extending and from simple systems to more complex systems without further calibration.
Beside static environments, mutualistic systems can also inhabit dynamic environments where they experience fluctuating physical and chemical cues or cohabitate with other populations. We verified that the theoretical criterion generally holds in both cases (Fig. S10a). However, the transition between collapse and coexistence does not strictly occur at 1, which further advocates for the necessity of the calibration process. With simulated data, we carried out the calibration procedure and verified that the applicability of our framework is well-maintained ( Fig. 4c-d, Fig. S10b, Supplementary Text section VI.C-D). The robustness of the framework suggests that it can be used to study microbial communities, of which advancements in both interpretation and prediction are in demand (50).

Discussion
The immense complexity and diversity of biological systems is intriguing and inspires the exploration of mechanistic details. However, these details can distract us from simple rules that emerge at a higher level. By abstracting away from low-level details, many simple rules for biological systems have been developed to enhance our understanding and provide predictive power (51,52). A classic example is the Hamilton's rule, which states that a cooperative trait will persist if / > , where is the relatedness of the recipient and the actor; is the benefit gained by the recipient; and is the cost to the actor. More recent examples include linear correlations underlying cell-size homeostasis in bacteria (53)(54)(55), ranking of quorum sensing modules according to their sensing potential (56,57), and the growth laws resulting from dynamic partitioning of intracellular resources (58,59). Beyond establishing another simple rule, we also demonstrated that one can purposefully seek an appropriate abstraction level where a simple unifying rule emerges over system diversity. If this rule anchors in the basic definition of a type of system, it can then be applied to diverse systems of the same type. Beyond microbial systems that we tested, our criterion in principle can also be applied to other systems of larger or smaller scales that share the same logic.
Although simple general rules in biology are powerful tools, their applicability to experimental systems can be limited by the difficulties in associating the abstracted parameters to lower-level mechanistic details and quantifying these details experimentally. This is evident in the application of Hamilton's rule to, for example, to experimental systems (43-45). For simple rules that are dictated by inequalities, our calibration procedure provides a tool to apply these rules directly to experimental systems. If one side of the inequality and some final outcomes can be measured, the other side can be calibrated as an empirical function. Although the procedure cannot further dissect the empirical function into specific mechanistic parameters, the function can serve as an overall summary of the underlying mechanistic details while bypassing the requirement of characterizing them individually. Our approach thus enables the downstream interpretation and prediction by these simple rules with readily-accessible measurements.

Benefit
Cost Stress   (b) Models originating from the basic logic of mutualism yield diverse coexistence criteria. Each line represents the generation of a model from the basic mutualism logic and branching represents different implementation details and system complexities. The circles represent the models and different colors represent diverse coexistence criteria derived from these models. This process aims to reflect the diversity of mutualistic systems in nature.
(c) A simple rule emerges at an appropriate level of abstraction. The lines represent the abstraction process that establishes ( ) > as the common structure shared by diverse criteria in panel b. ( ) represents effective benefit and is a complex function of model parameters , which include and . In general, increases with increasing and decreasing . The heatmap is generated using the following model: is the stress experienced by one population in the absence of its partner.
(d) Intuitive interpretation of the simple rule. The effective benefit must overcome stress for the system to coexist. Solid black line represents coexistence boundary and dashed black line represents baseline fitness level in the absence of partner. Blue represents a that is greater than (coexistence) and yellow represents a that is smaller than (collapse).
(e) ⁄ can predict various system outcomes. If the two features and are known, many downstream predictions, both qualitative and quantitative, can be made.
(f) Quantitative outcomes versus ⁄ . Simulation results show when ⁄ > 1, it is predictive of total density.
Note that the points do not necessarily lie on a single curve, but a positive trend is well-maintained. Other quantitative outcomes also follow similar positive trends when plotted against ⁄ . (a) The rationale behind the calibration procedure. Conventional approaches (denoted by dashed gray arrows) require quantifications of mechanistic parameters as functions of contextual variables ( ) and finding the appropriate structure of ( ) to construct ( ( )). However, both steps are challenging and require case-bycase applications due to the diversity of mutualistic interactions. Instead, using qualitative outcomes of the system, we can distill an empirical function ( ) to approximate the true ( ( )). ( )/ can then predict qualitative and quantitative outcomes. Dark blue indicates the data that are relatively easy to measure without requiring mechanistic understanding of the interaction.
(b) A schematic demonstrating the mathematical basis of the calibration procedure. 1 and 2 represent two system variables. A circle represents an observation at a particular = ( 1 , 2 ). 5 observations are shown. contains qualitative system outcomes for each observation. Closed circles indicate coexistence and open circles indicate collapse; the same notation scheme is used for all following figures. contains the measurement of stress for each observation (lighter colors indicate higher values). Using , and , a boundary that separates the two types of outcomes can be established (the red curve). According to our simple rule, = is true at the boundary; > is true for coexistence and < is true for collapse. Using these data and our simple rule, we can calibrate for a ( ) which then enables the interpretation and prediction of system outcomes. Refer to Movie S1 for a 3D visualization of the calibration.

(l)
( )/ is predictive of the normalized fold change of final total density relative to initial density.  (c) A system that is modulated by an oscillatory signal. The oscillatory signal S is described by I (intensity) and L (duration). The signal temporally modulates and . I and L are the two system variables used in calibration. The procedure achieves a prediction accuracy of 97.3% for new data.

Materials and Methods Supplementary Text
Tables S1 to S2 Figs. S1 to S10 Captions for Movies S1 to S4

Other Supplementary Materials for this manuscript include the following:
Movies S1 to S4 MATLAB code package

Model development.
We built mutualism models based on four key assumptions: a) Benefit shall increase growth rate or carrying capacity and is positively dependent on partner density. b) Cost shall decrease growth rate or carrying capacity. c) Stress shall produce negative growth of populations at some parameter combinations. d) Negative growth of a population shall be potentially counteracted by benefit provided by a partner, but further strengthened by cost.
See SI section II for detailed reasoning and implementation of each assumption.

Criteria derivation.
We calculated the analytical solutions of fixed points of the 52 models using MATLAB R2017a.
Then we identified the fixed points that represent stable coexistence. The coexistence criteria are derived by ensuring the fixed points are real positive numbers. We can then rearrange the inequality to have on one side. The other side of the inequality is then an expression of other parameters, which is expressed as ( ). More details are presented in SI section III. The MATLAB code of the models and the derivation and testing process is included in the Supplementary Materials.

Calibration of benefit landscape using SVM.
We used support vector machine (SVM) algorithms in MATLAB to implement the calibration. The input data are formulated into the following format: We designed kernels that have additive separability between and , which can be expressed in a general form: where is the kernel that dictates the shape of the empirical function of and is a constant that varies the weight of . The predictor trained using SVM can then be expressed as: According to our criterion, we know that = when ([ , ]) = 0. We can then derive from Eq. 6 a function of in terms of system variables : With each set of inputs, we used linear, quadratic, cubic, and sigmoidal kernels with different kernel parameters to train ( ). We then find the ( ) with lowest cross-validation classification loss and bootstrapped variance. The top ranked ( ) is then used along with measurements for interpretation and prediction. See section V of SI for detailed calibration method and the Supplementary Materials for the MATLAB code of the calibration process and sample data sets.

QS-based mutualism strains.
The two strains were constructed based on circuit components from a synthetic predator-prey system (66, 67). Both populations carry two plasmids. Briefly, M1 carries plasmids identical to the predator plasmids, denoted A1 for the module carrying ccdA (tet promoter (68) driving luxR and lasI followed by lux promoter driving ccdA) and B1 for the module carrying ccdB (Lac promoter (68) upstream of ccdB followed by tet promoter upstream of gfp). To construct M2, A1 was used as backbone. To obtain orthogonal communication, KpnI and NotI restriction digest cloning was used to replace luxR/lasI genes from A1 with lasR/luxI genes from the previously published prey plasmid (consisting of pLac lasRluxI CcdB (Kan R , p15A ori)). Reporter plasmid B1 is from (66). To construct B2, enzymes XhoI and KpnI were used to replace the tet promoter on prey plasmid with the ccdB module from B1. All M1 and M2 plasmids were verified using restriction digest and sequencing.

Growth conditions of QS-based synthetic system.
The experiments of QS-based mutualistic system were done in 96-well microtiter plates. PHbuffered M9 medium (M9 salt supplemented with 1mM thiamine, 0.2% casamino acid, 0.4% glucose, 2mM MgSO4, 0.1mM CaCl2 and buffered with 100mM MOPS with PH adjusted to 7.0) was used. 50 g/mL kanamycin and 100 g/mL chloramphenicol were added to the culture to maintain plasmids.
To measure circuit function, 4 ml LB media in a 14 ml culture tube was inoculated from single colony and incubated overnight at 37 o C at 250 r.p.m. The optical density is adjusted to 0.5 in M9 media

I. Previous mutualism models
Here we briefly summarize some limitations of previously published models in studying the transition between mutualism coexistence and collapse.

A. The Lotka-Volterra model
Consider the following model that follows the basic form of Lotka-Volterra model of mutualism. This formulation is a nondimensionalized form of previous models (69)(70)(71)(72)(73)(74). The parameters are renamed according to our parameter assignments.
Although the model formulation captures the logic of mutualism, it can generate unbounded growth of the two partners, which is not a biologically relevant state. For example, when = 1, 1 = 2 = 2 and 1 = A coexistence criterion for mutualism can be derived using L-V model formulation and is previously demonstrated (69)(70)(71)(72)(73): However, this criterion captures the transition between stable coexistence and unbounded growth. It also suggests that mutualism is destabilizing, since increasing the strength of mutualistic interaction (β 1 β 2 ) tends to violate the above condition (I. 3). The interpretation of this criterion can be contradictory to empirical observations that mutualism stabilizes community (75)(76)(77). Although Lotka-Volterra models are sufficient to answer many questions related to mutualistic systems, it has been proposed that this discrepancy between model dynamics and empirical observations can be attributed to its unrealistic assumptions (78-81).

B. Other variants
We found that general mutualism models that generate unbounded growth usually do not capture the transition between coexistence and collapse. For example, the following three models were established previously to find structurally stable mutualistic models (78). In contrast to the L-V model which implements the interaction as a linear term, the following models present three alternative ways of adding the interaction to the basic logistic growth equation.
To avoid generating unbounded growth, one strategy is to introduce saturating benefit (73). However, although preventing unbounded growth, these models may still not generate negative growth which can be potentially countered by the increase of partner density. For example, the following model is adapted from a previous work (81) which falls into this category. In addition, if no decreasing density is captured, this model will stabilize at its coexistence state with any positive initial density for both populations.

II. Building various mutualism models
The diversity of mutualistic systems impedes the formulation of a general mutualism model and the general proof of a single criterion. We generated various mutualism models to reflect the diversity of mutualistic systems in nature and examine how diverse the coexistence criteria are. We also aim to investigate whether there exists an invariant form that is preserved regardless of specific model implementations. If such an invariant form exists, it will then reveal a fundamental coexistence criterion that is originated from the basic logic of mutualism.

A. Incorporation of stress
We explicitly define stress as a reduction of growth rate or productivity of biological systems, which is consistent with previous works (85)(86)(87)(88). Stress is universal in biology because it is present whenever growth rate is lower than the optimal growth rate. In mutualistic systems, many studies have shown that stress factors alter the basal fitness of individual mutualists. These factors include nutrient limitation (89), rising temperature (90,91), rising CO2 levels (92), invasive species (93), etc. Beyond benefit and cost, it is known that these stress factors also determine mutualistic outcomes (94).
To study quantitatively how stress affects mutualistic outcomes, we include stress as a model parameter that reduces the growth rate or carrying capacity of individual mutualists in various biotic/abiotic contexts.
In previous models, stress has been included as linear turn-over rate (66,89,(95)(96)(97) or intraspecific competition (69,70). Although these studies have more specific terms describing the downward pressure they capture, we use "stress" as an appropriate umbrella term.
Stress, thus defined, plays an essential qualitative role in mutualism. First, an absence of stress would mean that populations operate at maximum fitness. Such populations could not benefit from mutualism, because they would already be operating at their optimal level. Second, in mutualistic models the downward pressure can resolve the unrealistic exponential growth of mutualistic partners by imposing an upper limit to mutual benefit (78). Third, using stress can dissect the baseline fitness of individual mutualists caused by biotic/abiotic factors from the effect of mutualistic interaction on fitness.

B. Model assumptions
The key step of generating various mutualism models is to establish the set of assumptions the models should follow. We first start off with the most apparent aspects of mutualism: benefit and cost. These two aspects lead to two assumptions: e) Benefit shall increase growth rate or carrying capacity and is positively dependent on partner density. f) Cost shall decrease growth rate or carrying capacity.
To study the transition between coexistence and collapse in mutualism systems, the models must be able to simulate collapse, leading to the third assumption where we explicitly introduce stress to achieve negative growth. In addition to generating negative growth, stress has a physical meaning which is the difference between maximal fitness and baseline fitness in the absence of the partner: g) Stress shall produce negative growth of populations with some parameter combinations.
The fourth assumption follows assumption c) to reinforce the effects of benefit and cost in mutualistic interaction: h) Negative growth of a population shall be potentially counteracted by benefit provided by a partner, but is further strengthened by cost.
Even when all the above assumptions are satisfied, we still need to verify that the models do not generate unbounded growth with any parameter set.

C. Modifications of the logistic growth equation
After establishing a minimal set of assumptions, we then need to establish a systematic way of generating diverse models.
Mathematically, there are infinite possible implementations of a mutualistic system. We attempt to cover the different common and plausible forms of kinetic models that have been adopted in previous studies. Previous models have captured benefit, cost and stress in many ways. The following is a short summary of how benefit, stress and cost are modeled previously that serve as building blocks for our own models. Many models use Hill equations to capture saturating benefit (66, 81,89,98,99). Cost is also an aspect that is widely modeled, which can be implemented many ways (81,100), such as independent or dependent of partner density. Although stress is not often generally discussed, one possible form of stress is selfregulation, also called density dependence or inter-population competition, which often appears in largescale models (69)(70)(71). Linear death rate (often imposed by dilution) is also another common form of stress (66, 89, [95][96][97].
Inspired by these observations, we first examined different methods to modify a logistic growth equation.
Consider a basic logistic growth equation: This equation only has two parameters, growth rate and carrying capacity . If we can derive a nondimensionalized logistic growth equation, it can be rewritten as: dX dτ = X(1 − X) (τ = t, Χ = N N m ⁄ ). (II. 2) The followings are common modifications of the above equation: 1. The growth rate can be modified by multiplying the right-hand-side with a constant: This equation modifies the growth rate from 1 to , where increasing increases growth rate and leaves the carrying capacity the same. can take any real number in this case. 2. The carrying capacity can also be modulated, leaving the growth rate unchanged: The carrying capacity becomes 1⁄ and > 0. This formulation does not generate negative growth. 3. The "1" in the logistic growth equation can also be modified: In this case, both growth rate and carrying capacity are scaled by a factor of . can take any real number because if the carrying capacity is negative, the growth rate will also be negative and thus the model will only generate bounded growth. 4. A first order term can be added to the equation: Although this modification is equivalent to the previous one, it is worth noting since this form is commonly used to represent death rate or turnover rate.
The above analysis demonstrates that 1, 3, and 4 are robust modifications that provide consistent model modulations with any real value of parameter. Following this analysis, we will then model stress, benefit and cost at these three locations in the logistic growth equation to capture the assumptions presented above for mutualistic interactions.

D. Specific implementations of benefit, stress and cost
To capture the model assumptions of a), b) and c), we used the following formulations to modify simple logistic growth equations at location 1, 3 or 4 mentioned above to capture the logic of mutualistic interaction.

a) Benefit shall be positively dependent on partner density.
This assumption can be modeled at all three locations by where 2 represents partner density (same as below) and and ′ are measures of strength of benefit. In both cases, the benefit function is bounded, which facilitates the stabilization of population densities (73).
b) Cost shall decrease growth rate or carrying capacity. We use as a measure of level of cost. Specifically, this is implemented at location 1 or 3 by multiplication 1 ( ≥ 1) (density-independent) or 1 2 + 1 Cost can also be implemented as a turnover rate at location 4 − ( ≥ 0) (density-independent) or − 2 ( ≥ 0) (density-dependent). Stress can also serve as a turnover rate at location 4 − .
(II. 12) Two or more of the parameters can share one location, so we permutated the three factors in the three locations and systematically generated 81 models (Table S1). The 81 models are then checked against assumption d) and verify that the models do not generate unbounded growth. See Supplementary Materials for the MATLAB code that investigates the model assumption d) and unbounded growth behavior. Out of the 81 models, 48 satisfy all 4 assumptions and do not generate unbounded growth with positive parameter values and positive initial densities.

III. Derivation of diverse criteria
This section demonstrates the diversity of criteria with various model formulations. We first use an example model to serve as the base model to demonstrate a detailed process to derive coexistence criteria. We then show the process of building various mutualism models with more complexity and deriving or approximating the corresponding criterion.

A. An example
Consider the following model as our base model for the analytical analyses: The non-dimensionalized version of the model is (model 21 in Table S1): Where τ = t, Χ 1 = Χ 2 = N N m , δ = d , β = K . The ranges of the parameters are: ≥ 1, δ ≥ 0, β ≥ 0 This system has 5 fixed points: By calculating the Jacobian, we found that the 4 th fixed point (III.8) represents stable coexistence. The 5 th fixed point is a saddle point and is unstable. For the 4 th fixed point to be in the first quadrant of the Real domain, 2 conditions must be satisfied at the same time: The first condition (III.10) derives the coexistence criterion (β+1) 2 4β ≥ and when this criterion is satisfied, the second condition is automatically satisfied if β > 1. When β ≤ 1, only facultative mutualism can be captured. Bifurcation analysis of the system shows Allee effect and how final density increases with increasing β, decreasing or decreasing .

B. Summary of criteria derived from symmetric models
Symmetric models in general generate interpretable analytical solutions. In Table S1, all the models are assumed to be symmetric between two populations. We found that these symmetric models already generate a diverse set of coexistence criteria.
Most of these criteria are derived using the same logic used for the base model presented above, except for the criteria of model 1, 4, 28, 31, 55 and 58 in Table S1, which are derived from setting the analytical solution of the saddle point lower than the initial density Χ 0 . This is because the steady states that represent coexistence for these models are constants, whereas the saddle points are modulated by model parameters.
For a detailed example of the derivation, refer to section III.C.2.
In natural mutualism systems, cost can also scale with its partner's density (100). In general, we observed that inclusion of density-dependent cost increases the complexity of the criteria. In addition, mutualists can also compete in within the same niche (101). In this case, the effective benefit term in general decreases from the model without competition, indicating the criteria also lumps the effect of competition.

C. Criteria from models with additional complexities
To relax some additional assumptions, such as complete symmetry, we also included additional levels of complexities to the base model by including asymmetry and turnover rate: dN 1 dt =

Criterion with two populations competing for resources
Our base model assumes the two populations having separate carrying capacities. However, in natural settings, mutualists often share resources (83). We can impose competition between them by adding − Χ 2 or − Χ 1 to location 3 in logistic growth equations of Χ 1 and Χ 2 , respectively (number 82 in Table S2): Using the same method described in section III.A, we derived the criterion of this model:

Criterion considering initial density
For this analysis, we used a symmetric model (number 83 in Table S2): Because most mutualism models generate Allee effect, coexistence can also be affected by initial density ( 0 ). However, we can derive a benefit term that depends on the initial density 0 (for both Χ 1 and Χ 2 ) to predict deterministically whether the two populations coexist or not. To coexist, 0 needs to be greater than the saddle point:

Criterion with asymmetric growth rate, cost, stress, and benefit
Asymmetric parameters of the mutualism are more realistic in capturing real-world mutualism systems, so we relaxed the assumptions that cost ( ), rescue strength (β) and stress (δ) are symmetric for both populations. The asymmetry of growth rate is captured by . In the following model (number 84 in Table  S2), we assumed Χ 1 and Χ 2 share the same carrying capacity. We found that separated carrying capacity also yields similar results.
Use the same logic presented in section III.A, the following criterion is a necessary condition for the 4 th fixed point to be present in the positive Real domain: ρ( 1 2 + 1 + 2 ) 2 4 1 2 ≥ 1 1 1 ρ + 2 2 2 .

(III. 29)
This criterion can also be rewritten as . Note that we specifically used Χ 1 as the reference population.

Criterion with turnover rate and asymmetric growth rate
If we consider both the asymmetry in growth rate ( ) and a turnover rate 0 , the model becomes (number 85 in Table S2): The analytical solution of this model is complex and involves solving 4 th order polynomials. However, we can introduce the concept of correction terms to approximate the criterion. When we only add in the model and assume 0 = 0, we get the following criterion: In addition, when we only add the 0 term in the model and assume = 1, we get (β(1 − 0 ) + 2) 2 4β ≥ .
Thus, we know that the criterion is more accurate when The approximated criterion is accurate when = 1.When → +∞ and 0 → 0, the criterion will also provide good approximations.

The criterion for N-mutualist systems
The following model is used to determine the criterion with N-mutualist model: The non-trivial steady state can be solved by solving the following equation for X * : n ε X * (1 − nX * ) − nδ β(n − 1)X * + 1 X * = 0 (III. 40) Using the same strategy in III.A, we derive the criterion for coexistence: where benefit becomes a function of both β and n. This result suggests that mutualism system can tolerate higher stress levels with increasing number of mutualists.

Mutualism model including cheater exploitation
The mutualism model including cheaters is adapted from III.22-23: We assume 1) cheaters and cooperators share the same carrying capacity, 2) cheaters accept benefit produced by cooperators but do not provide benefit, thus do not experience cost, and 3) there is a constant transfer rate ( ) from cooperator to cheater, representing mutation from the cooperator phenotype to cheater phenotype. Although this model does not generate stable coexistence of the cooperators, the time it takes for cheaters to take over the cooperators can serve as a metric for how stable the mutualistic system is.

Model structures that generate lower boundary for stress
We notice that to coexist, some model formulations not only require an upper boundary of , but also require a lower boundary as well. If is lower than a threshold, the system can be dominated by the population that is more fit. However, this lower boundary occurs due to the system dynamic shifting from a mutualism-dominated mode to a competition-dominated mode.
There are two types of system dynamics that can lead to a loss of partner. One is due to high stress where the weaker partner will go extinct and the stronger partner will suffer from the loss of its partner. The other is due to lack of stress where the system shifts to a competition-dominated interaction where the fitter partner survives better by excluding its partner. In our study, we only focus on the first type of partner loss since the second type is a dynamic of competitive systems instead of mutualistic systems.
In general, we observed that a model needs to be both asymmetric and potentially competitive to create a scenario where the fitter population excludes the weaker population. The lower boundary for decreases when the weaker population gives out more benefit to its partner or reduction of competition. This observation can potentially explain why some mutualistic systems collapse when external stress is reduced (89,102).

IV. Theoretical generality of the simple metric
To establish the generality of the theoretical criterion, we verified that it is applicable to both symmetric and asymmetric mutualistic systems. We also want to test that the predictive power of the predictor is maintained when the partners are obligatory or facultative. Because mutualists can compete for the same resources in nature, we also verified that coexistence of partners that compete for resources also can be depicted by the criterion. We show some typical results of the predictive power of the criterion in these cases with Fig. S2. Fig. 1f is generated using the model presented by equation III. 16-17, where ε ∈ [1, 1.2], δ ∈ [0, 2], β ∈ [2, 5] and = 1. The positive trend holds even when the mutualism is facultative (δ < 1). S2a is also generated using the model presented by equation III. 16-17. Both are generated with β ∈ [2,10], δ ∈ [1.2, 2], ε ∈ [1, 1.2] and = 1. The panel on the left assumes Χ 1 and Χ 2 have the same initial density which is a uniformly distributed between 0 and 0.5. The panel on the right assumes changing of individual initial densities in a linear fashion while maintaining the sum of initial density of Χ 1 and Χ 2 at 1.

Fig.
Where

B. The predictive accuracy is maintained in the presence of noise
In addition to investigating the probability of coexistence when initial density is randomly distributed (Fig.  S2a, e), we also explicitly modeled noise to test the effect of noise on the prediction accuracy of the criterion. We added Gaussian noise to our base model ( ): In Fig. S2g, we used the model above and the criterion tested is III.19. Multiple sets of initial densities are tested. Although we observed a decreasing accuracy with increasing noise, the prediction accuracy is robustly maintained (above 90%) even when the standard deviation of noise reaches 50% of signal strength.
The high maintenance of prediction accuracy is due to the stability of the mutualistic model. Because the two stable steady states exist on opposite sides of the separatrix, only when noise pushes densities across the separatrix (against the vector field), the outcome of the system will be altered. Otherwise, noise will not influence the final steady state.

V. Calibration using SVM
We chose SVM because it is a well-established algorithm and it only requires information of the support vectors which are usually data points along the boundary to obtain a classification boundary. Note that other algorithms can also be implemented for the same purpose. For example, when predicting probability of coexistence, logistic regression can be more suitable to directly predict probability of coexistence. See the Supplementary Materials for our MATLAB program which performs the calibration procedure using SVM.

A. Kernels
We used four types of kernels to cover different general shapes of the landscape: Linear kernel: Sigmoidal kernel: These four kernels follow the same structure: 〈[ , 1 ], [ , 2 ]〉 = 〈 , 〉 + ( 1 • 2 ), Our procedure allows any kernel structure that is supplied by the user, so other customized kernels can also be used for the calibration.

B. Standardizing input data
Standardizing input data is essential for robust learning of ( ). Before training the model with SVM, we first standardize the system variables and stress to mean 0 and variance 1.
C. From SVM output to calibrated ( ) SVM will output a predictor that separates coexistence and collapse indicated in the vector. Our program requires that the two classes are represented by 1 and -1 in the vector. Using kernels V.1-4 , we can write the predictor in the general form: where is a kernel parameter; represents coefficients associated with each input observation that are optimized by the SVM algorithm; 0 is the bias that is optimized by the algorithm. If we impose = at ([ , ]) = 0, we can obtain a function of in terms of input data: Using this function, we can obtain a value of with any system variable dictated by based on the limited number of observations used as inputs.
( ) can sometimes have wrong directionality, meaning that / > 1 is associated with collapse and / < 1 is associated with coexistence. We identify these cases by calculating If the above expression is negative (when the trained boundary have a learning accuracy greater than 50%), the calibrated ( ) has a wrong directionality. Flipping of the landscape is then required, and it is done by assuming: We name the output landscape with properly adjusted directionality ′. This landscape ′ describes the shape of the final calibrated ( ) but still needs to be rescaled for prediction purpose. The final ( ) is calculated by rescaling it back according to the mean and variance of the measurements.

D. Cross validation and bootstrapping
All the cross validations in this work are 10-fold cross validations. The cross-validation accuracy is represented by the average values. Bootstrapping is used to evaluate the degree of variation of quantified . The same number of data points as the input data is randomly sampled from the input data with replacement. We performed bootstrapping for 500 times. The variance and relative standard deviation (RSD) are then calculated based on the 500 bootstrapped ( ) quantified with 500 sets of sampled training data. The mean cross validation accuracy, the bootstrapped variance and relative standard deviation are then used to evaluate the accuracy of the ( ) outputs.

E. Twenty sets of simulations to establish and test the calibration procedure
To establish the calibration process, we first developed the method with simulated data where the true ( ) are known. This allows us to evaluate calibration results against the true function. We used the model presented in equation III.22-23 to conduct this set of analyses. Using 1 as the reference population, the ( ) can be expressed as: ( 1 + 1 1 + 1) However, depending on how system variables change the model parameters ( ( )) in the system variable space, the true ( ( )) can exhibit diverse shapes. With 1 , 2 ∈ [0,1], we constructed 20 arbitrary set of equations with different underlying functions. We also made sure the generated data have around 50:50 split of coexistence and collapse in the simulated results to reduce bias in the input data. This 50:50 split constraint also applies to the all other simulated data. Note that these functions are arbitrary and only serve the purpose of generating true ( ) that has various underlying functions. The data generated using these 20 models are included in the MATLAB code in the Supplementary Materials.
The 20 sets of ( ) are grouped into 5 types with each type containing 4 different examples. The 5 types of functions that describe ( ) include linear, quadratic, cubic, and Hill equation and another type of functions that are a mixture of the previous four types. As examples, the following are four sets of ( ) used to generate simulation results shown in Fig. S4b.

F. Calibration with simulated data with unknown ( ): an example
We want to test whether we can apply the calibration procedure to data that are generated by an arbitrary mutualism model. The specific model structure we used to generate data in Fig. 2c is: The analytical solution of this model cannot be expressed in a simple and explicit form. However, using simulated data, we can obtain an empirical function of on the system variable space that allow further prediction in the system variable space. The simulations are done using initial densities of [0.01, 0.01] for 10 unit-time.
To get the probability of coexistence, we ran the model 100 times with varying ratios of initial density while keeping the total initial density constant as 0.02. The density for a population is varied in a linear fashion from 0 to 0.02 and we terminated the simulation after 10 unit-time.
To calculate how well the systems can resist cheater exploitation, we modified the equations V. 16 − 17 to account for the emergence of cheaters: The simulations were done with = 10 −3 . The initial densities are [0.01, 0.01, 0, 0] for Χ 1 , Χ 2 , Χ 1c , and Χ 2c , respectively. The simulations were terminated after 1300 unit-time.

VI. Framework generality verified by complex mutualistic systems A. A 5-mutualist system
The data presented in Fig. 4a and Fig. S9a  All the parameters are linear combinations of 1 , 2 ∈ [0,1] and the coefficients are randomly generated.

B. The experimental 3-member mutualistic systems
The calibration shown in Fig. 4b is done on all 384 systems with triplicates that alternatively use each strain in a system as the reference strain. The δ measurements used in this calibration is the same as the measurements presented in Fig. 3j and Fig. S8b with the pairwise systems.

C. A mutualistic system in an oscillatory environment
The oscillatory signal is implemented as square pulses, where is time; is the duration of the pulses and the duration between the end of one pulse and the start of the next pulse; is the intensity of the signal. Throughout the simulation, δ 1 , δ 2 , β 1 , and β 2 are affected by this oscillating signal, where δ 1 = δ 2 = 0.8 + (VI. 7) β 1 = β 2 = 2(1 + ) (VI. 8) To simulate the data, we used a logarithm scale that vary from 10 −2 to 10 2 and vary from 10 −1 to 10 1 . The model used to generate the set of data in Fig. 4c and the left column of Fig. S10 is presented in equation III.22-III.23, where = 1.5. The initial densities used are [0.2, 0.2] for X 1 and X 2 respectively. The simulations are performed for 500 unit-time. The theoretical ( ( )) is calculated using III.32 with δ 1 = δ 2 = 0.8 + and β 1 = β 2 = 2(1 + ).

D. A mutualistic system cohabiting with other populations
To generate simulated data presented in Fig. 4d and the right column of Fig. S10 The initial densities are [1, 1, 1, 1, 1, 1, 1] for X 1 to X 7 respectively. The simulations were performed for 5000 unit-time. The theoretical ( ( ))/ in the right panel of Fig. S10a is directly calculated using III.32.
Violates assumption d) (β + 1) 2 4β > δ Violates assumption d)   Figure S1: Typical bifurcation diagrams of models that create unbounded growth. Previously-developed general coexistence criteria describe the boundary between stable coexistence (white regions) and unbounded growth (grey regions). Solid lines represent stable steady state for coexistence. Grey dashed lines represent steady states of a mutualist in the absence of its partner. This diagram is generated using model: The left panel is generated with β = 1 and the right panel is generated with δ = 1. Figure S2: The predictor δ ⁄ is a general metric for predicting performance of symmetric (panel a, b) and asymmetric mutualistic systems (panel c-f). All the traces are comprised of individual dots. Each dot represents result from one set of model parameter. a. Mutualistic systems have initial-density-dependent coexistence. In the left panel, we assume the same initial density of both populations that is a uniform random variable. In the right panel, we varied the ratio of the two initial densities and kept the total initial density constant. In both cases, ⁄ is predictive of coexistence probability. b. When cheaters can arise in a system, ⁄ is predictive of the time duration the system can persist before cheater exploitation. c. An asymmetric mutualistic system can be either obligatory (left panel), where both populations extinct when ⁄ < 1 or facultative (right panel), where one population can persist even when ⁄ < 1. The black dots represent total density and the two colors represent the total density of the two partners. The same representations are used in panel d. d. The predictive power of ⁄ for total density also holds for asymmetric systems that share carrying capacity. e. ⁄ is predictive of probability of coexistence for asymmetric systems. f. The predictive power of ⁄ for resistance to cheater exploitation also holds for asymmetric systems. Please refer to SI for models and parameter values used in generating this figure. g. The accuracy of criterion is robustly maintained after addition of extrinsic Gaussian noise. X axis indicates the standard deviation of the Gaussian noise. Note that the noise is significant considering the maximum total density of the system is 1. Please see Supplementary Text section IV.B for the simulation details. Accuracy Figure S3: A general method for quantifying . Regardless of different model formulations, can always be measured using either growth rate or final population size of an individual population. a. We simulated the growth of a population that is modulated by a death term . Stress is simply normalized by growth rate . In the simulations, the true ( ) depends on an independent variable . b. One type of model structure has modulating both growth rate ( ) and yield ( , final density). is the population size and is the carrying capacity. Based on the growth curve or final population size, measurements of growth rate ( ) and yield ( ) can be obtained when varying . The measurement ( ) can be calculated by 1 − max ( ) or 1 − max ( ) . Note that since is always non-negative, 1 − max ( ) is always positive, but can be greater than 1 when the population size decreases. To address this discrepancy, extrapolation of is needed to capture death rate (the black trace in the versus plot). c. In another type of model formulation, only modulates growth rate. Simulation shows that either equals to carrying capacity or 0. Thus, in this case, only can be used to calculate . The right-most panels show the comparison between the measured ( ) and true ( ). In all cases, is a good approximation of (the solid black line represents = ).   Figure S5: Application of the simple rule to a complex model with no explicit form of . a. 1 , 2 , 1 , 2 , 1 , and 2 vary with 1 and 2 . b. The calibrated ( ) along with is predictive of probability of coexistence. To get the probability of coexistence, at each ( 1 , 2 ), we ran 100 simulations with 100 different ratios of initial densities while keeping the total initial density the same. c. ( )/ is also predictive of how well the system can resist cheater exploitation. The y axis indicates the time the system can persist before the exploitation of cheater populations. d. The high pairwise R 2 values of the top 5 ( ) indicate high confidence of the calibration results. e. The top 5 ( ) from the calibration procedure. The title for each ( ) indicates its index corresponding to panel d and its corresponding training accuracy. Each ( ) can also be evaluated based on their bootstrapped relative standard deviation ( ).  The normalized equilibrium density of ∆ and ∆ were used to quantify for each strain. We only fitted density that is above the detection limit (10,000 cells/well). Since ∆ reaches below detection limit when [Leu] is relatively low, extrapolation is done to estimate the death rate at lower [Leu] (see Fig. S3a for the reasoning). As a result, of ∆ has a higher dynamic range than that of ∆ . b. 5 top ( ) show both consistency and discrepancy. All ( ) indicate an optimal initial density at an intermediate level. However, when [Trp] increase to 16μM, ( ) can either increase or decrease. c. Our approach can also predict probability of coexistence. We excluded the ratio of initial density, leaving [Trp] to construct a one dimensional . The boundary between coexistence and collapse is between 0.1 μM tryptophan, 0.8 μM leucine and 0.5μM tryptophan, 4.0 μM leucine. To calibrate a ( ) for this system, we assume supplementing the two amino acids does not change the cooperation capability. This calibration demonstrates that our procedure can also predict probability of coexistence for experimental systems. Figure S10: Mutualistic systems in dynamic environments. a. Two systems were analyzed: one system is a pairwise mutualism system that is modulated by oscillatory signals (the left panel) and the other system contains a pair of mutualists and 5 bystander populations that interact with the mutualists (the right panel). We verified that theoretically calculated / roughly holds as a predictor for total density. The transition between coexistence and collapse, however, do not occur strictly at 1. Each dot represents a simulation result from one set of model parameters. b. The data used for the calibrations and the calibrated ( ). The two ( )'s are used for the prediction plots shown in Fig. 4c, d. Movie S1. This video demonstrates the relative position of the boundary surface ( , ) = 0 and 5 observations in the schematic shown in Fig. 2b. Input data in Fig. 2b can be represented in a 3D space. The color of the dots and the z axis positions both indicate values. Coexistence is indicated by closed circles and collapse is indicated by open circles. The gray boundary surface separates coexistence and collapse. The surface is above observations that represent coexistence ( > ) and below observations that represent collapse ( < ). This surface can be directly interpreted as ( ) since ( , ) = 0 ⟹ ( , ) = 0 ⟹ ( ) .
Movie S2 -S4. These 3 videos show the calibrated surface of ( ) relative to input data in 3D space. The blue surfaces are the boundary ( , ) = 0 that separate coexistence and collapse. The surface is also equivalent to ( ). The indices on x and y axes in video 4 correspond to the same strain orders in Fig. 3j.