A unifying framework for interpreting and predicting mutualistic systems

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. Furthermore, even when a rule is established, it is often challenging to map it to mechanistic details and to quantify these details. Here we report a framework that addresses these challenges for mutualistic systems. We first deduce a general rule that predicts the various outcomes of mutualistic systems, including coexistence and productivity. We further develop a standardized machine-learning-based calibration procedure to use the rule 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.

By analyzing fixed points of several differential equation based models of mutualisms, Wu et al. have proposed a simple rule for coexistence of mutualistic species, based on generalized notions of benefit and cost associated with mutualism and the stress experienced by populations in the absence of mutualistic interactions. The rule states that coexistence of mutualistic partners requires the effective benefit to be larger than the stress. The authors have further devised a machine learning based calibration method to estimate the effective benefit based on measurements of a set of "context" variables. By applying this method to data obtained from simulations, previous experiments as well as a new experiment on synthetic mutualistic toxinantitoxin producing strains of E. coli the authors show that the effective benefit B can be used to make quantitative predictions of population size in mutualistic systems, without prior knowledge of mechanistic details. From a conceptual standpoint, and given the class of models studied by the authors, I am not sure if I find the "rule" postulated here (namely B(theta)>delta) to be surprising. However, I find the fact that the effective benefit B can be extracted from a set of independent observables without any knowledge of the actual benefit and cost associated with mutualism to be interesting and useful. Given that mutualisms are indeed observed in a diverse array of ecosystems, I think there is sufficient merit in devising predictive schemes such as the one developed in this paper. I therefore feel that the paper could be suitable for publication in Nature Communications.
Major comments: 1) My primary concern with the paper is how difficult I found it to understand the procedure from the main text, despite the fact that I am very much part of the target audience. Might it be possible to use a very simple model when describing the initial motivation (stress, benefit, etc) and then to later show that the approach can be generalized to the wide class of models? 2) Is it fair to say that the fact that the authors' approach works means that all of these mutualisms are effectively one-dimensional? Does this mean that there is a separation of timescales, or is that not necessary? Some more discussion of when the approach will fail would perhaps clarify the limits as well as why the procedure works at all. 3) My primary concern with the procedure is that it entirely based on the dichotomy of coexistence and collapse, and appears to rely on either of these two states being globally stable fixed points of the system for a given set of parameters. As the authors themselves have pointed out, in their analysis of the yeast cross-feeding mutualism data, their method fails in the limit of high amino acid concentrations, where the global fixed point corresponds to competitive exclusion. A similar problem might also occur if multiple stable attractors coexist, because in that case, the outcome of the mutualism would depend strongly on the initial conditions. 4) It is not clear to me whether the method works if the attractor is a limit cycle instead of a fixed point, as is the case with the bacterial cross-protection mutualism (Yurtsev, Conwill and Gore, PNAS, 2016). It may still be the case that the rule makes qualitatively correct predictions as to whether the mutualism survives or not, but I would expect the quantitative predictions to be more erratic. It would be helpful if the authors could elaborate on this.
Minor comments: 1) It would be nice to include a small discussion on how the calibration method could potentially be extended to include spatial structure. Is it possible to incorporate space as an additional variable in the set of "context" variables, or would it require the supervised learning algorithm to be run on a class of models that explicitly include spatial structure?
While my primary expertise is in the evolutionary analysis (and modelling) of mutualism, I very much enjoyed reading this ambitious paper on ecological dynamics of mutualistic interactions. It has a big aim: to find general rules predicting population dynamics across (almost) all mutualistic systems, particularly coexistence and collapse conditions, as well as quantitative predictions (density, coexistence probability, time to cheating take-over).
In ecology and evolution, we can sometimes be too affected by what Queller called the 'tyranny of detail' (Am Nat. 2017 Apr;189(4):345-353), to even attempt finding such general rules. Abstract, higher-level (coarse-grained, as the authors call them) rules, can be very useful in guiding research, however, and I really like how the authors try to nevertheless find them.
To do so, they analyse a range of 81 different mutualisms models, based on a large set of ODEs, and aiming to reflect (most of the) diversity of mutualistic systems out there (more later on this). Deriving coexistence criteria these models, they come up with a general rule which predicts mutualism coexistence as a function of the effective benefit and stress experienced by the mutualistic partners (Eq 1). They then proceed to show how these parameters, or approximations of them, could be empirically measured and do so using a new experiment, simulations and a number of previously published (microbial) datasets.
The general rule the authors derive seems novel to me but makes a lot intuitive sense, and the application to simulations and particularly real-life datasets are both informative and convincing. Together they make for a convincing analysis, and a very useful contribution, with many potential applications in the study of mutualisms. The paper is also well written and overall a pleasure to read, and the code to replicate all theoretical analyses and simulations has all been made available.
However, I have one major concern, which I feel the authors need to address because it has the potential to substantially affect the claimed generality of these results.
As I understand it, the various models analysed by the authors, all have in common that benefit is 'positively dependent on partner density' (Third line methods section, Section IIB Supplementary), although this benefit is bounded (section D Sup info, near equation II.8). However, while there are certainly situations in which this is undoubtedly true (for instance, in many of the lab microbial interactions studied empirically in this work), I am not necessarily convinced that such positive dependence on partner density hold generally among mutualistic interactions. After all, there is lots of evidence from natural systems of saturating benefits from increasing partner densities, and indeed also of shifts to negative effects beyond certain densities (e.g. Morris, Vázquez and Chacoff, 2010;Palmer and Brody, 2013 -> refs at the bottom of this review; the DeAngelis, Holland & Bronstein 2002 paper cited by the authors also gives a few examples). Such saturation of benefits makes a lot of intuitive sense: at some points all the flowers in a plant population are pollinated, or all the herbivores chased away by protective ants. Additional ants are then at best neutral, and potentially a fitness cost (see . The authors do mention previous work on saturating benefits in their supplementary information (in the section where they discuss previous models and again when they discuss their own modelling approach), for instance discussing some of the Holland & DeAngelis work on this, but as I understand it they don't incorporate any potential saturating effects (or even shifts to negative effects a high-densities) in their models. Of course, with some exceptions, we don't generally know to what extent real-life mutualisms typically experience (bounded) positive benefits from increased partner-densities as in the current work versus when saturating or even negative density effects start to appear. However, given that we have good evidence (and theoretical reasons) to think that they may commonly exist, this could undermine the claim of generality in this paper. Consequentially, I would like to see the authors either (i) include such effects in their models and analyse their impact, or (ii) convincingly explain and show why this concern is not relevant and their conclusions would hold even under saturating benefits, or (iii) scale back the claim of generality and clearly indicate that their results may not (all) apply to mutualisms with saturating effects, and are only directly applicable to the (much?) smaller subset of mutualisms with strict positive density effects.
I also have some more smaller remarks: -Please include line numbers throughout for a potential resubmission, including for the supplementary. Referencing sections of the manuscript is very cumbersome without them.
-I think it would be important to highlight more explicitly in the main text that the models analysed here are all ecological models, and do not include potential evolutionary responses. It's fine to not include evolutionary dynamics, but it's important for the reader to be aware of this limitation, particularly given that eco-evolutionary dynamics could actually be important, particularly in many of the microbial mutualistic systems considered here.
-Page 8 -Experimental Application of the metrics: for the first experimental application of the metric, I didn't fully understand what, if anything, the model equivalent of aTc is. The stress is imposed/modulated by IPTG, if I understand correctly, so is this just some external trigger of Quorum-sensing, independent of stress, that for whatever reason is required in this system? Would we not want QS to be triggered by the stress/stress-inducer itself? Some clarification for the reader would be helpful here.
-Main text table: I don't very much like cancer as an example of mutualism here. I know that the authors are using a slightly wider definition of mutualism here of two populations providing reciprocal benefits (first sentence of the ms), but most typically the term is used for situations where these populations are also of different species. I would suggest instead using an interspecific protection mutualism (e.g. ant-plant). -SI page 3 "In addition, although population collapse [..]the fitness of the benefit-receiver decreases with increasing partner density, which is contradictory to the basic logic of mutualism." -> I appreciate that given the definition of 'mutualism', at the stage where effects are negative the interaction is no longer strictly speaking a mutualism, but following my above general remark I wouldn't describe this effect as being contradictory to the basic logic of mutualism. The interaction at some density having a negative effect and then no longer being a mutualism is no way contradictory to it being a mutualist in other conditions, which I think what the models discussed here were trying to analyse.

Response to reviewers' comments
Reviewer #1 (Remarks to the Author): By analyzing fixed points of several differential equation based models of mutualisms, Wu et al. have proposed a simple rule for coexistence of mutualistic species, based on generalized notions of benefit and cost associated with mutualism and the stress experienced by populations in the absence of mutualistic interactions. The rule states that coexistence of mutualistic partners requires the effective benefit to be larger than the stress. The authors have further devised a machine learning based calibration method to estimate the effective benefit based on measurements of a set of "context" variables. By applying this method to data obtained from simulations, previous experiments as well as a new experiment on synthetic mutualistic toxin-antitoxin producing strains of E. coli the authors show that the effective benefit B can be used to make quantitative predictions of population size in mutualistic systems, without prior knowledge of mechanistic details. From a conceptual standpoint, and given the class of models studied by the authors, I am not sure if I find the "rule" postulated here (namely B(theta)>delta) to be surprising. However, I find the fact that the effective benefit B can be extracted from a set of independent observables without any knowledge of the actual benefit and cost associated with mutualism to be interesting and useful. Given that mutualisms are indeed observed in a diverse array of ecosystems, I think there is sufficient merit in devising predictive schemes such as the one developed in this paper. I therefore feel that the paper could be suitable for publication in Nature Communications.
We thank the reviewer for recognizing the relevance and significance of our work and for constructive suggestions and comments.
Major comments: 1) My primary concern with the paper is how difficult I found it to understand the procedure from the main text, despite the fact that I am very much part of the target audience. Might it be possible to use a very simple model when describing the initial motivation (stress, benefit, etc) and then to later show that the approach can be generalized to the wide class of models?
We thank the reviewer for this constructive comment. We have included in the main text a simple mutualism model to demonstrate the derivation process (page 5 line 78-84). In addition, to further clarify the calibration procedure, we have included a new Extended Data Figure 5, which provides a step-by-step illustration of the calibration procedure using the simulated data of a complex mutualism model. This figure is shown here as Fig. R1. and . For an experimental system, and would correspond to experimentally controllable parameters, which could affect multiple mechanistic parameters simultaneously. c. Simulated time courses with different and values. The group on the left simulates cocultures (blue: and red: ). The group on the right simulates monoculture, where . The growth rates of the monoculture can then be used to calculate . We directly used the parameter value of in our procedure. d. Detailed procedure of the calibration process: 1. Data collection and preprocessing: Collect the measurements of coexistence vs. collapse ( ) and stress ( ) on the domain of . Above a certain threshold of total final density, assign a value of 1 and assign -1 if total final density is below this threshold. Standardize the values of and values to have mean of and standard deviation of and name the normalized vectors and . For simplicity of presentation, the following and are standardized. In this example, we have 100 observations ( , 10 by 10) and two dimensions of ( : We empirically examined the top five SVM models, which were select by the lowest values. The weights of and were determined empirically using simulated data (see Extended Data Fig. 4c). 3. Train the top SVM models: Train each of the five SVM models with normalized input data to obtain the function that describes the boundary between the two classes. The functions can be expressed as: is the weight of observation , and is the bias term for the SVM model. Both and are optimized by the SVM algorithm.
is the kernel function; is a kernel parameter; , , and are input values for observation . and are independent variables of the function. A positive value of for a new set of and predicts coexistence and a negative value predicts collapse. 4. Quantify the top and assess its reliability: Quantify from each SVM model by imposing at . If has the correct directionality, get the function of by calculating and are calculated using the original measurements. See SI section V.C. for how to adjust for the directionality of . 1) Quantify relative standard deviations (RSD) of top : The sampling process is the same as step 2.2). Iterate 10 times to get 10 bootstrapped to quantify an RSD value on each pair and construct a variability landscape of . The title of each is the cross-validation accuracy of the model. 2) Quantify pair-wise consistencies of the top 5 models: To get R 2 of any two , the two scales of are unified by first using linear fitting to adjust one to conform to the scale of the other . R 2 is then calculated using the two that have adjusted scales. 3) Pick the best : The first is picked in this case since all top 5 have similar RSD and are all highly consistent. 5. Use in downstream predictions: Use the best to construct the metric . is indeed positively related to final total density in this example.

