# Estimating the probability of coexistence in cross-feeding communities

###### Abstract

The dynamics of many microbial ecosystems are driven by cross-feeding interactions, in which metabolites excreted by some species are metabolised further by others. The population dynamics of such ecosystems are governed by frequency-dependent selection, which allows for stable coexistence of two or more species. We have analysed a model of cross-feeding based on the replicator equation, with the aim of establishing criteria for coexistence in ecosystems containing three species, given the information of the three species’ ability to coexist in their three separate pairs, i.e. the long term dynamics in the three two-species component systems. The triple-system is studied statistically and the probability of coexistence in the species triplet is computed for two models of species interactions. The interaction parameters are modelled either as stochastically independent or organised in a hierarchy where any derived metabolite carries less energy than previous nutrients in the metabolic chain. We differentiate between different modes of coexistence with respect to the pair-wise dynamics of the species, and find that the probability of coexistence is close to for triplet systems with three pair-wise coexistent pairs and for so-called intransitive systems. Systems with two and one pair-wise coexistent pairs are more likely to exist for random interaction parameters, but are on the other hand much less likely to exhibit triplet coexistence. Hence we conclude that certain species triplets are, from a statistical point of view, rare, but if allowed to interact are likely to coexist. This knowledge might be helpful when constructing synthetic microbial communities for industrial purposes.

###### keywords:

Cross-feeding, Population dynamics, Replicator equation, Permanence, Coexistence, Stability^{†}

^{†}journal: Journal of Theoretical Biology

## 1 Introduction

In recent years it has become increasingly clear that microbial species form complex communities, and rarely exist in isolation from each other (Zelezniak et al., 2015). A common form of interaction occurs via the exchange of nutrients that are released by one species and absorbed and further metabolised by other species in the community. This phenomenon is known as cross-feeding or synthropy and has been observed in a wide range of systems such as the human gut flora (Belenguer et al., 2006), the interactions of sulfate-reducers and methane oxidisers in the deep sea (Hallam et al., 2004; Pernthaler et al., 2008), the degradation of pesticides (Katsuyama et al., 2009), methanogenic environments (Stams, 1994), and in soil nitrification (Costa et al., 2006).

Note that cross-feeding can come in different degrees of complexity and interdependence. For example, when studying a system of E. coli strains feeding off an inflow of glucose, one strain is known to partially degrade the glucose to acetate, which would then be consumed by a second strain. Thus, the second strain will be affected by a negative frequency-dependent selection, as it needs the primary degrader. Furthermore, it has been put forward that the primary strain is dependent on the second one, as the secondary metabolite could be toxic at high concentrations, see for example (Pelz et al., 1999).

Since the growth rate of a species within the community depends on the metabolites produced by other species, it indirectly depends on the frequency of other species. This implies that systems where cross-feeding is dominant are governed by frequency-dependent selection, which allows for both dominance and coexistence depending on the strength and sign of the interactions between the species (Bomze, 1983).

Frequency-dependent selection (together with other mechanisms such as spatial structure (Kerr et al., 2002)) is the most likely explanation for the stability of natural microbial ecosystems, such as the gut microbiota. Although we have mapped out a large number of microbial communities, we still lack a definite understanding of their dynamics and cannot explain why they are stable (Harcombe et al., 2014). This lack of knowledge becomes evident when it comes to assembling or constructing artificial communities that are stably maintained (Jagmann and Philipp, 2014), and the need for understanding has become even more pressing now that the potential for engineering microbial communities for specific industrial purposes has been unravelled (Großkopf and Soyer, 2014).

In order to build efficient microbial communities we could like to know if a given collection of species can form a stable ecosystem in which no species are outcompeted and driven to extinction. This question can be approached either from a top-down perspective using flux-balance analysis and co-occurrence data (Zelezniak et al., 2015; van Hoek and Merks, 2016), or from a detailed understanding of the population dynamics of cross-feeding systems. We take the latter approach and address the following straightforward and concrete question: if we have qualitative information about the pair-wise dynamics of three species, what can be said about the likelihood of coexistence in the three-species community?

Thus, this paper strives to outline what configurations of species triplets that are likely to form coexistent populations. More specifically, we would like to know the coexistence properties of a triplet based on known pair-wise interactions between the constituent species. Ideally we would like to have quantitative information about the interactions of the different species, but for many practical purposes this is too much to ask. We therefore settle for qualitative information, and assume that for each pair of species we know if they coexist or if one species outcompetes the other. However, this poses a problem since it is known that systems that are identical on the pair-wise level (in the above qualitative sense), might behave differently when all three species are present (Bomze, 1983). This implies that triplet coexistence cannot be determined from the three pairs in isolation, but is a property of the interactions in the complete triplet. But all is not lost, since we may still be able to say something about the probability of coexistence. With this in mind we set out to study cross-feeding systems in which one, two or three pairs of species (out of the three) co-exist in isolation, and also intransitive systems where no dominant species exists (like rock-scissors-paper), and we do this from a statistical point of view to estimate the probability of coexistence in triplets of species.

