# Delays in fitness adjustment can lead to coexistence

of hierarchically interacting species

###### Abstract

Organisms that exploit different environments may experience a stochastic delay in adjusting their fitness when they switch habitats. We study two such organisms whose fitness is determined by the species composition of the local environment, as they interact through a public good. We show that a delay in fitness adjustment can lead to coexistence of the two species in a metapopulation, although the faster growing species always wins in well-mixed competition experiments. Coexistence is favored over wide parameter ranges, and is independent of spatial clustering. It arises when individuals have different fitness values that can keep each other balanced.

How biodiversity – for example, the surprising coexistence of more than species in a single gram of soil Torsvik et al. (2002); Daniel (2005); Gans et al. (2005) – is stabilised is one of the most fundamental questions in ecology Levin (1992). Known stabilising factors derived from resource competition models include metabolic trade-offs Posfai et al. (2017) and reciprocal oscillations in population sizes Huisman and Weissing (1999). Cyclic competition models and their derivatives are also frequently employed to model biodiversity Szolnoki et al. (2014); Xue and Goldenfeld (2017); Shih and Goldenfeld (2014). However, recent experiments that have investigated interactions between a variety of different soil species in well-mixed pair-wise competition experiments found that these species could not be represented by cyclic competition models; instead, a few species outcompeted the others Higgins et al. (2017). Thus, microbial diversity in soil is thought to be supported by the highly porous and fragmented structure of this habitat Young and Crawford (2004), and the fluctuating environmental conditions that individual bacteria experience there Higgins et al. (2017).

Here, we ask how the coexistence of species is affected by spatial structure, characterized by intrinsic variations between local environments, using a simple model in which species adjust to these local environmental changes with a delay. Indeed, it is well-known that such delays in physiological responses (here, adjustment of fitness or growth rate) occur in microbes, following externally imposed changes in the environment Ackermann (2015); Lachmann and Jablonka (1996); Thattai and Van Oudenaarden (2004); Kussell and Leibler (2005), such as nutrient composition, or antibiotic stress Monod (1949); Lin and Kussell (2016); Fridman et al. (2014); Bollenbach et al. (2009); Yurtsev et al. (2013); Wu et al. (2013); Skanata and Kussell (2016). Here, we focus on changes that occur because species move between different habitats. We consider a system with two species, in which the dominant strain (fast grower) depends on the weaker strain (slow grower) for its fitness. Thus, the fitness of both species depends on the population structure of the local habitat. We assume that the fitness of an individual is not instantaneously reset when the environment changes, but is initially retained from its previous environment, as we discuss below.

We show that coexistence can arise in minimal two-species models, as a result of delays in fitness adjustment after a change of habitat, and nonlinearity of fitness functions. In particular, we conclude that the combination of underlying spatial structure – not just spatial segregation – and delay in adapting to environmental change can support biodiversity, by enabling individuals of the same species that have different fitness values to keep each other balanced. A key result of our study is that, because of this delay in fitness change, the outcomes of direct competition experiments between species in well-mixed systems may differ from those observed in spatially structured habitats.

Model. We study a metapopulation of locally well-mixed patches. Individuals can move or hop between patches with rate , which we refer to as the mixing rate, representing processes which couple individual patches (Fig. 1a). The local species composition on a patch with individuals determines their reference fitness on that patch. A well-known implementation of such an eco-evolutionary model is, for example, a public goods dilemma Velicer (2003). There, a ‘producer’-species produces a public good, which bestows a fitness benefit on the entire population Rainey and Rainey (2003); Greig and Travisano (2004); West et al. (2006); Diggle et al. (2007); Gore et al. (2009); Cordero et al. (2012); Drescher et al. (2014); Nadell et al. (2016); Becker et al. (2018), while imposing a fitness cost on its producer. Hence, the producer is at a disadvantage relative to a faster growing ‘non-producer’ species Hardin (1968); Axelrod and Hamilton (1981).

We study species that interact via such public goods as an examplary interaction topology that results in hierarchical population structure. As the precise dependence of a species’ fitness on species composition can vary between different experimental setups Damore and Gore (2012), it makes sense to study a conceptual model, which minimizes the number of tunable variables. Here, we assume that the growth rate or fitness of individuals depends on the fraction of slow growers on a patch, ,

(1) |