4
The same procedure is applied to the analysis of other simulated data and experimental data. The different sets of data only differ in terms of how they are generated and what and correspond to. e. The calibrated along with is predictive of probability of coexistence. To get the probability of coexistence at each ( ), we ran 100 simulations with 100 different ratios of initial densities while keeping the total initial density the same. f.
is also predictive of how well the system can resist cheater exploitation. The y axis indicates the time the system can persist before the cheater populations compete out the cooperators.
2) Is it fair to say that the fact that the authors' approach works means that all of these mutualisms are effectively one-dimnsional? Does this mean that there is a separation of timescales, or is that not necessary? Some more discussion of when the approach will fail would perhaps clarify the limits as well as why the procedure works at all.
We thank the reviewer for these insightful questions and observations. Our brief response is that yes, our claim is that, at a level of abstraction that is extremely general but still meaningful, it is possible to interpret and predict mutualism in a simple way with and .
Response to "one-dimensional": Although by formulation, the underlying models are at least two dimensions (as summarized in SI Table 1 and 2), one can indeed interpret our criterion as a "onedimensional" simplification and generalization of mutualistic systems. The emergence of our criterion indicates that, along a properly chosen trajectory (dictated by the shape of estimated B), certain coarsegrain outcomes of a mutualistic system (coexistence versus collapse, total density, probability of coexistence and persistence to cheater exploitation) approximately align with each other and with in a monotonic manner.
Response to "separation of time scale": From our derivation, the separation of time scales does not appear to be critical. The mechanism-based criteria derived from specific models (SI Table 1 and 2) emerge from the steady-state solutions of these models. As such, the potential separation of time scales in these models is masked in the corresponding criteria.
Response to "when the approach will fail": As demonstrated in the paper, our criterion and procedure work for the coarse-grained interpretation and prediction of the transition between mutualistic coexistence and collapse, total productivity, probability of coexistence, and time it takes for cheater exploitation. The criterion emerges from the analysis of steady-state solutions of two-population models. However, as indicated in the analysis in Figure 4, the criterion and procedure also work for higher-dimensional systems and for systems exposed to certain dynamic perturbations. These results further demonstrate that the steady-state assumption is not a critical assumption.
As noted below, our approach also works if a mutualistic system exhibits limit-cycle oscillations.
In general, it remains an open question whether and to what extent this approach would be applicable beyond the scope of mutualism and if the mutualistic system becomes much more complicated than what we have tested. We have included a sentence on page 13 (line 295-297) of the main text to discuss this limitation.

