MIIA predicts pairwise interactions in multi-species communities through the following two steps: (1) estimation of interaction coefficients in binary cultures, and (2) prediction of the shifts in interactions by additional members based on a minimal adjustment hypothesis.

Species that have strong dependencies on partners may not grow independently. Axenic population densities in this case can be unobtainable, which limits the originally described MIIA approach. In Equation (1), the values of ai,jB may explode for a very low axenic density of species i (i.e., xiA≪1) and also for a low population of its binary partner, i.e., xjB≪1; becomes unidentifiable if the values of xiA and xjB are unavailable. In this section, we provide an idea of handling the case where the value of xiA is extremely low or unavailable; in the next section, we will consider the case where the value of xjB is extremely low or unavailable.

Without compromising generality, we reformulated the model equations for binary and complex communities by multiplying xiA on both sides of Equations (1) and (2), which results in

where bi,jB and bi,jC are interaction coefficients scaled by xiA, i.e.,

The scaling above is always possible for any non-zero value of xiA. Thus, it implies that we translate ‘no growth’ as a minor presence – i.e., below the limit of detection – in place of absolute absence. The minimal adjustment rule of MIIA still applies to Equations (3) and (4) because the orthogonal relationship between binary and complex interaction coefficients established in the original space, remains valid on the new coordinates rescaled with a constant factor xiA (Figure 1B). We termed this extension scaled MIIA (s-MIIA).

The quantitative values of interaction coefficients predicted by the s-MIIA are different from those by the original MIIA. That is, what is determined with the s-MIIA is ‘scaled’ interaction coefficients (i.e., bi,jB and bi,jC), rather than ‘absolute’ values (i.e., ai,jB and ai,jC). However, it should be noted that the resulting scaled coefficients can also provide sufficient information required to predict how interaction changes in response to the addition of new neighbors, a primary question the MIIA aims to address. Thus, when the scope of prediction is confined to ‘relative’ interaction changes, both the original version and the s-MIIA generate the same result, i.e.,

or

where Δ⁢ai,jr⁢e⁢l and Δ⁢bi,jr⁢e⁢l denote the relative changes in interaction coefficients a_i,j and b_i,j that are predicted by the original and scaled MIIA, respectively. The superscripts B→C and C→B denote the interaction changes from binary to complex cultures and from complex to binary cultures, respectively. From Equations (7) and (8), therefore, it is clear that the predictions by the original formulation and the s-MIIA are exactly the same with respect to the relative changes in interaction coefficients.

While the scaling method above resolves the issue associated with extremely low or non-measurable axenic populations, identifying binary interaction coefficients is also dependent on the availability of binary co-culture data. We considered eight different growth scenarios in axenic and binary cultures and showed when binary interaction coefficients are identifiable by the original and scaled MIIA and when not (Table 1). For Cases I and II, binary interaction coefficients can be estimated either by the original formulation or the scaling method. Cases V and VI highlight the situations that could not be handled by the original formulation, but only by the scaling method. The remaining cases (i.e., Cases III, IV, VII, and VIII) are challenging scenarios, to which neither the original nor scaling method can be naively applied. However, we exclude Cases III and VIII from our consideration because these events might be rare if at all possible; therefore, our focus for demonstration purposes is placed on Cases IV and VII.

In order to estimate binary interaction coefficients for Cases IV and VII, we take bi,jB as an adjustable model parameter with an aim to determining its rational value or range based on the sensitivity analysis and the comparison with experimental evidence. The success of this parameterization strategy depends on several factors: (1) what binary interaction coefficients are unidentifiable, (2) what the complex community model [i.e., Equation (2)] looks like, and (3) the availability of experimental evidence to determine ranges of parameters. Figures 1C,D illustrate how the first and second factors affect predictions by considering the following two scenarios: (1) b1,3B were determined but not b1,2B, and (2) b1,2B was determined but not b1,3B. In the respective case, b1,3B and b1,2B are chosen as adjustable parameters. In the hypothetical scenario in Figures 1C,D, the prediction of interaction coefficients in the ternary community b1,2T and b1,3T is relatively less sensitive to the variation of the parameter b1,2B in comparison to b1,3B. This distinct sensitivity with respect to b1,2B and b1,3B of course depends on the slope of the linear ternary community model (l^T). As such, the sensitivity analysis evaluates the robustness of predictions and specifically shows which predictions are sensitive and which need additional data to reduce uncertainty.