### 1.1 Mathematical modelling of cross-feeding

The dynamics of large — consisting of some billion cells or more — and well-mixed populations of cross-feeding bacteria can be described by a system of coupled non-linear autonomous ordinary differential equations known as the replicator system of equations (Lundh and Gerlee, 2013). The assumption of well-mixedness allows us to disregard spatial effects in the system and in a large population, we may safely discard any stochastic individual interactions. The replicator system of equations has its origin in game theory (Hofbauer and Sigmund, 2002) where it describes an evolutionary game of strategies and players (Gokhale and Traulsen, 2010), which corresponds to species and steps in the metabolic process in the cross-feeding framework. This correspondence is due to the fact that each metabolic step considered introduces an a coupled interaction between the species that take part in the metabolic chain.

Replicator systems have been studied extensively (Bomze, 1983; Hofbauer and Sigmund, 2002; Gokhale and Traulsen, 2010; Lundh and Gerlee, 2013) in relation to game theory and population dynamics, and we apply this theory to the present problem in order to find how pair-wise dynamics influence triplet coexistence.

In the present setting, the fitness of a species is given by the amount of energy that the species can extract from the available nutrients excreted by another species. Here we interpret “energy” in a somewhat loose and abstract meaning. In this general setting the total energy is modelled as a sum over all possible metabolic interactions, represented as a series expansion of the fitness function. The model has been studied for two species by Lundh and Gerlee (2013), under the assumption that metabolites are only utilised by at most two species. In the same paper, the authors derived conditions for coexistence in a two-species population and the case of intransitive three-species populations.

The dynamics of a cross-feeding ecosystem need not be modelled in the game-theoretical framework of the replicator system described by Lundh and Gerlee (2013). Other ODE-systems have described cross-feeding as a direct interaction between the involved species (Bull and Harcombe, 2009; Estrela and Gudelj, 2010), by explicitly modelling nutrient uptake and mortality (Katsuyama et al., 2009), and by using adaptive dynamics (Doebeli, 2002). Agent-based models have also been used by Gerlee and Lundh (2010) and Crombach and Hogeweg (2009), whereas Pfeiffer and Bonhoeffer (2004) have studied the evolution of cross-feeding as a result of optimal ATP energy production in cells. In a recent study by Gedeon and Murphy (2015) cross-feeding in a chemostat environment was analysed. Under fairly general assumptions on the structure of the cross-feeding network they could show that there is a unique stable equilibrium that corresponds to the largest community of species that can be supported by the available resources, and that biomass production is maximised at this equilibrium point.

In this paper we take as a starting point the work of Lundh and Gerlee (2013) and derive conditions for coexistence. Although we are able to derive analytical expressions for the conditions of coexistence, the large number of parameters in the model makes it difficult to draw any direct conclusions. Instead we approach the problem from a statistical point of view and randomly generate a large number of three species systems. The interaction parameters that model the energy uptake of a species are modelled in two ways: either as independent random variables or according to a hierarchical model where energy gains further down in the metabolic chain are lower than energy gains from primary metabolites. For a given parameter model, interaction parameters are drawn randomly from certain probability distributions to estimate the likelihood of permanence from the coexistence criteria. Then, the relevant statistics of the sampled systems are computed and compared to the derived coexistence criteria.

## 2 Preliminaries

The replicator system of equations for a population of species with individual frequencies is defined as

(1) |

where denotes the derivative with respect to time of a species frequency , is the species fitness function and is the average fitness in the population. Intuitively, a species that is fitter than the population average will increase in proportion to its current frequency and a species less fit than the average will decrease correspondingly.

A fixed point of a system of ordinary differential equations is a point in the domain of such that . For the replicator system, the domain of definition for is the simplex

(2) |

due to the requirement that the species frequencies are positive and defined as fractions of the whole population. The stability of the fixed points is determined (Hofbauer and Sigmund, 2002) by the eigenvalues of the Jacobian

(3) |

Coexistence of species is related to the existence and location of the fixed points, but the system need not have stable fixed points in order to exhibit coexistence, and conversely the existence of stable fixed points does not imply coexistence. Rather, the property we are looking for is that of permanence, which is defined as:

###### Definition 1.

A replicator system (1) is considered permanent if for all initial states , we have that for all species and all .