Response to "why the procedure works at all":
The calibration procedure (based on SVM) we have developed works when we know a priori that a system is dictated by an inequality: in a general term, X > Y dictates qualitative outcomes and X/Y is positively correlated with quantitative outcomes.
The theoretical foundation of our procedure is the class of mutualism models we have constructed and analyzed. The class of models share two key properties: a) a general criterion that underlies coexistence of two populations; b) is positively correlated with total density, probability of coexistence and resistance to cheater. These positive correlations between and quantitative outcomes are predicted by our class of models.
The reason why this class of models behave similarly, we think, resides in the assumptions of our models. Our assumptions (see Method section and SI section II.B) strive to find the minimal description of mutualism: two populations help each other to reduce the stress that each receives and impose a cost to itself. Stripping away the specificities of criteria derived from each individual model, our criterion can be considered as a general quantitative description of the assumptions that are shared by all models.
We have included our clarification and discussion of the limitations of our procedure in main text (page 6 line 116-118 and page 7 line 141-142).
3) My primary concern with the procedure is that it entirely based on the dichotomy of coexistence and collapse, and appears to rely on either of these two states being globally stable fixed points of the system for a given set of parameters. As the authors themselves have pointed out, in their analysis of the yeast cross-feeding mutualism data, their method fails in the limit of high amino acid concentrations, where the global fixed point corresponds to competitive exclusion. A similar problem might also occur if multiple stable attractors coexist, because in that case, the outcome of the mutualism would depend strongly on the initial conditions. We thank the reviewer for these thoughtful comments.
Response to "global fixed points": We regret the lack of clarity on this point in our original manuscript. Our procedure does not require the collapse and coexistence states to be global fixed points. For example, a simple mutualism model can have both the collapse state (0,0) and the coexistence state with the same parameter set. In our SI tables where we calculated the criterion for coexistence, we analyzed the conditions where the fixed points corresponding to coexistence are real and positive. This means that if the initial density is sufficiently high (above the separatrix), the two populations coexist. The populations collapse when the coexistence stable steady state does not exist. Using the following model (eq. III.22-23 in Extended Data), we can generate vector fields that have multiple fixed points (Fig. R2):  Fig. R2. Vector fields of different parameter sets using the model presented in SI eq. III.22-23. The two vector fields in red circles both have two stable steady states. The criterion for the symmetric case is presented in Extended Data Table 1 model #48; the criterion incorporating initial density is in Extended Data Table 2  In addition, we have demonstrated that initial densities can be included in the criterion as a part of , beyond benefit ( ) and cost ( ). An example criterion that includes initial densities is shown in Extended Data Table 2, criterion #83 (derivation of the criterion can be found in SI. Section III.C.2). Thus, if initial density is a variable, it can be incorporated into the term of our criterion . If initial density is a random variable, we have also found that the probability of coexistence can be predicted by (Extended Data Figure 2a, e and Extended Data Figure 7c).

