key: cord-0429001-gakq6my0
authors: Jimenez, Paul; Scheuring, István
title: Density-dependent private benefit leads to bacterial mutualism
date: 2020-07-29
journal: bioRxiv
DOI: 10.1101/2020.07.28.224550
sha: 24a5faf780f0f1d25b06f708f0f0b4401a9b3244
doc_id: 429001
cord_uid: gakq6my0

Microorganisms produce materials leaked from the cell which are beneficial for themselves and their neighbors. We modeled the situation when cells can produce different costly secretions which increase the carrying capacity of the population. Strains that lose the function of producing one or more secretions avoid the cost of production and can exhaust the producers. However, secreting substances provides a private benefit for the producers in a density-dependent way. We developed a model to examine the outcome of the selection among different type of producer strains from the non-producer strain to the partial producers, to the full producer one. We were particularly interested in circumstances under which selection maintains partners that produce complementary secreted materials thus forming an interdependent mutualistic interaction. We show that interdependent mutualism is selected under broad range of conditions if private benefit decreases with density. Selection frequently causes the coexistence of more and less generalist cooperative strains, thus cooperation and exploitation co-occur. Interdependent mutual-ism is evolved under more specific circumstances if private benefit increases with density and these general observations are valid in a well-mixed and in a structured deme model. We show that the applied population structure supports cooperation in general, which, depending on the level of private benefit and intensity of mixing helps either the specialist or the generalist cooperators.

Microorganisms produce materials leaked from the cell which are beneficial for themselves and their neighbors. We modeled the situation when cells can produce different costly secretions which increase the carrying capacity of the population. Strains that lose the function of producing one or more secretions avoid the cost of production and can exhaust the producers. However, secreting substances provides a private benefit for the producers in a density-dependent way. We developed a model to examine the outcome of the selection among different type of producer strains from the non-producer strain to the partial producers, to the full producer one. We were particularly interested in circumstances under which selection maintains partners that produce complementary secreted materials thus forming an interdependent mutualistic interaction.

We show that interdependent mutualism is selected under broad range of conditions if private benefit decreases with density. Selection frequently causes the coexistence of more and less generalist cooperative strains, thus cooperation and exploitation co-occur. Interdependent mutualism is evolved under more specific circumstances if private benefit increases with density and these general observations are valid in a well-mixed and in a structured deme model. We show that the applied population structure sup-3 (Harcombe, 2010; Rivière et al.; Smith et al., 2019a) , which can explain unculturability of most of the microbial species (Kaeberlein et al., 2002; Stewart, 2012) . Some studies suggest that competitive interactions dominate over cooperative ones in microbial communities, i.e. mutualism is relatively rare in natural microbiota (Foster and Bell, 2012; Coyte et al., 2015) . Here, we reveal circumstances under which costly interdependent mutualism is stable from ecological and evolutionary points of view in microbial communities.

Previously, Oliviera et al. (Oliveira et al., 2014) studied the evolutionary stability of interdependent mutualism in a model of microbial communities. In this paper, we revisit and extend their model, thus we present their assumptions and results in more detail. They considered a situation where cells can produce different leaky materials that are essential for a higher carrying capacity. Variants may range from non-producers (complete cheaters) to generalist strains producing all types of secretions. Oliveira et al. assumed that numerous local populations are founded initially by some cells. Populations grow within local habitats following a dynamics depending on the proportions of different producer strains. After the growth phase, cells are mixed together and establish new local populations which then start to grow. This life cycle continues iteratively until a stationary state is reached. They showed that specialist interdependent producer strains are selected only if the cost of secretion is high and a medium number of individuals found a new habitat (genotype mixing or heterogeneity is at intermediate level). A low level of heterogeneity favors generalist cooperators, while high levels of heterogeneity prefers complete cheaters. Additionally, the coexistence of complete cheaters and generalist cooperators was possible under at very specific circumstances in their model.

Since the cost of producing leaky material is not high and a low level of heterogeneity can be found only in specific natural microbiota, the authors concluded that cooperation among generalists should be more common in nature than interdependent mutualism among specialists.

The dynamical model used by Oliveira et al. contains only one crucial parameter, the cost of producing secretions.

This simplifies analysis and the interpretation of results, but naturally neglects important characteristics of the problem being examined. One of the most important simplifications is that the model does not consider the private benefit of the producer. Moreover, it neglects that this private benefit is density dependent.

Neglecting the producer's private benefit assumes that the leaky materials are public goods available for all individuals equally, which is not necessarily true. In reality, there is inhomogeneity in the concentration of the secreted substances. Producer cells themselves and other cells in the producer's neighborhood typically detect higher concentrations of secretions than cells farther away from the producers (Gore et al., 2009; Morris et al., 2014; Lindsay et al., 2018; Smith et al., 2019b) .

Depending on the cell's mobility and the diffusion rate of secretions, either high or low density is more beneficial for producer strains as compared to non-producer ones. If diffusion of the secretions is faster than the cells movement, a concentration gradient is present around the producers. Therefore, increasing the density increases the proximity of cheaters to the cooperators which helps cheaters to exhaust cooperators intensively (Ross-Gillespie et al., 2009; Morris et al., 2014) . On the other hand, diffusion of the secretions can sometimes be slow compared to the cell's movement. Then, increasing density enhances the amount of public goods that can be captured by the cooperators themselves (Dobay et al., 2014) . Lindsay and his colleagues recently proved experimentally that both cases are realized (Lindsay et al., 2018) . Thus, we extended the model used by Oliveira et al. (2014) to include density-dependent private benefit of cooperators. Besides this extension, we also modified the basic dynamical system to describe the competition among the strains more adequately.