The number of fixed points, their location and stability properties form the basis for the classification of solutions to the system (1). Whether a fixed point is located on the boundary or in the interior of the domain of is of special interest, since a permanent system is characterised by either a stable interior fixed point or a non-edge cyclic trajectory around a center fixed point. In the permanent case, we have for all species and thus that the interior fixed point must satisfy for all , which means that the fitness of all species are equal

(4) |

The dynamics of the general three-species replicator system is studied and outlined in Bomze (1983).The characterisation is based on the number of fixed points and their locations in the interior and boundary of the simplex. The fitness function used in that paper has the form

(5) |

where the elements of a general normal-form payoff matrix are

(6) |

If a payoff matrix is not given in this zero-diagonal form, it may be transformed as such since the dynamics of the replicator system (1) does not change under column-wise addition of constants to the replicator system (Bomze, 1983).

### 2.1 Replicator system for cross-feeding

For the cross-feeding model at hand, we recall here the derivation fitness function introduced by Lundh and Gerlee (2013), where it is assumed that the fitness of a species depends on its capacity to extract energy from a primary nutrient and from nutrients produced by other species (including itself). Three main assumptions were made in order to simplify the model: Firstly, it was assumed that all uptake are equally efficient and linear with respect to the medium concentration, implicitly assuming that metabolites are scarce so that no saturation effects are present. The second assumption was that the amount of energy that can be extracted after the original metabolite being digested twice can be ignored. Thirdly, it was assumed that there is a separation in time scale between the dynamics of the metabolites and the bacterial population dynamics.

For an illustration of the hierarchy of metabolites and energy extraction, see Figure 1.

Species will gain an amount of energy when degrading the primary nutrient and when degrading metabolites excreted by species . Since these parameters determine how the species interact with one another we term them interaction parameters. To specify the interactions between three species we need parameters. This number will later be revised, so that a payoff matrix on normal form only needs 9 interaction parameters. As by assumption two, any interactions of higher order than two are assumed negligible, so that the fitness of a species is approximated by the sum of two terms:

(7) |

where is a constant which describes the the effectiveness of uptake of the energy source. Thus the product in the model corresponds to the total uptake of the metabolite of bacteria from species , and is the flow rate of metabolites into and out of the system.

This means that the concentration of the primary metabolite can be described by the following differential equation:

(8) |

The concentration of a first order metabolite is then expressed as

(9) |

Since we are looking for a steady state, we let the right hand sides in equations (8) and (9) be equal to equal to zero due to the third assumption above. This leads to a steady state of in primary resource at the level of , from which we define

(10) |

Similarly, the steady-state level for is

(11) |

The steady state fitness of (7) can thus be expressed as

(12) |

Note that this equation is a simplification of (Lundh and Gerlee, 2013, Eq. (9)) with the higher order terms ignored.

In a three-species replicator system, there is the possibility of stable interior fixed points as well as stable and unstable fixed points anywhere on the boundary of the system. Bomze (1983, 1995) characterises no less than 49 different types of phase portraits for a three-species replicator system, with the main division being the number of fixed points in the interior of the simplex (Bomze, 1983)

(13) |

### 2.2 Models for interactions parameters

The interactions parameters , are considered as random variables and modelled in two distinct ways. In the first scheme, and are assumed to be independent in the stochastic sense, and drawn from a uniform distribution on the unit interval, i.e. and for all species and . The uniform distribution is chosen due to its simplicity and can be seen to represent no prior knowledge on the interaction parameters. Please note that the positivity of the interaction terms (as they represent an energy gain) restricts us from using a normal distribution, which otherwise would have been a natural choice (Gokhale and Traulsen, 2010).

In the second scheme we constrain the amount of energy extracted at higher levels of cross-feeding by assuming that it is necessarily smaller than the amount extracted from the primary resource. This implies that for a fixed species , we have

(14) |

for all species . We implement this by letting and defining the second order terms as

(15) |

where the scaling factor is drawn from a -distribution.

In order to test the generality of our results we also consider the two schemes with exponentially distributed parameters with intensity (having mean 1/2, the same as the -parameters).

## 3 Analytical results

We are now ready to discuss how the results of Bomze (1983) and Lundh and Gerlee (2013) can be used to determine the probability of permanence for triplet systems. First, we will outline which of the previous results that are of interest and how they relate to triplet coexistence, as well as define a necessary condition for permanence in Section 3.1. Then, in Sections 3.2–3.4 we will describe the types of pair-wise interactions in the ecosystem that are of interest and how to define each of them.

In order to compare the replicator systems proposed by Lundh and Gerlee (2013) to that of Bomze (1983), we may use a technique described by Gerstung et al. (2011) and Stadler and Schuster (1990), namely that we define an alternative payoff matrix with elements

(16) |

so that the fitness function (12) may be written as

(17) |

which is equivalent to the linear fitness function (5). The proof of equivalence is straightforward and relies on the fact that (Gerstung et al., 2011).

