The validity of pairwise models in predicting community dynamics

Pairwise models are commonly used to describe many-species communities. In these models, a focal species receives additive fitness effects from pairwise interactions with other species in the community (“pairwise additivity assumption”), and all pairwise interactions are represented by a single canonical equation form (“universality assumption”). Here, we analyze the validity of pairwise modeling. We build mechanistic reference models for chemical-mediated interactions in microbial communities, and attempt to derive corresponding pairwise models. Even when one species affects another via a single chemical mediator, different forms of pairwise models are appropriate for consumable versus reusable mediators, with the wrong model producing qualitatively wrong predictions. For multi-mediator interactions, a canonical model becomes even less tenable. These results, combined with potential violation of the pairwise additivity assumption in communities of more than two species, suggest that although pairwise modeling can be useful, we should examine its validity before employing it.


Introduction
Multispecies microbial communities are ubiquitous. Microbial communities are important for industrial applications such as cheese and wine fermentation (van Hijum, Vaughan, and Vogel 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) or cause diseases (Kelly 1980). Community-level properties (e.g. species composition and biochemical activities) cannot be achieved, or achieved to the same extent, by summing the contributions of individual member species. Community-level properties are influenced by interactions wherein individuals alter the physiology of other individuals. To understand and predict properties of communities, choosing the appropriate mathematical model to describe species interactions is critical.
Two commonly-used modeling approaches are mechanistic modeling and pairwise modeling, each with its pros and cons. In mechanistic modeling, interaction mechanisms are explicitly modeled (Fig 1A and B, left panels). Thus, a mechanistic model requires discovering and quantifying interaction mechanisms (Fig 1 Table, "parameter" rows under "Mech." column).
2 Such a mechanistic model can in principle quantitatively predict community dynamics when species evolution is negligible. However, the complexity of microbial interactions and the difficulty in identifying and quantifying interactions have made it challenging to construct mechanistic models.
In contrast to mechanistic modeling, pairwise modeling considers only the fitness effects of pairwise species interactions (Figs 1A and B,right panels). Pairwise models have two central assumptions. First, the "universality" assumption: Regardless of interaction mechanisms, how one species affects another can be abstracted into a single canonical equation form so that only parameters can vary among interactions. Second, the "pairwise additivity" assumption: a focal species receives additive fitness effects from pairwise interactions with other species in the community. 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.
Pairwise models are easy to construct because they do not require knowledge of interaction mechanisms and need fewer parameters than mechanistic models (Fig 1 table). Parameters are relatively easy to estimate using community dynamics (Stein et al. 2013), or more systematically, using dynamics of monocultures and pairwise cocultures (Fig 2).
Not surprisingly, pairwise modeling has been commonly applied to communities (Wootton and Emmerson 2005). Pairwise models are often justified by their success in predicting ecological dynamics of two-species communities of prey-predation (Fig 1-FS1) (Volterra 1926;Wangersky 1978; "BiologyEOC -PopulationChanges" 2016) and competition (Gause 1934a;Gause 1934b). Pairwise modeling has been extended to model communities of more than two species (defined as "multispecies communities"), with empirical support from, for example, an artificial community of four competing protozoa species (Vandermeer 1969). Multispecies pairwise models have been extensively used to predict how perturbations to steady-state species composition exacerbate or decline over time (May 1972;Cohen and Newman 1984;Pimm 1982;Thébault and Fontaine 2010;Mougi and Kondoh 2012;Allesina and Tang 2012;Suweis et al. 2013;Coyte, Schluter, and Foster 2015).
However, pairwise modeling has known limitations. For instance, in a multispecies community, an interaction between two species can be altered by a third species (Levine 1976;Tilman 1987;Wootton 2002;Werner and Peacor 2003;Stanton 2003). Indirect interactions via a third species fall under two categories (Wootton 1993), which can be illustrated using the example of carnivore, herbivore, and plant. In an "interaction chain" (also known as "density-mediated indirect interactions"), a carnivore affects the density of an herbivore which in turn affects the density of plants. In "interaction modification" (also known as "trait-mediated indirect interactions" or "higher order interactions" (Vandermeer 1969;Wootton 1994;Billick and Case 1994;Wootton 2002)), a carnivore affects how often an herbivore forages plants. Interaction modification (but not interaction chain) violates the pairwise additivity assumption (Methods-Section 1). Interaction modification is thought to be common in ecological communities (Werner and Peacor 2003;Schmitz, Krivan, and Ovadia 2004). Limitation of pairwise modeling has also been studied experimentally (Dormann and Roxburgh 2005). However, empirically-observed failure of multispecies pairwise models could be due to limitations in data collection and analysis (Case and Bender 1981;Billick and Case 1994).
Given the benefits, limitations, and intellectual influence of pairwise modeling, we examine conditions under which pairwise models produce realistic predictions. Instead of investigating natural communities where interaction mechanisms can be difficult to identify, we start with in silico communities where species engage in predefined chemical interactions of the types commonly encountered in microbial communities. Based on these interactions, we construct mechanistic models, 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 tools; and the flexibility in creating different reference models allows us to explore a variety of conditions. This has allowed us to examine the domain of validity for pairwise modeling.