We show that density-dependent private benefit of cooperation significantly widens those circumstances where selection prefers specialist cooperators over generalist cooperators and complete cheaters. Further, we show that specialist cooperators can typically coexist with generalists or complete cheaters, if the private benefit of cooperation decreases with density.

In the following, we introduce the dynamical model systems from the simplest model with one secretion to those involving any number of secretions, and then present the structured population model. The first part of the Results section contains the invasion analysis and numerical studies in the local populations. Then we investigate the dynamics of the structured model numerically. The Discussion interprets and summarizes the main results.

First, we assumed that only one leaky material can be produced. There are two types of cells; type [0] which does not produce the leaky material and type [1] which does. The density of cheater (type [0]) and producer (type [1]) cells is denoted by g 0 and g 1 , respectively. Cheaters replicate with a replication rate r and follow a logistic growth with carrying capacity K base . Thus the dynamics of cheaters without the presence of cooperators is

With the help of the secreted materials, cells can utilize new sources of nutrients,consequently, we assume that these leaky materials increase the carrying capacity of the population. Since the efficiency of the secreted material depends on its concentration, the cooperative term of the carrying capacity K coop depends on the frequency of producers. We use a simple linear relationship between the frequency of producers and carrying capacity K coop (Oliveira et al., 2014) given by

where K max is a constant that describes the maximum carrying capacity which can be reached if all cells produce the secreted material. Thus, the presence of producers modifies the dynamics of the cheaters to

Producing leaky materials is costly, which means that the replication rate decreases or the death rate increases for the producer strain. Let r 1 = r − c denote the net growth rate of producer cells, where c (r > c > 0) is the cost of producing the leaky material. We assume that producer and non-producer strains differ only in the production of the leaky material, that is, all cells are genetic variants of closely related strains and thus, they compete for the same limiting resources with the same efficiency. Consequently, the strength of competition among the individuals of both strains should be the same. This term is r /(K base + K coop ) for cheaters (see eq. 3), which should be the same for cooperator strains. In this way, the dynamical model of the cooperator-cheater system will be the following 5 We note here that Oliveira et al. (2014) neglected the r 1 /r term in eq. (5) which cannot be supported by biological arguments and results in a structurally unstable model (for more details, see Discussion).

However, microbial populations are rarely well-mixed. Cells move in a limited and correlated manner even in liquid environments (Kümmerli et al., 2014; Sretenovic et al., 2017) and leaky materials diffuse with a limited speed or a fraction of them can bind to the surface of producers. A consequence of the nonrandom distribution of cells and the secreted material is that the relative fitness of cheaters and cooperators is density-dependent. Evidence that the relative fitness of cooperators can increase or decrease with density (see the Introduction) is supported by theoretical and experimental results. In this way, taking this effect into consideration we modify the availability of the secreted material in the function of density and trait-type in eqs. (4, 5). We introduce the ∆(g ) function which describes the density-dependent private benefit of producers, where g = g 0 + g 1 is the total density of the population. The cooperators' carrying capacity increases with ∆(g ), while the cheaters' carrying capacity decreases with this term.

Thus we obtain the following dynamics for the cheaters and cooperators.

In line with experimental observations (Lindsay et al., 2018) we used a saturating function for ∆(g ) if the private benefit of cooperators increases with density

and a similar decreasing function if the benefit of cooperators decreases with density

In both cases δ > 0 refers to the maximum intensity of the private benefit and K ∆ determines the density where this effect reaches half of the maximum value (∆(K ∆ ) = δ/2) ( Fig. 1 ).

We consider the case when microbes can extract two different materials which determines four different genotypes:

the cheater, producing neither material (denoted by [0,0]), the two types of specialist cooperators, producing one of the materials (denoted by [1,0] and [0,1]), and the generalist cooperator, producing both materials ([1,1]). With the same assumptions applied above, the dynamics of traits are described by

The private benefit of cooperators as a function of the population density. Private benefit either increases (blue) or decreases (green) with density. The functions are characterized by parameters K ∆ and δ.

where g i j represents the density and r i j = r − [ (i + j ) ] ν c is the replication rate of genotype [i , j ] and g = g i j = g 00 + g 10 + g 01 + g 11 equals to the total density.

We introduce a new parameter ν > 0 to characterize the cost function of producing more different secretions. ν refers to accelerating cost with the number of produced materials. ν = 1 is selected in most cases, that is the cost is simply additive, but we also studied the effect of decelerating and accelerating costs as well. The two secretions can act in an additive, multiplicative, or synergistic (super-additive) ways on the carrying capacity (K coop ). In the multiplicative case, both secretions can act only together. Then, coexisting specialist cooperators are in an interdependent mutualistic interaction in this case. Since we were interested in the possibility of interdependent mutualism between bacterial strains, we assumed that secretions act multiplicatively on K coop (Oliveira et al., 2014) . Thus we use the following multiplicative K coop function K coop = K max g 11 + g 10 g i j g 11 + g 01 g i j .

Similarly, we can define a dynamic model of strains producing any number of different leaky materials (see the Supplementary Material). Later we will study the behavior of 1-secretion, 2-secretions and 3-secretions models in detail.