By the property mentioned in Section 2 that a payoff matrix may be transformed by column-wise addition and subtraction, we will use an alternate form of (16) to be able to directly use the results of Bomze (1983), namely:

(18) |

where we have the definition of total energy uptake

(19) |

### 3.1 Stable interior fixed points

We are now to derive general conditions for existence and stability of fixed points (FPs) in the interior of the state space , defined as (13). A stable interior fixed point is not a sufficient condition for permanence, as there are replicator systems with a stable fixed point in the interior of the simplex that cannot be reached from all initial states. These systems will be called conditionally permanent, and an example of trajectories in permanent vs. conditionally permanent systems is shown in Figure 2. In biological terms, conditional permanence means that a species triplet may be stable at certain initial frequencies but not at others. Dynamical systems sometimes exhibit stable trajectories in so-called limit cycles, where solutions tend towards a cyclic trajectory where all frequencies are non-zero. For a replicator system of less than four species, no such limit cycles are possible, as proven by Hofbauer (1981).

The method we use is due to Bomze (1983) and also used by Stadler and Schuster (1990) and is based on finding two coordinates that define a unique fixed point

(20) |

that lies in the interior of when are positive. The coordinates are given by

(21) | |||

(22) |

where are the co-factors of the payoff matrix (18), defined as

(23) | ||||

(24) | ||||

(25) |

where .

The criterion for existence of the fixed point (20) is that are real and positive, which occurs when all of the have the same sign,

(26) |

where denotes the sign function. The stability properties of the fixed point is determined by the Jacobian’s (3) eigenvalues

(27) |

where we have defined, for notational convenience,

(28) | ||||

(29) |

The system is stable when both eigenvalues are real and negative or when the real part of the complex eigenvalues are negative, i.e., when

(30) |

In conclusion the conditions for existence and stability of fixed points in the interior of the simplex are (i) that the co-factors have the same sign (26) and (ii) that the real part of the eigenvalues are negative (30) (see Table 1).

### 3.2 Centre fixed points and cyclic trajectories

We noted before that permanence might be the result of cyclic trajectories in the interior of the state space. It is therefore valuable to know for which systems this occurs. If the eigenvalues to the payoff matrix of the replicator system are strictly imaginary, i.e. with a zero real part, in a neighbourhood of an inner fixed point, then all trajectories in the neighbourhood will be cyclic and the fixed point is a center. An example of a center system is shown in Figure 3.

The center fixed point is the limiting case between complex eigenvalues with negative and positive real parts. However, from a probabilistic point of view, we may note that strictly imaginary eigenvalues require an exactly zero real part of (27). Since all terms depend on the random variables and , we have that equality constitutes a zero-measure set and hence we expect a zero probability of finding such systems from randomly generated parameters.

### 3.3 Pair-wise coexistence

We have now set out the general criteria for triplet coexistence and will now investigate how pair-wise behaviour affects coexistence. In a permanent three-species system, it may be the case that the involved species are pair-wise coexistent when isolated from the third species. In terms of dynamical systems, the criterion for pair-wise coexistence means that there exists one or more stable fixed points on the non-corner edges of the simplex. If the fixed point on a given edge is stable, any trajectories along that edge will converge to the fixed point as .

In the interior of the simplex near a stable boundary fixed point there are however two possible behaviours: a semi-stable fixed point will repel the trajectories in the part of its neighborhood that lies in the interior of the simplex, whereas a stable fixed point will attract any trajectory in its neighborhood.

We now analyse the situation where we have pair-wise coexistence of one or more of the three species pairs, with a method described by Stadler and Schuster (1990). These boundary fixed points must be semi-stable, and hence attract trajectories on the boundary and repel trajectories in the interior of the simplex. To analyse this situation we use the payoff matrix (16) on its normal form

(31) |

that corresponds the homogenous form of the replicator system. We find that the fixed points on the non-corner edge of the simplex are

(32) | |||

(33) | |||

(34) |

under the conditions

(35) | |||

(36) | |||

(37) |

which ensure that each pair of payoff elements, for example and of (32), have the same sign so that the coordinates of are properly defined on the intervals when normalised.

In order to have permanence, we require that the edge fixed points are semi-stable fixed points which attract trajectories on the edge and repel trajectories in the interior of the simplex, so-called saddle points. The previous conditions hold for both stable and unstable fixed points, and to single out the fixed points which are stable along the edges, we modify the conditions (35)-(37) to require positiveness for each element in the pairs of payoff elements

(38) | |||

(39) | |||

(40) |

so that edge-bound trajectories on both sides of the fixed point will tend to the fixed point (Stadler and Schuster, 1990). We will now discuss the cases where three, two and one fixed points exists on the boundary of the system. Biologically this corresponds to situations where three, two and one of the three pairs exhibit stable coexistence in isolation.