Response to "competitive exclusion":
We apologize that we have not made it clear in our main text that our criterion does not describe the transition between coexistence and competitive exclusion. The winner in competitive exclusion increases in fitness when the loser is excluded. In contrast, collapse corresponds to the transition where mutualistic partners experience a reduction in fitness when a population is extinct (main text page 3, line 34). In simulations, we found that when there is no competition, there is no competitive exclusion. Thus, the transition from coexistence to competitive exclusion is a property of competition, and not of mutualism.
The various criteria we presented in SI Table 1-2 were not derived in a way to predict the transition between coexistence and competitive exclusion. However, by imposing the fixed point or to be stable, future studies can use a similar approach to derive a general criterion that describes the transition between coexistence and competitive conclusion. Suppose the criterion is , where is a lumped term describing effective competition, we can still use our calibration procedure to calibrate for an empirical .
The following model demonstrates the relationships between these two transitions (SI eq. III.30-31), which incorporates both asymmetry and competition. If we do not include competition in the model ( ), competitive exclusion is not present: Using , and , we generated the following phase diagram of system behaviors (Fig. R3).

Fig. R3.
Both upper bound (blue) and lower bound (black) exist for when there is competition. The region below the black dots represents competitive exclusion. The region above the blue dots represents collapse. This diagram shows that transition into competitive exclusion and mutualistic collapse are two distinct boundaries. When there is no competition, the transition between coexistence and competitive exclusion is not present.
Since our criterion does not apply to the transition from coexistence to competitive exclusion, our prediction does not apply to the case when is lowered to an extent that competition starts to dominate mutualistic interaction. This applies to the yeast auxotroph example where becomes low when the two amino acids are supplied at high doses, as also noted by the reviewer.
In our initial submission, we included section III.C.7 in SI to discuss the distinctions between these two types of transitions. We have updated this section to include Fig. R3 to show the difference between the two types of boundaries and how the lower boundary of only occurs with competition. In addition, we have clarified this point in multiple places of our main text, including page 3 line 34, page 10 line 217-218 and added a new paragraph on page 6 line 102-109 to clarify this point.