Estimation of binary interaction coefficients in MIIA (both original and new formulations) (i.e., ai,jB’s and bi,jB’s) include the ratio of variables as shown in Equations (1) and (3). A caution is needed in calculating them from multiplicate data because the mean value of ratio can be different from the ratio of means. This so-called bias issue is fundamentally associated with the small sample size (Cochran, 1977) and often arises in biological experiments where the sample size is typically limited to 3 to 5. To evaluate the impact of the variability in small-size samples on the ratio estimation for the datasets used in our case studies, we first accounted for all possible combinations of three replicates of each variable, e.g., 27 combinations from three variables (xiA, xiB, and xjB) and then compared estimated values of ai,jB’s and bi,jB between the following two cases: (1) taking the mean of ratio and (2) taking the ratio of means. For the Wang et al.’ (2017) data used in the first case study (that shows relatively more significant variation than the data of the second case study), we observed only negligible differences between the two cases, i.e., <1.4% for bi,jB; <2.8% for ai,jB.

The study published from Wang et al. (2017) provides population data for three soil organisms that are known to grow in all combinations of axenic, binary and ternary cultures. The three bacterial species analyzed in their work include: Leuconostoc lactis (LL), Janthinobacterium lividum (JL), and Lactococcus piscium (LP). The copy number of the 16S rRNA genes – as determined by quantitative PCR – has been used as a proxy for the density of each species (see Supplementary Table S1 for the raw data retrieved from the original paper). We used this consortium data as a technical proof-of-concept example to illustrate how the scaled MIIA works and how to interpret its predictions.

The original MIIA predicted absolute interaction coefficients in binary and ternary communities (ai,jB and ai,jT) (Figure 2A), while the s-MIIA estimated their scaled values (ai,jB⁢xiA and ai,jT⁢xiA, i.e., bi,jB and bi,jT) (Figure 2B). Both methods predicted sign changes for some cases: interaction coefficients a_LL,JL and a_LL,LP were positive (i.e., promotive) in binary cultures, but became negative (i.e., inhibitory) in the ternary culture. This prediction was consistent with the experimental observations in the original paper by Wang et al. (2017).

It should be noted that the s-MIIA still allows us to assess the changes in interactions in the ternary community because the same constant (i.e., xiA) was multiplied on both ai,jB and ai,jT. As a critical difference from the original method, the prediction of the s-MIIA is limited to relative comparison of interaction parameters only for the same organism that is influenced by others. That is, we can compare b_i,j values across different j’s along each row in the matrix of interaction coefficients (Figure 2B), but not across different i’s because the values of the scaling constant xiA are different among species. As expected, predicted relative changes of interactions coefficients – e.g., defined in Equation (7) – were identical between the original and scaled methods (Figure 2C). The plus and minus signs in the table of Figure 2C denote the change of interactions in the ternary community in positive and negative directions relative to interactions in binary communities. In Figure 2D, we illustrated how the relationships between two species can be changed in binary and ternary communities. Wang et al.’s (2017) data showed reasonably small standard deviations of population densities for both axenic and binary cultures (from 0.006 to 0.069) and the ternary culture (from 0.022 to 0.094) where standard deviations were calculated using a multiplicative lognormal noise function as described in Supplementary Table S1 as well as in the previous paper by Song et al. (2019). The measurement error in this range did not deteriorate the predictive power of s-MIIA (Supplementary Figure S1; see also Song et al., 2019).