In the case of pair-wise coexistence of all three species we have existence of unique fixed points as of (32)-(34) on the edges of the simplex when the conditions (35)-(37) are fulfilled. Furthermore, the fixed points are semi-stable when there exists a stable interior fixed point and the conditions (38)-(40) hold (Table 1). We note that the criteria for semi-stability of the boundary fixed points imply existence of the interior fixed point. The conditions are collected in Table 2, where all conditions are necessary for a system with three pair-wise fixed points and one stable triplet fixed point in the interior of the state space.

Fixed point | Existence | Stability | |
---|---|---|---|

(32) | (35) | ||

(33) | (36) | ||

(34) | (37) | ||

(20) | (26) | (30) |

For a system with two pairs that are coexistent in isolation from the third species, we require that any two of the conditions of edge fixed points in Table 2 hold. This corresponds to the case where one of the edge fixed points have migrated to a corner of the simplex when payoff entries have different signs. Recall that a stable fixed point in the interior of the simplex is not in itself a sufficient condition for permanence, as there are initial states that do not converge to the fixed point.

### 3.4 Intransitivity and permanence

A special type of pair-wise relation, which is known to promote coexistence (Kerr et al., 2002), is that of an intransitive species triplet. In this case the species dominate each other in a circular fashion, just as in the game rock-scissors-paper. This can be expressed as

(41) |

where is the -th row, -th column element of the payoff matrix (16), and we consider all indices modulo . This condition follows from considering the ordered pairs such that the species outcompetes species when no other species are present. In terms of the phase portrait of the replicator system, this means that the fixed point is unstable and that the fixed point is stable.

The criterion for permanence for an intransitive three-species replicator system is

(42) |

as stated in Theorem 2 of Lundh and Gerlee (2013). The permanence factors are defined as

(43) |

and are positive as long as the system is permanent per condition (41). In conclusion, the three-species replicator system (1) is pair-wise intransitive and permanent if condition (42) holds together with (41). The conditions are collected in Table 3.

## 4 Numerical evaluation of permanence criteria

For a general three-species system, the derived criteria for existence and stability of fixed points in the 2-simplex are hard to interpret due to the high dimensionality of the parameters (9 interaction parameters in total). We therefore investigate them from a probabilistic point of view rather than study them analytically. The statistical properties of three-species systems are evaluated numerically for random interactions parameters drawn according to the Uni(0,1)-distribution in the independent and hierarchical parameter model (see section 2.2 for details). We also extend the analysis to the case where the interaction parameters are drawn from an exponential distribution (with mean ). The results were averaged across different realisations (independent draws of the interaction parameters) and we used parameter values and (Lundh and Gerlee, 2013). The empirical probabilities of finding a system with a certain property are estimated as the fraction of systems that satisfy the conditions (equalities and/or inequalities) that correspond to the property.

### 4.1 Stable interior fixed point

Existence of the interior fixed point is determined by the signs of the sub-determinants (23)-(25), which gives the criterion

(44) |

The necessary criterion for stability of the fixed point is negative real parts of the eigenvalues (27) to the interior fixed point, which is ensured by

(45) |

Table 4 collects the results from the independent and hierarchical model of the interaction parameters for uniformly distributed random variables. We note that for both models the probability of existence of a fixed point is on the order of 1% and the probability that this fixed point is also stable is roughly one order of magnitude smaller.

### 4.2 Pair-wise coexistence

For a pair-wise coexistent triplet, we require one, two or three semi-stable fixed points on the non-corner boundary of the simplex (in addition to a stable interior fixed point). The criteria for semi-stability that are found in Table 2 are sufficient also for existence of the boundary fixed points.

Pr( semi-stable FPs on boundary) | ||||
---|---|---|---|---|

model | Uni, Indep. | 12 | 6.0 | 2.2 |

Uni, Hier. | 2.2 | 1.2 | 1.4 |

We note that the probability of coexistence is decreasing with the number of required fixed points on the boundary. This is not surprising since more inequalities need to hold for two and three pairs.

### 4.3 Intransitivity and permanence

Existence of intransitive systems is determined by the conditions (41), that need to hold for all pairs where we as usual consider the indices modulo 3. The results are collected in Table 6, where we see that the uniform model is one order of magnitude more likely to be permanent compared to the tree hierarchy model. Intransitive systems contain a fixed point in the interior of the simplex which is either stable or unstable, as the probability of a cycle fixed point is a zero-measure set. Also, the conditions that determine stability of the fixed point are symmetric, which shows that the probability of permanence should be of the probability of intransitivity. The reported numerical results are in agreement with this theory.

### 4.4 Comparison of coexistent systems