where we set the basal level of fitness to ensure that the fast grower can grow by itself (i.e. in the absence of the public good), and discuss the nonlinearity coefficient in the next paragraph. We distinguish between the fast-growing and slow-growing strain by introducing a fitness difference between them. Thus, the reference fitness of the fast and slow growers in the same local environment are and , respectively.

Our choice of fitness function means that the local amount of the public good (such as invertase for yeast or glutamine synthetase for B. subtilis Gore et al. (2009); Liu et al. (2015)) on a given patch increases with increasing producer fraction on that patch, which in turn enhances the fitness of the species on the patch. A coefficient implies that the impact of the public good on an individual’s fitness saturates when a large amount of public good is present (see e. g. Refs. Gore et al. (2009); Archetti et al. (2015); Heilmann et al. (2015)). By contrast, means that a significant fraction of producers is required on a patch in order for the fitness to change noticeably Cornforth et al. (2012). We note that this type of interaction, where the fitness depends on a group of players, is referred to as an -player game in the context of game theory Mitteldorf and Wilson (2000); Taylor (1992); van Baalen and Rand (1998); Archetti and Scheuring (2012).

Because of the different numbers of slow growers that can inhabit a patch, individuals of the same species may have different fitness values (Fig. 1). For example, for a local patch with individuals, there are different possible compositions and thus different fitness types for each species (Fig. 1b). These different fitness types have fitness values for for the slow growers and for for the fast growers, where we have shifted the index so that the fitness of the least fit individual can be denoted by in both cases. Since the number of individuals may differ between patches, there is a large variety of different fitness types in the metapopulation at any one time.

The local species composition thus determines the reference fitness of an individual, which in turn determines its growth rate: Individuals reproduce at rates proportional to their fitness values by replacing another individual on the same patch by an identical copy of itself. This implementation of reproduction via the so-called Moran process Moran (1962) is a technicality and ensures that the number of individuals on a patch does not change (as would be the case in a chemostat Novick and Szilard (1950)), thereby avoiding stabilisation of the slow growers by known effects which we do not consider here, such as by disproportional growth of producer patches Chuang et al. (2009); Cremer et al. (2012).

Delay. So far, we have specifically used the term ‘reference fitness’ in order to distinguish between the fitness prescribed by the local environment, and the actual fitness value of an individual. Figure 1c shows the fitness of a slow-grower (black line) during a short span of its lifetime: at and , the composition of its local patch, and thus its reference fitness (dashed grey line), changes. We assume that there is a delay between the change in species composition of a patch and the time at which the individual adjusts to the reference fitness. This assumption seems appropriate, the exoproduct concentration on the patch will equilibrate over a certain timescale owing to, for example, the slow production of exoproducts as a response to detected changes in patch composition (which may in turn depend on the internal proteome of an individual at a certain time Wolf et al. (2008); Snijder and Pelkmans (2011); Raser and O’Shea (2006)), or the slow diffusion of the produced public good locally on a patch Allen et al. (2013); Borenstein et al. (2013). We absorb all these various biological processes into one stochastic fitness adjustment rate for each individual. Thus, Fig. 1c shows that the slow grower retains its fitness from the previous environment for a period that scales on average as .

Coexistence. We first discuss how the presence of a delay in fitness adjustment delay alters the species composition of the entire metapopulation. Without this delay (), the slow grower would consistently die out for all finite values of mixing rate and fitness difference , and also for all fitness nonlinearities . Extinction occurs because the fitness of the slow grower is always lower than that the fitness of the fast grower, and because we choose the system size to be large enough such that stochastic fluctuations do not influence the outcome Shnerb et al. (2000); Traulsen et al. (2005, 2006); Reichenbach et al. (2006); Dobramysl et al. (2018).

We explore the effect of delay with a Gillespie simulation Gillespie (1976), and show the resulting phase diagram for fitness functions with in Fig. 2a. Strikingly, a phase of coexistence of both species extends over a broad regime of fitness difference and mixing rates . The coexistence phase begins at mixing rates less than the basal fitness , and broadens with increasing mixing rate. As it is most pronounced for high mixing rates, we explain its presence first by discussing the well-mixed limit for all different values of . In doing so, we will also show why coexistence only occurs for .

Well-mixed limit. Without the delay in fitness adjustments, the slow grower would die out in the well-mixed limit, as reproduction takes place on fully updated patches and the fast grower is fitter for every given patch. The delay in attaining the reference fitness means that different fitness types within a species (with fitnesses adjusted to previous environments) are present on a local patch. Reproduction then occurs within a random sample of the fitness types present in the well-mixed metapopulation; thus, one can imagine that reproduction occurs across the entire population with all different possible fitness types present.