| The structured population model

The above defined dynamical model describes the dynamics of cooperators and cheaters following a simple life cycle and living in an intensively mixed population. Density-dependent private benefit is the consequence of the limited diffusion of leaky materials and the limited motion of cells, therefore, some spatial structure is assumed in the background. This model reflects microbiomes living in aquatic or viscous but homogeneous habitats. However, numerous microbial communities follow a life cycle which can be modelled better by the so-called structured deme model (Wilson, 1975; Cremer et al., 2012; Oliveira et al., 2014; Chuang et al., 2009; Nadell et al., 2016) . Specialists can produce only certain secretions depending on their genotype (green stars and blue squares) while generalists produce both secretions. Each local population grows according to eq. (10) reaching a steady state, and then they merge, establishing new local populations randomly. This cycle repeats until the genotype frequencies reach a steady state for the whole community.

First, we studied the dynamical behavior of the 1-secretion model through analytical and numerical methods. Let us start the analysis with the case when δ = 0, where the system is well-mixed, and there is no private benefit. It is easy to show that cheaters always exclude cooperators from the population in this case, since dg 0 /d t > dg 1 /d t for every g 0 , g 1 > 0. This result is not surprising since this is a dynamic model of the public goods game in a well-mixed environment. However, if δ > 0 then the system's dynamical behavior becomes more complex. Assume that the cheater is the resident genotype and the cooperator is the invader. That is g 0 → K base , while g 1 → 0 and K coop → 0.

Thus

Since r > r 1 , dg 1 /d t < 0, this means that rare cooperators can never invade the population of cheaters. Similarly, rare cheaters cannot invade the population of cooperators if dg 0 /d t < 0 in the g 0 → 0 limit. Since K coop → K max , g 1 → g * 1 , (where g * 1 is the equilibrium density of the cooperators is the positive solution of g 1 = r 1

, then by comparing eqs. (6) and (7) we have

which supports dg 0 /d t < 0. After rearranging eq. (13), we find that cooperators are resistant against the invasion of rare cheaters if

Consequently, if eq. (14) is valid then the monomorph cheater and monomorph cooperator populations are the alternative locally asymptotic states of the system. If eq. (14) is not true and ∆(g ) is an increasing function then

is always true, since g 0 + g 1 decreases due to the invasion of cheaters and ∆(g ) also decreases. In this way, the cheater strategy is globally stable. The situation is different if ∆(g ) is a decreasing function and eq. (14) is not true.

Thus rare cheaters invade cooperators. However, cheaters' invasion causes g 0 + g 1 to decrease and ∆(g ) to increase.

Consequently, the fitness of cheaters decreases while the fitness of cooperators increases. If there is a g * * 0 , g * * 1 > 0 density pair where dg * * 0 d t > dg * * 1 d t = 0 then this is a stable equilibrium point, in which cooperators and cheaters will be in a stable coexistence.

Based on the above results, we can give a sufficient condition when the cheater strategy invades and excludes the cooperator strategy. Since δ ≥ ∆(g ), it follows from eq. (14) that if

then cheaters win against cooperators. Without loss of generality, we might assume that r = 1. (A different r only rescales the time unit.) Assuming a low cost (c r ) and a high cooperative carrying capacity compared to the basic (K base K max ) the above relation simplifies to δ/c < 1/2.

The qualitative analysis described above is confirmed and visualized by numerical simulations. Assuming that all parameters are fixed except δ and c, and assuming that c << r then eqs. (14) and (15) are determined only by δ/c, the maximal relative private benefit of cooperation. Fig. 3 presents the stable equilibria of eqs. (6) and (7) Density F I G U R E 3 Stable fixed points of the 1-secretion model as a function of δ/c. Red lines denotes the equilibrium density of cheaters, black lines denotes the equilibrium density of cooperators. The cheater state is globally stable for δ/c below a critical level and locally stable above it (horizontal red line at density 1). a) Private benefit decreases with density. b) Private benefit increases with density. c = 0.01, r = 1, K base = 1, K max = 10 5 , K ∆ = K max /2 δ/c 1/2. There is a polymorphic stable equilibrium of cheaters and cooperators between 1/2 δ/c 3/2 , while the cheater strategy becomes a locally stable state if δ/c 1/2. If δ/c 3/2 the system will be bi-stable with alternative monomorphic stable states, the dynamics converges either to the cooperative or to the cheater state ( Fig.   3a ). The situation is simpler if the private benefit increases with density. Then again the cheater strategy is the only stable state below a critical δ/c, and the system becomes bi-stable with the pure cooperative and pure cheater state above this δ/c level (Fig. 3b) . Naturally, there is an unstable polymorphic fixed point between the alternative stable states, which is not depicted in the figure. For the analysis of the 2-secretions model, again it can be shown easily that cheaters win over cooperators if there is no private benefit (δ = 0). Using a similar argument from the 1-secretion model, we can give a sufficient condition for the competitive dominance of cheaters over generalist cooperators ( g 00 d t > g 00 d t , for every g 00 , g 11 > 0). Knowing that ∆(g ) ≤ δ cheaters are dominant over generalist cooperators if

Assuming that δ 1 and thus δ 2 δ, cheaters exclude the generalist cooperators if

Similarly, cheater is dominant over the specialist cooperator [1, 0] 

where r i j = r 10 or r 01 . We assume here that both specialist cooperators are present in the population, thus K coop > 0.