We have so far investigated how likely the system is to exhibit certain properties on the level of pairs of species, given certain assumptions on the interaction terms. We now return to the original question, which was: If we have qualitative information about the pair-wise dynamics of three species, what can be said about the likelihood of coexistence in the three species community?

This question can be answered by looking at the probability of permanence conditioned on existence of the system. For example, if we know that among the three species two of them coexist in pairs, how likely is the system to exhibit coexistence of all three species? We investigated this for pair-wise coexistence and intransitive triplets and the results are collected in Table 7. First, we note that there is a pattern to the pair-wise coexistent systems in that the conditional probability of permanence is increasing with the number of coexistent pairs. This pattern is the same for both the independent and hierarchical model for the interactions parameters. We also note that the conditional probability of permanence in triplets with three coexistent pairs is by far larger than the probability for systems with two or one coexistent pairs. Finally, although intransitive systems are unlikely to exist when the interactions parameters are random (as shown in Table 6), close to half of the intransitive systems are permanent. This means that the intransitive property, and not the permanence, is the limiting factor. The converse is true for systems with one or two coexistent pairs on the boundary, as these systems are fairly likely to exist (Table 5) but have a low conditional probability of permanence.

In order to test the robustness of these results, we have calculated the corresponding probabilities with exponentially distributed parameters. This scenario describes a situation where any (finite) energy uptake is possible, but where the probability decreases exponentially with the amount of energy extracted. The rate parameter was set to 2 so as to ensure the same mean value of the random parameters as in the case with Uni(0,1)-distributed parameters. In Table 7, we see that the trend in the pair-wise coexistent systems is similar: the conditional probability of permanence is increasing with the number of fixed points on the boundary. Also, nearly half of the intransitive systems are permanent, and we note that the probability is highest for the tree hierarchy model.

## 5 Discussion

We have investigated the statistical properties of the replicator system (1) via the parameters , that represent the energy uptake of a species (Lundh and Gerlee, 2013). These parameters are modeled as random variables and we use two example distributions for the parameters. The Uni(0,1) distribution is chosen to investigate the properties of the system under the assumption that any energy uptake is equally likely when normalised onto the [0,1]-interval, and the Exp(2) distribution is chosen to test the robustness of the results. Other possible distributions are log-normal and the gamma distributions, where a certain interval for the extracted energy may be specified.

Two schemes are used for the random interactions parameters the independent where any species is allowed to extract any amount of energy from any metabolite and, in particular, may extract more energy from a derived metabolite than from the primary nutrient. This is the most general case of cross-feeding, where different species may be specialised on different nutrients. To investigate the case where a species is not able to extract more energy further down in the metabolic chain, the tree hierarchy model which satisfies the condition is used as an alternative. In general the results are similar across both the different distributions and the two models, which suggests a certain robustness. However one trend seems to be that for both distributions the tree hierarchy model exhibits a lower probability of permanence. The reasons for this is that in the tree hierarchy model the first-order interactions () are larger than the second-order interactions (), which tends to suppress coexistence. This can be understood intuitively by considering the limiting case where the second-order terms tend to zero. In that case coexistence is improbable, and the species with the largest first-order term dominates no matter how small the relative advantage compared to the other species.

We have shown that triplet coexistence is more or less likely to occur in the three-species systems, depending on the pair-wise interactions of the systems. We find that nearly half of the systems with three coexistent pairs and intransitive systems are permanent, where the corresponding numbers for systems with one and two coexistent pairs are closer to one percent and ten percent, respectively. This is an important point when considering the possibilities for engineering permanent ecosystems: if an intransitive or pair-wise coexistent triplet is found, then it is likely to also be permanent.

The connection to normal form games suggests that a similar statistical analysis could be readily applied to standard three strategy games. In that case one would have to decide upon a reasonable stochastic model for the elements of the payoff matrix. Here a normal distribution could be used since we no longer have any restriction on the signs of the payoffs. A natural extension would be to look at larger number of species (more than 3 strategies in a game theory context) and also include more levels of interaction in the fitness function (more than 2 players in the game). Such an analysis has been carried out by Gokhale and Traulsen (2010), and they could show that the maximum number of internal equilibria grows as , where is the number of players and is the number of strategies. However, numerical results showed that the fraction of games (generated at random) that allow for coexistence of all strategies/species rapidly approaches zero as and grow. Here our methodology could be applied to derive criteria that could aid in generating and finding games that exhibit coexistence and permanence.

One can compare the results in Table 4 with (Han et al., 2012, Theorems 1 and 2) where one can find that their random game with two players and three strategies has a probability of for the existence of an internal fixed point and a probability for a stable internal fixed point. So, at least for these cases, we see that the cross-feeding set up makes it less probable for both existence and stability of fixed points than for the random matrix set up in (Han et al., 2012). The explanation for this can be found in the generation of the cross-feeding random matrix, see Equation (31) which consists of elements, , that are not independent since and share one stochastic variable in one of the two terms of . This positive correlation makes it harder to find a system where the second order metabolites are compensating enough with the first order pay-off.