We consider the average number of individuals in the following, even though the number of individuals per patch may vary in principle due to mixing. This approach simplifies the description of our system, because it means that there are then only distinct fitness types per species in the metapopulation (Fig. 1 b).

The dynamics of the system is described by differential equations, one for each fitness type. We denote the density of a slow grower of type by with (Fig. 2b). The same applies for fast growers (). As the total number of individuals does not change, . How the density of each type changes due to fitness adjustments is a combinatorial problem that depends on the proportions of slow and fast growers in the entire system, or, more precisely, on the different possible combinations of groups of players that can be drawn from the metapopulation. The time evolution for the slow grower fitness type of density is described by

(2) | ||||

where the sum over runs over all fitness types unless otherwise specified. The first term in brackets corresponds to reproduction; the second corresponds to the delayed fitness adjustments to local environments of players (we discuss this in detail for the exemplary value in the supplement).

Before analyzing the fixed point structure of Eq.2, we use arguments from game theory to explain how coexistence can arise Maynard Smith (1982); Drossel (2001); Nowak and Sigmund (2004); Frey (2010). To do so, we investigate the evolutionary stability of the absorbing states (all fast or all slow growers) by asking whether patches with configurations corresponding to these states are invasible. This involves comparing the fitness values of the absorbing states with those of potentially invading states, corresponding to the other possible patch configurations. The panels in Fig. 2c show the fitness of slow and fast growers within a given patch configuration (shown in the supplement or for in Fig 1b), with increasing numbers of slow growers. The configuration of the fittest slow-grower contains slow growers, while the configuration for the fittest fast-grower contains only slow growers. Hence, for zero fitness difference , the fittest fast grower is slightly less fit than the fittest slow grower (left panel Fig. 2c). Thus for , the absorbing state of all slow growers is evolutionarily stable.

An increasing fitness difference reduces the fitness values of the slow grower. The fitness of the slow grower in its strongest configuration and the strongest fast grower are equal at (Fig. 2c, middle panel). Hence the state of all slow growers becomes invasible by a patch configuration with a single fast-grower for . If that single fast grower then reproduces and adjusts its fitness to that of the invaded patch, this patch configuration may in turn become liable to invasion by slow growers. This intricate interplay indicates that both species may survive. Indeed, this state of coexisting different fitness types must survive, because the absorbing state of all fast growers itself also is not evolutionarily stable (as can be shown analogously to the instability of the absorbing state of slow growers). This coexistence is stable up to a fitness difference (Fig. 2c, right panel). For the state of all fast growers becomes evolutionary stable. This game theoretical analysis is corroborated by a linear stability analysis of the fixed points of the set of Eq. 2.

As the Jacobian evaluated at the absorbing fixed points is sparse (cf. Supplement), its eigenvalues are easily accessible. For as considered here, all eigenvalues apart from one are always negative. The eigenvalue that can change sign is given by for the fast, and for the slow grower absorbing state. Thus, the fast (slow) grower fixed point is stable for () .

We plot and and thus the regions of fixed point stabilities in Fig. 2d for . For and fitness differences the outcome is bistable, meaning that the survivor depends on the initial conditions. For our initial conditions of equal amounts of fast and slow grower, the slow grower survives up to the grey dashed line. For , both absorbing fixed points are unstable in the shaded region, and a third stable coexistence fixed point emerges. This fixed point is present also for , but it is unstable. The shading in the region where coexistence can occur for shows the species composition at the fixed point, as obtained from Eq. 2. As expected, the fraction of slow growers decreases with increasing fitness difference .

This fixed point structure is robust: Even though the extent of the coexistence region depends on , it will also occur for different , as we discuss in the supplement. Importantly, coexistence does not depend on the length of delay, as the eigenvalue have the structure discussed above for . This condition on ensures that the fitness adjustment is slow enough for different fitness types to be generated repeatedly between reproduction events.

Before returning to the full phase diagram, we note that the fact that nonlinear fitness functions can yield coexistence in -player games has been discussed previously, albeit in different setups Motro (1991); Hauert et al. (2006); Archetti and Scheuring (2011); Assaf et al. (2013); Archetti and Scheuring (2012); Li and Thirumalai (2018). In our model, the -player structure alone does not lead to coexistence: only the additional timescale Levin (2000); Goldenfeld and Woese (2011); Roca et al. (2009); Wu et al. (2009) due to delay allows that reproduction to occur in a well-mixed system consisting of all possible fitness types. These fitness types can balance each other during reproduction.