We extend our analysis to a more challenging case where the full estimation of binary interaction coefficients by the original MIIA is not possible due to the ineffective growth of some of the member species in axenic and binary cultures. These datasets were obtained from defined mixed cultures using four bacterial strains previously studied by Kato et al. (2008): a cellulose-degrading anaerobe (Clostridium straminisolvens CSK1), a saccharide-utilizing anaerobe (Clostridium thermosuccinogenes FG4), a peptide- and acetate-utilizing aerobe (Pseudoxanthomonas taiwanensis M1-3) and a peptide-, glucose-, and ethanol-utilizing aerobe (Brevibacillus agri M1-5), which are hereafter denoted by CS, CT, PT and BA, respectively. A number of different types of interactions such as trophic interactions, competition, and lethal inhibition as well as growth promotion have been detected in these bacterial communities (Kato et al., 2008; Yamamoto et al., 2010). When cultured together, these four organisms formed a stable consortium. With simpler memberships that do not contain all four members, however, some of the organisms could not grow as mentioned above. For instance, the growth of anaerobe CS requires the presence of aerobes (such as PT or BA) for the removal of oxygen; the anaerobe CT not only depends on aerobes, but also on CS that can degrade cellulose into saccharides. We summarized all interaction features below and provided the raw data retrieved from Kato et al. (2008) in Supplementary Table S2:

Strong growth dependency of CS and CT on other member species led to the absence of population data derived from axenic cultures and certain binary parings. Consequently, the original MIIA identified only 33% (4 out of 12) of the binary interaction coefficients (Supplementary Table S3). With this partial identification of interaction coefficients in binary communities, the original MIIA could not provide any predictions for ternary and quaternary communities except the CS-PT-BA consortium, the estimation of which was limited to 33% of interactions. The s-MIIA generated improved results by estimating 66% of the binary interaction (8 out of 12) and predicting 100% of interaction coefficients in the CS-PT-BA consortium. Prediction for other multi-species communities was limited however: only 16.7% in the CT-PT-BA consortium and no predictions for all other communities (Supplementary Table S3). This result shows that despite improvement, the scaling method alone is not sufficient to handle the cases where organisms cannot co-grow in binary communities (i.e., Cases IV and VII in Table 1), therefore requiring additional analyses to overcome this limitation.

While the scaling method provided expanded estimates of binary interaction coefficients in comparison to the original method, the datasets from Kato et al. (2008) still contain several binary interaction coefficients that remain unknown. These unknown binary coefficients (bC⁢S,C⁢TB, bC⁢T,C⁢SB, bP⁢T,C⁢TB, and bB⁢A,C⁢TB) are all associated with CT, which shows no growth in the axenic or even binary cultures. Our strategy to overcome this limitation is to take these four coefficients as adjustable parameters to examine how MIIA predictions would vary as their functions (see section Materials and Methods) and to determine their most plausible values or ranges based on the consistency between final model predictions with any experimental observations.

No direct experimental evidences were available from Kato et al. (2008) that can be used to decisively determine specific values of the four unknown binary parameters. However, they provided an interspecies interaction network for the four-member consortium that could be used as a basis to inform some of the unidentifiable binary interaction coefficients. This interaction network was derived from a combined analysis of three complementary datasets including those reported elsewhere (Kato et al., 2004, 2005; Yamamoto et al., 2010). To summarize their integrative analysis of experimental data: (1) they analyzed substrate utilization profiles and metabolites of each member under axenic culture conditions to propose exchange scenarios among the members; (2) they examined how the growth of a member in pure cultures can be promoted or suppressed when a cell-free culture filtrate of another member was added to the growth media; (3) they compared populations of a member (e.g., CS) between in the presence and the absence of another member (e.g., PT) to examine the effects of one on others (in the above example, effects of PT on the population of CS).