By neglecting the δ 2 term as before the above relation leads to

Since eq. (17) follows from eq. (19), eq. (19) guarantees the competitive superiority of cheaters over any cooperative strategy. Again, if K base K max and r = 1, then δ/c < 1/2 gives a sufficient condition for δ/c below which the cheater is the only stable state. The consequence of the above analysis is that if δ/c is greater than a critical value then the competitive superiority of the cheaters is not supported. We examine the local asymptotic stability of the single strategies in this case by studying their resistance against the rare invader strategies. Following the same way of thinking as in the 1-secretion model (see eq. 12) it can be shown that neither the rare generalist cooperator nor the pair of rare specialists can invade the population of cheaters.

Let us now consider the case when the generalist cooperator is the resident strategy. 

where r i j refers to either r 10 or r 01 and g * 11 (= r 11

is the equilibrium density of the generalist strategy if it is alone in the population. Assuming again that δ 1, thus neglecting the ∆ 2 term and if

then the generalist cooperator strategy is stable against the invasion of rare specialists. Applying a very similar calculation we determine that the generalist cooperator is resistant against the invasion of the rare cheater strategy if ∆(g * 11 ) > 1 2 r − r 11 r + r 11 1 + K base K max .

Consequently, if both eqs. (21) and (22) are valid, then the generalist cooperator strategy is stable against all possible rare invaders. Whether eq. (21) or (22) is the stricter condition depends on ν and c. The linear case (ν = 1) in eq. (22) follows from eq. (21), that is if the generalist strategy is resistant against the invasion of specialist cooperators then it is stable against cheaters as well.

Since we have already shown that rare cooperators never invade the cheater population, the remaining step is to study the stability of specialist cooperators against the invading generalist cooperators and cheaters. We know that the dynamics of strategy [0,1] and [1,0] are determined by the same equations, consequently g * 01 = g * 10 and thus K coop = K max /4 in the steady state, if the cheater and generalist cooperator strategies are missing from the system.

After the same analysis as before, we conclude that specialist cooperator [0,1] is resistant against the invasion of the rare cheater strategy if

and resistant against the invasion of the rare generalist cooperator strategy if

The same relations are valid for genotype [1, 0] , only the indices are changed. The above conditions can be satisfied simultaneously only if (r − r 01 /(r + r 01 ) < (r 01 − r 11 )/2r 11 , which can be true if ν ≥ 1, that is when the cost is additive or accelerating.

The last option is that the specialist and generalist cooperators invade the population each others mutually. Generalist cooperators invade specialists when eq. (24) is not valid, and specialists invade generalists if eq. (21) is not true. These conditions can be satisfied if ∆(g * 11 ) < ∆(g * 01 ). Since g * 11 > g * 01 (except if c and ν are unrealistically large), mutual invasion is possible if ∆(g ) is a decreasing function. The consequence of mutual invasion is that generalist and specialist strategies could not exclude each other, but they coexist in this case.

The common nature of the above relations is that if all parameters are fixed except c and δ and it is assumed that c r then δ/c is the only parameter which determines the dynamics of the system. To demonstrate this, let us consider the eq. (22). After some rearrangement we have δ/c > 1

where G * 11 is either g * 11 /K ∆ or K ∆ /g * 11 depending on whether the private benefit decreases or increases with density. Since c r the right hand side of the relation is approximately independent of c. The relation is determined only by δ/c.

Besides the qualitative analysis presented above, we performed numerical simulations to reveal the behavior of the system in a more comprehensive way and to show how the parameters affect the dynamics. The most important results are summarized in Fig. 4 and Fig. 5 . We selected a basic parameter set and showed the stable states of the dynamics as a function of δ/c with this parameter set ( Fig. 4 b and Fig. 5 b) . In line with the mathematical analysis, the cheater strategy is the only stable state if δ/c is below a critical level. Above this level, the cheater strategy becomes only locally stable (a rare alternative strategy cannot invade) and alternative stable states emerge as well. When the private benefit decreases with density, first the cheater and the specialist cooperators will be in coexistence. We note that this stable state remains hidden in the local analysis since rare specialists never invades the resident cheater strategy. As δ/c increases further cheater strategy is dismissed by the generalist cooperator in the polymorphic equilibrium, and increasing δ/c even further results in that the generalist cooperator and the cheater will be the two alternative stable states of the dynamics (Fig. 4) .

Again in harmony with the local mathematical analysis, the situation is totally different if the private benefit increases with density. The system remains bistable, that is, the equilibrium states depend on the initial densities, but there is no coexistence of different strategies. In Fig. 5 , we depict the stable states starting the dynamics from three qualitatively different initial densities which correspond to either cheater, the generalist cooperator, or the specialist cooperators being the dominant strategies initially and the others are rare. As Fig. 5b shows if the cheater strategy is the dominant initially, then always the cheater wins the selection (line ch on Fig. 5 ). If generalist is the dominant initially then above a critical δ/c this strategy will be the winner, while below cheater wins (ge lines on The stable densities at the basic parameters (r = 1, c = 0.01, ν = 1, K base = 1, K max = 10 5 , K ∆ = K max /2). c) Stable densities at a higher cost (c = 0.05). d) Stable densities when ν is smaller than one (ν = 0.5). e) Stable densities when ν is bigger than one (ν = 1.5). f) Stable densities at smaller K max (K max = 10 3 ).