Results
Establishing a mechanistic reference model In our mechanistic models (Fig 1A, left), we focus on chemical interactions which are widespread in microbial communities (Fig1-FS2) (Stams 1994;Czárán, Hoekstra, and Pagie 2002;Duan et al. 2009). A mechanistic model includes a set of species as well as chemicals that mediate interactions among species. Species Si could release or consume chemical Cj, and chemical Cj could increase or decrease the growth rate of species Sk.
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 modeling. Additive fitness effects have been observed for certain "homologous" metabolites. For example, in multi-substrate carbon-limited chemostats of E. coli, the fitness effects from glucose and galactose were additive (Lendenmann and Egli 1998). "Heterologous" metabolites (e.g. carbon and nitrogen sources) likely affect cell fitness in a multiplicative fashion. However, if released mediators cause small changes to the concentrations already in the environment, then additivity approximation may still be acceptable. For example, suppose that the fitness influences of released carbon and nitrogen with respect to those already in the environment are wc and wn, respectively. If wc, wn<<1, the additional relative fitness influence will be (1+wc)(1+wn)-1 ≈ wc+wn. However, even among homologous metabolites, fitness effects may not be additive (Hermsen et al. 2015). "Sequential" metabolites (e.g. diauxic shift) provide such examples.
Si and Ci are state variables representing the concentrations of species Si and chemical Ci, respectively. ri0 is the basal fitness of an individual of species Si (the net growth rate of a single individual in the absence of any intra-species or inter-species interactions  (Gopalsamy 1992), and saturable Lotka-Volterra, where the fitness effect of Sj on Si saturates at high density of Sj (Thébault and Fontaine 2010). Since we model continuous growth which does not impose carrying capacity and since chemical influence from one species to another is likely saturable, we have adopted the saturable Lotka-Volterra as our canonical pairwise model: Here, rij is the maximal positive or negative fitness effect of Sj on Si, and Kij is the Sj exerting half maximal fitness influence on Si (parameter definition in Fig 1 table). When j=i, nonzero rii and Kii reflect density-dependent growth effect in Si.
From a mechanistic model, we derive a pairwise model either analytically or numerically ( Fig  2A). In the latter case ( Fig 2B-C, Methods-Section 2), we should already have a pre-specified pairwise model (e.g. the canonical pairwise model) in mind. We then use the dynamics of monocultures and pairwise cocultures obtained from the mechanistic model to find parameters that minimize the difference between the two models within a training time window T. Specifically, we define a distance measure D as the fold-difference between the dynamics from the two models, averaged over any time interval  and species number N ( Fig 2C): Here Si,pair and Si,mech are Si calculated using pairwise and mechanistic models, respectively. Since species with densities below a set extinction limit, Sext, are assumed to have gone extinct in the model, we set all densities below the extinction limit to Sext in calculating D to avoid singularities. Within the training window T, minimizing () DT using a nonlinear least square routine yields parameters of the best-matching pairwise model. We then use D outside the training window to quantify how well the best-matching pairwise model predicts the mechanistic model.

Reusable and consumable mediators are best represented by different forms of pairwise models
To build a pairwise model, we must accurately represent the fitness influence of one species on another (rij and Kij). Even though this basic process seems straightforward as outlined in Fig 2, in practice, challenges may arise. For example, identifying the set of best-matching parameters for nonlinear functions may not be straightforward, and measurement errors further hamper parameter estimation. Partly due to these challenges, studies on deriving pairwise model parameters for a given community are scarce (Pascual and Kareiva 1996;Stein et al. 2013), despite the popularity of pairwise models. In this section, we analytically derive pairwise models from mechanistic models of two-species 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;N.S. Jakubovics 2010) or consumable such as metabolites (Stams 6 1994;Freilich et al. 2011). To facilitate mathematical analysis, we consider community dynamics within a dilution cycle. We show that a single canonical pairwise model may not encompass these different interaction mechanisms.
Consider a commensal community where species S1 stimulates the growth of species S2 by producing a reusable ( Fig 3A) or a consumable chemical C1 (Fig 3B). When C1 is reusable, the mechanistic model can be transformed into a pairwise model (Fig 3A), provided that the concentration of the mediator (which is initially zero) has acclimated to be proportional to the producer population size (Fig 3A; Fig 3-FS1). This pairwise model takes the canonical form (compare with Fig 1B right). Thus, the canonical pairwise model is appropriate, regardless of whether the producer coexists with the consumer, outcompetes the consumer, or is outcompeted by the consumer.
If C1 is consumable, different scenarios are possible within a dilution cycle (Methods-Section 3). Case I: When supplier S1 always grows faster than consumer S2 ( 21 10 20 SC r r r  ), C1 will accumulate within each dilution cycle (proportional to S1) without plateauing to a steady state (Fig 3-FS2  ( , ) 1 0 (the "f-zero-isocline", blue line in Fig 3-FS4 where ω and ψ are constants (Fig 3B-ii). As expected, parameter estimation for the alternative pairwise model is more accurate after the initial period of time (Fig 3-FS6 bottom panels).
Case III: When supplier S1 always grows slower than consumer S2 ( 10 20 rr  ), Eqs. S3-2 and S3-5 above are still valid. Thus, consumable C1 declines to zero concentration as it is consumed by S2 whose relative abundance over S1 eventually exponentially increases at a rate of is required for the alternative pairwise model to work (Fig3-FS8A).
Otherwise, we can get qualitatively wrong results (Fig 3-FS5B and D).
If any of the requirements on initial conditions or time scale (Methods-Section 3; Fig3-FigS8) is violated, the alternative pairwise model may fail to represent the mechanistic model. For example, if dilution cycles are such that the community can never approach the f-zero-isocline, then C1 could accumulate proportionally to S1 even if 21 10 20 0 SC r r r    . In this case, the canonical but not the alternative pairwise model is appropriate (Fig 3-FS7). Similarly, a gentler dilution scheme, which allows the community to remain near the f-zero-isocline, leads to better parameter estimation for a pairwise model (Fig 3-FS6). if additionally, the half-saturation constant K for C1 consumption ( 12 CS K ) is identical to that for C1's influence on the growth of consumer ( 21 SC K ) (the "K assumption"), and if S2 has not gone extinct. This equation form has precedence in the literature (e.g. (Mougi and Kondoh 2012)), where the interaction strength r21 reflects the fact that the consumable mediator from S1 is divided among consumer S2. 8 The canonical and the alternative models are not interchangeable. The alternative pairwise model is not predictive of community dynamics where C1 accumulates without reaching a steady state (Fig 3-FS2, compare dashed and solid lines). Similarly, an interaction mediated by a consumable mediator that eventually reaches steady state can often be described by the alternative pairwise model (Fig 3-FS3, dashed lines). However, if we were to use a canonical model, model accuracy becomes uncertain: Parameters estimated at the steady state can predict community dynamics at steady state (Fig 3-FS3A), but not community dynamics when initial species ratios differ from the steady state ratio (Fig 3-FS3B and C, compare dotted with solid lines).
We have shown here that even when one species affects another species via a single mediator, either a canonical or an alternative form of pairwise model may be more appropriate in different situations. Choosing the appropriate model depends on whether the mediator is reusable or consumable and if consumable, how the fitness of two species compare. Additionally, a pairwise model may approximate the mechanistic model properly only under appropriate initial conditions and after the estimated time scale (Fig 3-FS5 ,B and D;. Considering that reusable and consumable mediators are both common, our results call for revisiting the universality assumption of pairwise modeling.

Multi-mediator interactions require pairwise models different from single-mediator interactions
A species often affects another species via multiple mediators (Kato et al. 2008;Traxler et al. 2013;Kim, Lee, and Ryu 2013). For example, a subpopulation from one species might die and release numerous chemicals that can affect another species in different ways. Here we examine cases where S1 releases two chemicals C1 and C2 which additively affect the growth of S2 (Fig  4). We ask when two mediators can mathematically be regarded as one mediator (to facilitate further abstraction into a pairwise model), and how multi-mediator interactions affect pairwise modeling.
When both mediators are reusable (Methods-Section 4), their combined effect ( ) generally cannot be modeled as a single mediator except under special conditions (Fig 4). These special conditions (Methods-Section 4) include: (1) mediators share similar "potency" (Fig 4C, diagonal), or (2) one mediator has much stronger "potency" than the other (i.e. one mediator dominates the interaction; note the log scale in Fig  4C).
When both mediators are consumable as in Cases II and III, the interaction term becomes ). All these forms deviate from the canonical form.
In summary, when S1 influences S2 through multiple mediators, rarely can we approximate them as a single mediator. Multiple mediators generally make equations of pairwise modeling more complex than single mediators, casting further doubt on the usefulness of a universal form for all community interactions.

A multispecies pairwise model can work for interaction chains but generally not for interaction modifications
For a community with more than two species, can we construct a multispecies pairwise model from two-species pairwise models? The answer is yes for an interaction chain mediated by chemicals ( Fig 5A), so long as mediators between different species pairs are independent and each species pair can be represented by a pairwise model. The equation form of the multispecies pairwise model can vary (Methods-Section 5), as discussed in previous sections.
Consistent with previous work (Methods-Section 1), interaction modification can cause a multispecies pairwise model to fail. For example, S1 releases C1 which stimulates S2 growth; C1 is consumed by S3 and stimulates S3 growth ( Fig 5C). Here, the presence of S3 changes the strength of interaction between S1 and S2, an example of interaction modification. Viewing this differently, S1 changes the nature of interactions between S2 and S3: S2 and S3 do not interact in the absence of S1, but S3 inhibits S2 in the presence of S1. This causes the three-species pairwise model to make qualitatively wrong conclusions about species persistence even though each species pair can be described by a pairwise model ( Fig 5D). As expected, if S3 does not remove C1, the three-species pairwise model works ( Fig 5-FS1, A-B).
Interaction modification can occur even in communities where no species changes "the nature of interactions" between any other two species ( Fig 5E). Here, both S1 and S3 contribute reusable C1 to stimulate S2. S1 promotes S2 regardless of S3; S3 promotes S2 regardless of S1; S1 and S3 do not interact regardless of S2. However, a multispecies pairwise model assumes that the fitness effects from the two producers on S2 will be additive, whereas in reality, the fitness effect on S2 saturates at high C1. As a result, even though the dynamics of each species pair can be represented by a pairwise model (Fig 5F right, purple), the three-species pairwise model fails to capture community dynamics ( Fig 5F). Thus, the nonlinearity in how a mediator affects a species can also violate the additivity assumption of a pairwise model. As expected, if C1 affects S2 in a linear fashion, the community dynamics is accurately captured in the multispecies pairwise model (Fig 5-FS1, C-D).
In summary, for chemical-mediated indirect interactions, a multispecies pairwise model can work for the interaction chain category but generally not the interaction modification category.

Discussion
Multispecies pairwise models are widely used in theoretical research due to their simplicity. In two-species interactions such as prey-predation based on contact-dependent inhibition (instead of diffusible chemical mediators), Lotka-Volterra pairwise models can in fact be the mechanistic representation of interactions and thus predictive of community dynamics (Fig 1-FS1). The inadequacy of multispecies pairwise models has been discussed theoretically (Wootton 2002;Wootton and Emmerson 2005) and empirically (Case and Bender 1981;Dormann and Roxburgh 2005;Aschehoug and Callaway 2015), although the reasons for model failures in explaining particular experimental results are often unclear (Billick and Case 1994).
Here, we have considered the validity of pairwise models in well-mixed two-and three-species communities where all species interactions are known and thus community dynamics can be described by a mechanistic reference model. We have focused on chemical-mediated interactions commonly encountered in microbial communities (Fig 1-FS2) (Kato et al. 2005;Gause 1934a;Ghuysen 1991;Nicholas S Jakubovics et al. 2008;Chen et al. 2004;D'Onofrio et al. 2010;Johnson et al. 1982;Hamilton and Ng 1983). To favor the odds of successful pairwise modeling, we have also assumed that different chemical mediators exert additive fitness effects on a target species.
What are the conditions under which the influence of one species on another can be represented by a canonical two-species pairwise model (the universality assumption)? When an interaction employs a single mediator, then a canonical saturable pairwise model ( Fig 1B) will work for a reusable mediator (Fig 3A). For a consumable mediator, depending on the fitness relationship between the two species, either a canonical ( Fig 3A) or an alternative ( Otherwise, qualitatively wrong predictions can ensue (Fig 3-FS5, B and D). In all cases, acclimation time is required for a pairwise model to converge to the mechanistic model (Fig3-FS8B; Fig3A). If one species influences another through multiple mediators, then in general, these mediators may not be regarded as a single mediator (Fig 4 and Methods-Section 4). Thus, conditions for a working canonical pairwise model become even more restrictive.
In communities of more than two species, indirect interactions via a third species can occur. For indirect interactions in the form of interaction chains, as long as each two-species segment of the chain engages in independent interactions and can be represented by a pairwise model, then multispecies pairwise models will generally work (Figs 5A-B, Methods-Section 5). However, depending on whether each mediator is reusable or not, equation forms will vary. For indirect interactions in the form of interaction modification (higher-order interactions), even if each species pair can be accurately represented by a pairwise model, a multispecies pairwise model may fail (Fig 5, C-F). Interaction modification includes trait modification (Wootton 2002;Werner and Peacor 2003;Schmitz, Krivan, and Ovadia 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, Kelsic, and Kishony 2016). Moreover, interaction mediators are often 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 by-products are also consumed by multiple species (Kolenbrander 2000).
Pairwise modeling (or variations of it) still has its uses when simulating a particular community phenomenologically. One can even imagine that an extended pairwise model (e.g. can serve as a general-purpose model for pairwise interactions via a single mediator. Even the effects of indirect interactions may be quantified and included in the model by incorporating higher-order interaction terms (Case and Bender 1981;Worthen and Moore 1991), although many challenges will need to be overcome (Wootton 2002). In the end, although these strategies may lead to a sufficiently accurate phenomenological model for specific cases, "one-formfitting-all" may generate erroneous predictions when modeling different communities.
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 (Levins 1966;Durrett and Levin 1994; Damore and Gore 2012)? Tilman argued that if a small number of mechanisms (e.g. the "axes of trade-offs" in species' physiological, morphological, or behavioral 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 S1 grows exponentially (not explicitly depicted in equations in Fig 6) and the growth rate of S2, which is normally zero, is promoted by mediator C from S1 in a linear fashion (Fig 6). If we do not know how S1 stimulates S2, we can still construct a pairwise model ( Fig  6A). If we know the identity of mediator C and realize that C is consumable, then we can instead construct a mechanistic model incorporating C (Fig 6B). However, if C is produced from a precursor via an enzyme E released by S1, then we get a different form of mechanistic model ( Fig 6C). If, on the other hand, E is anchored on the membrane of S1 and each cell expresses a similar amount of E, then equations in Fig 6D are mathematically equivalent to Fig 6B. 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 interaction mechanisms fall within the domain of validity for a pairwise model. In these cases, pairwise modeling provides the appropriate level of abstraction, and constructing a pairwise model can be far easier than measuring the many parameters required by a mechanistic model. However, if we do not know whether pairwise modeling is valid, we will need to be cautious about indiscriminative use of pairwise models since they can fail to even qualitatively capture community dynamics (e.g. in We will need to be equally careful in extrapolating and generalizing conclusions obtained from pairwise models. Considering recent advances in identifying and quantifying interactions, we advocate a transition to models that incorporate interaction mechanisms at the appropriate level of abstraction.

Section 1. Interaction modification but not interaction chain violates the pairwise additivity assumption
In a pairwise model, the fitness of a focal species Si is the sum of its "basal fitness" (ri0, the net growth rate of a single individual in the absence of any intra-species or inter-species 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, () ij j fS describes how Sj, the density of species Sj, positively or negatively affects the fitness of Si, and is a linear or nonlinear function of only Sj and not of a third species.
Indirect interactions via a third species fall under two categories (Wootton 1993). The first type is known as "interaction chain" or "density-mediated indirect interactions". For example, the consumption of plant S1 by herbivore S2 is reduced when the density of herbivore is reduced by carnivore S3. In this case, the three-species pairwise model does not violate the pairwise additivity assumption (compare with Eq. S1-1) (Case and Bender 1981;Wootton 1994).
The second type of indirect interactions is known as "interaction modification" or "trait-mediated 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, Krivan, and Ovadia 2004). For example, when carnivore is present, herbivore will spend less time foraging and consequently plant density increases. In this case, f12 in Eq S1-2 is a function of both S2 and S3, violating the pairwise additivity assumption.

Section 2. Summary of simulation files
Simulations are based on Matlab® and executed on an ordinary PC. Steps are: Step 1: Identify monoculture parameters ri0, rii, and Kii ( Fig 2C, Row 1 and Row 2).
Step 3: Calculate distance D between population dynamics of the reference mechanistic model and the approximate pairwise model over several generations outside of the training window to assess if the pairwise model is predictive.
Fitting is performed using nonlinear least square (lsqnonlin routine) with default optimization parameters. The following list describes the m-files used for different steps of the analysis:

File name Function
FitCost_BasalFitness. Estimates alternative pairwise model interaction parameter (r21) in cases where we know that S2 is only affected by S1 and that KS2C1=KC1S2 to accelerate optimization Section 3. 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 Fig3-FS7 for an example where dilution cycles violate requirements for deriving a pairwise model). We further assume 10 0 r  and 20 0 r  so that species cannot go extinct in the absence of dilution. See Fig3-FS8 for a summary of this section.
When S1 releases a consumable mediator which stimulates the growth of S2, the mechanistic model as per Fig 3B, (Eq. S3-2).
Case I:

21
10 20 SC r r r  Since producer S1 always grows faster than consumer S2, 0 S R  as → ∞. Define Since However, if C1 has not yet reached steady state, imposing steady state assumption would falsely predict Rs at steady state and thus remaining at its initial value (Fig 3-FS3, dotted lines). Since 1 dC dt in Eq. S3-1 is the difference between two exponentially growing terms, we factor out the exponential term S1 to obtain 1 1 1 2 1 1 12 where ω and ψ are constants (Fig 3B, ii).
For certain conditions (which will be discussed at the end of this section, Fig3-FS8A), this alternative model can make reasonable predictions of community dynamics even before the community reaches the steady state (Fig 3-FS3  Equations S3-7 and S3-10 allow us to construct a phase portrait where the x axis is 1 C and the y axis is ˆS R (Fig 3-FS4A-D). Note that at steady state, 1( , ) S CR= (1, 1). Setting Eq. S3-9 to zero: C  . 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-zero-isocline during which the alternative pairwise model can be derived (Fig 3-FS4A-D).
It is difficult to solve Eq. S3-7 and S3-8 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 f t , the time scale for the community to approach the f-zero-isocline. Note that f t 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 f t is that 10 S is sufficiently high (Fig3-FS8B) to avoid the long lag phase that is otherwise required for the mediator to accumulate to a high enough concentration.
From Eq. S3-11, the asymptotic value for the f-zero-isocline is   12 1( This is plotted as a black dotted line in Fig 3-FS4

C S C S S
The approximation in the last step is due to the very definition of Case II-1: . The initial steepness of the community dynamics trajectory (Eq. S3-13) will be much smaller than that of the f-zero-isocline (Eq. S3-14) if     a typical trajectory moves toward the f-zero-isocline almost horizontally (Fig 3-FS4A). The unscaled form of Eq. S3-19 is Thus,  

CS s s s S CS C S C S C S S C S C s s S s C S C S SC C S C S S C S C s CS SC
Taking the derivative of both sides, and using Eq. S3-21 and Eq. S3-22, we have    

C S C S S C S C C S C S s s S C S C S SC C S C S C S S C S C CS SC
The solution to the above equation is: The second term approaches zero much faster compared to the first term due to the negative exponent with an exponential term. Thus,

S C S C C S S C C S C S C S S C S C C S
Similar to Case II-2, if Eq. S3-19 is satisfied, a typical trajectory is illustrated in Fig 3-FS4C where the trajectory is almost horizontal and 1 C increases to much greater than 1 before the system reaches the f-zero-isocline.

Case III:
10 20 rr  In this case, supplier S1 always grows slower than S2. As t , 21 S R S S    and 1 0 C  . The phase portrait is separated into two parts by the f-zero-isocline (Fig 3-FS4E), where, as in Eq. S3-5,   S is small, there is a lagging phase during which the trajectory rises steeply before leveling off (Fig 3-FS4H). Although the alternative pairwise model can be derived once the trajectory hits the f-zero-isocline, f t takes a complicated form. Here we consider two cases where 10 S is large enough so that we can approximate the trajectory as a straight line going through (0, RS(t=0)) (Fig3-FS4, F, G). Graphically, 10 S is large enough so 27 that the green segment in Fig3-FS4F, whose length is  S R , is much smaller than (0) S R . In other words, From Eqs. S3-1 and S3-2   12 12 1 2 1 1 1 20 10 Case III-2.

S C S C S R 
If 10 S is large enough so that Eq. S3-28 is satisfied, a typical example is displayed in Fig3-FS4G. The trajectory moves with a small positive slope so that the intersection of the community 28 dynamics trajectory with the f-zero-isocline is near the black dotted line

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 C1 cannot be specified in pairwise model, problems could occur. To illustrate, we examine the phase portrait of the pairwise equation Below, we plot Eq. S3-31 under different parameters (Fig 3-FS5) to reveal conditions for convergence between mechanistic and pairwise models.  If 0   (Fig 3-FS5B ): RS increases exponentially in mechanistic model (Eq. S3-2). Thus, consumable C1 will decline toward zero as C1 is consumed by S2 whose relative abundance over S1 exponentially increases. Thus, from mechanistic model (Eq. S3-2), RS eventually increases exponentially at a rate of 20 10 rr  .

Section 4. Conditions under which a pairwise model can represent one species influencing another via two reusable mediators.
Here, we examine a simple case where S1 releases reusable C1 and C2, and C1 and C2 additively affect the growth of S2 (see example in Fig 4).
D can be close to zero when (i)

Section 5. Multi-species pairwise model for an interaction chain.
Without loss of generality, we consider an example where each step of an interaction chain is best represented by a different form of pairwise model. Suppose that S1 releases a reusable mediator C1 that promotes S2 and that S2 releases a consumable mediator C2 that promotes S3. The interaction mediated by C2 can be described by the alterative pairwise model. The corresponding community-pairwise model then will be: To see this, note that after a transient period of time, the mechanistic model of the three-species community is: Equations S5-2 and S5-1 are equivalent. Figures   Fig 1. The abstraction of interaction mechanisms in a pairwise model compared to a mechanistic model. (A) The mechanistic model (left) considers a bipartite network of species and chemical interaction mediators. A species can produce or consume chemicals (open arrowheads pointing towards and away from the chemical, respectively). A chemical mediator can positively or negatively influence the fitness of its target species (filled arrowhead and bar, respectively). The corresponding pairwise model (right) includes only the fitness effects of species interactions, which can be positive (filled arrowhead), negative (bar), or zero (line terminus). (B) In the example here, species S1 releases chemical C1, and C1 is consumed by species S2 and promotes S2's fitness. In the mechanistic model, the three equations respectively state that 1) S1 grows exponentially at a rate r10, 2) C1 is released by S1 at a rate 11 CS  and consumed by S2 with saturable kinetics, and 3) S2's growth (basal fitness r20) is influenced by C1 in a saturable fashion. In the pairwise model here, the first equation is identical to that of the mechanistic model. The second equation is similar to the last equation of the mechanistic model except that r21 and K21 together reflect how the density of S1 (S1) affects the fitness of S2 in a saturable fashion. For all parameters with double subscripts, the first subscript denotes the focal species or chemical, and the second subscript denotes the influencer. Note that unlike in mechanistic models, we have omitted "S" from subscripts in pairwise models (e.g. r21 instead of 21 SS r ) for simplicity. Analytically deriving a pairwise model from a mechanistic model allows us to uncover approximations required for such a transformation (top). Alternatively (bottom), through a "training window" of the mechanistic model population dynamics, we can numerically derive parameters for a preselected pairwise model that best fits the mechanistic model. We then quantify how well such a pairwise model matches the mechanistic model under conditions different from those of the training window. (B) A mechanistic model of three species interacting via two chemicals (left) can be translated into a pairwise model of three interacting species (center). S1 inhibits S1 and promotes S2 (via C1). S2 promotes S2 and S3 (via C2) as well as S1 (via removal of C1). S3 promotes S1 and inhibits S2 (via removal of C1 and C2, respectively). Take interactions between S2 and S3 for example: the saturable Lotka-Volterra pairwise model will require estimating ten parameters (colored, right), some of which (e.g. r33 in this case) may be zero. (C) In the numerical method, the six monoculture parameters (ri0, rii, and Kii, i=2, 3; green and red) are first estimated from training window T (within a dilution cycle) of monoculture mechanistic models. Subsequently, the four interaction parameters (rij and Kij, i  j, olive) can be estimated from the training window T of the S2 + S3 coculture mechanistic model. Parameter definitions are described in Fig 1. To estimate parameters, we use an optimization routine to minimize D , the fold-difference (shaded area) between dynamics from a pairwise model (dotted lines) and the mechanistic model (solid lines) averaged over T. In all simulations, to ensure that resources not involved in interactions are never limiting, a community is diluted back to its inoculation density when total population increases to a high-density threshold. Too frequent dilutions will allow only small changes in population dynamics within a dilution cycle or T, which is not suitable for estimating pairwise models. Too infrequent and large dilutions will cause large fluctuations in dynamics, which can sometimes violate conditions for deriving pairwise models (Fig 3-FS6; Fig 3-FS7). Under most cases we have tested, small variations in dilution frequency do not affect our conclusions. See Methods-Section 2 for relevant Matlab codes. Interactions mediated via a single reusable or consumable mediator are best represented by different forms of pairwise models. S1 stimulates the growth of S2 via a reusable (A) or a consumable (B) chemical C1. In mechanistic models of the two cases, equations for S1 and S2 are identical but equations for C1 are different. In (A), C1 can be solved to yield 1 1 1 1 1 1 1 1 1 10 10 10 10 10 10 1 10 10 assuming zero initial C1. We have approximated C1 by omitting the second term (valid after the initial transient response has passed so that C1 has become proportional to S1). This approximation allows an exact match between the canonical pairwise model (Fig 1B, right) and the mechanistic model (ii), and thus justifies the pairwise model. In (B), depending on the relative growth rates of the two species, and if additional requirements are satisfied (Methods-Section 3; Fig 3-FS7; Fig 3-FS8), canonical or alternative pairwise model should be used. Thus, depending on whether the mediator is consumed or reused, the most appropriate pairwise model (colored) takes different forms. of Ci by S1 or low Ci required to achieve half-maximal influence on S2). (B) Under what conditions can an interaction via two reusable mediators be approximated by a pairwise model? (C) Under restricted conditions, two reusable mediators can be consolidated into a single mediator. We can directly compute the best-fitting pairwise model parameters over a training window of T by minimizing D (Methods-Section 4, Eq. S4-2). Here, the difference D between the two models over T =10 generations is plotted over a range of potencies KC1 and KC2. The canonical pairwise model is valid (blue regions indicating small difference) when KC1 ≈ KC2 or when one interaction is orders of magnitude stronger than the other interaction (Methods-Section 4). (D) A community where the canonical pairwise model is not valid. Here, KC1=10 3 and KC2=10 6 . We estimate the best-fitting pairwise model by minimizing D (Methods-Section 4, Eq. S4-2) in three training windows (spanning 10 generations of growth for S1). At various S1, we calculate the fitness effect of S1 on S2 using the pairwise model and the mechanistic model (B). In two of the three training windows, the two models fail to match. In the training window with the lowest S1, the two models match because the effect of C2 is negligible in this range (KC2>>S1, condition iib in Methods-Section 4). These mismatches mean that a pairwise model cannot consistently capture reference dynamics. Simulation parameters are listed in Fig 4-SD1.  Fig 5. Interaction chain but not interaction modification may be represented by multispecies pairwise model. We examine three-species communities engaging in indirect interactions. Each species pair is representable by a two-species pairwise model (purple in the right columns of B, D, and F). We then use these two-species pairwise models to construct a three-species pairwise model, and test how well it predicts the dynamics from mechanistic model. In B, D, and F, left panels show dynamics from the mechanistic models (solid lines) and three-species pairwise models (dotted lines). Right panels show the difference metric D calculated over population densities after taking dilution into consideration. (A-B) Interaction chain: S1 affects S2, and S2 affects S3. The two interactions employ independent mediators C1 and C2, and both interactions can be represented by the canonical pairwise model. The three-species pairwise model matches the mechanistic model in this case. Simulation parameters are provided in Fig 5-SD1. (C-F) Interaction modification. (C-D) S3 consumes C1, a mediator by which S1 stimulates S2. Parameters are listed in Fig 5-SD2. (E-F) S1 and S3 both supply C1 which stimulates S2. Simulation parameters are listed in Fig 5-SD3. In both interaction modification cases, the three-species pairwise model fails to predict reference dynamics. How one species (S1) may influence another (S2) can be mechanistically modeled at different levels of abstraction. For simplicity, here we assume that interaction strength scales in a linear (instead of saturable) fashion with respect to mediator concentration or species density. The basal fitness of S2 is zero. (A) In the simplest form, S1 stimulates S2 in a pairwise model. (B) In a mechanistic model, we may realize that S1 stimulates S2 via a mediator C which is consumed by S2. The corresponding mechanistic model is given. (C) Upon probing more deeply, it may become clear that S1 stimulates S2 via an enzyme E, where E degrades an abundant precursor (such as cellulose) to generate mediator C (such as glucose).
In the corresponding mechanistic model, we may assume that E is released by S1 at a rate 1 ES  and that E liberates C at a rate CE  . (D) If instead E is anchored on the cell surface (e.g. in cellulose degradation via cellulosome), then E is proportional to S1. If we substitute E into the second equation, then (B) and (D) become equivalent. Thus, when enzyme is anchored on cell surface but not when enzyme is released, the mechanistic knowledge of enzyme can be neglected.   We use the mechanistic model for a reusable mediator to generate reference dynamics of S1, S2, and C1 over 150 generations of community growth. The basal fitness of S1 and S2 in pairwise models are identical to those in mechanistic models, and here rii or Kii (i = 1, 2) are irrelevant due to the lack of intrapopulation interactions. We use every 10 community doublings of reference dynamics as training windows to numerically estimate best-matching canonical pairwise model parameters r21 and K21. Dashed and solid rectangles represent training windows before and after acclimation, respectively. Note that population fractions (instead of population densities) are plotted, which fluctuate less during dilutions compared to mediator concentration. (B) Pairwise model parameters estimated after acclimation (solid rectangle) match their analytically-derived counterparts (black dotted lines) better than those estimated before acclimation (dashed rectangle). (C) A pairwise model generated from population dynamics before acclimation (top) predicts future reference dynamics less accurately than that generated after acclimation (bottom). (D) Quantification of the difference between pairwise and mechanistic models before (dashed) or after (solid) acclimation. All parameters are listed in Fig 3-SD1.  Fig 3-FS2. A canonical pairwise model, but not the alternative pairwise model, is suitable for a consumable mediator that accumulates without reaching a steady state within each dilution cycle. In a commensal community, the consumable mediator C1 accumulates as the consumer S2 gradually goes extinct. Pairwise model parameters are estimated from the mechanistic model dynamics in the training window (magenta, between 50 and 60 generations). The canonical model ( Fig 3A) shows dynamics (dotted) that match those of the mechanistic model (solid). As expected, the alternative pairwise model (dashed) fails. Thus, accumulating C1 can be regarded as a reusable mediator. Note that population fractions (instead of population densities) are plotted, which fluctuate less during dilutions compared to mediator concentration. All parameters are listed in Fig 3-SD2.  Fig 3-FS4. Community trajectory approaching the f-zero-isocline allows us to use the alternative pairwise model approximation. S1 releases a consumable metabolite C1 which stimulates S2 growth. In all panels, brown circles indicate the C1 and RS (=S2/S1) of a community, and are separated by 1/4 of community doubling time. In the vicinity of the f-zero-isocline (   We consider commensalism through a consumable mediator, where the producer (blue) and the consumer (green) eventually reach a steady state. We compare 1000-fold (A) and 10-fold (B) dilution steps to examine how fluctuations caused by dilutions affect parameter estimation. We use every 10 community doublings of reference dynamics as training windows to numerically estimate best-matching alternative pairwise model parameters. In (A), compared to (B), we see larger errors in estimating the interaction strength r21 compared to the true value (calculated from Fig 3B). In both (A) and (B), parameter estimations are less accurate if estimated before C1 has stabilized (rigorously speaking, before the community dynamics has approached the f-zero-isocline). Note that population fractions (instead of population densities) are plotted, which fluctuate less during dilutions compared to mediator concentration. All parameters are listed in Fig 3-SD6.   Fig 5C, if S3 does not remove the mediator of interaction between S1 and S2, a three-species pairwise model accurately matches the mechanistic model. Simulation parameters are provided in Fig 5-SD4. (C-D) As a control for Fig 5E, we ensured that fitness effects from multiple species are additive. In this case, a three-species pairwise model can represent the mechanistic model. To ensure the linearity of fitness effects, we have used a larger value of half saturation concentration ( 21 SC K = 10 9 instead of 10 5 in Fig 5E-F). We have adjusted the interaction coefficient accordingly such that the overall interaction strength exerted by S1 and S3 on S2 is comparable to that in Fig 5E-F (as evident by comparable population compositions). Since the interaction influences under these conditions remain in the linear range, the three-species pairwise model accurately predicts the reference dynamics. Simulation parameters are provided in Fig 5-SD5.