An overall interspecies interaction network among the four species depicted by Kato et al. (2008) show various types of interspecies interactions among member species (Figure 3A). For demonstration purpose, we took this network as a reliable specific interaction scenario, while it may not represent true interactions in the quaternary community. Figure 3B is the prediction of the s-MIIA obtained based on default parameters, bC⁢S,C⁢TB=bC⁢T,C⁢SB=bP⁢T,C⁢TB=bB⁢A,C⁢TB=0, i.e., by assuming neutral interactions for unknown interactions in binary communities. Surprisingly, this default prediction showed no contradiction with the interaction network experimentally derived by Kato et al. with respect to the positive and negative effects among four species although providing two additional interactions (dotted lines in Figure 3B) that were not seen in Figure 3A. As to the magnitude of interactions, however, there were some discrepancies between the two networks. In particular, Kato et al. identified two stronger interactions, i.e., the promotive effect of CS on CT and inhibitive effect of PT on BA (Figure 3A), the former of which was however not captured in the default network (Figure 3B). As addressed in the previous section, the comparison of interactions predicted by the s-MIIA is relative and therefore limited to the effects of different species for the same species that is influenced. With this limitation in mind, we tried to adjust unknown binary interaction parameters so that interaction coefficients in the quaternary community are consistent with the experimentally derived network. In this case, we were able to obtain a stronger effect of CS on CT in the quaternary network by increasing the value of bC⁢T,C⁢SB to 2.5 based on the parameter sensitivity analysis. The resulting interaction network and sensitivity profiles were shown in Figures 3C,D.

The sensitivity analysis was extended to all parameters (Figure 4). Interaction coefficients in ternary and quaternary cultures that show weak dependency on binary interaction parameters indicate that these insensitive coefficients can be robustly predicted regardless of assumed binary parameter values. Similarly, this implies that the process for determining interaction coefficients is sensitive to the assumed binary parameters (including those that may change their signs) and that additional experimental understanding on interactions in ternary and/or quaternary cultures can be required.

While an appropriate value for bC⁢T,C⁢SB was determined by integrating the parameter sensitivity analysis result and the experimental understanding from Kato et al., we lack further data or experimental observations that can be used to determine other parameters, i.e., bC⁢S,C⁢TB, bP⁢T,C⁢TB, and bB⁢A,C⁢TB. In order to examine the effects of the choice of these parameters, we considered the following three cases:

As mentioned previously, Case 1 denotes default setting that assumes neutral interactions; Case 2 is a simple adjustment of a minimum number of parameters (i.e., bC⁢T,C⁢SB) in order to match with experimental understanding. Case 3 represents an example of alternative parameter setting that leads to the same quaternary network as in Case 2, but different predictions in ternary networks.

We provided graphical representations of binary interactions for the three cases described above (Figure 5A). These three scenarios equally predict the structure of the interaction network in the quaternary community experimentally determined by Kato et al. (Figure 5B). Despite variations of assumed binary coefficients among three cases, predicted interactions in ternary communities were fairly robust. Particularly, Cases 1 and 2 generated exactly the same ternary interaction networks (left panel of Figure 5C). Overall, we found that out of 24 possible interactions in ternary communities, 75% of interactions were consistently predicted across three different parameter settings, while 25% predictions were not. When binary and quaternary communities were included, the portion of robust predictions increased to 79% (i.e., 38 out of the 48 interactions in total of binary, ternary, and quaternary communities). Therefore, these results (including both case-independent and case-specific predictions) can be used to inform the design of new experiments for further validation.

We finally examined how interspecies interactions could be modulated by the addition of new members or the loss of existing members. Interaction networks predicted across binary, ternary and quaternary cultures (Figure 5) showed that interspecies interactions can be membership dependent, while specific predictions were case-dependent. For example, predictions for Cases 1 and 2 showed that CS and CT, who could not coexist in the binary culture (Figure 5A), exhibit mutualism in the presence of PT, but show antagonism in the presence of BA (the left panel of Figure 5C) or both PT and BA (Figure 5B). By contrast, Case 3 predicted no such shifts, i.e., the relationship of CS and CT remained antagonistic regardless of who are their neighbors (the right panels of Figures 5A–C). Similar results were predicted for the influence of CT on BA, i.e., Cases 1 and 2 showed strong neighbor dependence of their interaction, but Case 3 did not. Interestingly, the influence of CT on PT was predicted to be membership dependent for all cases: the effect of CT on PT was negative in the quaternary culture, but the same in binary and ternary cultures was neutral in Cases 1 and 2, while positive in Case 3. We quantified modulation of interactions across binary and complex communities for Cases 1 to 3 (Figure 6). In all cases, we were able to predict considerable shifts in interspecies interactions across binary, ternary and quaternary communities.