The effect of parameters on the dynamical behavior of the 2-secretions system is presented in the other panels of Fig. 4 and Fig. 5 . We depict the same bifurcation diagrams, however K ∆ , c, ν and K max were modified separately.

As we have shown analytically, increasing the parameter c modifies the dynamical behavior of the system only very moderately at the small c limit (Fig. 4 c and Fig. 5 c) . Similarly, changing the K max value causes only minor changes in the behavior, still the K max K base limit remains valid (Fig. 4 f and Fig. 5 f) . As our mathematical analysis predicts, the generalist cooperator strategy is supported if ν < 1 and suppressed if ν > 1 (Fig. 4d ,e and Fig. 5d ,e). The effect of changing K ∆ is presented in (Fig. 4a and Fig. 5a ). If the private benefit decreases with density, then decreasing K ∆ causes a faster loss of private benefits with density, this generalist strategy can be invaded by specialists strategies more easily (see eq. 21). Therefore, specialist strategies are present in stable coexistence with generalist strategies in a wider δ/c interval at higher K ∆ (compare Fig. 4a and b) . In contrast, if the private benefit increases with density, then decreasing K ∆ increases the private benefit faster at smaller densities, and thus the generalist strategy is more successful against the invasion of the cheater strategy.

Our analysis also included the 3-secretions model. Since analytical derivation becomes even more complex for three traits and thus for eight genotypes, and the results remain qualitatively the same, we focus only on the numerical (ch) indicates the bifurcation diagram when cheaters are dominant initially (g 00 (0) = 1, g 10 (0) = g 01 (0) = g 11 (0) = 0.1). In this case, cheaters always win the selection (horizontal red line). Arrow (ge) and dashed line indicate the bifurcation diagram when generalist is common initially (g 00 (0) = 0.1, g 10 (0) = g 01 (0) = 0.1, g 11 (0) = 10 4 ). Arrow (sp) indicates the steady states when specialist cooperators are common initially (g 00 (0) = g 11(0) = 0.1, g 10 (0) = g 01(0) = 10 4 ). There can be an interval of δ/c where specialists (blue) are resistant against the invasion of other strategies. Below that interval the cheaters and above that the generalist cooperators excludes the specialists. b) The stable densities at the basic parameters (r = 1, c = 0.01, ν = 1, K base = 1, K max = 10 5 , K ∆ = K max /2.) a) Stable densities at a smaller K ∆ (K ∆ = K max /4). c) Stable densities at a higher c (c = 0.05). d) Stable densities when ν is smaller than one (ν = 0.5). e) Stable densities when ν is bigger than one (ν = 1.5). f) Stable densities at smaller K max (K max = 10 3 ).

of the C h and C o1 emerges as an alternative stable state. As δ/c increases further,the C h strategy is excluded from this stable polymorphic state, and the C o1 and C o2 strategies will coexist. For even larger δ/c values, C o1 is excluded by the C o3 and C o2 strategies which are in coexistence. Further increase of δ/c lead to the exclusion of the C o2 and the generalist cooperator (C o3) and cheater (C h) strategies are the two alternative stable states (Fig. S1a) . Similar to the 2-secretions model, the coexistence of strategies is not possible if the private benefit increases with population density. Below a critical level of δ/c the cheater strategy is the only stable state, while above this value one of the cooperator strategies (typically the generalist cooperator) becomes an alternative stable state (Fig. S1b ).

In the applied structured population model it is common that strategies that increase the average fitness of the local population, can spread and even fixate in the whole population even though these can be competitively inferior to other strategies locally (Smith, 1964; Wilson, 1975; Chuang et al., 2009 ). In our model system, the cooperative genotypes are probably supported by the population structure. For the 1-secretion model the intuitive explanation is simple: consider a case when the cheater excludes the cooperator (δ/c < 1/2) in a local population. This will happen in every local population where there is at least one cheater initially. In those local populations which contain only cooperators initially, populations grow to a much higher population density (K max K base ) than in those "infected" by the cheaters. Consequently, the local cooperative populations send a lot of cooperators to the next generation, balancing or even overriding the weakness of cooperators in competition with cheaters locally. However, here we studied a more complex situation where besides the complete cheaters, the specialist and generalist cooperators are in competition with each other where the private benefit of the producers is density dependent.

To compare the dynamical behavior of the local population to the structured population, we used the same parameter set for the simulations as before. Further, we selected the δ/c values to scan the qualitatively different regimes of the dynamics. The total population is divided into m = 10 3 local populations and the average number of founder individuals, n varies from 2 to 50 through the simulations (low and high heterogeneity). Initial densities are the same for all strategies during the simulations, however, the random selection of individuals into the local populations can form many different initial configurations and the local populations always start to grow from a different composition of genotypes.

According to the detailed numerical simulations in the 1-secretion model, cooperators exclude cheaters or are in a stable coexistence with them if a few individuals form a local population initially (roughly n < 10), while the equilibria are similar to those that follow from local population dynamics if there are numerous founder individuals (n = 50) (Fig.   S2 ). The population structure clearly supports cooperation.

As we have shown above, there are qualitatively different regions in dynamical behavior of 2-secretions models.