Response to "multiple attractors":
Indeed, it is unclear if our criterion is directly applicable when a system consists of multiple attractors, all corresponding to coexistence. This will likely depend on specific systems. For instance, for such a complex mutualistic system, it may still be possible to have a single boundary separating collapse and coexistence. In such an example, however, the quantitative predictive power of the estimated function will likely be limited. Further studies are needed to formulate such a model (in a relevant manner) and to test the applicability of our criterion and approach.
We included a discussion on page 13 line 295-297 of our manuscript to clarify this raised point.
4) It is not clear to me whether the method works if the attractor is a limit cycle instead of a fixed point, as is the case with the bacterial cross-protection mutualism (Yurtsev, Conwill and Gore, PNAS, 2016). It may still be the case that the rule makes qualitatively correct predictions as to whether the mutualism survives or not, but I would expect the quantitative predictions to be more erratic. It would be helpful if the authors could elaborate on this.
We thank the reviewer's constructive comments and pointing us to this highly relevant paper. We have recaptured the simulation data presented in (Yurtsev, Conwill and Gore, PNAS, 2016), where a limit cycle exists for a pair of mutualistic bacterial strains. We have found that our calibration procedure still applies and provides high prediction accuracy.
By examining the paper suggested by reviewer, we found that the density of the two populations oscillate with a phase difference, where one increases, the other decreases. This complementary growth dynamic creates a rather stable total density overtime. Using the final timepoint of the final cycle for each condition, we classified the simulations into coexistence and collapse ( Figure R4d). Using data in Figure R4, we conducted the calibration procedure and our procedure provides an estimation of effective benefit and an average cross-validation accuracy of 96.8%. In this case, the overall densities do not have a wide distribution, so the quantitative predictive power of were not tested ( Figure R4e). This example further illustrates the applicability of our criterion and calibration procedure to diverse mutualistic systems. We have included Figure R4 as Extended Data Figure 11 in and included a description of the system on page 12 (line 263-271) in the main text.  Calibration results show that our calibration procedure still provides high prediction accuracy. The total densities corresponding to coexistence do not have a wide distribution, so the quantitative predictive power of were not tested Minor comments: 1) It would be nice to include a small discussion on how the calibration method could potentially be extended to include spatial structure. Is it possible to incorporate space as an additional variable in the set of "context" variables, or would it require the supervised learning algorithm to be run on a class of models that explicitly include spatial structure?
We thank the reviewer for this insightful suggestion. We have not attempted to incorporate spatial structure into our criterion. Indeed, as the reviewer mentioned, a possible way to do so is to use different metrics of spatial structure as one dimension in the context variable. For example, the context variable can be the distance between the seedings of two partners or the degree of intermixing of the initial seedings. The calibrated will then also be a function of these spatial contexts. We would expect decreasing seeding distance or increasing the extent of intermixing increases the benefit level and thus increases .
We also agree with the reviewer that an alternative is to explicitly include spatial structure in the models. When doing so, the criterion can be applied to local segments where the homogeneity assumption is appropriate.
In light of the reviewer's suggestion, we have revised our main text on page 7 line 135 and page 13 line 290-295 to include discussions on the possibility of incorporating spatial component in system contexts.