The understanding of the balance between these types is important: It explains why coexistence occurs not just in the well-mixed limit, but also for more realistic intermediate mixing rates. The extension of this coexistence phase to lower mixing rates in the phase diagram in Fig. 2a indicates that the presence of a well-mixed coexistence fixed point impacts the population structure even in not well-mixed situations.

Phase diagram.
The phase diagram in Fig. 2a shows the average extinction times and phase boundaries for , and a shorter delay of for comparison for nonlinear fitness functions with ^{1}^{1}1We aborted calculations after and verified that aborted simulations yield slow- and fast-grower coexistence towards the end of their simulation runs..
The white triangles mark the fitness differences between which coexistence is expected from the well-mixed limit (cf. Fig. 2d).
Our phase boundaries do not exactly match the well-mixed expectation due to the finite size effects:
When the proportions of slow and fast growers at the coexistence fixed point become strongly biased, stochastic fluctuations can lead to extinction even when both species should coexist.
Increasing the size of the metapopulation would progressively extend the phase boundaries, which we marked with dashes where we expect system size effects to dominate.

In the weak mixing limit, the phase boundary marks a transition from a fast growing to a slow growing phase without a coexistence region. The phase boundary in this limit can be understood using invasion probabilities of single individuals on locally fixated patches Altrock and Traulsen (2009); Antal and Scheuring (2006); Blythe and McKane (2007); Bauer and Frey (2018b). For , the transition between the low and high mixing limit occurs approximately at Bauer and Frey (2018b). For , the two phase boundaries in Fig. 2a indicate that the lowest mixing rate at which the coexistence occurs also increases with ; we verified this general trend, but a precise investigation of these values would go beyond the scope of this work, as the phase diagram remains qualitatively the same.

We find it remarkable that coexistence can occur over such a wide parameter range, as this feature indicates that the stabilisation of populations afforded by delayed fitness adjustment may be a realistic effect. We note that our results do not depend on the formation of spatial patterns Nowak and May (1992); Reichenbach et al. (2007); Pigolotti et al. (2017). The behaviour of our model relies on the fact that the fittest weak individual can be fitter than the fittest fast grower; while we do not include private components of the public good Li and Thirumalai (2018); Gore et al. (2009); Zhang and Rainey (2013) explicitly, these would practically have a similar effect for larger groups of microbes. The advantage of our model is that it is mathematically more tractable, while a realistic model of private public goods includes additional parameters, but of course applies more directly to individual experimental situations.

Conclusion. We have shown that species can coexist in spatially fragmented systems with delays in fitness adjustment even when one species always dominates over the other in a locally well-mixed environment. This coexistence occurs because the combination of delay with the interaction via a public good means that species compete via a variety of different fitness types (cf. Fig. 2a) which can balance each other for some fitness functions. Our work presents a conceptual study of a two-species population in a fragmented environment, motivated by a study of soil bacteria. The local pore structure and differing mixing rates, depending on e. g. water content, nutrients and connectivity within the sample Coyte et al. (2017); Peaudecerf et al. (2018); Pérez-Reche et al. (2012); Dechesne et al. (2010), and additional mircobial interactions and community structures would need to be considered to model the situation in soil accurately. However, the robustness of coexistence over large ranges of fitness difference and mixing rate suggests that delayed responses to changes of habitat may be a significant factor in the maintenance of biodiversity, even though we only study one specific interaction. Thus, our study indicates that it is important to mimic the spatial structure when biodiversity in spatially structured habitats is investigated (see e.g. Gandhi et al. (2016); Limdi et al. (2017); Peaudecerf et al. (2018)), even in experimental configurations where spatial assortment Van Dyken et al. (2013); Pande et al. (2016) does not occur.

Acknowledgement. We are grateful to Jonas Denk, Heinrich Jung and David Muramatsu for discussions, and acknowledge funding from the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 660363, an LMU Research Fellowship, and the German Excellence Initiative via the programme ‘NanoSystems Initiative Munich’.

## References