In a recent study by Huang et al. (2012) the emergence of polymorphism was studied in an evolutionary variant of the replicator equation. With a small mutation rate the payoff matrix was expanded with new random entries corresponding to the introduction of a new strategy/species into the ecosystem. When analysing coexisting triplets that emerged in this evolutionary process they found that almost all of them exhibited pair-wise coexistence, while only small fraction shwed Prisoner’s dilemma like dominance when considered as pairs. That result is in line with what one would expect from our analysis. Pair-wise coexistence might be an unlikely scenario, but when it occurs then it is highly likely (compared to the other cases) to yield permanence of the triplet.

The replicator system of equations can be shown to be equivalent to the Lotka–Volterra system of ODEs (Bomze, 1983). Studies of a such a Lotka–Volterra system with normally distributed random interaction parameters showed (Coyte et al., 2015) that cooperative interactions in a consortium are not likely to form stable permanence for a large number of species. In the model, the dependency between cooperating species that allowed cross-feeding or otherwise mutualistic groups to thrive also caused instability - when one of the cooperating species declined in numbers, others were likely to follow.

We hope that the results derived in this paper, and future extensions of it, will be useful for experimentalists trying assemble stable microbial communities. Our results could aid in the construction of artificial communities that perform well-defined industrial functions. With the current advances in synthetic biology our predictions could also be tested by using genetic cross-feeding switches that are tuneable to achieve interaction strengths within a given range (Kerner et al., 2012). By generating a large number of different mutants one could investigate if indeed three pair-wise coexistent species are most likely to show triplet coexistence, and if the fraction of such system is approximately 50%.

In conclusion, we think that this and similar bottom-up studies (Gedeon and Murphy, 2015), that strive for an understanding of basic mechanisms in cross-feeding, have a lot to offer when it comes to understanding complex microbial ecosystems.

## 6 Acknowledgements

The authors would like to thank the both anonymous reviewers for their careful reading, hard questions and valuable suggestions that significantly improved this paper. PG would like to acknowledge funding from the Swedish Foundation for Strategic Research (Grant no. AM13-0046) and Vetenskapsrådet (Grant no. 2014-6095)

## 7 References

## References