Thus, to receive a comprehensive picture about the dynamical behavior of the structured population, we selected δ/c = 0.3, 0.6, 1, and 2, when the private benefit decreases (see Fig. 4 ) and δ/c=0.5, 1, and 1.5, when the private benefit increases with density (see Fig. 5 ). Fig. 6 summarizes the dynamical behavior of the structured populations. One of the most unambiguous characteristics of the selection dynamics in a structured population is that specialist cooperators will be the winners of the selection or a member of a coexisting community even under wider circumstances than in unstructured populations.

Interestingly, specialist cooperators exclude the alternative genotypes even for a low maximal relative benefit when n is at the intermediate level. A very low n favours the generalist cooperators, while the cheater strategy excludes the others at a high n as in the local populations (Fig. 6 a,b) . The explanation is the following: although generalist cooperators are competitively inferior to cheaters and specialist cooperators at a small δ/c, there is some probability that an initial local population is formed only by one genotype at very low n (n < 4). If only the generalist cooperators inhabit a local population initially then this population reaches a very high density (nearly K max ) at the end of the growth phase, thus transferring a high number of generalists to the next generation. Naturally, it is highly probable that local populations will be founded by specialist cooperators and cheaters initially, however, specialist cooperators have a significantly smaller carrying capacity (nearly K max /4) than the population of generalists, not to mention that the cheaters' carrying capacity is K base . As a consequence, generalist cooperators exclude their competitors. As n increases, the probability of having a founder population with a single strategy decreases. This is detrimental to the generalist cooperators since they are competitively inferior to the other strategies while the specialist cooperators suffer less. Furthermore, higher n increases the probability that interdependent mutualistic couples be formed initially. Consequently, specialist cooperators win the competition for an intermediate n.

As n becomes large (n ≥ 50)

almost every local population contains at least one cheater strategy initially, which excludes the cooperators in the local population. Thus, the cheaters win as in the well-mixed local population.

Similarly to the single population model, specialist cooperators can coexist with either cheaters (Fig. 6c) or generalist cooperators (Fig. 6e) when the private benefit decreases with density and δ/c is not too small or not too high.

According to the above argument, very low n is favorable for the generalist cooperators, while an intermediate level of n decreases the success of cheaters in competition with specialist cooperators (Fig. 6c) . Population structure does not change selection in populations at high δ/c (Fig. 6g) , that is the generalist cooperator wins the competition. When the private benefit increases with density and the δ/c ratio is above a critical level, then the generalist cooperators are the typical winners of the selection as in the single population model (Fig. 6d,f) .

We analyzed how the dynamical behavior of the structured population depends on the other parameters of the model. Similar to the local population, decreasing K ∆ enhances the parameter space where specialist cooperators are favored (Fig. 7a,b) . Likewise, the accelerating cost of producing multiple materials, that is if u > 1, favors the specialist cooperators (Fig. 7c, d) . Naturally, changes of these parameters in the opposite directions weaken the success of specialists (not shown).

Knowing that the dynamical behavior of the local population for the 3-secretions model is similar to the 2-

The effect of model parameters on selection. a), c) private benefit decreases with density, b), d) private benefit increases with density. a) A smaller K ∆ supports the interdependent mutualists if the private benefit decreases with density, c), d) an accelerating production cost (ν > 1) favor interdependent mutualism. Changing the parameters in the opposite direction favors the generalist producer. Other parameters are the same as in Fig 6. secretions model, it is not surprising that this qualitative similarity continues in the structured population. Similarly to the 2-secretions model we found that the C h strategy wins the competition only if the relative private benefit of cooperation is small and the mean number of founder cells is high (Fig. S3 ), while C o3 always wins the competition if there are only two or three founder cells in the local population (Fig. S3 ). Further, if the relative private benefit is high, then C o3 will be the winner of the selection, independently of the number of founder cells (Fig. S3 g, h) . C o3 is competitively inferior to any other strategy if δ/c is small or not big ( Fig. S1 and Fig. S4 ), so this strategy is excluded in every local population where there are individuals from the other strategies initially. This happens more and more frequently as n increases. On the contrary, C o2 cooperators are successful when there are only a few founder individuals locally and δ/c is small or moderate (Fig. S3a, b,c) . In these cases, there is a high probability that two different genotypes of C o2 are placed together in the same local population. They mutually help each other; thus this local population can reach a high population density. Although C h or any member of C o1 excludes the mutualistic C o2 pairs (see Fig. S1 and Fig. S4 ), the chance of these competitors being together with the C o2 pairs in the founder populations is low if n is not too large. Increasing the number of founders further makes co-residence of C h and C o1 with C o2 more probable initially. Consequently, the relative benefit of C o 2 decreases. In parallel, the complete co-residence of the interdependent C o1 community occurs with higher probability in local populations, increasing the relative benefit of this strategy. Thus C o1 will be the winner of the selection or will be dominant in the whole population at a medium n (Fig. S3 a,b, c) . As in the 2-secretions model, different cooperative strategies are in coexistence only if the private benefit decreases with density, δ/c is roughly 0.5 to 1.5 and n is not too small to prevent the emergence of cooperative coalitions initially (Fig. S3 c, e) . According to the numerical simulations (not shown), the parameter dependence of the results remains the same as in the 2-secretions model: smaller K ∆ and ν > 1 are favorable for the specialist cooperators, while ν < 1 increases the parameter space where the generalist win the competition.

We studied a model of microbial mutualism where the competing cells can produce different types of leaky materials.