- Torsvik et al. (2002) V. Torsvik, L. Øvreås, and T. F. Thingstad, Science 296, 1064 (2002).
- Daniel (2005) R. Daniel, Nat. Rev. Microbiol. 3, 470 (2005).
- Gans et al. (2005) J. Gans, M. Wolinsky, and J. Dunbar, Science 309, 1387 (2005).
- Levin (1992) S. A. Levin, Ecology 73, 1943 (1992).
- Posfai et al. (2017) A. Posfai, T. Taillefumier, and N. S. Wingreen, Physical Review Letters 118, 028103 (2017).
- Huisman and Weissing (1999) J. Huisman and F. J. Weissing, Nature 402, 407 (1999).
- Szolnoki et al. (2014) A. Szolnoki, M. Mobilia, L.-L. Jiang, B. Szczesny, A. M. Rucklidge, and M. Perc (2014).
- Xue and Goldenfeld (2017) C. Xue and N. Goldenfeld, Phys. Rev. Lett. 119, 268101 (2017).
- Shih and Goldenfeld (2014) H. Y. Shih and N. Goldenfeld, Phys. Rev. E. 90, 050702 (2014).
- Higgins et al. (2017) L. M. Higgins, J. Friedman, H. Shen, and J. Gore, bioRxiv p. 175737 (2017).
- Young and Crawford (2004) I. M. Young and J. W. Crawford, Science 304, 1634 (2004).
- Ackermann (2015) M. Ackermann, Nat. Rev. Microbiol. 13, 497 (2015).
- Lachmann and Jablonka (1996) M. Lachmann and E. Jablonka, J. Theor. Biol. 181, 1 (1996).
- Thattai and Van Oudenaarden (2004) M. Thattai and A. Van Oudenaarden, Genetics 167, 523 (2004).
- Kussell and Leibler (2005) E. Kussell and S. Leibler, Science (New York, N.Y.) 309, 2075 (2005).
- Monod (1949) J. Monod, Ann. Rev. Microbiol. 3, 371 (1949).
- Lin and Kussell (2016) W.-H. Lin and E. Kussell, Curr. Biol. 26, 1486 (2016).
- Fridman et al. (2014) O. Fridman, A. Goldberg, I. Ronin, N. Shoresh, and N. Q. Balaban, Nature 513, 418 (2014), ISSN 0028-0836.
- Bollenbach et al. (2009) T. Bollenbach, S. Quan, R. Chait, and R. Kishony, Cell 139, 707 (2009).
- Yurtsev et al. (2013) E. A. Yurtsev, H. X. Chao, M. S. Datta, T. Artemova, and J. Gore, Mol. Sys. Biol. 9, 683 (2013).
- Wu et al. (2013) A. Wu, K. Loutherback, G. Lambert, L. Estevez-Salmeron, T. D. Tlsty, R. H. Austin, and J. C. Sturm, Proc. Natl. Acad. Sci. USA 110, 16103 (2013).
- Skanata and Kussell (2016) A. Skanata and E. Kussell, Phys. Rev. Lett. 117, 038104 (2016).
- Velicer (2003) G. J. Velicer, Trends Microbiol. 11, 330 (2003).
- Rainey and Rainey (2003) P. B. Rainey and K. Rainey, Nature 425, 72 (2003).
- Greig and Travisano (2004) D. Greig and M. Travisano, Proc. R. Soc. Lond. B (Suppl.) 271, S25 (2004).
- West et al. (2006) S. A. West, A. S. Griffin, A. Gardner, and S. P. Diggle, Nat. Rev. Microbiol. 4, 597 (2006).
- Diggle et al. (2007) S. P. Diggle, A. S. Griffin, G. S. Campbell, and S. a. West, Nature 450, 411 (2007).
- Gore et al. (2009) J. Gore, H. Youk, and A. van Oudenaarden, Nature 459, 253 (2009).
- Cordero et al. (2012) O. X. Cordero, L.-A. Ventouras, E. F. DeLong, and M. F. Polz, Proc. Nat. Acad. Sci. USA 109, 20059 (2012).
- Drescher et al. (2014) K. Drescher, C. D. Nadell, H. A. Stone, N. S. Wingreen, and B. L. Bassler, Curr. Biol. 24, 50 (2014).
- Nadell et al. (2016) C. D. Nadell, K. Drescher, and K. R. Foster, Nat. Rev. Microbiol. 14, 589 (2016).
- Becker et al. (2018) F. Becker, K. Wienand, M. Lechner, E. Frey, and H. Jung, Sci. Rep. 8, 4093 (2018).
- Hardin (1968) G. Hardin, Science 169, 1243 (1968).
- Axelrod and Hamilton (1981) R. Axelrod and W. D. Hamilton, Science 211, 1390 (1981).
- Damore and Gore (2012) J. A. Damore and J. Gore, J. Theor. Biol. 299, 31 (2012).
- Liu et al. (2015) J. Liu, A. Prindle, J. Humphries, M. Gabalda-Sagarra, M. Asally, D.-Y. D. Lee, S. Ly, J. Garcia-Ojalvo, and G. M. Süel, Nature 523, 550 (2015).
- Archetti et al. (2015) M. Archetti, D. A. Ferraro, and G. Christofori, Proc. Nat. Acad. Sci. USA 112, 1833 (2015).
- Heilmann et al. (2015) S. Heilmann, S. Krishna, and B. Kerr, Frontiers in Microbiology 6, 1 (2015).
- Cornforth et al. (2012) D. M. Cornforth, D. J. T. Sumpter, S. P. Brown, and Å. Brännström, The American Naturalist 180, 296 (2012).
- Mitteldorf and Wilson (2000) J. Mitteldorf and D. S. Wilson, J. Theor. Biol. 204, 481 (2000).
- Taylor (1992) P. D. Taylor, Evolutionary Ecology 6, 352 (1992).
- van Baalen and Rand (1998) M. van Baalen and D. A. Rand, J. Theor. Biol. 193, 631 (1998).
- Archetti and Scheuring (2012) M. Archetti and I. Scheuring, J. Theor. Biol. 299, 9 (2012).
- Moran (1962) P. A. P. Moran, The Statistical Process of Evolutionary Theory (Clarendon Press, 1962).
- Novick and Szilard (1950) A. Novick and L. Szilard, Science 3, 715 (1950).
- Chuang et al. (2009) J. S. Chuang, O. Rivoire, and S. Leibler, Science 323, 272 (2009).
- Cremer et al. (2012) J. Cremer, A. Melbinger, and E. Frey, Sci. Rep. 2, 281 (2012).
- Wolf et al. (2008) D. M. Wolf, L. Fontaine-Bodin, I. Bischofs, G. Price, J. Keasling, and A. P. Arkin, PLoS ONE 3 (2008).
- Snijder and Pelkmans (2011) B. Snijder and L. Pelkmans, Nature Reviews Molecular Cell Biology 12, 119 (2011).
- Raser and O’Shea (2006) J. M. Raser and E. K. O’Shea, Science 304, 1811 (2006).
- Allen et al. (2013) B. Allen, J. Gore, and M. A. Nowak, eLife 2013, 1 (2013), ISSN 2050084X.
- Borenstein et al. (2013) D. B. Borenstein, Y. Meir, J. W. Shaevitz, and N. S. Wingreen, PLoS ONE 8 (2013).
- Shnerb et al. (2000) N. M. Shnerb, Y. Louzoun, E. Bettelheim, and S. Solomon, Proc. Nat. Acad. Sci. USA 97, 10322 (2000).
- Traulsen et al. (2005) A. Traulsen, J. C. Claussen, and C. Hauert, Phys. Rev. Lett. 95, 238701 (2005).
- Traulsen et al. (2006) A. Traulsen, J. C. Claussen, and C. Hauert, Phys. Rev. E 74, 011901 (2006).
- Reichenbach et al. (2006) T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. E 74, 051907 (2006).
- Dobramysl et al. (2018) U. Dobramysl, M. Mobilia, M. Pleimling, and U. C. Täuber, J. Phys. A: Math. Theor. 51, 063001 (2018).
- Gillespie (1976) D. T. Gillespie, J. Comp. Phys. 22, 493 (1976).
- Maynard Smith (1982) J. Maynard Smith, Evolution and the Theory of Games (Cambridge University Press, 1982).
- Drossel (2001) B. Drossel, Adv. Phys. 50, 209 (2001).
- Nowak and Sigmund (2004) M. Nowak and K. Sigmund, Science (New York, N.Y.) 303, 793 (2004).
- Frey (2010) E. Frey, Physica A 389, 4265 (2010).
- Motro (1991) U. Motro, J. Theor. Biol. 151, 145 (1991).
- Hauert et al. (2006) C. Hauert, M. Holmes, and M. Doebeli, Proc. R. Soc. Lond. B 273, 2565 (2006).
- Archetti and Scheuring (2011) M. Archetti and I. Scheuring, Evolution 65, 1140 (2011).
- Assaf et al. (2013) M. Assaf, M. Mobilia, and E. Roberts, Phys. Rev. Lett. 111, 238101 (2013).
- Li and Thirumalai (2018) X. Li and D. Thirumalai, bioRxiv p. 288670 (2018).
- Levin (2000) S. A. Levin, Ecosystems 3, 498 (2000).
- Goldenfeld and Woese (2011) N. Goldenfeld and C. Woese, Annu. Rev. Condens. Matter 2, 375 (2011).
- Roca et al. (2009) C. P. Roca, J. A. Cuesta, and A. Sánchez, Phys. Rev. E 80, 046106 (2009).
- Wu et al. (2009) Z. X. Wu, Z. Rong, and P. Holme, Phys. Rev. E 80, 036106 (2009).
- Altrock and Traulsen (2009) P. M. Altrock and A. Traulsen, Phys. Rev. E 80, 011909 (2009).
- Antal and Scheuring (2006) T. Antal and I. Scheuring, Bull. Math. Biol. 68, 1923 (2006), eprint 0509008.
- Blythe and McKane (2007) R. A. Blythe and A. J. McKane, J. Stat. Mech. p. P07018 (2007).
- Bauer and Frey (2018b) M. Bauer and E. Frey, Phys. Rev. E 97, 042307 (2018b).
- Nowak and May (1992) M. A. Nowak and R. M. May, Nature 359, 826 (1992).
- Reichenbach et al. (2007) T. Reichenbach, M. Mobilia, and E. Frey, Nature 448, 1046 (2007).
- Pigolotti et al. (2017) S. Pigolotti, M. Cencini, D. Molina, and M. A. Muñoz, J. Stat. Phys. pp. 1–30 (2017).
- Zhang and Rainey (2013) X. X. Zhang and P. B. Rainey, Evolution 67, 3161 (2013).
- Coyte et al. (2017) K. Z. Coyte, H. Tabuteau, E. A. Gaffney, K. R. Foster, and W. M. Durham, Proc. Nat. Acad. Sci. USA 114, E161 (2017).
- Peaudecerf et al. (2018) F. J. Peaudecerf, F. Bunbury, V. Bhardwaj, M. A. Bees, A. G. Smith, R. E. Goldstein, and O. A. Croze, Phys. Rev. E 97, 022411 (2018).
- Pérez-Reche et al. (2012) F. J. Pérez-Reche, S. N. Taraskin, W. Otten, M. P. Viana, L. da F. Costa, and C. A. Gilligan, Phys. Rev. Lett. 109, 098102 (2012).
- Dechesne et al. (2010) A. Dechesne, G. Wang, G. Gulez, D. Or, and B. F. Smets, Proc. Nat. Acad. Sci. USA 107, 14369 (2010).
- Gandhi et al. (2016) S. R. Gandhi, E. A. Yurtsev, K. S. Korolev, and J. Gore, Proc. Natl. Acad. Sci. USA 113, 6922 (2016).
- Limdi et al. (2017) A. Limdi, A. Perez-Escudero, A. Li, and J. Gore, bioRxiv p. 201723 (2017).
- Van Dyken et al. (2013) J. D. Van Dyken, M. J. I. Müller, K. M. L. Mack, and M. M. Desai, Curr. Biol. 23, 919 (2013).
- Pande et al. (2016) S. Pande, F. Kaftan, S. Lang, A. Svato, S. Germerodt, and C. Kost, ISME J. 10, 1413 (2016).