11
Reviewer #2 (Remarks to the Author): While my primary expertise is in the evolutionary analysis (and modelling) of mutualism, I very much enjoyed reading this ambitious paper on ecological dynamics of mutualistic interactions. It has a big aim: to find general rules predicting population dynamics across (almost) all mutualistic systems, particularly coexistence and collapse conditions, as well as quantitative predictions (density, coexistence probability, time to cheating take-over).
In ecology and evolution, we can sometimes be too affected by what Queller called the 'tyranny of detail' (Am Nat. 2017 Apr;189(4):345-353), to even attempt finding such general rules. Abstract, higher-level (coarse-grained, as the authors call them) rules, can be very useful in guiding research, however, and I really like how the authors try to nevertheless find them.
To do so, they analyse a range of 81 different mutualisms models, based on a large set of ODEs, and aiming to reflect (most of the) diversity of mutualistic systems out there (more later on this). Deriving coexistence criteria these models, they come up with a general rule which predicts mutualism coexistence as a function of the effective benefit and stress experienced by the mutualistic partners (Eq 1). They then proceed to show how these parameters, or approximations of them, could be empirically measured and do so using a new experiment, simulations and a number of previously published (microbial) datasets.
The general rule the authors derive seems novel to me but makes a lot intuitive sense, and the application to simulations and particularly real-life datasets are both informative and convincing. Together they make for a convincing analysis, and a very useful contribution, with many potential applications in the study of mutualisms. The paper is also well written and overall a pleasure to read, and the code to replicate all theoretical analyses and simulations has all been made available.
We thank the reviewer for finding this work significant and enjoyable to read, and for making insightful and constructive suggestions and comments.
However, I have one major concern, which I feel the authors need to address because it has the potential to substantially affect the claimed generality of these results.
As I understand it, the various models analysed by the authors, all have in common that benefit is 'positively dependent on partner density' (Third line methods section, Section IIB Supplementary), although this benefit is bounded (section D Sup info, near equation II.8). However, while there are certainly situations in which this is undoubtedly true (for instance, in many of the lab microbial interactions studied empirically in this work), I am not necessarily convinced that such positive dependence on partner density hold generally among mutualistic interactions. After all, there is lots of evidence from natural systems of saturating benefits from increasing partner densities, and indeed also of shifts to negative effects beyond certain densities (e.g. Morris, Vázquez and Chacoff, 2010;Palmer and Brody, 2013 -> refs  The authors do mention previous work on saturating benefits in their supplementary information (in the section where they discuss previous models and again when they discuss their own modelling approach), for instance discussing some of the Holland & DeAngelis work on this, but as I understand it they don't incorporate any potential saturating effects (or even shifts to negative effects a highdensities) in their models. Of course, with some exceptions, we don't generally know to what extent real-life mutualisms typically experience (bounded) positive benefits from increased partner-densities as in the current work versus when saturating or even negative density effects start to appear. However, given that we have good evidence (and theoretical reasons) to think that they may commonly exist, this could undermine the claim of generality in this paper. Consequentially, I would like to see the authors either (i) include such effects in their models and analyse their impact, or (ii) convincingly explain and show why this concern is not relevant and their conclusions would hold even under saturating benefits, or (iii) scale back the claim of generality and clearly indicate that their results may not (all) apply to mutualisms with saturating effects, and are only directly applicable to the (much?) smaller subset of mutualisms with strict positive density effects.
We thank the reviewer for the insightful comments and suggested references. We understand the concerns about generality the reviewer raised, and sincerely appreciate having our attention drawn to this apparent weakness in our presentation. We clarify here how our work does take into account the effects the reviewer has mentioned, and also describe how we have worked to make this clearer in the manuscript.
We completely agree that accounting for saturating benefits in our models is important, and they are therefore implied in all our models. As asymptotic behaviors, we modeled the effect of benefit saturating with increasing partner density. In SI section II.D. we state that our models use Hill equations to capture the effect of benefit ( ) on partner fitness: ( or Fig. R5. For specific values, the benefit that one population receives from the other can saturate. In Fig. R5, when is not present ( ) the fitness of is 0 and -1 for the left and right panel respectively. When the density of increases, the fitness of increases and then asymptotically 13 plateaus. We have included this figure in SI section II. D, and also mention it in SI II.C line 146-148 to clarify the importance of saturating effects.
Similarly, in Extended Data Table 1, we have included models where cost has a constant effect, or the effect increases either linearly or in a saturating fashion with partner density (refer to SI section II.D.). For the dependent cases, cost is implemented by either: We have included this figure in SI section II. D to provide a visualization of how modulates partner fitness.
We can examine the overall effects of and by visualizing how our models capture the relationships between the density of and fitness. For example, model #79: captures the net effect as first saturating and then decreasing. Assuming , the instant growth rate of becomes: If we choose specific values for this function ( , , , ), we can easily see the overall fitness effect of both benefit and cost which can increase at low partner density and decrease at high partner density: 14 Using model #25, which differs from model #79 only in its constant effect of : Again, assuming , the instant growth rate of becomes: If we use the same set of parameters as in Fig. R7, the fitness of is: Fig. R8. When the effect of does not increase with partner density, fitness of is a monotonically increasing function of .
As the third example, a different shape can be generated using model #73: The instant growth rate of is: Using , , , , we get the following relationship:  In summary, we believe that the way we model benefit and cost is consistent with the suggestions made by the reviewer. To clarify this in the manuscript, we have edited the description of the model formulation and the rationale in the main text (page 4 line 73). As we described above, we have also edited the supplemental material section II.C, expanded SI section II.D to include Fig. R5 and Fig.  R6, and added a new section II.E to include Fig. R7-R9 to showcase different types of density vs. net effect curves our models capture.
I also have some more smaller remarks: -Please include line numbers throughout for a potential resubmission, including for the supplementary. Referencing sections of the manuscript is very cumbersome without them.
We have included line numbers on our main text and supplementary information.
-I think it would be important to highlight more explicitly in the main text that the models analysed here are all ecological models, and do not include potential evolutionary responses. It's fine to not include evolutionary dynamics, but it's important for the reader to be aware of this limitation, particularly given that eco-evolutionary dynamics could actually be important, particularly in many of the microbial mutualistic systems considered here.
We thank the reviewer for the suggestion. We now note in the main text (page 4, lines 64-65) that our models do not address evolutionary dynamics.
-Page 8 -Experimental Application of the metrics: for the first experimental application of the metric, I didn't fully understand what, if anything, the model equivalent of aTc is. The stress is imposed/modulated by IPTG, if I understand correctly, so is this just some external trigger of Quorum-sensing, independent of stress, that for whatever reason is required in this system? Would we not want QS to be triggered by the stress/stress-inducer itself? Some clarification for the reader would be helpful here.
We apologize for the lack of clarity in our presentation. We have now included additional clarification for our first experimental system. The reviewer did understand it correctly that QS is triggered by aTc and is independent of stress, which is induced by IPTG. We have clarified this point on page 9 line 188 and 190-191 of the main text.
-Main text table: I don't very much like cancer as an example of mutualism here. I know that the authors are using a slightly wider definition of mutualism here of two populations providing reciprocal benefits (first sentence of the ms), but most typically the term is used for situations where these populations are also of different species. I would suggest instead using an interspecific protection mutualism (e.g. ant-plant).
We thank the reviewer for the suggestion. We recognize and agree that mutualism has been traditionally used to describe interspecific interactions. As general background examples, we agree that the ant-plant mutualism is more fitting than cancer. As such, we have replaced the latter with the former in Table 1. -SI page 3 "In addition, although population collapse [..]the fitness of the benefit-receiver decreases with increasing partner density, which is contradictory to the basic logic of mutualism." -> I appreciate that given the definition of 'mutualism', at the stage where effects are negative the interaction is no longer strictly speaking a mutualism, but following my above general remark I wouldn't describe this effect as being contradictory to the basic logic of mutualism. The interaction at some density having a negative effect and then no longer being a mutualism is no way contradictory to it being a mutualist in other conditions, which I think what the models discussed here were trying to analyse.