The production of leaky materials was modeled to be costly but producing them serves as a private benefit for the producer cell, and this private benefit depends on population density. Further, we considered a situation when all of the different leaky materials are necessary to cause a positive effect and the specialist producers are completely interdependent on each other. Different analyses with competition between cell types, either the one's non-producing leaky materials to the types which produce all different secretions were done, in an intensively mixed population with simple life cycle and in a structured population with a more detailed life cycle.

Naturally, if the private benefit is missing or weak compared to the production cost, then the situation is identical to the public goods game and then the cheater strategy is the only stable state in the intensively mixed population. However, as we showed with mathematical and numerical methods, mutually dependent cooperative strains are present in a steady state in a well mixed population for a broad range of conditions, if the private benefit decreases with density ( Fig. 4, Fig S1a) . Depending on the maximal relative private benefit (δ/c), the mutually dependent cooperative strains coexist either with more or less generalist cooperative strains ( Fig. 4 and Fig. S1 a) . These mutualist strains are either exploited by strains producing a lower number of different materials including the cheater strain, or these strains exploit the more generalist strain in the coexisting community. In those cases, when the private benefit increases with density (and δ/c is above a critical level) the population has two equilibrium states. One of them is the cheater strategy and the other is typically the generalist cooperator strategy. Mutually dependent cooperative strains win the competition only under very specific circumstances ( Fig. 5 and Fig. S1 ).

In agreement with the results of previous work (Oliveira et al., 2014; Estrela et al., 2016) , if producing more secretions is more costly, that is if there is a tradeoff between the synthesis of the secretions, then the generalist strategy is suppressed and interdependent mutualism is supported. Contrary, if producing more secretions is less costly, that is if the synthesis of different secretions is present in multiple genotypes then the generalist strategy is more supported, and specialists are suppressed (Fig. 4, Fig. 5, Fig. 7) .

The speed of change of density dependence when density is small is determined by K ∆ in our model. A quicker decrease (smaller K ∆ ) of density dependence of the private benefit helps the interdependent mutualists to coexist with the generalist (Fig. 4a,b and Fig. 7a ). The situation is the opposite if the private benefit increases with density: faster saturation to the maximum private benefit increases the success of the generalists and suppresses the specialists (Fig. 5a,b and Fig. 7b ). The intuitive explanation is the following: specialists have a competitive benefit over the generalists at lower densities since their density independent growth rate is larger than the generalist's growth rate ( r 10 = r 01 > r 11 ), but the generalists can have a fitness benefit over specialists at higher densities because they can reach a higher equilibrium density. The relative dominance of these opposite effects is affected by the density dependence of the private benefit. If it is strong at lower densities and negligible at high densities, this supports the specialists, if it is significant at high densities this promotes the generalist strategy.

In similar work, Estrela et al. (Estrela et al., 2016) importantly because no private benefit was included in their model. Furthermore, they used a population dynamical model in which all genotypes coexist in a neutrally stable equilibria in the local populations, which makes the model structurally unstable (Arnold, 1988) .

Theoretical and experimental studies pointed out earlier that cooperators are supported in the spatially structured populations (Nowak, 2006; Chuang et al., 2009) . However, these studies considered only the selection of a cooperative and a cheater strategy, different common goods and the mutual dependence of specialized cooperators, in other words, the presence of partial cheaters was not possible. One of the few exceptions is the work by Oliveira et al. (2014) who they focused on the possibility of mutualistic interactions between specialized cooperators in a structured population. They have found that microbial mutualism is possible only if the number of founder individuals is not extremely large (n = 4 − 8) and the cost of producing the leaky material is high (c > 0.1). This result led them to the conclusion that interdependent mutualism among the partners competing for the same niche should be rare in microbial communities. We used in our model the same population structure as Oliveira et al (2014), but our conclusion is significantly different because density dependent private benefit is taken into account. Comparing the dynamical behavior in the structured model with the well-mixed model (see e.g. Fig. 4, Fig. 5 and Fig. 6) , it is apparent that dynamics in the local population are strongly retained in the structured population, although the population structure further increases the possibility of microbial cooperation. Naturally, if n is small then the generalist genotype will be the winner of the competition, but specialist mutualists win or coexist with other genotypes for a larger number of founder cells under very different circumstances (Fig. 6 Fig. S3 ). Further, the mutually dependent cooperators are selected in the structured population even in those cases when the cheater strategy is the winner of the competition in the well-mixed population. Our results suggest that interdependent mutualism should not be rare in real microbial communities, especially in those in which their private benefit decreases with density. The selection of interdependent cooperation and coexistence of traits with different levels of cooperation, in general, is in harmony with the Black Queen Hypothesis which explains the loss of essential leaky functions of microbes as an adaptive process (Morris et al., 2014; Morris, 2015) . The consequence of the adaptive loss of the production of leaky materials is that the effectiveness of communities with specialist cooperators is sub-optimal since maximal population density can be attained if generalist cooperators dominate the community (Fig. 4) .

The above conclusion is modulated by the fact that we studied a worst case scenario from the point of view of cross-feeding, since we assumed that all secretions are needed for the enhanced carrying capacity. Further, we assumed that genotypes use the same niche, in this way, competition is maximal between the different strategies.

Ruling out these assumptions, the possibility of coexistence of different genotypes increases and we hypothesize that cross-feeding, the complex networks of more and less specialized cooperators, and complete cheaters will be even more widespread in these communities including those where private benefit increases with density. These questions could be explored in a further study.