## Appendix A Supplemental Information

### a.1 Structure of differential equations in the well-mixed limit

We discuss the structure of the differential equations describing the well-mixed limit in more detail for the exemplary value of .

We consider an average number of individuals on a patch, even though in principle the number of individuals per patch can vary when individuals move between patches. For , which we discuss here, we then need to consider different fitness types of slow growers with fraction for , and different types of fast growers with total population fraction for . The total density of slow growers is ( for fast-growers). The fitness values of the different slow growers are given by , and the fitness values of the different fast growers are .

The temporal change for each fitness type , can then be split into a term corresponding to reproduction, , and a term for fitness adjustments, , which we explain separately in the following.

Reproduction can occur between any two individuals in the well-mixed limit and occurs proportionally to an individual’s fitness, such that part of the differential equation for reproduction, , reads

The first term corresponds to individuals of type reproducing by replacing any other individual in the well-mixed metapopulation, while the other two terms symbolise the reproduction of another species by replacement of . This type of reproduction means that , i.e. the number of individuals in the metapopulation is conserved.

Fitness adjustments only occur within a species (i.e. a slow grower can adjust its fitness to a different value, but will still remain a slow grower). Each individual adjusts its fitness in a group of players, so that one can effectively consider -player reaction terms. The fraction of individuals of the unfittest slow-growing type of fitness increases when fitness updates of a slow-grower occur in groups that contain precisely one slow grower: for this term, one thus has to consider all different combinations of players with only one slow grower and five fast growers, i.e. . The term is weighed by all different possibilities with which this combination can be drawn out of the entire metapopulation, which in this case adds up to . Similarly, there are fitness adjustment events which lead to a decrease in the fraction of this fitness type, for example when this unfittest slow grower has adjusts to an environment with one other slow grower (of any type) and four fast growers. All other terms can be calculated analogously. Thus, for an average of players per patch, the differential equations for fitness adjustment (here referred to as ) for the slow growers read