- Belenguer et al. (2006) Belenguer, A., et al., 2006. Two routes of metabolic cross-feeding between Bifidobacterium adolescentis and butyrate-producing anaerobes from the human gut. Applied and Environmental Microbiology 72, 3593–3599.
- Bomze (1983) Bomze, I. M., 1983. Lotka-volterra equation and replicator dynamics: A two-dimensional classification. Biological Cybernetics 48, 201–211.
- Bomze (1995) Bomze, I. M., 1995. Lotka-volterra equation and replicator dynamics: New issues in classification. Biological Cybernetics 72, 447–453.
- Bull and Harcombe (2009) Bull, J. J., Harcombe, W. R., 2009. Population dynamics constrain the cooperative evolution of cross-feeding. PLoS ONE 4 (e4115).
- Costa et al. (2006) Costa, E., Perez, J., Kreft, J.-U., 2006. Why is metabolic labour divided in nitrification? Trends in Microbiology 14, 213–219.
- Coyte et al. (2015) Coyte, K. Z., Schluter, J., Foster, K. R., Nov. 2015. The ecology of the microbiome: Networks, competition and stability. Science 350 (6261), 663–666.
- Crombach and Hogeweg (2009) Crombach, A., Hogeweg, P., 2009. Evolution of resource cycling in ecosystems and individuals. BMC Evolutionary Biology 9 (1), 122.
- Doebeli (2002) Doebeli, M., 2002. A model for the evolutionary dynamics of cross-feeding polymorphisms in microorganisms. Population Ecology 44, 59–70.
- Estrela and Gudelj (2010) Estrela, S., Gudelj, I., 2010. Evolution of cooperative cross-feeding could be less challenging than originally thought. PLoS ONE 5 (e14121).
- Gedeon and Murphy (2015) Gedeon, T., Murphy, P., Sep. 2015. Dynamics of Simple Food Webs. Bulletin of Mathematical Biology, 1–21.
- Gerlee and Lundh (2010) Gerlee, P., Lundh, T., 2010. Productivity and diversity in a cross-feeding population of artificial organisms. Evolution 64, 2716–2730.
- Gerstung et al. (2011) Gerstung, M., Nakhoul, H., Beerenwinkel, N., 2011. Evolutionary games with affine fitness functions: Applications to cancer. Dynamical Games and Applications 1, 370–385.
- Gokhale and Traulsen (2010) Gokhale, C., Traulsen, A., 2010. Evolutionary games in the multiverse. Proceedings of the National Academy of Sciences 107, 5500–5504.
- Großkopf and Soyer (2014) Großkopf, T., Soyer, O. S., Apr. 2014. Synthetic microbial communities. Current Opinion in Microbiology 18, 72–77.
- Hallam et al. (2004) Hallam, S. J., et al., 2004. Reverse methanogenesis: testing the hypothesis with environmental genomics. Science 305, 1457–1462.
- Han et al. (2012) Han, T. A., Traulsen, A., Gokhale, C. S., 2012. On equilibrium properties of evolutionary multi-player games with random payoff matrices. Theoretical Population Biology 81 (4), 264 – 272.
- Harcombe et al. (2014) Harcombe, W. R., Riehl, W. J., Dukovski, I., Granger, B. R., Betts, A., Lang, A. H., Bonilla, G., Kar, A., Leiby, N., Mehta, P., Marx, C. J., Segrè, D., May 2014. Metabolic Resource Allocation in Individual Microbes Determines Ecosystem Interactions and Spatial Dynamics. Cell Reports 7 (4), 1104–1115.
- Hofbauer (1981) Hofbauer, J., 1981. On the occurrence of limit cycles in the volterra–lotka equation. Nonlinear Analysis, Theory, Methods and Applications 5, 1003–1007.
- Hofbauer and Sigmund (2002) Hofbauer, J., Sigmund, K., 2002. Evolutionary Games and Population Dynamics. Cambridge University Press, Cambridge.
- Huang et al. (2012) Huang, W., Haubold, B., Hauert, C., Traulsen, A., Jun. 2012. Emergence of stable polymorphisms driven by evolutionary games between mutants. Nature Communications 3, 919–7.
- Jagmann and Philipp (2014) Jagmann, N., Philipp, B., Aug. 2014. Design of synthetic microbial communities for biotechnological production processes. Journal of Biotechnology 184, 209–218.
- Katsuyama et al. (2009) Katsuyama, C., et al., 2009. Complementary cooperation between two syntrophic bacteria in pesticide degradation. Journal of Theoretical Biology 256, 644–654.
- Kerner et al. (2012) Kerner, A., Park, J., Williams, A., Lin, X. N., 2012. A programmable Escherichia coli consortium via tunable symbiosis. PLoS ONE 7 (3), e34032.
- Kerr et al. (2002) Kerr, B., Riley, M. A., Feldman, M. W., Bohannan, B. J. M., Jul. 2002. Local dispersal promotes biodiversity in a real-life game of rock–paper–scissors. Nature 418 (6894), 171–174.
- Lundh and Gerlee (2013) Lundh, T., Gerlee, P., 2013. Cross-feeding dynamics described by a series expansion of the replicator equation. Bulletin of Mathematical Biology 75, 709–724.
- Pelz et al. (1999) Pelz, O., Tesar, M., Wittich, R., Moore, E. R. B., Timmis, K. N., Abraham, W. R., 1999. Towards elucidation of microbial community metabolic pathways: unravelling the network of carbon sharing in a pollutant-degrading bacterial consortium by immunocapture and isotopic ratio mass spectrometry. Environmental Microbiology 1 (2), 167–174.
- Pernthaler et al. (2008) Pernthaler, A., et al., 2008. Diverse syntrophic parterships from deep-sea methane vents revealed by direct cell capture and metagenomics. Proceedings of the National Academy of Sciences USA 105, 7052–7057.
- Pfeiffer and Bonhoeffer (2004) Pfeiffer, T., Bonhoeffer, S., 2004. Evolution of cross-feeding in microbial populations. American Society of Naturalists 163, E126–E135.
- Stadler and Schuster (1990) Stadler, P. F., Schuster, P., 1990. Dynamics of small autocatalytic networks – 1. bifurcations, permanence and exclusion. Bulletin of Mathematical Biology 52, 485–508.
- Stams (1994) Stams, A. J., 1994. Metabolic interactions between anaerobic bacteria in methanogenic environments. Antonie Van Leeuwenhoek 66, 271–294.
- van Hoek and Merks (2016) van Hoek, M. J., Merks, R. M., June 2016. Emergence of microbial diversity due to cross-feeding interactions in a spatial model of gut microbial metabolism, bioRxiv pre-print.
- Zelezniak et al. (2015) Zelezniak, A., Andrejev, S., Ponomarova, O., Mende, D. R., Bork, P., Patil, K. R., May 2015. Metabolic dependencies drive species co-occurrence in diverse microbial communities. Proceedings of the National Academy of Sciences of the United States of America 112 (20), 6449–6454.