In summary, the character of density dependence of the private benefit determines essentially how easily interdependent mutualism can evolve. If the private benefit decreases with density, then interdependent mutualist genotypes typically coexist with cooperative strains while if the private benefit increases with density then interde- δ/c Density F I G U R E S 4 Competition between the cooperators if C h is missing. a) The private benefit decreases with density. All cooperative strategies are in competition. Initial densities are 0.1 for all strategies. b), c), d) The private benefit increases with density and pairwise competition. b) Competition between C o1 and C o2. c) Competition between C o1 and C o3, d) Competition between C o2 and C o3. Initial density is 1 for the more generalist strategy and 0.1 for less generalist strategy in all cases while other strategies are not included. Color code: C o1 (blue), C o2 (green), C o3 (black).

Game theory of public goods in one-shot social dilemmas without assortment

Evolution of optimal Hill coefficients in nonlinear public goods games

Geometrical Methods In The Theory Of Ordinary Differential Equations

Microbial contributions to climate change through carbon cycle feedbacks. The ISME journal

Comparative analysis of dna extraction methods to study the body surface microbiota of insects: A case study with ant cuticular bacteria

Insights into the coral microbiome: underpinning the health and resilience of reef ecosystems

Cooperation in microbial communities and their biotechnological applications

The human microbiome: at the interface of health and disease

Simpsons paradox in a synthetic microbial system

Public good dynamics drive evolution of iron acquisition strategies in natural bacterioplankton populations

The ecology of the microbiome: Networks, competition, and stability

Growth dynamics and the evolution of cooperation in microbial populations

Microbial communication, cooperation and cheating: quorum sensing drives the evolution of cooperation in bacteria

Interaction effects of cell diffusion, cell density and public goods properties on the evolution of cooperation in digital microbes

Solutions to the public goods dilemma in bacterial biofilms

Private benefits and metabolic conflicts shape the emergence of microbial interdependencies

Diminishing returns in social evolution: the not-so tragic commons

Competition, not cooperation, dominates interactions among culturable microbial species

Snowdrift game dynamics and facultative cheating in yeast

Novel cooperation experimentally evolved between species

Metabolic resource allocation in individual microbes determines ecosystem interactions and spatial dynamics

Evolutionary establishment of moral and double moral standards through spatial interactions

Microbial ecology along the gastrointestinal tract

The role of soil microorganisms in plant mineral nutritioncurrent knowledge and future directions

Isolating "uncultivable" microorganisms in pure culture in a simulated natural environment

The ecogenetic link between demography and evolution: can we bridge the gap between theory and data?

Habitat structure and the evolution of diffusible siderophores in bacteria

When increasing population density can promote the evolution of metabolic cooperation

Symbiotic bacteria on the cuticle of the leaf-cutting ant acromyrmex subterraneus subterraneus protect workers from attack by entomopathogenic fungi

Syntrophic exchange in synthetic microbial communities

Black queen evolution: the role of leakiness in structuring microbial communities

Coexistence of evolving bacteria stabilized by a shared black queen function

Spatial structure, cooperation and competition in biofilms

Five rules for the evolution of cooperation

Evolutionary limits to cooperation in microbial communities

Fitness and stability of obligate cross-feeding interactions that emerge upon gene loss in bacteria

Regulation of bacteriocin production in streptococcus mutans by the quorum-sensing system required for development of genetic competence

Resolving the tragedy of the commons: the feedback between population density and intraspecific conflict

Mutual cross-feeding interactions between bifidobacterium longum subsp. longum ncc2705 and eubacterium rectale atcc 33656 explain the bifidogenic and butyrogenic effects of arabinoxylan oligosaccharides

Density dependence and cooperation: theory and a test with bacteria

Microbiota in social brain

Group selection and kin selection

The classification and evolution of bacterial cross-feeding

The public and private benefit of an impure public good determines the sensitivity of bacteria to population collapse in a snowdrift game

An early mechanical coupling of planktonic bacteria in dilute suspensions

Growing unculturable bacteria

Phase transitions and volunteering in spatial public goods games

The ecology and evolution of social behavior in microbes

The microbiome and mental health: Hope or hype?

Microbes in the coral holobiont: partners through evolution, development, and ecological interactions

Role of the gut microbiota in nutrition and health

A theory of group selection

Bacterial cheating drives the population dynamics of cooperative antibiotic resistance plasmids

] is denoted by black. a) stable states occur when private benefit decreases with density. b) Stable states occur when private benefit increases with density. Arrows (C h, C o1) indicate the stable state when the cheater or C o1 strategies are common initially (g i j k = 0.1, (i , j , k ) {0, 1}, except g 000 (0) = 1, or g i j k (0) = 0.001, (i , j , k ) {0, 1}, except g 110 = g 101 = g 011 = 10 3 ). Arrow (C o3) and dashed lines indicate the stable states when generalists are common (g i j k (0) = 0.001, (i , j , k ) {0, 1}, except g 111 (0) = 10 3 ) and arrow (C o2) denotes the steady states when

Let us assume that strains can produce n ≥ 1 number of different secretions. Let us denote a genotype of a strain with [i 1 , i 2 , ...i n ], where i k ∈ {0, 1} denotes whether the actual material is secreted (i k = 1) or is not secreted (i k = 0) by the strain. Thus the general dynamical system is 

Here are the figures of 1-secretion and 3-secretions models, which we refer in the main text.