These terms should be multiplied by fitness adjustment rate . As discussed for the most unfit slow grower, there is always one positive term for each (the term with slow growers is positive), since fitness adjustment in this group would increase the number of in the system, while all other terms are negative. The binomial term in front of each term indicates how many possibilities there are to distribute the number of fast growers among the individuals, while the number before the binomial counts the possibilities to order the specific slow grower among the other slow growers.

The full differential equation for any slow grower is then

The concise form for these differential equations for general is shown in the main text of the paper. The differential equations for different fitness fractions for the fast growers can be expressed analogously. The constraint can be used to eliminate one of these differential equations.

It is clear that (more precisely ) and (more precisely ) are both fixed points of the systems, corresponding to survival of the weakest fast and strongest slow grower. We next explore whether these fixed points are stable for general .

The Jacobian of the differential equations of eq. 1 in the main text is a dimensional matrix. However, at these fixed points, the matrix is sparse: for the fast growing fixed point , the derivative of the update terms of the differential equation for fitness types with respect to all fast-growing fitness types are zero because they all involve terms of . Similarly, derivatives of the reproduction terms with respect to fitness types retain terms containing , and are thus also all zero. The eigenvalues of the Jacobian are thus the eigenvalues of the block corresponding to fast growers and slow growers separately. Analogous arguments show that even these submatrices are highly sparse, and that the diagonal elements are actually the eigenvalues of these Jacobians. All of these eigenvalues apart from two contain a negative term , which for is negative; of the remaining two eigenvalues, one is also always negative, and the other determines the stability of the fixed point and is discussed and analysed in the main text of the paper. We conclude this section by noting that this last eigenvalue corresponds to the fitness difference between the unfittest and fittest individuals of both types (depending on whether the absorbing state of all fast or all slow growers is studied). Figure S1 shows the fitness values of all different different fitness types. This figure is a reproduction of the left panel of Fig. 2c in the main text, with the only difference that here we show the full patch configuration representing the different fitness types for both species on the top and bottom -axis (slow and fast grower, respectively), for clarity.

