LotkaVolterra pairwise modeling fails to capture diverse pairwise microbial interactions
 Cited 36
 Views 2,849
 Annotations
Abstract
Pairwise models are commonly used to describe manyspecies communities. In these models, an individual receives additive fitness effects from pairwise interactions with each species in the community ('additivity assumption'). All pairwise interactions are typically represented by a single equation where parameters reflect signs and strengths of fitness effects ('universality assumption'). Here, we show that a single equation fails to qualitatively capture diverse pairwise microbial interactions. We build mechanistic reference models for two microbial species engaging in commonlyfound chemicalmediated interactions, and attempt to derive pairwise models. Different equations are appropriate depending on whether a mediator is consumable or reusable, whether an interaction is mediated by one or more mediators, and sometimes even on quantitative details of the community (e.g. relative fitness of the two species, initial conditions). Our results, combined with potential violation of the additivity assumption in manyspecies communities, suggest that pairwise modeling will often fail to predict microbial dynamics.
https://doi.org/10.7554/eLife.25051.001eLife digest
From the soil to our body, microbes, such as bacteria, are everywhere and affect us in many ways. Many microbes perform important roles in natural environments and for our health, but some of them can cause harm and lead to diseases. Often, microbes affect and interact with each other within large groups or communities. Because of their widespread ramifications, it is important to understand how microbial communities work.
In addition to experiments, mathematical modeling offers one way to gain insight into the dynamics of microbial communities. A model commonly used to describe the interactions between organisms is the socalled ‘pairwise model’. Pairwise models can be useful to predict the dynamics of a community in which two species physically interact, such as a predatorprey community. However, it was unknown if this model was suitable to adequately predict the dynamics of microbial species in communities. Microbes often interact via chemicals that diffuse in the environment. For example, one microbe might provide food for another microbe or release toxins to kill it. However, a pairwise model does not consider food or toxins, but only how one microbe stimulates or inhibits the growth of another.
Momeni et al. simulated different scenarios commonly found in microbial communities to test whether a pairwise model could capture how, for example, chemicals released by one bacterial species would either help others to grow or stop them from growing. The results showed that for many scenarios, pairwise models cannot qualitatively represent the dynamics of a microbial community.
A next step will be to work on the limitations of current experimental technologies and mathematical models to improve the understanding of microbial communities. This knowledge could be used to develop new strategies for ecosystem engineering, such as for example making soils more fertile to improve crop yields, or tackling antibiotic resistance of bacteria.
https://doi.org/10.7554/eLife.25051.002Introduction
Multispecies microbial communities are ubiquitous. Microbial communities are important for industrial applications such as cheese and wine fermentation (van Hijum et al., 2013) and municipal waste treatment (Seghezzo et al., 1998). Microbial communities are also important for human health: they can modulate immune responses and food digestion (Round and Mazmanian, 2009; Kau et al., 2011) during health and disease. Properties of the entire community (‘community properties’, e.g. species dynamics, ability to survive internal or external perturbations, and biochemical activities of the entire community) are influenced by interactions wherein individuals alter the physiology of other individuals (Widder et al., 2016). To understand and predict community properties, choosing the appropriate mathematical model to describe species interactions is critical.
A mathematical model ideally focuses only on details that are essential to community properties of interest. However, it is often unclear a priori what the minimal essential details are. We define ‘mechanistic models’ as models that explicitly consider interaction mediators as state variables. For example, if species S_{1} releases a compound C_{1} which stimulates species S_{2} growth upon consumption by S_{2}, then a mechanistic model tracks concentrations of S_{1}, C_{1}, and S_{2} (Figure 1A and B, left panels). Note that mechanistic models used here still omit molecular details such as how chemical mediators are received and processed by recipients and how mediators subsequently act on recipients. In contrast, LotkaVolterra (‘LV’) pairwise models only consider the fitness effects of interactions. Specifically, LV models assume that the fitness of an individual is the sum of its basal fitness (the net growth rate of an individual in isolation) and fitness influences from pairwise interactions with individuals of the same species and of every other species in the community (‘additivity’ assumption). Furthermore, regardless of interaction mechanisms or quantitative details of a community, all fitness influences are typically expressed using a single equation form wherein parameters can vary to reflect the signs and magnitudes of fitness influences (‘universality’ assumption). Thus in the example above, a pairwise model only describes how S_{1} increases the fitness of S_{2} (Figure 1A and B, right panels).
LV pairwise models are popular. LV pairwise modeling has successfully explained the oscillatory dynamics of hare and its predator lynx (Figure 1—figure supplement 1) (Volterra, 1926; Wangersky, 1978; BiologyEOC, 2016). Pairwise models have also been instrumental in delineating conditions for multiple carnivores to coexist when competing for herbivores (MacArthur, 1970; Chesson, 1990). In both cases, mechanistic models and pairwise models happen to be mathematically equivalent for the following reasons. In the harelynx example, both species are also interaction mediators, and therefore pairwise and mechanistic models are identical. In the second example, if herbivores (mediators of competitive interactions between carnivores) rapidly reach steady state, herbivores can be mathematically eliminated from the mechanistic model to yield a pairwise model of competing carnivores (MacArthur, 1970; Chesson, 1990). Pairwise models are often used to predict how perturbations to steadystate species composition exacerbate or decline over time (May, 1972; Thébault and Fontaine, 2010; Mougi and Kondoh, 2012; Allesina and Tang, 2012; Suweis et al., 2013; Coyte et al., 2015). Although most work are motivated by contactdependent preypredation (e.g. harelynx) or mutualisms (e.g. plantpollinator) where LV models could be identical to mechanistic models, these work do not explicitly exclude chemicalmediated interactions where species are distinct from interaction mediators.
The temptation of using pairwise models is indeed high, including in microbial communities where many interactions are mediated by chemicals (Mounier et al., 2008; Faust and Raes, 2012; Stein et al., 2013; Marino et al., 2014; Coyte et al., 2015). Even though pairwise models do not capture the dynamics of chemical mediators, predicting species dynamics is still highly desirable in, for example, forecasting species diversity and compositional stability. For chemicalmediated interactions, LV pairwise models are far easier to construct than mechanistic models for the following reasons. Mechanistic models would require knowledge of chemical mediators, which are often challenging to identify. Since chemical mediators are explicitly modeled, mechanistic models require more equations and parameters than their cognate pairwise models (Figure 1, Table). Pairwise model parameters are relatively easy to estimate using community dynamics or dynamics of monocultures and pairwise cocultures (Mounier et al., 2008; Stein et al., 2013; Guo and Boedicker, 2016). Consequently, pairwise modeling has been liberally applied to microbial communities.
LV pairwise models have been criticized when applied to communities of more than two species (referred to as ‘multispecies communities’) (Levine, 1976; Tilman, 1987; Wootton, 1993, 2002; Werner and Peacor, 2003; Stanton, 2003; Schmitz et al., 2004). This is because a third species can influence interactions between a species pair (‘indirect interactions’), which sometimes violates the additivity assumption of pairwise models. For example, a carnivore can indirectly increase the density of a plant by decreasing the density of an herbivore (‘interaction chain’; ‘densitymediated indirect interactions’). A carnivore can also decrease how often an herbivore forages plants (‘interaction modification’, ‘traitmediated indirect interactions’, or ‘higher order interactions’) (Vandermeer, 1969; Wootton, 1994; Billick and Case, 1994; Wootton, 2002). In interaction modification, foraging per herbivore decreases, whereas in interaction chain, the density of herbivores decreases. Interaction modification (but not interaction chain) violates the additivity assumption (MethodsInteraction modification but not interaction chain violates the additivity assumption) (Tilman, 1987; Wootton, 1994; Schmitz et al., 2004) and can cause the pairwise model to generate qualitatively wrong predictions. Indeed, pairwise models largely failed to predict biomass and species coexistence in threespecies and sevenspecies plant communities (Dormann and Roxburgh, 2005), although reported failures of pairwise models could be due to limitations in data collection and analysis (Case and Bender, 1981; Billick and Case, 1994).
Here, we examine the universality assumption of pairwise models when applied to microbial communities (or any community that employs diverse chemicalmediated interactions). Microbes often influence other microbes in a myriad of fashions, via consumable metabolites, reusable signaling molecules, or a combination of chemicals (Figure 2). Can a single equation form, traditionally employed in pairwise models, qualitatively describe diverse interactions between two microbial species? The answer is unclear. On the one hand, pairwise models have been applied successfully to diverse microbial communities. For example, an LV pairwise model and a mechanistic model both correctly predicted ratio stabilization and spatial intermixing between two stronglycooperating populations exchanging diffusible essential metabolites (Momeni et al., 2013). In other examples, pairwise models largely captured competition outcomes and metabolic activities of threespecies and fourspecies artificial microbial communities (Vandermeer, 1969; Guo and Boedicker, 2016; Friedman et al., 2017). On the other hand, pairwise models often failed to predict species coexistence in sevenspecies microbial communities (Friedman et al., 2017), although this could be due to interaction modification discussed above.
Instead of investigating natural communities where interaction mechanisms can be difficult to identify, we use in silico communities. In these communities, two species interact via mechanisms commonly encountered in microbial communities, including growthpromoting and growthinhibiting interactions mediated by reusable and consumable compounds (Figure 2) (Stams, 1994; Czárán et al., 2002; Duan et al., 2009). We construct mechanistic models for these twospecies communities and attempt to derive from them pairwise models. A mechanistic reference model offers several advantages: community dynamics is deterministically known; deriving a pairwise model is not limited by inaccuracy of experimental methods; and the flexibility in creating different reference models allows us to explore a variety of interaction mechanisms. We demonstrate that a single pairwise equation form often fails for commonlyencountered diverse pairwise microbial interactions. We conclude by discussing when pairwise models might or might not be useful, in light of our findings.
Results
Throughout this work, we consider communities grown in a wellmixed environment where all individuals interact with each other with an equal chance. A wellmixed environment can be found in industrial fermenters. Moreover, at a sufficiently small spatial scale, a spatiallystructured environment can be approximated as a wellmixed environment, as chemicals are uniformlydistributed locally. Motile organisms also reduce the degree of spatial structure. A wellmixed environment allows us to use ordinary differential equations (ODEs), which are more tractable than partial differential equations demanded by a spatiallystructured environment. This in turn allows us to sometimes analytically demonstrate failures of pairwise models.
Mechanistic model versus pairwise model
A mechanistic model describes how species release or consume chemicals and how chemicals stimulate or inhibit species growth (Figure 1A left). In contrast, in pairwise models, interation mediators are not explicitly considered (Figure 1A right). Instead, the growth rate of an individual of species S_{i} is the sum of its basal fitness ($r}_{i0$, net growth rate of the individual in the absence of any intraspecies or interspecies interactions) and fitness effects from intraspecies and interspecies interactions. The fitness effect from species S_{j} to species S_{i} is represented by ${f}_{ij}\left({S}_{j}\right)$, where $S}_{j$ is the density of species S_{j}. ${f}_{ij}\left({S}_{j}\right)$ is a linear or nonlinear function of only $S}_{j$ and not of another species. When $j=i$, ${f}_{ii}\left({S}_{i}\right)$ represents densitydependent fitness effect within S_{i} (e.g. densitydependent growth inhibition or stimulation).
In a multispecies pairwise model, a single form of ${f}_{ij}$ is used for all pairwise species interactions. For example, the most popular LV model is linear LV:
Here, ${r}_{i0}$ is the basal fitness of an individual of S_{i}, and can be positive, negative, or zero; ${r}_{ij}$ is the fitness effect per S_{j} individual on S_{i}. Positive, negative, or zero ${r}_{ij}$ represents growth stimulation, inhibition, or no effect, respectively. An example of linear LV is the logistic LV pairwise model traditionally used for competitive communities:
Here, nonnegative ${r}_{i0}$ is the basal fitness of S_{i}; positive ${\Lambda}_{ij}$ is the carrying capacity imposed by limiting shared resource (e.g. space or carbon source) such that a single S_{i} individual will have a zero net growth rate when competing with a total of ${\Lambda}_{ij}$ individuals of S_{j}.
Alternative forms of fitness effect ${f}_{ij}$ (Wangersky, 1978) include LV with delayed influence, where the fitness influence of one species on another may lag in time (Gopalsamy, 1992), or saturable LV (Thébault and Fontaine, 2010) where
Here, ${r}_{i0}$ is the basal fitness of an individual of S_{i}, ${r}_{ij}$ is the maximal fitness effect species S_{j} can exert on S_{i}, and $K}_{ij$ (>0) is the $S}_{j$ at which half maximal fitness effect on S_{i} is achieved. ${r}_{i0}$ and ${r}_{ij}$ can be positive, negative, or zero. Note that at a low concentration of influencer, the saturable form can be converted to a linear form.
Our goal is to test whether a single equation form of pairwise model can qualitatively predict dynamics of species pairs engaging in various types of interactions commonly found in microbial communities (e.g. Figure 2). To do so, we use a combination of analytical and numerical approaches (Figure 1—figure supplement 2). Analytically deriving a pairwise model from a mechanistic model not only reveals assumptions required to generate the pairwise model, but also alleviates any concern that we may have failed to identify the optimal pairwise model parameters. When interactions become more complex (e.g. involving multiple mediators), we take the numerical approach, which is typically used to infer pairwise models from experimental results (Pascual and Kareiva, 1996). In the numerical approach, we mimic experimentalists by first deciding on a pairwise model to be used, and then employing a nonlinear least squares routine to numerically identify model parameters that minimize the average difference $\overline{D}$ between pairwise and mechanistic model dynamics within a training time window $T$ (Figure 1—figure supplement 2; MethodsSummary of simulation files). To evaluate how well a pairwise model predicts longterm mechanistic model dynamics, we ‘buy time’ by introducing 'dilutions' in numerical simulations of both models and quantify their difference $\overline{D}$.
Reusable versus consumable mediators require different pairwise models
In this section, we analytically derive pairwise models from mechanistic models of twospecies communities where one species affects the other species through a single mediator. The mediator is either reusable such as signaling molecules in quorum sensing (Duan et al., 2009; Jakubovics, 2010) or consumable such as metabolites (Stams, 1994; Freilich et al., 2011) (Figure 2). We show that a single pairwise model may not encompass these different interaction mechanisms and that for consumable mediator, the choice of pairwise model also depends on details such as the relative fitness and initial densities of the two species.
Consider a commensal community where species S_{1} stimulates the growth of species S_{2} by producing a reusable (Figure 3A) or a consumable (Figure 3B) chemical C_{1}. We consider community dynamics where species are not limited by any abiotic resources, such as within a dilution cycle of a turbidostat experiment where all other metabolites are in excess.
When C_{1} is reusable, the mechanistic model (Figure 3A,i) can be transformed into a saturable LV pairwise model (compare Figure 3A,ii with Equation 3), especially after the concentration of the mediator (which is initially zero) has acclimated to be proportional to the producer population size (Figure 3A legend; Figure 3—figure supplement 1). This saturable LV pairwise model is valid regardless of whether the producer coexists with the consumer, outcompetes the consumer, or is outcompeted by the consumer.
If C_{1} is consumable, different scenarios are possible (Figure 3B; Methods).
Case I: When supplier S_{1} always grows faster than consumer S_{2} (the basal fitness of S_{1} is higher than the maximal fitness of S_{2}), $C}_{1$ will eventually accumulate proportionally to $S}_{1$ (Figure 4A left; MethodsDeriving a pairwise model for interactions mediated by a single consumable mediator Case I). In this case, C_{1} may be approximated as a reusable mediator and can be predicted by the saturable LV pairwise model (Figure 4A right, compare dotted and solid lines).
Case II: When S_{1} and S_{2} can coexist (the basal fitness of S_{1} is higher than the basal fitness of S_{2} but less than the maximal fitness of S_{2}), a steady state solution for ${C}_{1}$ and species ratio ${R}_{S}={S}_{2}/{S}_{1}$ exists (Figure 4B; MethodsDeriving a pairwise model for interactions mediated by a single consumable mediator Case II, Equation 11). To arrive at a pairwise model, we will need to eliminate ${C}_{1}$ which is mathematically possible (i.e. after community dynamics converges to the ‘$f$zeroisocline’ on the phase plane of mediator ${C}_{1}$ and species ratio $R}_{S$, as depicted by blue lines in Figure 3—figure supplement 2A–D). However, the derived pairwise model differs from the saturable LV model:
where constants ${r}_{20}$, ${r}_{{S}_{2}{C}_{1}}$, and $\omega =1{K}_{{S}_{2}{C}_{1}}/{K}_{{C}_{1}{S}_{2}}$ can be positive, negative, or zero, and ψ =$\left({K}_{{S}_{2}{C}_{1}}{\alpha}_{{C}_{1}{S}_{2}}\right)/\left({K}_{{C}_{1}{S}_{2}}{\beta}_{{C}_{1}{S}_{1}}\right)$ is positive (see Figure 1 table for parameter definitions and see Equation 13 in Methods). We will refer to this equation as ‘alternative pairwise model’, although the fitness influence term is a function of both $S}_{1$ and $S}_{2$ instead of the influencer $S}_{1$ alone as defined in the traditional LV pairwise model.
Case III: When supplier S_{1} always grows slower than consumer S_{2}, i.e. when the basal fitness of S_{1} ($r}_{10$) is less than the basal fitness of S_{2} ($r}_{20$), consumable $C}_{1$ declines to zero concentration. This is because C_{1} is consumed by S_{2} whose relative abundance over S_{1} eventually exponentially increases at a rate of ${r}_{20}{r}_{10}$. Similar to Case II, under certain conditions (i.e. after community dynamics converges to the $f$zeroisocline as seen in Figure 3—figure supplement 2E–H), the alternative pairwise model (Equation 4) can be derived (MethodsDeriving a pairwise model for interactions mediated by a single consumable mediator, Case III).
For both Case II and Case III, we analytically demonstrate that in the absence of dilutions, alternative pairwise model dynamics can converge to mechanistic model dynamics (see Figure 3—figure supplements 3 and 5 for initial condition requirement and time scale of convergence). However, if initial $S}_{1$ and $S}_{2$ are such that the time scale of convergence is long compared to the duration of one dilution cycle (e.g. Figure 3—figure supplement 2C and G), then we will have to perform dilutions and the saturable LV model can sometimes be more appropriate than the alternative model (Figure 3—figure supplement 4). Thus in these cases, whether a saturable LV or an alternative model is more appropriate also depends on initial conditions.
The alternative model (Equation 4) can be further simplified to
if additionally, the halfsaturation constant $K$ for C_{1} consumption (${K}_{{C}_{1}{S}_{2}}$) is identical to that for C_{1}’s influence on the growth of consumer (${K}_{{S}_{2}{C}_{1}}$), and if S_{2} has not gone extinct. This equation form has precedence in the literature (e.g. [Mougi and Kondoh, 2012]), where the interaction strength $r}_{21$ reflects the fact that the consumable mediator from $S}_{1$ is divided among consumer $S}_{2$. Thus, we can regard the alternative model (Equation 4) or its simplified version (Equation 5) as a ‘divided influence’ model.
The saturable LV model and the alternative model are not interchangeable (Figure 4). When a consumable mediator accumulates without reaching a steady state within each dilution cycle (Figure 4A left; inset: $C}_{1$ eventually becomes proportional to $S}_{1$), the saturable LV model is predictive of community dynamics (Figure 4A right, compare dotted and solid lines). In contrast, predictions from the alternative pairwise model are qualitatively wrong (Figure 4A right, compare dashed and solid lines). When a consumable mediator eventually reaches a nonzero steady state within each dilution cycle (Figure 4B, black), could a saturable LV model still work? The saturable LV model derived from training window i (initial 10 generations) fails to predict species coexistence regardless of initial species ratios (Figure 4B left magenta box, compare solid with dotted). In comparison, the saturable LV model derived from training window ii (at steadystate species ratio) performs better, especially if the starting species ratio is identical to that of the training dynamics (Figure 4B, top panel in right magenta box). However, at a different starting species ratio, the saturable LV model fails to predict which species dominates the community (Figure 4B, bottom panel in right magenta box). In contrast, community dynamics can be described by the alternative pairwise model derived from either window i or ii (Figure 4B, compare dashed and solid lines in left and right magenta boxes).
We have shown here that even when one species affects another species via a single mediator, either a saturable LV model or an alternative pairwise model may be appropriate. The appropriate model depends on whether the mediator is reusable or consumable, how fitness of the two species compare, and initial species densities (Figure 3; Figure 3—figure supplements 2–5). Choosing the wrong pairwise model generates qualitatively flawed predictions (Figure 4). Considering that reusable and consumable mediators are both common in microbial interactions, our results call for revisiting the universality assumption of pairwise modeling.
Twomediator interactions require pairwise models different from singlemediator interactions
A species often affects another species via multiple mediators (Kato et al., 2008; Yang et al., 2009; Traxler et al., 2013; Kim et al., 2013). For example, a fraction of a population might die and release numerous chemicals, and some of these chemicals can simultaneously affect another individual. Here we examine the case where S_{1} releases two reusable chemicals C_{1} and C_{2}, both affecting the growth of S_{2} (Figure 5A). Since the effect of each mediator can be described by a saturable LV pairwise model (Figure 3A), we ask when the two mediators can be mathematically regarded as one mediator and described by a saturable LV pairwise model (Figure 5B).
We assume that fitness effects from different chemical mediators on a focal species are additive. Not making this assumption will likely violate the additivity assumption essential to pairwise models. Additive fitness effects have been observed for certain ‘homologous’ metabolites. For example, in multisubstrate carbonlimited chemostats of E. coli, the fitness effects from glucose and galactose were additive (Lendenmann and Egli, 1998). ‘Heterologous’ metabolites such as carbon and nitrogen sources likely affect cell fitness in a multiplicative fashion. However, if $W}_{C$ and $W}_{N$, the fitness influences of released carbon and nitrogen with respect to those already in the environment, are both small (i.e. $W}_{C$, $W}_{N$< < 1), the additional relative fitness influence will be additive: $(1+{W}_{C})(1+{W}_{N})1\approx {W}_{C}+{W}_{N}$. However, we need to keep in mind that even among homologous metabolites, fitness effects may not be additive (Hermsen et al., 2015). ‘Sequential’ metabolites (e.g. diauxic shift) provide another example of nonadditivity. Similar to the previous section, we assume that all abiotic resources are unlimited.
For the two reusable mediators, depending on their relative ‘potency’ (defined in Figure 5A legend), their combined effect generally cannot be modeled as a single mediator except under special conditions (MethodsConditions under which a saturable LV pairwise model can represent one species influencing another via two reusable mediators). These special conditions include: (1) mediators share similar potency (Figure 5—figure supplement 1B), or (2) one mediator dominates the interaction (Figure 5—figure supplement 1C). If these conditions are not satisfied, we can easily find examples where saturable LV pairwise models derived from a lowdensity community and from a highdensity community have qualitatively different parameters (Figure 5—figure supplement 1D). Consequently, the future dynamics of a lowdensity community can be predicted by a saturable LV model derived from a lowdensity community but not by a model derived from a highdensity community (Figure 5C). Thus, even though each mediator can be modeled by saturable LV, their joint effects may or may not be modeled by saturable LV depending on the relative potencies of the two mediators and sometimes even on initial conditions (high or low initial $S}_{1$).
Similarly, when both mediators are consumable and do not accumulate (as in Cases II and III of Figure 3B), the fitness effect term becomes $\frac{{r}_{{S}_{2}{C}_{1}}{S}_{1}}{{\omega}_{{C}_{1}}{S}_{1}+{\psi}_{{C}_{1}}{S}_{2}}+\frac{{r}_{{S}_{2}{C}_{2}}{S}_{1}}{{\omega}_{{C}_{2}}{S}_{1}+{\psi}_{{C}_{2}}{S}_{2}}$. Except under special conditions (e.g. when ${\omega}_{{C}_{1}}$ and ${\omega}_{{C}_{2}}$ are zero, or when ${\omega}_{{C}_{1}}/{\omega}_{{C}_{2}}={\psi}_{{C}_{1}}/{\psi}_{{C}_{2}}$, or when one mediator dominates the interaction), the two mediators may not be regarded as one. By the same token, when one mediator is a steadystate consumable and the other is reusable, they generally may not be regarded as a single mediator and would require yet a different pairwise model (i.e. with the fitness effect term $\frac{{r}_{{S}_{2}{C}_{1}}{S}_{1}}{{\omega}_{{C}_{1}}{S}_{1}+{\psi}_{{C}_{1}}{S}_{2}}+\frac{{r}_{{S}_{2}{C}_{2}}{S}_{1}}{{S}_{1}+{K}_{{S}_{2}{C}_{2}}{r}_{10}/{\beta}_{{C}_{2}{S}_{1}}}$).
In summary, when S_{1} influences S_{2} through multiple mediators, rarely can we approximate them as a single mediator. Sometimes, a pairwise model derived from one community may not apply to communities initiated at different densities (Figure 5C; Figure 5—figure supplement 1D). This casts further doubt on the usefulness of a single pairwise model for all pairwise microbial interactions.
LV competition model can fail if two competing species engage in an additional interaction
So far, by assuming that abiotic resources are always present in excess (e.g. in turbidostats), we have not considered species competition for abiotic resources. In this section, we consider a competitive commensal community in a batch environment where S_{1} and S_{2} compete for an essential shared resource C_{1} supplied by the environment at a constant rate (e.g. constant light), and S_{1} supplies an essential consumable metabolite C_{2} to promote S_{2} growth (Figure 6A, left). We show that an LV pairwise model works for some but not all communities even though these communities qualitatively share the same interaction mechanism.
In our mechanistic model (MethodsCompetitive commensal interaction, Equation 47), the fitness of S_{2} is multiplicatively affected by C_{1} and C_{2} (Mankad and Bungay, 1988). We choose parameters such that the effect from C_{2} to S_{2} is far from saturation (e.g. linear with respect to $C}_{2$ and $S}_{1$) to simplify the problem. In our LV pairwise model (Figure 6A, right; MethodsCompetitive commensal interaction, Equation 48), intra and interspecies competition is represented by the traditional logistic LV model (Equation 2; Gause, 1934; Thébault and Fontaine, 2010; Mougi and Kondoh, 2012). We then introduce a linear term ($r}_{21}{S}_{1$) to describe the fitness effect of commensal interaction.
We tested various sets of mechanistic model parameters where the two species coexist in a steady fashion (Figure 6B), or one species goes extinct (Figure 6C), or species composition fluctuates (Figure 6D). LV pairwise models deduced from a fixed period of training time could predict future dynamics in the first two cases, but failed to do so in the third case. Thus, depending on dynamic details of communities, a pairwise model sometimes works and sometimes fails.
To summarize our work, even for pairwise microbial interactions, depending on interaction mechanisms (reusable versus consumable mediator, single mediator versus multiple mediators), we will need to use a plethora of pairwise models to avoid qualitative failures in predicting which species dominates a community or whether species coexist (Figures 3, 4 and 5). Sometimes, even when different communities share identical interaction mechanisms, depending on details such as relative species fitness, interaction strength, and initial conditions, the bestfitting pairwise model may or may not predict future dynamics (Figure 3B, Figure 3—figure supplement 4, Figure 4, Figure 5, Figure 5—figure supplement 1, and Figure 6). This defeats the very purpose of pairwise modeling – using a single equation form to capture fitness effects of all pairwise species interactions regardless of interaction mechanisms or quantitative details. In a community of more than two microbial species, interaction modification can cause pairwise models to fail (Figure 7). Even if species interact in an interaction chain and thus interaction modification does not occur, various chain segments may require different forms of pairwise models. Taken together, a pairwise model is unlikely to be effective for predicting community dynamics especially if interaction mechanisms are diverse.
Discussions
Multispecies pairwise models are widely used in theoretical ecology due to their simplicity. These models assume that all pairwise species interactions can be captured by a single pairwise model regardless of interaction mechanisms or quantitative details of a community (universality assumption). This assumption may be satisfied if, for example, interaction mediators are always species themselves (e.g. preypredation in a food web) so that pairwise models are equivalent to mechanistic models. However, interactions in microbial communities are diverse and often mediated by chemicals (Figure 2). Here, we consider the validity of universality assumption of pairwise models in wellmixed, twospecies microbial communities. We have focused on various types of chemicalmediated interactions commonly encountered in microbial communities (Figure 2) (Kato et al., 2005; Gause, 1934; Ghuysen, 1991; Jakubovics et al., 2008; Chen et al., 2004; D'Onofrio et al., 2010; Johnson et al., 1982; Hamilton and Ng, 1983). For each type of species interaction, we construct a mechanistic model to generate reference community dynamics (akin to experimental results). We then attempt to derive the bestmatching pairwise model and ask how predictive it is.
We first consider cases where abiotic resources are in excess. When one species affects another species via a single chemical mediator, either the saturable LV or the alternative pairwise model is appropriate, depending on the interaction mechanism (consumable versus reusable mediator), relative fitness of the two species, and initial conditions (Figure 3; Figure 3—figure supplement 2 to Figure 3—figure supplement 5). These two models are not interchangeable (Figure 4). If one species influences another species through multiple mediators, then in general, these mediators may not be regarded as a single mediator nor representable by a single pairwise model. For example, for two reusable mediators, unless their potencies are similar or one mediator is much more potent than the other, saturable LV model parameters can qualitatively differ depending on initial community density (Figure 5—figure supplement 1D). Consequently, a pairwise model derived from a highdensity community generates false predictions for lowdensity communities (Figure 5C), limiting the usefulness of pairwise models. We then consider a community where two species compete for a shared resource while engaging in commensalism via a single chemical mediator. We find that the bestfitting LV pairwise model can predict future dynamics in some but not all communities, depending on parameters used in the mechanistic model (Figure 6). Thus, although a single equation form can work in many cases, it generates qualitatively wrong predictions in many other cases.
In communities of more than two microbial species, indirect interactions via a third species can occur. When indirect interactions take the form of interaction chains, if each chain segment of two species engages in an independent interaction and can be represented by a pairwise model, then multispecies pairwise models can work (Figure 7AB). However, as discussed above, pairwise equation forms may vary among chain segments depending on interaction mechanisms and quantitative details of a community. When indirect interactions take the form of interaction modification, even if each species pair can be accurately represented by a pairwise model, a multispecies pairwise model may fail (Figure 7C–F, ). Interaction modification includes trait modification (Wootton, 2002; Werner and Peacor, 2003; Schmitz et al., 2004), or, in our cases, mediator modification. Mediator modification is very common in microbial communities. For example, antibiotic released by one species to inhibit another species may be inactivated by a third species, and this type of indirect interactions can stabilize microbial communities (Kelsic et al., 2015; Bairey et al., 2016). As another example, interaction mediators are often generated by and shared among multiple species. For example in oral biofilms, organic acids such as lactic acid are generated from carbohydrate fermentation by many species (Bradshaw et al., 1994; Marsh and Bradshaw, 1997; Kuramitsu et al., 2007). Such byproducts are also consumed by multiple species (Kolenbrander, 2000).
One can argue that an extended pairwise model (e.g. $\frac{d{S}_{2}}{dt}={r}_{20}{S}_{2}+\frac{{r}_{{S}_{2}C}{S}_{1}}{\varsigma +\omega {S}_{1}+\psi {S}_{2}}{S}_{2}$) embodying both the saturable form and the alternative form can serve as a generalpurpose model at least for pairwise interactions via a single mediator. In fact, even the effects of indirect interactions may be quantified and included in the model by incorporating higherorder interaction terms (Case and Bender, 1981; Worthen and Moore, 1991), although with many challenges (Wootton, 2002). In the end, although these strategies may lead to a sufficiently accurate phenomenological model especially within the training window, they may fail to predict future dynamics.
When might a pairwise model be useful? First, pairwise models have been instrumental in understanding ecological phenomena such as preypredator oscillatory dynamics and coexistence of competing predator species (Volterra, 1926; MacArthur, 1970; Case and Casten, 1979; Chesson, 1990). In these cases, mechanistic models are either identical to pairwise models or can be transformed into pairwise models under simplifying assumptions. Second, pairwise models of pairwise species interactions can provide a bird’seye view of strong or weak stimulatory or inhibitory interactions in a community. For example, Vetsigian et al., 2011 found that interactions between soilisolated Streptomyces strains are enriched for reciprocity – if A inhibits or promotes B, it is likely that B also inhibits or promotes A (Vetsigian et al., 2011). Third, pairwise models have been useful in qualitatively understanding species assembly rules in small communities (Friedman et al., 2017). That is, qualitative information regarding species survival in competitions among a small number of species may be used to predict survival in more diverse communities within a similar time window. Fourth, a pairwise model can serve as a starting point for generating hypotheses on species interactions (e.g. Li et al., 2015). Note that when applied to microbial communities (Mounier et al., 2008; Stein et al., 2013; Marino et al., 2014), a fitting pairwise model means that the training dynamics of the community under investigation can be approximated by a theoretical community where species interactions satisfy the additivity and universality assumptions of pairwise models. Even though the theoretical community is likely different from the real community, hypothesis formulation can still be valuable. Finally, pairwise models can be useful in making predictions of limited scales. For example, Stein et al. used 2/3 of community dynamics data as a training set to derive a multispecies pairwise model, and in the bestcase scenario, the model generated reasonable predictions on the remaining 1/3 of data (Stein et al., 2013). However, as we have shown, pairwise models can generate qualitatively wrong predictions (Figures 4–7), especially if interaction mechanisms are diverse such as in microbial communities. Not surprisingly, predicting qualitative consequences of species removal or addition using a pairwise model has encountered difficulties, especially in communities of more than three species (Mounier et al., 2008; Friedman et al., 2017).
An alternative to a pairwise model is a mechanistic model. How much information about interaction mechanisms do we need to construct a mechanistic model? That is, what is the proper level of abstraction which captures the phenomena of interest, yet avoids unnecessary details (Li et al., 2015; Durrett and Levin, 1994)? For example, Tilman argued that if a small number of mechanisms (e.g. the ‘axes of tradeoffs’ in species traits) could explain much of the observed pattern (e.g. species coexistence), then this abstraction would be highly revealing (Tilman, 1987). However, the choice of abstraction is often not obvious. Consider for example a commensal community where S_{1} grows exponentially (not explicitly depicted in equations in Figure 8) and the net growth rate of S_{2}, which is normally zero, is promoted by mediator C from S_{1} in a linear fashion (Figure 8). If we do not know how S_{1} stimulates S_{2}, we can still construct an LV pairwise model (Figure 8A). If we know the identity of mediator C and realize that C is consumable, then we can instead construct a mechanistic model incorporating $C$ (Figure 8B). However, if C is produced from a precursor via an enzyme E released by S_{1}, then we get a different form of mechanistic model (Figure 8C). If, on the other hand, E is anchored on the membrane of S_{1} and each cell expresses a similar amount of $E$, then equations in Figure 8D are mathematically equivalent to Figure 8B. This simple example, inspired by extracellular breakdown of cellulose into a consumable sugar C (Bayer and Lamed, 1986; Felix and Ljungdahl, 1993; Schwarz, 2001), illustrates how knowledge of mechanisms may eventually help us determine the right level of abstraction.
In summary, under certain circumstances, we may already know that microbial interaction mechanisms fall within the domain of validity for a particular pairwise model. In these cases, a pairwise model provides the appropriate level of abstraction, and constructing such a pairwise model is much easier than a mechanistic model (Figure 1). However, if we do not know whether a pairwise model is valid, we will need to be cautious since pairwise models can fail to even qualitatively capture pairwise microbial interactions. We need to be equally careful when extrapolating and generalizing conclusions obtained from pairwise models, especially for communities where species interaction mechanisms are diverse. Considering recent advances in identifying and quantifying interactions, we advocate a transition to models that incorporate interaction mechanisms at the appropriate level of abstraction.
Materials and methods
Interaction modification but not interaction chain violates the additivity assumption
In a pairwise model, the fitness of a focal species S_{i} is the sum of its ‘basal fitness’ ($r}_{i0$, the net growth rate of a single S_{i} individual in the absence of any intraspecies or interspecies interactions) and the additive fitness effects exerted by pairwise interactions with other members of the community. Mathematically, an $N$species pairwise model is often formulated as
Here, ${f}_{ij}\left({S}_{j}\right)$ describes how $S}_{j$, the density of species S_{j}, positively or negatively affects the fitness of S_{i}, and is a linear or nonlinear function of only $S}_{j$.
Indirect interactions via a third species fall under two categories (Wootton, 1993). The first type is known as ‘interaction chain’ or ‘densitymediated indirect interactions’. For example, the consumption of plant S_{1} by herbivore S_{2} is reduced when the density of herbivore is reduced by carnivore S_{3}. In this case, the threespecies pairwise model
does not violate the additivity assumption (compare with Equation 6 (Case and Bender, 1981; Wootton, 1994).
The second type of indirect interactions is known as ‘interaction modification’ or ‘traitmediated indirect interactions’ or ‘higher order interactions’ (Vandermeer, 1969; Wootton, 1994; Billick and Case, 1994; Wootton, 2002), where a third species modifies the ‘nature of interaction’ from one species to another (Wootton, 2002; Werner and Peacor, 2003; Schmitz et al., 2004). For example, when carnivore is present, herbivore will spend less time foraging and consequently plant density increases. In this case, $f}_{12$ in Equation 7 is a function of both $S}_{2$ and $S}_{3$, violating the additivity assumption.
Summary of simulation files
Simulations are based on Matlab and executed on an ordinary PC. Steps are:
Step 1: Identify monoculture parameters $r}_{i0$, $r}_{ii$, and $K}_{ii$ (Figure 1—figure supplement 2C, Row 1 and Row 2).
Step 2: Identify interaction parameters $r}_{ij$, $r}_{ji$, $K}_{ij$, and $K}_{ji$ where $i\ne j$ (Figure 1—figure supplement 2C, Row 3).
Step 3: Calculate distance $\overline{D}$ between population dynamics of the reference mechanistic model and the approximate pairwise model over a period of time outside of the training window to assess if the pairwise model is predictive.
Fitting is performed using nonlinear least squares (lsqnonlin routine) with default optimization parameters. The following list describes the mfiles used for different steps of the analysis:
File name  Function 

FitCost_BasalFitness  Calculates the cost function for monocultures (i.e. the difference between the target mechanistic model dynamics and the dynamics obtained from the pairwise model) 
FitCost_BFSatLV.m  Calculates the cost function for communities (i.e. the difference between the target mechanistic model dynamics and the dynamics obtained from the saturable LV pairwise model) 
FitCost_BFSatLV_Dp.m  Calculates the cost function for communities (i.e. the difference between the target mechanistic model dynamics and the dynamics obtained from the alternative pairwise model) 
DynamicsMM_WM_MonocultureDpMM.m  Returns growth dynamics for monocultures, based on the mechanistic model 
DynamicsMMSS_WM_NetworkDpMM.m  Returns growth dynamics for communities of multiple species, based on the mechanistic model 
DynamicsWM_NetworkBFSatLV.m  Returns growth dynamics for communities of multiple species, based on the saturable LV pairwise model 
DynamicsWM_NetworkBFSatLV_Dp.m  Returns growth dynamics for communities of multiple species, based on the alternative pairwise model 
DeriveBasalFitnessMM_WM_DpMM.m  Estimates monoculture parameters of pairwise model (Step 1) 
DeriveBFSatLVMMSS_WM_DpMM.m  Estimates saturable LV pairwise model interaction parameters (Step 2) 
DeriveBFSatLVMMSS_WM_DpMM_Dp.m  Estimates alternative pairwise model interaction parameters (Step 2) 
DeriveBFSatLVMMSS_WM_DpMM_r21.m  Estimates saturable LV pairwise model interaction parameters ($r}_{21$ and $K}_{21$) in cases where we know that S_{2} is only affected by S_{1}, to accelerate optimization 
DeriveBFSatLVMM_WM_DpMM_Dp_r21.m  Estimates alternative pairwise model interaction parameter ($r}_{21$) in cases where we know that S_{2} is only affected by S_{1} and that $K}_{S2C1}={K}_{C1S2$ to accelerate optimization 
DynamicsWM_NetworkBFLogLV_DI.m  Returns growth dynamics for communities of two species competing for an environmental resource while engaging in an additional interaction, based on the logistic LV pairwise model (Figure 6) 
C2Sp2_ARCLi_NoSatDp_FitBFLogLV_DI.m  Estimates logistic LV pairwise model interaction parameters for communities of two species competing for an environmental resource while engaging in an additional interaction, and compares community dynamics from pairwise and mechanistic models (Figure 6) 
Dynamics_WM_NetworkDpMM_ODE23.m  Defines differential equations when using Matlab’s ODE23 solver to calculate community dynamics 
Case_C1Sp2_CmnsDp_ODE23.m  Example of using Matlab ODE23 solver for calculating community dynamics 
Deriving a pairwise model for interactions mediated by a single consumable mediator
To facilitate mathematical analysis, we assume that requirements calculated below are eventually satisfied within each dilution cycle (see Figure 3—figure supplement 4 for an example where dilution cycles necessitated by long convergence time violate requirements for a pairwise model to converge to the mechanistic model). We further assume that ${r}_{10}>0$ and ${r}_{20}>0$ so that species cannot go extinct in the absence of dilution. See Figure 3—figure supplement 5 for a summary of this section.
When S_{1} releases a consumable mediator which stimulates the growth of S_{2}, the mechanistic model as per Figure 3B, is
Let ${C}_{1}(t=0)={C}_{10}=0$; $S}_{1}(t=0)={S}_{10$; and $S}_{2}(t=0)={S}_{20$. Note that the initial condition ${C}_{10}=0$ can be easily imposed experimentally by prewashing cells. Under which conditions can we eliminate $C}_{1$ so that we can obtain a pairwise model of $S}_{1$ and $S}_{2$?
Define ${R}_{S}={S}_{2}/{S}_{1}$ as the ratio of the two populations.
Case I: ${r}_{10}{r}_{20}>{r}_{{S}_{2}{C}_{1}}$
Since producer S_{1} always grows faster than consumer S_{2}, ${R}_{S}\to 0$ as $t\to \infty $. Define $\stackrel{~}{C}}_{1}={C}_{1}/{S}_{1$ (‘~' indicating scaling against a function).
Since $R}_{S$ declines exponentially with a rate faster than $\left{r}_{20}+{r}_{{S}_{2}{C}_{1}}{r}_{10}\right$, we can ignore the third term of the right hand side of Equation 10 if it is much smaller than the first term. That is,
Thus for $t\gg \mathrm{ln}\left(\frac{{\alpha}_{{C}_{1}{S}_{2}}{R}_{S}\left(0\right)}{{\beta}_{{C}_{1}{S}_{1}}}\right)/{r}_{20}+{r}_{{S}_{2}{C}_{1}}{r}_{10}$, $\frac{d{\stackrel{~}{C}}_{1}}{dt}\approx {\beta}_{{C}_{1}{S}_{1}}{r}_{10}{\stackrel{~}{C}}_{1}$. When initial $\stackrel{~}{{C}_{1}}$ is 0, this equation can be solved to yield: $\stackrel{~}{{C}_{1}}\approx {\beta}_{{C}_{1}{S}_{1}}\left(1\mathrm{exp}\left({r}_{10}t\right)\right)/{r}_{10}$. After time of the order of $1/{r}_{10}$, the second term can be neglected. Thus, $\stackrel{~}{{C}_{1}}\approx {\beta}_{{C}_{1}{S}_{1}}/{r}_{10}$ after time of the order of $max(\mathrm{ln}\left(\frac{{\alpha}_{{C}_{1}{S}_{2}}{R}_{S}\left(0\right)}{{\beta}_{{C}_{1}{S}_{1}}}\right)/{r}_{20}+{r}_{{S}_{2}{C}_{1}}{r}_{10},1/{r}_{10})$. Then ${C}_{1}$ can be replaced by $\left({\beta}_{{C}_{1}{S}_{1}}/{r}_{10}\right){S}_{1}$ in Equation 8, and a saturable LV pairwise model can be derived.
Case II: ${r}_{{S}_{2}{C}_{1}}>{r}_{10}{r}_{20}>0$
For Equation 8, we find that a steady state solution for ${C}_{1}$ and ${R}_{S}$, denoted respectively as $C}_{1}^{\ast$ and $R}_{S}^{\ast$, exist. They can be easily found by setting the growth rates of S_{1} and S_{2} to be equal, and $d{C}_{1}/dt$ to zero.
However, if $C}_{1$ has not yet reached steady state, imposing steady state assumption would falsely predict $R}_{S$ at steady state and thus remaining at its initial value (Figure 4Bii , dotted lines). Since $d{C}_{1}/dt$ in Equation 8 is the difference between two exponentially growing terms, we factor out the exponential term $S}_{1$ to obtain
where $f\left({C}_{1},{R}_{S}\right)=1\frac{{\alpha}_{{C}_{1}{S}_{2}}}{{\beta}_{{C}_{1}{S}_{1}}}\frac{{C}_{1}}{{C}_{1}+{K}_{{C}_{1}{S}_{2}}}{R}_{S}$. When $f\approx 0$, we can eliminate $C}_{1$ and obtain an alternative pairwise model
Or
where ω and ψ are constants (Figure 3Bii).
For certain conditions (which will be discussed at the end of this section, Figure 3—figure supplement 5A), this alternative model can make reasonable predictions of community dynamics even before the community reaches the steady state (Figure 4Bii , compare dashed and solid lines). Below we discuss the general properties of community dynamics and show that there exists a time scale $t}_{f$ after which it is reasonable to assume $f\approx 0$ and the alternative model can be derived. We also estimate $t}_{f$ for several scenarios.
We first make $C}_{1$ and $R}_{S$ dimensionless by defining $\hat{C}}_{1}={C}_{1}/{C}_{1}^{\ast$ and $\hat{R}}_{S}={R}_{S}/{R}_{S}^{\ast$ (‘ ^ ' indicating scaling against steady state values). Equation 9 can then be rewritten as
where $\hat{K}}_{{S}_{2}{C}_{1}}={K}_{{S}_{2}{C}_{1}}/{C}_{1}^{\ast$.
From Equations 8 and 11, we obtain
or
where $\hat{\beta}}_{{C}_{1}{S}_{1}}={\beta}_{{C}_{1}{S}_{1}}/{C}_{1}^{\ast$ and $\hat{K}}_{{C}_{1}{S}_{2}}={K}_{{C}_{1}{S}_{2}}/{C}_{1}^{\ast$.
Using these scaled variables, $f$ (i.e. the square bracket in Equation 15) can be rewritten as
and
Equations 14 and 17 allow us to construct a phase portrait where the x axis is $\hat{C}}_{1$ and the y axis is $\hat{R}}_{S$ (Figure 3—figure supplement 2A–D). Note that at steady state, $\left({\hat{C}}_{1},\text{}{\hat{R}}_{S}\right)=(1,1)$. Setting Equation 16 to zero:
defines the $f$zeroisocline on the $\hat{C}}_{1}{\hat{R}}_{s$ phase plane (i.e. values of $\left({\hat{C}}_{1},\text{}{\hat{R}}_{S}\right)$ at which $f\left({\hat{C}}_{1},{\hat{R}}_{S}\right)=0$ and thus $\hat{C}}_{1$ can be eliminated to obtain a pairwise model; Figure 3—figure supplement 2A–D blue lines). As shown in Figure 3—figure supplement 2A, the phase portrait is divided into four regions by the $f$zeroisocline (blue) and the steady state ${\hat{C}}_{1}=1$ (vertical solid line), and grey arrows dictate the direction of the community dynamics trajectory ($\hat{C}}_{1$, $\hat{R}}_{S$). Starting from 'initial state' (${\hat{C}}_{1}\text{(}t=0\text{)}=0$, ${\hat{R}}_{S}\left(t=0\right)$), the trajectory moves downward right (brown circles and orange lines in Figure 3—figure supplement 2A–D) until it hits ${\hat{C}}_{1}=1$. Then, it moves upward right and eventually hits the $f$zeroisocline. Afterward, the trajectory moves toward the steady state (green circles) very closely along (and not superimposing) the $f$zeroisocline during which the alternative pairwise model can be derived (Figure 3—figure supplement 2A–D).
It is difficult to solve Equations 14 and 15 analytically because the detailed community dynamics depends on the parameters and the initial species composition in a complicated way. However, under certain initial conditions, we can estimate ${t}_{f}$, the time scale for the community to approach the $f$zeroisocline. Note that ${t}_{f}$ is not a precise value. Instead it estimates the acclimation time scale after which a pairwise model can be derived.
One assumption used when estimating all ${t}_{f}$ is that ${S}_{10}$ is sufficiently high (Figure 3—figure supplement 5B) to avoid the long lag phase that is otherwise required for the mediator to accumulate to a high enough concentration.
From Equation 18, the asymptotic value for the $f$zeroisocline is
This is plotted as a black dotted line in Figure 3—figure supplement 2A–D.
Below we consider three different initial conditions for ${\hat{R}}_{S}\left(t=0\right)$:
Case II1. ${\hat{R}}_{S}\left(t=0\right)\gg max(1,{\hat{K}}_{{S}_{2}{C}_{1}}^{1})$
From Equation 11, this becomes ${R}_{S}\left(0\right)/{R}_{S}^{\ast}\gg max(1,\frac{{r}_{10}{r}_{20}}{{r}_{20}+{r}_{{S}_{2}{C}_{1}}{r}_{10}})$.
A typical trajectory of the system is shown in Figure 3—figure supplement 2B: at time $t=0$, using Equations 14 and 15, the community dynamics trajectory (orange solid line in Figure 3—figure supplement 2B inset) has a slope of
From Equation 18, the slope of the $f$zeroisocline (blue line in Figure 3—figure supplement 2B inset) at ${\hat{R}}_{S}={\hat{R}}_{S}\left(0\right)$ is
The approximation in the last step is due to the very definition of Case II1: ${\hat{R}}_{S}\left(t=0\right)\gg 1$. The initial steepness of the community dynamics trajectory (Equation 20) will be much smaller than that of the $f$zeroisocline (Equation 21) if
If we do not scale, together with Equation 11, this becomes:
In this case, the community dynamics trajectory before getting close to the $f$zeroisocline can be approximated as a straight line (the orange dotted line) and the change in $\hat{R}}_{S$ can be approximated by the green segment in the inset of Figure 3—figure supplement 2B. Since the green segment, the orange dotted line and the red dashed line form a right angle triangle, the length of green segment can be calculated once we find the length of the red dashed line $\mathrm{\Delta}{\hat{C}}_{1}$, which is the horizontal distance between (${\hat{C}}_{1}\left(0\right)$, ${\hat{R}}_{S}\left(0\right)$) and the $f$zeroisocline and can be calculated from Equation 16:
which yields
The green segment $\mathrm{\Delta}{\hat{R}}_{S}$ is then the length of red dashed line ($\mathrm{\Delta}{\hat{C}}_{1}$, Equation 24 ) multiplied with $\frac{d{\hat{R}}_{S}}{d{\hat{C}}_{1}}}_{{\hat{C}}_{1}\left(0\right)=0$ (Equation 20), or
Note that if Equation 22 is satisfied, $\mathrm{\Delta}{\hat{R}}_{S}\ll {\hat{R}}_{S}\left(0\right)$. What is the time scale ${t}_{f}$ for the community to traverse the orange dotted line to be close to the $f$zeroisocline? Since from Equation 14 $\frac{d{\hat{R}}_{S}}{dt}=\left({r}_{20}+{r}_{{S}_{2}{C}_{1}}\frac{{\hat{C}}_{1}}{{\hat{C}}_{1}+{\hat{K}}_{{S}_{2}{C}_{1}}}{r}_{10}\right){\hat{R}}_{S}\approx \left({r}_{20}{r}_{10}\right){\hat{R}}_{S}$. In Equation 14, the second term in the parenthese can be dropped if
In case II1, before the system reaches the $f$zeroisocline, from Equation 24, ${\hat{C}}_{1}\le \mathrm{\Delta}{\hat{C}}_{1}<1/{\hat{R}}_{S}(0)$ thus
From the top portion of Equation 11,
thus
According to the condition, ${\widehat{R}}_{S}\left(0\right)\gg max\left(1,\underset{{S}_{2}{C}_{1}}{\overset{1}{\widehat{K}}}\right)$. If ${\widehat{K}}_{{S}_{2}{C}_{1}}^{1}>1$, then ${\widehat{R}}_{S}\left(0\right)\gg {\widehat{K}}_{{S}_{2}{C}_{1}}^{1}$, ${\widehat{R}}_{S}\left(0\right){\widehat{K}}_{{S}_{2}{C}_{1}}\gg 1$ and ${\hat{K}}_{{S}_{2}{C}_{1}}<1$.
If ${\hat{K}}_{{S}_{2}{C}_{1}}^{1}<1$, then ${\hat{R}}_{S}(0)\gg 1$
Thus, the above approximation of Equation 14 is valid, and we obtain
${t}_{f}\approx \mathrm{ln}\left(\frac{{\hat{R}}_{S}\left(0\right)+\mathrm{\Delta}{\hat{R}}_{S}}{{\hat{R}}_{S}\left(0\right)}\right)/\left({r}_{20}{r}_{10}\right)$.
Since here $\mathrm{\Delta}{\hat{R}}_{S}\ll {\hat{R}}_{S}\left(0\right)$ and $\mathrm{ln}\left(1+x\right)\sim x$ for small $x$, together with Equation 25, we have
If unscaled, using Equation 11, this becomes
Case II2. ${\hat{R}}_{S}\left(t=0\right)$ is comparable to 1
That is, ${R}_{S}\left(t=0\right)\approx {R}_{S}{}^{*}$. If ${S}_{10}$ is low, a typical example is shown in Figure 3—figure supplement 2D. Here, because it takes a while for $C}_{1$ to accumulate, during this lagging phase ${\hat{R}}_{S}\left(t\right)\approx {\hat{R}}_{S}\left(0\right)\mathrm{exp}\left({r}_{10}{r}_{20}t\right)$ and there is a sharp plunge in $\hat{R}}_{S$ before the trajectory levels off and climbs up. Although the trajectory eventually hits the $f$zeroisocline where the alternative pairwise model can be derived, estimating ${t}_{f}$ is more complicated. Here we consider a simpler case where ${S}_{10}$ is large enough so that the trajectory levels off immediately after $t=0$, and ${\hat{R}}_{S}\approx 1$ before the trajectory hits the $f$zeroisocline (Figure 3—figure supplement 2A). Since $\hat{R}}_{S$ decreases until ${\hat{C}}_{1}=1$ and from Equation 20, and similar to the reasoning in Case II1, if
or if
a typical trajectory moves toward the $f$zeroisocline almost horizontally (Figure 3—figure supplement 2A). The unscaled form of Equation 28 is
To calculate the time it takes for the trajectory to reach the $f$zeroisocline, let ${\mathrm{\Delta}}_{s}{\hat{C}}_{1}={\hat{C}}_{1}1$ and ${\mathrm{\Delta}}_{s}{\hat{R}}_{S}={\hat{R}}_{S}1$ at any time point $t$ respectively represent deviation of (${\hat{C}}_{1}\left(t\right)$, ${\hat{R}}_{S}\left(t\right)$) away from their steady state values of (1, 1). We can thus linearize Equation 14 and Equation 15 around the steady state. Note that since at the steady state $f$= 0, thus ${\Delta}_{s}f=f$.
Rewrite Equation 14 as
$\frac{d{\hat{R}}_{S}}{dt}=\left({r}_{20}+{r}_{{S}_{2}{C}_{1}}\frac{{\hat{C}}_{1}}{{\hat{C}}_{1}+{\hat{K}}_{{S}_{2}{C}_{1}}}{r}_{10}\right){\hat{R}}_{S}=h\left({\hat{C}}_{1},{\hat{R}}_{S}\right)$.
We linearize this equation around the steady state ${\hat{C}}_{1}=1,\text{}\text{}{\hat{R}}_{S}=1$
$\frac{d\left(1+{\mathrm{\Delta}}_{s}{\hat{R}}_{S}\right)}{dt}=h\left(1+{\mathrm{\Delta}}_{s}{\hat{C}}_{1},1+{\mathrm{\Delta}}_{s}{\hat{R}}_{S}\right)\approx h\left(1,1\right)+{\mathrm{\Delta}}_{s}{\hat{C}}_{1}\frac{\mathrm{\partial}h}{\mathrm{\partial}{\hat{C}}_{1}}{}_{\left({\hat{C}}_{1}=1,{\hat{R}}_{S}=1\right)}+{{\mathrm{\Delta}}_{s}{\hat{R}}_{S}\frac{\mathrm{\partial}h}{\mathrm{\partial}{\hat{R}}_{S}}}_{\left({\hat{C}}_{1}=1,{\hat{R}}_{S}=1\right)}$.
At steady state, $\frac{d{\hat{R}}_{S}}{dt}=h\left(1,1\right)=0$. Thus, $r}_{20}+\frac{{r}_{{S}_{2}{C}_{1}}}{1+{\hat{K}}_{{S}_{2}{C}_{1}}}{r}_{10$=0.
Thus,
Recall Equations 15 and 17 as
$\frac{d{\hat{C}}_{1}}{dt}={\hat{\beta}}_{{C}_{1}{S}_{1}}\left(1\frac{{\hat{C}}_{1}\left(1+{\hat{K}}_{{C}_{1}{S}_{2}}\right)}{{\hat{C}}_{1}+{\hat{K}}_{{C}_{1}{S}_{2}}}{\hat{R}}_{S}\right){S}_{1}={\hat{\beta}}_{{C}_{1}{S}_{1}}f\left({\hat{C}}_{1},{\hat{R}}_{S}\right){S}_{1}$.
Linearize around the steady state ${\hat{C}}_{1}=1,{\hat{R}}_{S}=1$ (note $f$(1,1)=0):
Thus,
Similar to the above calculation, we expand $f$ (Equation 16) around steady state 0,
Utilizing Equation 30, Equation 31, and Equation 32,
Taking the derivative of both sides, and using Equation 31 and Equation 32, we have
$\begin{array}{c}\frac{{d}^{2}f}{d{t}^{2}}=\frac{{\hat{K}}_{{C}_{1}{S}_{2}}{\hat{\beta}}_{{C}_{1}{S}_{1}}}{1+{\hat{K}}_{{C}_{1}{S}_{2}}}\frac{d\left({S}_{1}f\right)}{dt}+\frac{{r}_{{S}_{2}{C}_{1}}{\hat{K}}_{{S}_{2}{C}_{1}}}{{\left(1+{\hat{K}}_{{S}_{2}{C}_{1}}\right)}^{2}}{\hat{\beta}}_{{C}_{1}{S}_{1}}\left({\mathrm{\Delta}}_{s}{\hat{C}}_{1}\frac{{\hat{K}}_{{C}_{1}{S}_{2}}}{1+{\hat{K}}_{{C}_{1}{S}_{2}}}+{\mathrm{\Delta}}_{s}{\hat{R}}_{S}\right){S}_{1}\\ =\frac{{\hat{K}}_{{C}_{1}{S}_{2}}{\hat{\beta}}_{{C}_{1}{S}_{1}}}{1+{\hat{K}}_{{C}_{1}{S}_{2}}}\frac{d\left({S}_{1}f\right)}{dt}\frac{{\hat{\beta}}_{{C}_{1}{S}_{1}}{r}_{{S}_{2}{C}_{1}}{\hat{K}}_{{S}_{2}{C}_{1}}}{{\left(1+{\hat{K}}_{{S}_{2}{C}_{1}}\right)}^{2}}f{S}_{1}\end{array}$.
The solution to the above equation is:
where $a={r}_{{S}_{2}{C}_{1}}{\hat{K}}_{{S}_{2}{C}_{1}}{\hat{\beta}}_{{C}_{1}{S}_{1}}{S}_{10}/{\left(1+{\hat{K}}_{{S}_{2}{C}_{1}}\right)}^{2}$ and $b={\hat{\beta}}_{{C}_{1}{S}_{1}}{\hat{K}}_{{C}_{1}{S}_{2}}{S}_{10}/\left(1+{\hat{K}}_{{C}_{1}{S}_{2}}\right)$ are two positive constants. $D}_{1$ and $D}_{2$ are two constants that can be determined from the initial conditions of $\hat{R}}_{S$ and $\hat{C}}_{1$. $M(\kappa ,\mu ,z)$ and $W(\kappa ,\mu ,z)$ are Whittaker functions with argument z. As $z\to \infty $ (http://dlmf.nist.gov/13.14.E20 and http://dlmf.nist.gov/13.14.E21)
$\begin{array}{c}M\left(\kappa ,\mu ,z\right)\sim \mathrm{exp}\left(z/2\right){z}^{\kappa}\\ W\left(\kappa ,\mu ,z\right)\sim \mathrm{exp}\left(z/2\right){z}^{\kappa}\end{array}$.
Thus when ${e}^{{r}_{10}t}b/{r}_{10}\gg 1$,
$\begin{array}{c}f\sim {D}_{1}\mathrm{exp}\left[\frac{b}{2{r}_{10}}{e}^{{r}_{10}t}\frac{{r}_{10}t}{2}+\frac{{e}^{{r}_{10}t}b}{2{r}_{10}}\right]{\left(\frac{{e}^{{r}_{10}t}b}{{r}_{10}}\right)}^{\left(\frac{1}{2}+\frac{a}{b{r}_{10}}\right)}+{D}_{2}\mathrm{exp}\left[\frac{b}{2{r}_{10}}{e}^{{r}_{10}t}\frac{{r}_{10}t}{2}\frac{{e}^{{r}_{10}t}b}{2{r}_{10}}\right]{\left(\frac{{e}^{{r}_{10}t}b}{{r}_{10}}\right)}^{\left(\frac{1}{2}+\frac{a}{b{r}_{10}}\right)}\\ ={\left(\frac{b}{{r}_{10}}\right)}^{\left(\frac{1}{2}+\frac{a}{b{r}_{10}}\right)}{D}_{1}\mathrm{exp}\left(\left({r}_{10}+\frac{a}{b}\right)t\right)+{\left(\frac{b}{{r}_{10}}\right)}^{\left(\frac{1}{2}+\frac{a}{b{r}_{10}}\right)}{D}_{2}\mathrm{exp}\left(\frac{b}{{r}_{10}}{e}^{{r}_{10}t}+\frac{a}{b}t\right)\end{array}$.
The second term approaches zero much faster compared to the first term due to the negative exponent with an exponential term. Thus,
Thus, when ${e}^{{r}_{10}t}b/{r}_{10}\gg 1$, ${\Delta}_{s}f$ =$f$ approaches zero at a rate of $r}_{10}+\frac{{r}_{{S}_{2}{C}_{1}}{\hat{K}}_{{S}_{2}{C}_{1}}\left(1+{\hat{K}}_{{C}_{1}{S}_{2}}\right)}{{\hat{K}}_{{C}_{1}{S}_{2}}{\left(1+{\hat{K}}_{{S}_{2}{C}_{1}}\right)}^{2}$. Therefore, as a conservative estimation, for
the community is sufficiently close to fzeroisocline.
Case II3. ${\hat{R}}_{S}\left(t=0\right)\ll 1/\left(1+{\hat{K}}_{{C}_{1}{S}_{2}}\right)$ or ${R}_{S}\left(0\right)\ll {\beta}_{{C}_{1}{S}_{1}}/{\alpha}_{{C}_{1}{S}_{2}}$
Similar to Case II2, if Equation 28 is satisfied, a typical trajectory is illustrated in Figure 3—figure supplement 2C where the trajectory decreases slightly until ${\hat{C}}_{1}=1$. $\hat{C}}_{1$then increases to much greater than one before the system reaches the fzeroisocline. ${t}_{f}$ can then be estimated from $t}_{f1$, the time it takes for $\hat{C}}_{1$ to reach 1 and $t}_{f2$ , the time takes for $\hat{R}}_{S$ to increase to $1/\left(1+{\hat{K}}_{{C}_{1}{S}_{2}}\right)$. Using Equation 15, since $\hat{R}}_{S$ decreases very little, and ${\hat{R}}_{S}\left(t=0\right)\ll 1/\left(1+{\hat{K}}_{{C}_{1}{S}_{2}}\right)$,
Therefore, ${\hat{C}}_{1}\approx \frac{{\hat{\beta}}_{{C}_{1}{S}_{1}}{S}_{1}(0)}{{r}_{10}}({e}^{{r}_{10}t}1)$.
During $t}_{f1$, $\hat{C}}_{1$ increases from 0 to 1. Thus, $1\approx \frac{{\hat{\beta}}_{{C}_{1}{S}_{1}}{S}_{1}(0)}{{r}_{10}}({e}^{{r}_{10}{t}_{f1}}1)$. $t}_{f1}\approx \mathrm{ln}({r}_{10}/({\hat{\beta}}_{{C}_{1}{S}_{1}}{S}_{10})+1)/{r}_{10$.
If
$t}_{f1}\approx ({\hat{\beta}}_{{C}_{1}{S}_{1}}{S}_{10}{)}^{1}\ll {r}_{10}^{1$.
Using Equation 14 and since $\hat{C}}_{1$ is very large,
$\frac{d{\hat{R}}_{S}}{dt}=\left({r}_{20}+{r}_{{S}_{2}{C}_{1}}\frac{{\hat{C}}_{1}}{{\hat{C}}_{1}+{\hat{K}}_{{S}_{2}{C}_{1}}}{r}_{10}\right){\hat{R}}_{S}\approx \left({r}_{20}+{r}_{{S}_{2}{C}_{1}}{r}_{10}\right){\hat{R}}_{S}$.
This yields
${t}_{f2}\approx \mathrm{ln}\left(\frac{1}{{\hat{R}}_{S}\left(0\right)\left(1+{\hat{K}}_{{C}_{1}{S}_{2}}\right)}\right)/\left({r}_{20}+{r}_{{S}_{2}{C}_{1}}{r}_{10}\right)$,
and a conservative estimation of $t}_{f$ is
In the unscaled form, this becomes:
Case III: ${r}_{10}<{r}_{20}$
In this case, supplier S_{1} always grows slower than S_{2}. As $t\to \infty $, ${R}_{S}={S}_{2}/{S}_{1}\to \infty $ and ${C}_{1}\to 0$. The phase portrait is separated into two parts by the fzeroisocline (Figure 3—figure supplement 2E), where, as in Equation 12,
$f({C}_{1},{R}_{S})=1\frac{{\alpha}_{{C}_{1}{S}_{2}}}{{\beta}_{{C}_{1}{S}_{1}}}\frac{{C}_{1}}{{C}_{1}+{K}_{{C}_{1}{S}_{2}}}{R}_{S}=0$ or ${R}_{S}=\frac{{\beta}_{{C}_{1}{S}_{1}}}{{\alpha}_{{C}_{1}{S}_{2}}}\left(1+\frac{{K}_{{C}_{1}{S}_{2}}}{{C}_{1}}\right)$.
Note that the asymptotic value of ${R}_{S}$ (black dotted line, Figure 3—figure supplement 2E–H) is
From Equation 9, $d{R}_{S}/dt>0$. From Equation 8, below the $f$zeroisocline, $dC{}_{1}/dt>0$ and above the $f$zeroisocline, $dC{}_{1}/dt<0$. Thus if the system starts from $(0,{R}_{S}(0\left)\right)$, the phase portrait dictates that it moves with a positive slope until a time of a scale ${t}_{f}$ when it hits the $f$zeroisocline, after which it moves upward to the left closely along the $f$zeroisocline (Figure 3—figure supplement 2E). After ${t}_{f}$, the alternative pairwise model can be derived. Although ${t}_{f}$ is difficult to estimate in general, it is possible for the following cases.
Case III1. ${R}_{S}\left(0\right)\gg {\beta}_{{C}_{1}{S}_{1}}/{\alpha}_{{C}_{1}{S}_{2}}$
Similar to Case II2, if ${S}_{10}$ is small, there is a lagging phase during which the trajectory rises steeply before leveling off (Figure 3—figure supplement 2H). Although the alternative pairwise model can be derived once the trajectory hits the $f$zeroisocline, ${t}_{f}$ takes a complicated form. Here we consider two cases where ${S}_{10}$ is large enough so that we can approximate the trajectory as a straight line going through $(0,{R}_{s}(t=0))$ (Figure 3—figure supplement 2F). Graphically, ${S}_{10}$ is large enough so that the green segment in Figure 3—figure supplement 2F, whose length is $\Delta {R}_{S}$, is much smaller than ${R}_{S}\left(0\right)$. In other words,
From Equations 8 and 9
$\frac{d{R}_{S}}{d\left({C}_{1}/{K}_{{C}_{1}{S}_{2}}\right)}}_{{C}_{1}\left(0\right)=0}={\frac{d{R}_{S}/dt}{d{C}_{1}/dt}}_{t=0}{K}_{{C}_{1}{S}_{2}}=\frac{\left({r}_{20}{r}_{10}\right){R}_{S}\left(0\right){K}_{{C}_{1}{S}_{2}}}{{\beta}_{{C}_{1}{S}_{1}}{S}_{1}\left(0\right)$.
$\Delta ({C}_{1}/{K}_{{C}_{1}{S}_{2}})$, the red segment in Figure 3—figure supplement 2F, is the horizontal distance between $(0,{R}_{S}(0\left)\right)$ and the $f$zeroisocline and
Thus, if
or
then from Equation 9 and $r}_{20}>{r}_{10$, the upper bound of ${t}_{f}$ can be calculated as
Case III2. $R}_{S}\left(0\right)\ll {\beta}_{{C}_{1}{S}_{1}}/{\alpha}_{{C}_{1}{S}_{2}$
A typical example is displayed in Figure 3—figure supplement 2G. The trajectory moves with a small positive slope so that the intersection of the community dynamics trajectory with the fzeroisocline is near the black dotted line ${\beta}_{{C}_{1}{S}_{1}}/{\alpha}_{{C}_{1}{S}_{2}}$ (Equation 39) where ${C}_{1}/{K}_{{C}_{1}{S}_{2}}$ is large. The upper bound of ${t}_{f}$ can thus be estimated from Equation 9:
which yields a conservative estimate of
Conditions for the alternative pairwise model to approximate the mechanistic model
Cases II and III showed that population dynamics of the mechanistic model could be described by the alternative pairwise model. However, since the initial condition for $C}_{1$ cannot be specified in pairwise model, problems could occur. To illustrate, we examine the phase portrait of the pairwise equation
where $\omega =1\frac{{K}_{{S}_{2}{C}_{1}}}{{K}_{{C}_{1}{S}_{2}}}$, $\psi =\frac{{\alpha}_{{C}_{1}{S}_{2}}{K}_{{S}_{2}{C}_{1}}}{{\beta}_{{C}_{1}{S}_{1}}{K}_{{C}_{1}{S}_{2}}}$. From Equations 8 and 13,
Below, we plot Equation 43 under different parameters (Figure 3—figure supplement 3) to reveal conditions for convergence between mechanistic and pairwise models.
Case II (${r}_{{S}_{2}{C}_{1}}>{r}_{10}{r}_{20}>0$): steady state $R}_{S}^{\ast$ exists for mechanistic model.
If $\omega =1{K}_{{S}_{2}{C}_{1}}/{K}_{{C}_{1}{S}_{2}}\ge 0$ (Figure 3—figure supplement 3A): When $R}_{S}<{R}_{S}^{\ast$, $d{R}_{S}/dt$ is positive. When $R}_{S}>{R}_{S}^{\ast$, $d{R}_{S}/dt$ is negative. Thus, wherever the initial $R}_{S$, it will always converge toward the only steady state $R}_{S}^{\ast$ of the mechanistic model.
If $\omega <0$ (Figure 3—figure supplement 3B): $\omega +\psi {R}_{S}=0$ or ${R}_{S}=\omega /\psi $ creates singularity. Pairwise model $R}_{S$ will only converge toward the mechanistic model steady state if
Case III (${r}_{10}<{r}_{20}$): $R}_{S$ increases exponentially in mechanistic model (Equation 9). Thus, $C}_{1$ will decline toward zero as C_{1} is consumed by S_{2} whose relative abundance over S_{1} exponentially increases. Hence, according to Equation 9, $R}_{S$ eventually increases exponentially at a rate of ${r}_{20}{r}_{10}$.
If $\omega \ge 0$ (Figure 3—figure supplement 3C): Equation 43 $\frac{d{R}_{S}}{dt}=\left({r}_{20}+\frac{{r}_{{S}_{2}{C}_{1}}}{\omega +\psi {R}_{S}}{r}_{10}\right){R}_{S}$>0. Thus, Equation 43, which is based on alternative pairwise model, also predicts that $R}_{S$ will eventually increase exponentially at a rate of ${r}_{20}{r}_{10}$, similar to the mechanistic model.
If $\omega <0$ (Figure 3—figure supplement 3D): ${R}_{S}\left(0\right)>\omega /\psi $ (Equation 44) is required for unbounded increase in $R}_{S$ (similar to the mechanistic model). Otherwise, $R}_{S$ converges to an erroneous value instead.
Conditions under which a saturable LV pairwise model can represent one species influencing another via two reusable mediators
Here, we examine a simple case where S_{1} releases reusable C_{1} and C_{2}, and C_{1} and C_{2} additively affect the growth of S_{2} (see example in Figure 5). Similar to Figure 3A, the mechanistic model is:
Now the question is whether the saturable LV pairwise model
can be a good approximation.
For simplicity, let’s define ${K}_{C1}={K}_{{S}_{2}{C}_{1}}{r}_{10}/{\beta}_{{C}_{1}{S}_{1}}$ and ${K}_{C2}={K}_{{S}_{2}{C}_{2}}{r}_{10}/{\beta}_{{C}_{2}{S}_{1}}$. Small $K}_{Ci$ means large potency (e.g. small $K}_{C2$ can be caused by small ${K}_{{S}_{2}{C}_{2}}$ which means low $C}_{2$ required to achieve half maximal effect on S_{2}, and/or large synthesis rate ${\beta}_{{C}_{2}{S}_{1}}$). Since $S}_{1$ from pairwise and mechanistic models are identical, we have
$\overline{D}$ can be close to zero when (i) ${K}_{C1}\approx {K}_{C2}$ or (ii) $\frac{{r}_{{S}_{2}{C}_{1}}{S}_{1}}{{S}_{1}+{K}_{C1}}$ and $\frac{{r}_{{S}_{2}{C}_{2}}{S}_{1}}{{S}_{1}+{K}_{C2}}$(effects of C_{1} and C_{2} on S_{2}) differ dramatically in magnitude. For (ii), without loss of generality, suppose that the effect of C_{2} on S_{2} can be neglected. This can be achieved if (iia) ${r}_{{S}_{2}{C}_{2}}$ is much smaller than ${r}_{{S}_{2}{C}_{1}}$, or (iib) ${K}_{{C}_{2}}$ is large compared to S_{1}.
Competitive commensal interaction
For the community in Figure 6A, our mechanistic model is:
Here, $S}_{1$ and $S}_{2$ are the densities of the two species; $r}_{i0$ is the basal net growth rate of S_{i} (negative, representing death in the absence of the essential shared resource C_{1}); C_{1} is supplied at a constant rate $\beta}_{0$; ${\beta}_{{C}_{2}}{}_{{S}_{1}}$is the production rate of C_{2} by S_{1}; ${\alpha}_{{C}_{i}}{}_{{S}_{j}}$ is the amount of resource C_{i} consumed to produce a new S_{j} cell.
The growth of S_{2} is controlled by $C}_{1$ and $C}_{2$. When $C}_{1$ is limiting ($C}_{1}/{K}_{{S}_{2}{C}_{1}}\ll {C}_{2}/{K}_{{S}_{2}{C}_{2}$), the fitness influence of the two chemicals on S_{2} becomes:
which is the standard Monod equation. A similar argument holds for limiting $C}_{2$. We have intentionally chosen very large ${K}_{{S}_{2}}{}_{{C}_{2}}$ to ensure that the fitness effect of C_{2} on S_{2} is linear with respect to $C}_{2$. This way, we minimize the number of pairwise model parameters that need to be estimated.
For our LV pairwise model, to capture intraspecies competition, we use
where nonnegative $b}_{i0$ represents the maximal birth rate of S_{i} at nearly zero population density (no competition), and nonnegative $d}_{i$ represents the constant death rate of S_{i}. Positive ${\kappa}_{i}$ is the ‘carrying capacity’ imposed by the limiting resource, and is the $S}_{i$ at which birth rate becomes zero. This equation can be simplified to:
When ${\Lambda}_{i}$>0 (i.e. when ${b}_{i0}>{d}_{i}$), this resembles standard LV model traditionally used for competitive interactions (compare to Equation 2; Gause, 1934; Thébault and Fontaine, 2010; Mougi and Kondoh, 2012).
Thus, for the competitive commensal community, we have:
Here, birth rate of each species is reduced by competition from the two species, and ${\Lambda}_{ij}$ is the carrying capacity such that a single S_{i} individual will have a zero birth rate when encountering a total of ${\Lambda}_{ij}$ individuals of S_{j}. For S_{2}, We used $\left({b}_{20}+{r}_{21}{S}_{1}\right)\left(1\frac{{S}_{1}}{{\Lambda}_{21}}\frac{{S}_{2}}{{\Lambda}_{22}}\right){S}_{2}$ instead of ${b}_{20}\left(1\frac{{S}_{1}}{{\Lambda}_{21}}\frac{{S}_{2}}{{\Lambda}_{22}}\right){S}_{2}+{r}_{21}{S}_{1}{S}_{2}$ so that when the shared resource is exhausted (i.e. $1\frac{{S}_{1}}{{\Lambda}_{21}}\frac{{S}_{2}}{{\Lambda}_{22}}$=0), S_{2} does not keep growing due to the presence of S_{1}.
References
 1

2
Analysis of doublesubstrate limited growthBiotechnology and Bioengineering 20:183–202.https://doi.org/10.1002/bit.260200203

3
Highorder species interactions shape ecosystem diversityNature Communications 7:12285.https://doi.org/10.1038/ncomms12285
 4
 5

6
PopulationChangesAccessed June 2016.
 7

8
Testing for higher order interactionsThe American Naturalist 118:920–929.https://doi.org/10.1086/283885

9
Global stability and multiple domains of attraction in ecological systemsThe American Naturalist 113:705–714.https://doi.org/10.1086/283427
 10

11
MacArthur's consumerresource modelTheoretical Population Biology 37:26–38.https://doi.org/10.1016/00405809(90)90025Q
 12
 13
 14

15
Experimental evidence rejects pairwise modelling approach to coexistence in plant communitiesProceedings of the Royal Society B: Biological Sciences 272:1279–1285.https://doi.org/10.1098/rspb.2005.3066

16
Chemical interactions between organisms in microbial communitiesContributions to Microbiology 16:1–17.https://doi.org/10.1159/000219369

17
The importance of being discrete (and spatial)Theoretical Population Biology 46:363–394.https://doi.org/10.1006/tpbi.1994.1032

18
Microbial interactions: from networks to modelsNature Reviews Microbiology 10:538–550.https://doi.org/10.1038/nrmicro2832

19
The cellulosome: the exocellular organelle of ClostridiumAnnual Review of Microbiology 47:791–819.https://doi.org/10.1146/annurev.mi.47.100193.004043
 20

21
Community structure follows simple assembly rules in microbial microcosmsNature Ecology & Evolution 1:0109.https://doi.org/10.1038/s415590170109
 22

23
Serine betalactamases and penicillinbinding proteinsAnnual Review of Microbiology 45:37–67.https://doi.org/10.1146/annurev.mi.45.100191.000345

24
Stability and Oscillations in Delay Differential Equations of Population DynamicsSpringer Science & Business Media.

25
Modeling community population dynamics with the opensource language RMethods in Molecular Biology 1151:209–231.https://doi.org/10.1007/9781493905546_15
 26
 27

28
A growthrate composition formula for the growth of E.coli on coutilized carbon substratesMolecular Systems Biology 11:801.https://doi.org/10.15252/msb.20145537
 29

30
Talk of the town: interspecies communication in oral biofilmsMolecular Oral Microbiology 25:4–14.https://doi.org/10.1111/j.20411014.2009.00563.x

31
Saccharification of complex cellulosic substrates by the cellulase system from Clostridium thermocellumApplied and Environmental Microbiology 43:1125–1132.

32
Stable coexistence of five bacterial strains as a cellulosedegrading communityApplied and Environmental Microbiology 71:7099–7106.https://doi.org/10.1128/AEM.71.11.70997106.2005

33
Network relationships of bacteria in a stable mixed cultureMicrobial Ecology 56:403–411.https://doi.org/10.1007/s0024800793574
 34
 35
 36

37
Oral microbial communities: biofilms, interactions, and genetic systemsAnnual Review of Microbiology 54:413–437.https://doi.org/10.1146/annurev.micro.54.1.413

38
Interspecies interactions within oral microbial communitiesMicrobiology and Molecular Biology Reviews 71:653–670.https://doi.org/10.1128/MMBR.0002407
 39

40
Competitive interactions in ecosystemsThe American Naturalist 110:903–910.https://doi.org/10.1086/283116

41
American Scientist421–431, The Strategy of Model Building in Population Biology, American Scientist, 54.

42
Which games are growing bacterial populations playing?Journal of the Royal Society, Interface 12:20150121.

43
Species packing and competitive equilibrium for many speciesTheoretical Population Biology 1:1–11.https://doi.org/10.1016/00405809(70)900390

44
Model for microbial growth with more than one limiting nutrientJournal of Biotechnology 7:161–166.https://doi.org/10.1016/01681656(88)900624
 45

46
Physiological approaches to the control of oral biofilmsAdvances in Dental Research 11:176–185.https://doi.org/10.1177/08959374970110010901
 47
 48
 49

50
Microbial interactions within a cheese microbial communityApplied and Environmental Microbiology 74:172–181.https://doi.org/10.1128/AEM.0133807
 51
 52

53
The gut microbiota shapes intestinal immune responses during health and diseaseNature Reviews Immunology 9:313–323.https://doi.org/10.1038/nri2515
 54

55
The cellulosome and cellulose degradation by anaerobic bacteriaApplied Microbiology and Biotechnology 56:634–649.https://doi.org/10.1007/s002530100710

56
A review: the anaerobic treatment of sewage in UASB and EGSB reactorsBioresource Technology 65:175–190.https://doi.org/10.1016/S09608524(98)000467

57
Metabolic interactions between anaerobic bacteria in methanogenic environmentsAntonie Van Leeuwenhoek 66:271–294.https://doi.org/10.1007/BF00871644

58
Interacting guilds: moving beyond the pairwise perspective on mutualismsThe American Naturalist 162:S10–S23.https://doi.org/10.1086/378646
 59
 60
 61

62
The importance of the mechanisms of interspecific competitionThe American Naturalist 129:769–774.https://doi.org/10.1086/284672
 63

64
Application of stateofart sequencing technologies to indigenous food fermentationsCurrent Opinion in Biotechnology 24:178–186.https://doi.org/10.1016/j.copbio.2012.08.004
 65
 66

67
Variations and fluctuations of the number of individuals in animal species living togetherJournal Du Conseil Conseil Permanent International Pour Exploration De La Mer 1931:3–51.

68
LotkaVolterra population modelsAnnual Review of Ecology and Systematics 9:189–218.https://doi.org/10.1146/annurev.es.09.110178.001201
 69
 70
 71

72
The nature and consequences of indirect effects in ecological communitiesAnnual Review of Ecology and Systematics 25:443–466.https://doi.org/10.1146/annurev.es.25.110194.002303

73
Indirect effects in complex ecosystems: recent progress and future challengesJournal of Sea Research 48:157–172.https://doi.org/10.1016/S13851101(02)001491

74
HigherOrder interactions and indirect effects: a resolution using laboratory Drosophila communitiesThe American Naturalist 138:1092–1104.https://doi.org/10.1086/285271

75
Translating metabolic exchange with imaging mass spectrometryNature Chemical Biology 5:885–887.https://doi.org/10.1038/nchembio.252
Decision letter

Bruce LevinReviewing Editor; Emory University, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
[Editors’ note: a previous version of this study was rejected after peer review, but the authors submitted for reconsideration. The first decision letter after peer review is shown below.]
Thank you for submitting your work entitled "The validity of pairwise models in predicting community dynamics" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor, Bruce R. Levin, and a Senior Editor. The reviewers have opted to remain anonymous.
Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.
This review and editorial decision, has two elements. Part 1 has been put together by the Reviewing editor but is the product of a discussion between the editor and the reviewers. Part 2 are the separate comments and suggestions of the reviewers. There is a fair amount of overlap. BRL agrees with all of the critical comments and suggestions of the reviewers. Following our discussion, BRL sent the first part of this recommendation to the reviewers. Save for a typo, they approved of this collective review and recommendation.
Part 1.
Recommendation:
The subject of this report is right on. In addition to the considerable industrial and ecological and medical interest in bacterial communities, thanks to the "microbiome mania" the need to understand the processes that determine and maintain the structure of these communities has become increasingly important. Central to acquiring this understanding will doubtless be mathematical and/or computer simulations models of the population and evolutionary dynamics of bacteria and the interactions between these populations and the physical, chemical and biotic factors that determine their densities and relative frequencies in communities. How to construct and analyze the properties of these models is not at all clear at this time; the longstanding fissure between population and community ecology has yet to be breached. This report considers the problems of modeling bacterial populations and communities.
As trendy and important as this subject is, we don't consider this report to be of sufficient general interest to be published in ELife. Although it may be of interest to some theoreticians, we believe it would be of little interest and utility to population and evolutionary biologists and biometricians doing experimental or other empirical research with bacteria and bacterial communities. It is too abstract and doesn't address specific questions. Some readers may also see it as a complex, equation illustrated rant telling them what they already know about the limitations of pairwise models. We believe it would be more suitable for a fine but more specialized journal like PLoS Computational Biology, where the readers will be more attune to mathematical modeling of the sort considered in this report.
Collective comments and suggestions:
1) – We agree with the George Box adage, "All models are wrong, some are useful". However, we don't see how their pairwise models would be useful for understanding the processes that determine "distribution and abundance" of bacterial strains and species in communities. To be sure, even without understanding the mechanisms responsible, it would be useful to be able to predict the distribution and abundance of species from empirical estimates of the parameter of these "pairwise" models. But as we see this report as a cautionary tale saying that it is unlikely that pairwise models will be able to achieve these ends much beyond three species and if that. Moreover, even if it were possible to provide predictive pairwise models for specific simple communities, perhaps the macrobiotic of yogurt, there would be little or no generality, like to cheese or bread starters, much less natural communities.
2) Arguably (not to BRL) mechanistic models are most useful when they are wrong, when there are substantial qualitative rather than just small quantitative differences between the predictions of the model and observations in the empirical studies. That way one knows that there are fundamental errors in the biological, chemical and environmental assumptions upon which the model was based. In the best of cases, a wrong mechanistic model will point to the biological and other assumptions that have to be modified to obtain a better fit and thereby increase our understanding of the mechanisms responsible for the phenomenon under study. It is not clear how the pairwise models of the sort they are considering could achieve this end when they fit or don't fit.
3) What one means by a "mechanistic model" and mechanisms in general is in the eyes, mind and perspective of the beholder. The authors can do a better job of providing the readers with a clear distinction between what they mean by pairwise rather than mechanistic models. They present the logistic model and its extensions, like the Gause competition models, as pairwise rather than mechanistic models. We agree, but some may not. These models make the reasonable assumption that the rate of growth of populations decline as they approach the point of saturation of the environment, the parameter K. But the authors and most others interpret that parameter to reflect resource limitation, which seems mechanistic albeit not explicit in the model. Moreover, the ecological nature of the equilibrium in Logistic models and its extensions is not defined. Is that equilibrium, K, the density at stationary phase in batch culture? Or is the population at equilibrium in an environment in which resources are continually made available and bacteria and wastes are continually being removed, like a chemostat or turbidostat?
4) There is a real need to develop models that deal with the inconvenient reality of the physical structure of real habitats of bacteria. The ODE mechanistic and pairwise models considered in this report do not address this reality, which doubtless contributes to the distribution and abundance of species and strains of bacteria in natural communities.
Part 2.
Reviewer #1:
In this work, the authors study the validity of utilizing phenomenological models of pairwise interactions to describe the dynamics of ecological communities. To do so, they employ theory and simulations to compare the dynamics of detailed mechanistic models to the ones of appropriately parameterized pairwise models.
Improving our understanding of the capabilities and limitations of pairwise models is an important and timely goal, especially given their recent popularity in modeling microbial communities. I also appreciate the authors' approach of trying to delineate the types of mechanistic situations in which the pairwise approximation is valid.
However, the aspects of the problem that the authors focused on are not always the ones I feel are the most interesting and useful. Specifically, the authors focused on identifying situations in which a specific pairwise model can provide an exact description of the community dynamics. It is not surprising that phenomenological pairwise models often fail to capture exactly the dynamics of more complex mechanistic models. Nonetheless, there may be many situations where they are still useful, either by providing approximate description of the dynamics, or by capturing important qualitative features, such as the existence of oscillations, the set of coexisting species, presence of alternative stable states, etc. To me, understanding the conditions under which pairwise models are not even approximately correct, or make qualitatively wrong prediction is one of the important outstanding challenges in community modeling. I would encourage the authors to use their approach to provide insight into these questions, but, at the very least, they should clarify the difference between different failures of pairwise models. Additionally, they should be more explicit about cases in which no pairwise model would capture the community dynamics (e.g. in cases of interaction modification, or higher order interaction), versus ones in which the pairwise model that was considered wasn't adequate, but there may be a different one that is.
The authors also do not consider cases where interactions are mediated by externally supplied, abiotic mediators, rather than ones produced by the species themselves. The author's model is an extension of MacArthur's competition model, whose link to pairwise models have been extensively studied (e.g. Chesson, P. "MacArthur's Consumer Resource Model." Theoretical Population Biology 37, no. 1 (1990): 2638.).
What is not captured is competition of abiotic resources, such as the models common in David Tilman's work. The authors should make these distinctions clear, and put their work in the context of previous work. I also believe that extending the work to include competition for abiotic resources would add significant value to the work.
Reviewer #2:
The work of Momeni et al. "The validity of pairwise models in predicting community dynamics" explores conditions under which mechanistic models of species interacting via chemical mediators can be reduced to pairwise models. Pairwise models do not require a full mechanistic understanding of the nature of interactions within the community and thus use fewer parameters, which is why they are often used. The authors report that in many cases, pairwise models fail to predict, quantitatively or qualitatively, the dynamics of many species communities, which is why they should be used with caution. The authors do a good job in exploring the various scenarios under which pairwise models are not sufficient to capture these dynamics.
A point that deserves more attention however is how these results are immediately relevant for past and future studies of microbial communities? The authors provide references to several studies that employed pairwise models to simulate community dynamics. However, many of those do not model chemically mediated interactions, but rather direct interactions, such as the predator prey examples mentioned in the introduction and are thus not directly relevant to what is studied here. Can the authors provide more specific examples of works that employ pairwise models to model chemically mediated interactions?
The authors clearly describe several regimes under which pairwise models fail. However, to identify whether the community being modeled falls into these regimes, one often needs to have a quite detailed mechanistic understanding of the interactions within the community in the first place. What would then be the advantages of using a pairwise model, given that this information is available? Also, it is not clear if this is an exhaustive list of regimes that they look at or are there other regimes out there to be explored? Thus, we are left with a conundrum that the work does not address: why the results are useful for systems where we do not know all the details of the interactions?
In general, for dynamical systems it is not at all surprising that 'the devil is in the details' of interactions for many a system, as weak couplings can often have important roles. Especially for systems that are under evolutionary pressure I can see how small couplings can be essential. Thus, models with underdetermined number of parameters will often fail, especially as long as the main state variables of a system are unknown, or when these tend to shift in importance from system to system.
Thus, to me the strong 'punch' in the message of the paper is missing, especially for a wider audience outside of the ecology crowd.
As a reviewer I am part of this outside audience I am not an ecologist nor am I a theorist.
https://doi.org/10.7554/eLife.25051.050Author response
[Editors’ note: the author responses to the first round of peer review follow.]
This review and editorial decision, has two elements. Part 1 has been put together by the Reviewing editor but is the product of a discussion between the editor and the reviewers. Part 2 are the separate comments and suggestions of the reviewers. There is a fair amount of overlap. BRL agrees with all of the critical comments and suggestions of the reviewers. Following our discussion, BRL sent the first part of this recommendation to the reviewers. Save for a typo, they approved of this collective review and recommendation.
Part 1.
Recommendation:
The subject of this report is right on. In addition to the considerable industrial and ecological and medical interest in bacterial communities, thanks to the "microbiome mania" the need to understand the processes that determine and maintain the structure of these communities has become increasingly important. Central to acquiring this understanding will doubtless be mathematical and/or computer simulations models of the population and evolutionary dynamics of bacteria and the interactions between these populations and the physical, chemical and biotic factors that determine their densities and relative frequencies in communities. How to construct and analyze the properties of these models is not at all clear at this time; the longstanding fissure between population and community ecology has yet to be breached. This report considers the problems of modeling bacterial populations and communities.
As trendy and important as this subject is, we don't consider this report to be of sufficient general interest to be published in eLife. Although it may be of interest to some theoreticians, we believe it would be of little interest and utility to population and evolutionary biologists and biometricians doing experimental or other empirical research with bacteria and bacterial communities.
The famous quote of George Box (“All models are wrong, some are useful”) reflects the need to omit less important details and focus on important aspects of a system when constructing a model. This quote creates a dilemma for experimentalists: how would they decide which models are sufficiently useful? Their fears are legitimate: although all models are wrong, there are different reasons for being “wrong”. Certain models are “wrong” because details (e.g. the precise values of parameters) are off. Other models are wrong because the fundamental assumptions are wrong (e.g. blending inheritance). The first class can still yield insightful information that transcends details (e.g. the famous LV pairwise model capturing the oscillatory dynamics of predatorprey communities). We presume that this is the reason why theoretical papers are still being written and some of them are still published in highprofile journals. Given the complexity of microbial communities, the temptation of using simple pairwise modeling by “biologists and biometricians doing experimental or other empirical research with bacteria and bacterial communities” is in fact strong (see Introduction, fourth paragraph). Thus, our work is directly relevant to people who work on microbial communities.
It is too abstract and doesn't address specific questions. Some readers may also see it as a complex, equation illustrated rant telling them what they already know about the limitations of pairwise models. We believe it would be more suitable for a fine but more specialized journal like PLoS Computational Biology, where the readers will be more attune to mathematical modeling of the sort considered in this report.
The specific question we have addressed is “Can we use a single equation form to represent diverse pairwise microbial interactions, as done in multispecies pairwise modeling?” The answer is “No”, as elaborated in Figures 3–6. This fundamental flaw we have uncovered is orthogonal to the known limitation of pairwise modeling (i.e. modification of pairwise interactions by a third species). To our knowledge, our findings have not been discussed in the literature. This is a major contribution, given the popularity and intellectual impact of pairwise models.
From reading reviewers’ comments, we realize that our writing may have been too abstract. Thus, we have reduced the use of equations in the main text (only keeping the essential ones) and have relegated most of our mathematical derivations to Methods so that general readers don’t have to read them, but for those who might be interested, the information is readily available. To demonstrate success or failure of pairwise modeling, we have replaced abstract figures of parameter comparisons with easiertoread figures of population composition dynamics. Finally, since it is already known that interaction modification by a third species undermines multispecies pairwise modeling, we have cut down descriptions on this and moved it (Figure 7) from Results to Discussion.
Collective comments and suggestions:
1) We agree with the George Box adage, "All models are wrong, some are useful". However, we don't see how their pairwise models would be useful for understanding the processes that determine "distribution and abundance" of bacterial strains and species in communities.
We are arguing the opposite: pairwise models are not as useful as one might hope when applied to microbial communities and should be used and interpreted with great caution.
To be sure, even without understanding the mechanisms responsible, it would be useful to be able to predict the distribution and abundance of species from empirical estimates of the parameter of these "pairwise" models. But as we see this report as a cautionary tale saying that it is unlikely that pairwise models will be able to achieve these ends much beyond three species and if that.
As stated above, we argue that a single equation form, which is the standard of practice in pairwise modeling, often fails even for twospecies microbial communities.
Moreover, even if it were possible to provide predictive pairwise models for specific simple communities, perhaps the macrobiotic of yogurt, there would be little or no generality, like to cheese or bread starters, much less natural communities.
The possible lack of generality among different communities does not negate the importance of being able to understand communities of interest.
2) Arguably (not to BRL) mechanistic models are most useful when they are wrong, when there are substantial qualitative rather than just small quantitative differences between the predictions of the model and observations in the empirical studies. That way one knows that there are fundamental errors in the biological, chemical and environmental assumptions upon which the model was based. In the best of cases, a wrong mechanistic model will point to the biological and other assumptions that have to be modified to obtain a better fit and thereby increase our understanding of the mechanisms responsible for the phenomenon under study. It is not clear how the pairwise models of the sort they are considering could achieve this end when they fit or don't fit.
We agree with you that failure of pairwise models may not be as informative as failure of mechanistic models. However, to be fair, pairwise models can be useful in some aspects. We have added this paragraph to the Discussion:
“When might pairwise modeling be useful? First, pairwise modeling has been instrumental in understanding ecological phenomena such as preypredator oscillatory dynamics and coexistence of competing predator species (Volterra 1926; MacArthur 1970; Case and Casten 1979; Chesson 1990). […] Not surprisingly, predicting qualitative consequences of species removal or addition using pairwise modeling has encountered difficulties, especially in communities of more than three species (Mounier et al. 2008; Friedman, Higgins, and Gore 2016).”
3) What one means by a "mechanistic model" and mechanisms in general is in the eyes, mind and perspective of the beholder. The authors can do a better job of providing the readers with a clear distinction between what they mean by pairwise rather than mechanistic models.
This is a very good point. We have now provided a clear definition of what we mean by “mechanistic models” in the Introduction:
“We define “mechanistic models” as models that explicitly consider interaction mediators as state variables. For example, if species S1 releases a compound C1 which stimulates species S2 growth upon consumption by S2, then a mechanistic model tracks concentrations of S1, C1, and S2 (Figure 1A and B, left panels). Note that mechanistic models used here still omit molecular details such as how chemical mediators are received and processed to affect recipient.”
They present the logistic model and its extensions, like the Gause competition models, as pairwise rather than mechanistic models. We agree, but some may not. These models make the reasonable assumption that the rate of growth of populations decline as they approach the point of saturation of the environment, the parameter K. But the authors and most others interpret that parameter to reflect resource limitation, which seems mechanistic albeit not explicit in the model. Moreover, the ecological nature of the equilibrium in Logistic models and its extensions is not defined. Is that equilibrium, K, the density at stationary phase in batch culture? Or is the population at equilibrium in an environment in which resources are continually made available and bacteria and wastes are continually being removed, like a chemostat or turbidostat?
In our definition, the logistic model is not considered mechanistic because it does not model interaction mediators (the limiting metabolite). We agree that the logistic model can nicely capture resource competition, as shown by Gause. To clarify (Figure 1 table), in our mechanistic model, parameter K_{SC}is the concentration of mediator C where half maximal fitness effect on recipient species S is achieved. In our pairwise model, parameter K_{ij}is the density of the interactioninitiating species S_{j}where half maximal fitness effect on recipient species S_{i}is achieved. Communities in our original manuscript were grown in a turbidostat so that no metabolites other than interaction mediators were limiting. We intentionally did so to avoid introducing additional mediators and additional parameters. However, as per reviewer 1’s request, we have added an example where two species in a batch environment compete for an abiotic metabolite while simultaneously engaging in commensalism (Figure 6). We show that although a logistic LV pairwise model can work in some cases, it generates wrong predictions in other cases, depending on parameter choices of the mechanistic model.
4) There is a real need to develop models that deal with the inconvenient reality of the physical structure of real habitats of bacteria. The ODE mechanistic and pairwise models considered in this report do not address this reality, which doubtless contributes to the distribution and abundance of species and strains of bacteria in natural communities.
Spatial structure is undoubtedly important, as demonstrated in our earlier work (eLife 00230, eLife 00960). In this work, we consider a wellmixed environment. We have added a paragraph to justify our choice in Results:
“Throughout this work, we consider communities grown in a wellmixed environment where all individuals interact with each other at an equal chance. […] This in turn allows us to sometimes analytically demonstrate failures of pairwise modeling.”
If a model fails to capture microbial interactions in a wellmixed environment, a spatiallystructured environment is unlikely to rescue this failure.
Part 2.
Reviewer #1:
In this work, the authors study the validity of utilizing phenomenological models of pairwise interactions to describe the dynamics of ecological communities. To do so, they employ theory and simulations to compare the dynamics of detailed mechanistic models to the ones of appropriately parameterized pairwise models.
Improving our understanding of the capabilities and limitations of pairwise models is an important and timely goal, especially given their recent popularity in modeling microbial communities. I also appreciate the authors' approach of trying to delineate the types of mechanistic situations in which the pairwise approximation is valid.
However, the aspects of the problem that the authors focused on are not always the ones I feel are the most interesting and useful. Specifically, the authors focused on identifying situations in which a specific pairwise model can provide an exact description of the community dynamics. It is not surprising that phenomenological pairwise models often fail to capture exactly the dynamics of more complex mechanistic models.
Our writing did at times create the impression that we were only interested in exact matches, but this is not our intention. Figures 4, 5, 6 and 7, and Figure 3—figure supplement 4 now demonstrate qualitative failures of using a single pairwise equation form to describe diverse pairwise microbial interactions and diverse details of microbial communities. We have also revised our writing in Results to clarify our intention.
Our goal is to test whether a single form of pairwise model can qualitatively predict dynamics of species pairs engaging in various types of interactions commonly found in microbial communities (e.g. Figure 2).
We summarize our results in the Discussion:
“We first consider cases where abiotic resources are in excess. When one species affects another species via a single chemical mediator, either saturable LV or alternative pairwise model is appropriate, depending on interaction mechanism (consumable versus reusable mediator), relative fitness of the two species, and initial conditions (Figure 3; Figure 3—figure supplement 2 to Figure 3—figure supplement 5). […] Thus, although a single equation form can work in many cases, it generates qualitatively wrong predictions in many other cases.
Nonetheless, there may be many situations where they are still useful, either by providing approximate description of the dynamics, or by capturing important qualitative features, such as the existence of oscillations, the set of coexisting species, presence of alternative stable states, etc.
Yes, we agree. Please see the fifth paragraph of the Discussion. We have also added the following paragraph to the Introduction:
“LV pairwise models are popular. LV pairwise modeling has successfully explained the oscillatory dynamics of hare and its predator lynx (Figure 1—figure supplement 1) (Volterra 1926; Wangersky 1978; “BiologyEOC – PopulationChanges” 2016). […] In the second example, if herbivores (mediators of competitive interactions between carnivores) rapidly reach steady state, herbivores can be mathematically eliminated from the mechanistic model to yield a pairwise model of competing carnivores (MacArthur 1970; Chesson 1990).
To me, understanding the conditions under which pairwise models are not even approximately correct, or make qualitatively wrong prediction is one of the important outstanding challenges in community modeling. I would encourage the authors to use their approach to provide insight into these questions, but, at the very least, they should clarify the difference between different failures of pairwise models. Additionally, they should be more explicit about cases in which no pairwise model would capture the community dynamics (e.g. in cases of interaction modification, or higher order interaction), versus ones in which the pairwise model that was considered wasn't adequate, but there may be a different one that is.
We have reorganized our text so that in the Discussion, we recap the two failure reasons of pairwise modeling: i) Interaction modification/higher order interactions in communities of three or more species violate the additivity assumption (this is previously known; Figure 7). In this case, even if each pairwise species interaction can be described by a pairwise model, the multispecies pairwise model will fail. ii) A single equation form cannot describe diverse pairwise interactions between microbes, thus violating the universality assumption (this is the focus of our paper; Figures 3, 4, 5 and 6).
We would like to emphasize that we are not trying to argue that “no pairwise model would capture the community dynamics”. This is likely true for complex interactions mediated by chemicals, but a strong statement like this would require a mathematical proof, which is beyond the scope of this paper. Instead, we argue that (last paragraph of Results):
“To summarize our work, even for pairwise microbial interactions, depending on interaction mechanisms (reusable versus consumable mediator, single mediator versus multiple mediators), we will need to use a plethora of pairwise models to avoid qualitative failures in predicting which species dominate a community or whether species coexist (Figures 3, 4 and 5). […] Taken together, pairwise modeling is unlikely to be effective for predicting community dynamics especially if interaction mechanisms are diverse.”
The authors also do not consider cases where interactions are mediated by externally supplied, abiotic mediators, rather than ones produced by the species themselves. The author's model is an extension of MacArthur's competition model, whose link to pairwise models have been extensively studied (e.g. Chesson, P. "MacArthur's Consumer Resource Model." Theoretical Population Biology 37, no. 1 (1990): 2638.).
What is not captured is competition of abiotic resources, such as the models common in David Tilman's work. The authors should make these distinctions clear, and put their work in the context of previous work. I also believe that extending the work to include competition for abiotic resources would add significant value to the work.
Originally, we intentionally left out the influence of external nutrients (e.g. competition for abiotic resources) so that we could focus on interactions mediated by a single chemical. If failure of pairwise modeling is established in these simpler cases, then our conclusion on the limitation of pairwise modeling is already valid. Nevertheless, we have now added a section on two species competing for an abiotic resource while engaging in a commensal interaction. Again, we observed that depending on parameter choices of the mechanistic model, the bestfitting pairwise model would work in some communities, but fail in other communities (Figure 6).
Reviewer #2:
The work of Momeni et al. "The validity of pairwise models in predicting community dynamics" explores conditions under which mechanistic models of species interacting via chemical mediators can be reduced to pairwise models. Pairwise models do not require a full mechanistic understanding of the nature of interactions within the community and thus use fewer parameters, which is why they are often used. The authors report that in many cases, pairwise models fail to predict, quantitatively or qualitatively, the dynamics of many species communities, which is why they should be used with caution.
Our focus is not “manyspecies communities.” Instead, we show that a single pairwise equation form (as used in pairwise modeling) cannot capture diverse pairwise chemicalmediated microbial interactions. Please see the last paragraph of the Results.
The authors do a good job in exploring the various scenarios under which pairwise models are not sufficient to capture these dynamics.
A point that deserves more attention however is how these results are immediately relevant for past and future studies of microbial communities?
We have moved an originally supplementary figure to the main text. In this figure (the new Figure 2), we list many examples of microbial interactions mediated by reusable or consumable compounds, which cannot be modeled by a single form of pairwise model (Figures 3 and 4 and Figure 3—figure supplement 4). Multimediator interactions, which pose challenges for pairwise modeling (Figure 5), are also a norm rather than an exception:
“A species often affects another species via multiple mediators (Kato et al. 2008; Yang et al. 2009; Traxler et al. 2013; Kim, Lee, and Ryu 2013). For example, a fraction of a population might die and release numerous chemicals, and some of these chemicals can simultaneously affect another individual.”
As per reviewer 1’s request, we have also added an example where two species compete for an abiotic resource while engaging in a commensal interaction. Again, we find that a pairwise model sometimes works but sometimes fails (Figure 6). Given the ubiquity of various interaction mechanisms that we consider in our work, our results are immediately relevant for studies of microbial communities.
The authors provide references to several studies that employed pairwise models to simulate community dynamics. However, many of those do not model chemically mediated interactions, but rather direct interactions, such as the predator prey examples mentioned in the introduction and are thus not directly relevant to what is studied here.
We have added this to our Introduction:
“Pairwise models are often used to predict how perturbations of steadystate species composition exacerbate or decline over time (May 1972; Thébault and Fontaine 2010; Mougi and Kondoh 2012; Allesina and Tang 2012; Suweis et al. 2013; Coyte, Schluter, and Foster 2015). Although most work are motivated by contactdependent preypredation (e.g. harelynx) or mutualisms (e.g. plantpollinator) where LV models could be identical to mechanistic models, these work do not explicitly exclude chemicalmediated interactions where species are distinct from interaction mediators. “
Relevant details of these papers are below:
(May, Nature1972)
This paper models densitydependent selfinhibition for each species (a_{ii}< 0), and positive or negative interspecies interactions (a_{ij}). The work is abstract, with no mentioning of particular interaction mechanisms.
(Thébault and Fontaine, Science2010)
This work examines two types of networks: trophic (all interactions are +/) and mutualistic (all interactions are +/+). For both cases, a logistic LV form is used for intraspecies competition and a saturable LV form is used for the effect of interspecies interactions. Aside from pestplant and plantpollinator images to illustrate trophic and mutualistic interactions, there is no mentioning that such model may not be applied to interactions where interaction mediators are distinct from species themselves.
(Mougi and Kondoh, Science2012)
This work argues that studies on purely mutualistic or purely trophic networks may not be adequate and examines communities with both types of interactions. The model is fairly general, with one assumption that distinguishes it from previous models: "with a biologically feasible assumption that interaction strengths decrease with increasing resource species, due to an allocation of interacting effort". This makes it closer to our alternative model (Eq. 5). The work does not mention that such models can fail for microbial communities when interactions are mediated by chemicals.
(Allesina and Tang, Nature2012)
This works uses the same framework as May, but examines what would happen when interactions are mostly trophic or mostly mutualistic (or a mix), and when more realistic network structures are imposed. They broadly speak of ecological networks, without mentioning specific mechanisms that are valid or not valid.
(Suweis et al., Nature2013)
This work examines mutualistic communities (mentioning plantanimal interactions in their motivations). This work uses two types of equations: Holling Type I (linear LV, similar to May's equations) and Holling Type II (nonlinear, similar to (Thebault and Fontaine, 2010) or the saturable LV). Beyond occasionally mentioning plantanimals (as an example of cooperation), they don't delve into which interaction mechanisms are valid or invalid.
(Coyte et al., Science2015)
This paper specifically mentions human gut microbiome as their motivation and as examples throughout. They use equations similar to May's for their analysis.
We do not intend to undermine these works by suggesting that they should have made this distinction clear. However, without an explicit analysis such as offered by this work, pairwise models have been applied liberally to any, including microbial, communities where interaction mediators are chemicals and not species.
Can the authors provide more specific examples of works that employ pairwise models to model chemically mediated interactions?
Given the universality of chemicalmediated interactions in microbial communities (Figure 2), a sizable fraction of microbial interactions are mediated by chemicals. In fact, the field of “secreted metabolomics” is devoted to identifying chemical interaction mediators. We have added this to the Introduction:
“The temptation of using pairwise models is indeed high, including in microbial communities where many interactions are mediated by chemicals (Mounier et al. 2008; Faust and Raes 2012; Stein et al. 2013; Marino et al. 2014; Coyte, Schluter, and Foster 2015).”
Coyte, K.Z., Schluter, J., and Foster, K.R. (2015). The ecology of the microbiome: Networks, competition, and stability. Science 350, 663–666.
Faust, K., and Raes, J. (2012). Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538–550.
Marino, S., Baxter, N.T., Huffnagle, G.B., Petrosino, J.F., and Schloss, P.D. (2014). Mathematical modeling of primary succession of murine intestinal microbiota. Proc. Natl. Acad. Sci. 111, 439–444.
Mounier, J., Monnet, C., Vallaeys, T., Arditi, R., Sarthou, A.S., Hélias, A., and Irlinger, F. (2008). Microbial interactions within a cheese microbial community. Appl. Environ. Microbiol. 74, 172–181.
Stein, R.R., Bucci, V., Toussaint, N.C., Buffie, C.G., Rätsch, G., Pamer, E.G., Sander, C., and Xavier, J.B. (2013). Ecological Modeling from TimeSeries Inference: Insight into Dynamics and Stability of Intestinal Microbiota. PLoS Comput Biol 9, e1003388.
This list is by no means exhaustive.
The authors clearly describe several regimes under which pairwise models fail. However, to identify whether the community being modeled falls into these regimes, one often needs to have a quite detailed mechanistic understanding of the interactions within the community in the first place. What would then be the advantages of using a pairwise model, given that this information is available?
Please see the fifth paragraph of the Discussion. Also, in Discussion:
“In summary, under certain circumstances, we may already know that microbial interaction mechanisms fall within the domain of validity for a particular pairwise model. In these cases, pairwise modeling provides the appropriate level of abstraction, and constructing a pairwise model is much easier than a mechanistic model (Figure 1).”
Also, it is not clear if this is an exhaustive list of regimes that they look at or are there other regimes out there to be explored? Thus, we are left with a conundrum that the work does not address: why the results are useful for systems where we do not know all the details of the interactions?
The point of our work is to caution against indiscriminative use of pairwise modeling. Our work is by no means an exhaustive list of regimes that work or not work for pairwise modeling. Nor is this our goal. By demonstrating that a single form of pairwise equation, as traditionally used in pairwise modeling, cannot qualitatively capture diverse interaction mechanisms between two microbial species, we reveal an unknown/undiscussed limitation of pairwise modeling. We hope that our work will help putting pairwise modeling to proper uses (Discussion, fifth paragraph). Even though our work emphasizes the importance of delineating interaction mechanisms, we do not argue that all details of interactions must be known for a model to be useful. In fact, a good model abstracts away nonessential details. For example, interactions mediated by a single reusable mediator can be described by a saturable LV model without the knowledge of mediator identity. Our Figure 8 provides another example of omitting unnecessary details (Figure 8D is equivalent to Figure 8B). We also added to the Introduction:
“Can a single pairwise model traditionally employed in pairwise modeling qualitatively describe diverse interactions between two microbial species?[…] On the other hand, pairwise models often fail to predict species coexistence in sevenspecies microbial communities (Friedman, Higgins, and Gore 2016), but this could be due to interaction modification discussed above.”
In general, for dynamical systems it is not at all surprising that 'the devil is in the details' of interactions for many a system, as weak couplings can often have important roles. Especially for systems that are under evolutionary pressure I can see how small couplings can be essential. Thus, models with underdetermined number of parameters will often fail, especially as long as the main state variables of a system are unknown, or when these tend to shift in importance from system to system.
Thus, to me the strong 'punch' in the message of the paper is missing, especially for a wider audience outside of the ecology crowd.
As a reviewer I am part of this outside audience I am not an ecologist nor am I a theorist.
You are right in that models with underdetermined number of parameters and state variables unsurprisingly fail (although we provide as many cases where they work as cases where they fail to work, see Figures 3, 4, 5, 6 and 7). According to this view, then no pairwise modeling paper should be published. However, the reality is that these papers are published in highprofile generalaudience journals and that pairwise model remains one of the most popular tools in theoretical ecology. To our knowledge, our finding (a single pairwise model cannot describe diverse pairwise microbial interactions) has not been discussed in the literature. “One man’s trivial is another man’s fundamental”, although here we are not talking about “another man” but many other men! Our work points out a fundamental flaw in arguably one of the most influential models in ecology when applied to microbial communities, and thus deserves to be in a highprofile, generalaudience journal.
https://doi.org/10.7554/eLife.25051.051Article and author information
Author details
Funding
Boston College
 Babak Momeni
NIH Office of the Director
 Babak Momeni
 Wenying Shou
W. M. Keck Foundation
 Babak Momeni
 Li Xie
 Wenying Shou
Fred Hutchinson Cancer Research Center
 Li Xie
 Wenying Shou
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Reviewing Editor
 Bruce Levin, Emory University, United States
Publication history
 Received: January 11, 2017
 Accepted: March 18, 2017
 Accepted Manuscript published: March 28, 2017 (version 1)
 Accepted Manuscript updated: April 3, 2017 (version 2)
 Version of Record published: June 13, 2017 (version 3)
Copyright
© 2017, Momeni et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,849
 Page views

 806
 Downloads

 36
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Computational and Systems Biology
 Neuroscience

 Computational and Systems Biology
 Genetics and Genomics