### a.2 Impact of nonlinearity coefficient and

The phase boundary for in Fig. S2a (in black) shows that low-grower stability increases with mixing up to . For , the phase boundary in the high mobility limit is not well-defined (i.e. there is no phase transition in the system), as the system is bistable and thus depends strongly on initial conditions. The region where between and of all simulation runs yield slow growers in shown in light blue. Within this blue region, simulation runs sometimes yield to survival of slow growers or of fast growers, but not to coexistence. The average position of the phase boundary is shown in dark blue as a guide to the eyes. We note that the phase boundary in the small mixing limit changes quantitatively but not qualitatively with both and .

Figure S2b shows the boundaries for the fixed point stability and for (black, solid and long dashes) and (blue, short dashes and dashed dotted). Even for increased the coexistence region (between the lines) extends over a large range of fitness difference . This fact means that coexistence is also robust even for higher , especially when one considers that the fitness function may be even more nonlinear (decreasing ). In addition, when it comes to debating the extent of different phases, our explanation of how the different fitness types needs to balance each other makes clear that the importance of for the phase stability is that it sets the difference between the least and most fit individuals of both species. This fitness difference could (for larger populations) also be an effect observed due to private components of the public good, as discussed in the main text. Because different experimental cases would require different fitness functions and more specific modelling, it is reasonable to limit the discussion of the point of the importance of in our model to only providing evidence for the fact that coexistence also also arises and is robust for larger .