Coexistence and cooperation in structured habitats

Background Natural habitats are typically structured, imposing constraints on inhabiting populations and their interactions. Which conditions are important for coexistence of diverse communities, and how cooperative interaction stabilizes in such populations, have been important ecological and evolutionary questions. Results We investigate a minimal ecological framework of microbial population dynamics that exhibits crucial features to show coexistence: Populations repeatedly undergo cycles of separation into compartmentalized habitats and mixing with new resources. The characteristic time-scale is longer than that typical of individual growth. Using analytic approximations, averaging techniques and phase-plane methods of dynamical systems, we provide a framework for analyzing various types of microbial interactions. Population composition and population size are both dynamic variables of the model; they are found to be decoupled both in terms of time-scale and parameter dependence. We present specific results for two examples of cooperative interaction by public goods: collective antibiotics resistance, and enhanced iron-availability by pyoverdine. We find stable coexistence to be a likely outcome. Conclusions The two simple features of a long mixing time-scale and spatial compartmentalization are enough to enable coexisting strains. In particular, costly social traits are often stabilized in such an environment—and thus cooperation established.

First, we restate the dynamics already disclosed in the main text for reference: Within a single deme we model the growth dynamics of multiple strains N i (t) and the decrease of the resource concentration S(t): (S1.1a) where ∂ t is a time derivative, and different strains i are characterized by their growth rates α i and yields ϕ i . These growth rates α i (t) and yields ϕ i (t) will be separated into three terms, (S1.2b) Here, α is the average growth rate and ϕ the average yield, computed as arithmetic mean over all strains. The deviations from these average values for individual strains are denoted by δα i and δϕ i . With this definition, the average growth rate α can be used as unit of time αt, and the average yield acts as unit for substrate concentrations Sϕ, which then counts potentially growing cells. In addition, we collect all time-dependencies in the terms A(t) and Y (t). In the public good scenarios treated later, these time-dependent terms are used to couple to the public good dynamics. For this section that only treats resource consumption, we assume they are constant during the growth period, A(t) = Y (t) = 1. However, we set A(t) = 0 when resources are depleted within a deme, S(t) = 0 for t > T depl . Thus, the decreasing resource concentration due to growing populations in Eq. (S1.1b) acts as a timer to stop the dynamics. We transform all population sizes N i (t) to total population size N (t) = i N (t) and fractions X i (t) = N i (t)/N (t) of different strains, which turns Eqs. (S1.1) into ∂ t X i = αA(t) δα i X i 1 − j δα j δα i X j , (S1.3b) (S1.3c) Eq. (S1.3b) describes the frequency changes in a multispecies Lotka-Volterra model, also known as replicator dynamics, which has been studied extensively [2][3][4]. For us, however, the dynamics of the population size N is important as well, because interactions between microbial strains do not only change the frequency dynamics, but also how frequency dynamics couples to the population size dynamics. Throughout our work, we assume that different strains are seeded into new demes by independent Poisson distributions, each with parameter n i which is its mean size per deme in the pool: Averages over all droplets using this Poisson distribution are written as angular brackets, such that average population sizes N are N = n P n n N(T mix ; n) . (S1.5) These averages are always evaluated at the mixing time T mix , such that the global mixing step is represented exactly in this averaging. In some cases we also want to indicate the parameters of the seeding distribution directly: Then, we write the variables that we average over as arguments and the parameters after a vertical line, N(n) n (τ ) , which indicates that the average should be computed for the indicated variables with these parameters of the distribution. As already stated in the main text, the changes in inoculum size and composition between cycles are given by ∆n (τ +1) = d N − n (τ ) , (S1.6a) ∆x (τ +1) = ∆X + Cov X, N/ N , (S1.6b) which is used here for reference within this appendix.
In the following, we present solutions to the withindeme dynamics given by Eqs. (S1.3), insert them into this cycle mapping in Eqs. (S1.6), and then evaluate these expressions to obtain the isoclines, ∆n = 0 and ∆x i = 0.

B. Dynamics of a single strain
At first we solve the within-deme dynamics for a single strain, with constant growth rate α and constant yield ϕ, feeding on a single finite resource, and without any other environmental interactions. This is the most basic model that fits our framework, and has the advantage that we can obtain analytic solutions that allow insight about the behavior of the system. With only a single train we can drop any index i for this section.
As long as cells are growing we find N (t; n) = n exp αt as solution to Eq. (S1.1a). This expression can be used in the resource dynamics Eq. (S1.1b), which we integrate from S(0) = S 0 to S(T depl ) = 0: we obtain  The inoculum size for a single is averaged with a Poisson distribution over all possible inoculum sizes. This smears the sharp transition from starting with enough cells that deplete all nutrients to inoculum sizes that are still in the exponential growth phase at Tmix ( light colored thin lines). As explained in the text, existence of a fixed point depends on dN (Tmix; 1) > 1, and thus how many cells in the inoculum. If this condition is not met, the population will be washed out over multiple cycles.
S 0 = n ϕ exp(αT depl ) − 1 , which can be transformed to find the explicit depletion time T depl as a function of n and the other parameters. Assuming that growth ceases upon resource depletion, α(t > T depl ) = 0, the full solution for population size is In the dynamics over multiple cycles, this population sizes is evaluated at time T mix . Thus, population size becomes a function of inoculum size n only. Depending on this number resources are either used up and final population size is given by the second expression, or the population is still in its growth phase and we need to use the first expression. These time-dependent solutions are shown in Fig. S1.1A, while panel B shows the mapping for inoculum sizes over cycles, and includes these solutions as function of n as light colored lines. Explicitly, the mapping for the average inoculum size n (τ ) takes the form In these expressions, the Poisson distribution smoothly interpolates the two branches of the population size in Eq. (S1.7), which is shown in Fig. S1.1B as darker solid lines. This cycle mapping for a single strain has either one or two fixed points, but only one of them is stable: One fixed point always exists for empty inocula, n = 0, as cells cannot spontaneously grow out of nothing. If this is the only fixed point, extinction is a stable state of the dynamics. However, if a second fixed point exists, it will occur on the branch of the solution that indicates depleted resources, see Eq. (S1.7). This second fixed point with n > 0 is always stable when it exists, and extinction, n = 0, becomes an unstable fixed point upon appearance of this non-zero fixed point.
Stability of a fixed point in this one-dimensional mapping is determined by the slope of final population sizes with respect to the average inoculum size, ∂ n N n | n=n , evaluated at the fixed point n . If the absolute value of this expression is smaller than 1, the fixed point is stable, otherwise it is unstable. Before doing the calculation explicitly, however, we introduce an auxiliary step: taking a derivative with respect to the parameter of a Poisson distribution, we find which simply leads to a shift in this summation index, P n n F (n + 1) − F (n) , (S1.9) because we are only concerned with functions that exhibit the property F (0) = 0. This property just indicates that nothing grows in an empty inoculum. The expression (S1.9) generalizes straightforward to multidimensional variables, ∂ ni P n n , where the shift occurs only for strain i. With this calculation in place, the slope of the cycle mapping for a single strain at the origin, n = 0, is where we used the Poisson distribution for parameter n = 0 P n 0 = δ n0 . When nutrients are not yet depleted at mixing, the Eq. (S1.10) evaluates to d exp αT mix . Thus, large dilution rates and/or long mixing times enable a stable fixed point, which exists when d exp αT mix > 1. On the other hand, if nutrients are already depleted at mixing, the expression in Eq. (S1.10) simply evaluates to dS 0 ϕ. Thus, depending on parameters, this extinction fixed point can be either stable or unstable. Further analytic progress can be made with the assumption that mixing time is large enough that even a single cell uses up all nutrients. Then, the cycle mapping in Eq. (S1.8) only contains terms in the second sum, as n thres (T mix ) = 0. This sum can be computed, and we find the mapping From this expression, we can compute the fixed point n , where W is the Lambert-W function. For larger dS 0 ϕ, the argument in W becomes very small, and W becomes linear in its argument, W (−y exp(y)) ≈ −y exp(−y). This leads to an exponential correction for the single strain fixed point, n ≈ dS 0 ϕ 1 − exp(−dS 0 ϕ) , which can explain, for instance, why the population size isocline in Fig. 3 is shifted slightly off the value of the linear approximation in Eq. (S1.11). Stability of this non-zero fixed point can be checked by (S1.12) which is smaller than 1 for any realistic value of d, since the term involving the exponential is bounded from above by 1/e. Thus, the fixed point given by Eq. (S1.11) will be stable. Its existence can be determined by the instability of the extinction fixed point at the origin, n = 0, which occurs for dS 0 ϕ 1.

C. Dynamics of multiple strains
The main technical challenge for multiple strains is to determine the depletion time T depl (n) for arbitrary inoculum vector n. This was possible for a single strain, as shown in the previous section. In this section we will introduce the expansion factor ξ, that lets us circumvent this problem and allows progress with approximations.
First, we define this average expansion factor of a population as (S1.13) In the special case of constant growth rate, α(t) = α, this reduces to ξ(t) = exp(αt). Using ξ, the solutions to total size N and population fractions X i in Eqs. (S1.3) are with n and x the initial conditions for total population and fractions in a given deme. A time-dependence in these solutions enters only via the integration limit in ξ(t). In Eqs. (S1.14) we also see the justification for the name "expansion factor": populations expand from initially n cells to the final numbers of roughly nξ cells, with corrections due to different growth rates. One step towards obtaining solutions is to compute the expansion factor at the time of depletion, ξ(T depl ). Its derivation involves integrating resource use, Eq. (S1.3c), from S(0) = S 0 to S(T depl ) = 0, which gives For time-constant yields, Y (t) = 1, and a sufficiently large amount of resources, such that ξ 1, this relation simplifies: The integral with its explicit time-dependence drops, and T depl only enters via the final value of the expansion factor ξ(T depl ). Consequently, Eq. (S1.15) reduces to which we can approximately solve as For small differences between the strains, δα i , δϕ i 1, we have G ≈ 1 and thus the leading order is simply ξ ≈ S 0 ϕ/n.
The expansion factor ξ(T depl ) can be interpreted as a (non-linear, but monotonic) transformation of the depletion time T depl . With this viewpoint demes with an increased yield relative to the average, i x i δϕ i > 0, strains will grow longer (larger ξ), while in demes with growth faster than average, i x i δα i > 0, the time to depletion decreases (smaller ξ). The weighted growth rate differences also feature a coefficient log(S 0 ϕ/n), which reflects the fact that growth rate differences usually have a larger impact on the dynamics and can be easier (and faster) selected, compared to population level traits like yield and its differences between strains.
Since we assume that growth stops when resources are depleted, also the value of ξ stays constant after the depletion time T depl . For the single strain case, we showed that a strain that does not deplete all nutrients when seeded alone, will ultimately go extinct. Thus, the relevant case is faster depletion than mixing, T depl < T mix , such that ξ(T mix ) = ξ(T depl ), and we will use the Eq. (S1.16) as ξ(T mix ).
The expansion factor will also be used later to mediate the coupling to additional public good dynamics. Influences of varying growth growth rates are already contained in the definition (S1.13) itself. Any changes in yield can be treated in the full resource-use equation (S1.15). The two scenarios with collective antibiotic resistance and iron-extraction via pyoverdine will be examples for either of these extensions.

D. Computing isoclines
To develop approximate expressions for the isoclines, we use the solutions of the within-deme dynamics, Eq. (S1.14), and insert them into the changes over cycles, Eq. (S1.6). Using a weak selection approximation (δα i 1), we expand to first order in growth rate differences, O(δα). Specifically, this will involve the approximation ξ δαi ≈ 1 + δα i log ξ. For the change in total inoculum size, we obtain This is a general expression in terms of ξ that will also hold for public goods dynamics below. Inserting the approximation for ξ obtained specifically for resource consumption, Eq. (S1.16), is exactly the next step. the change in inoculum size is given by Thus, ∆n (τ +1) becomes independent of growth rate differences (up to first order), as all terms in δα i cancel with the dependencies contained in ξ, see Eq. (S1.13).
To obtain the isocline equation, we set ∆n = 0 in (S1.18), and find which is Eq. (10) in the main text. We see that the isocline position of the isocline is mostly determined by dS 0 ϕ. It interpolates between the single strain fixed points on the boundaries, given by Eq. (S1.11). The difference in inoculum size between the two fixed points is 2dS 0 ϕδϕ 1 ; thus smaller differences in yield make it actually more unlikely to find a stable fixed point. This can be observed directly in Fig. 4, where the coexistence regions for small δϕ 1 becomes a narrow region where δα 1 and δϕ 1 are almost proportional to each other. The composition isocline is obtained by similar methods. The two terms in its change over cycles, ∆x i = ∆X i + Cov X i , N/ N , can be expanded as with W = nξ/ nξ denoting relative deviations from expected total sizes. These general expression hold for any interaction between strains encoded in the expansion factor ξ, and we will also use them for the two public good scenarios with antibiotics and pyoverdine. After inserting ξ into the change of composition surprisingly many terms cancel up to first order of strain differences δα i and δϕ i . The resulting expressions are For two strains, they simplify to The two independent Poisson distributions for inoculum sizes n 1 and n 2 can be transformed into a Poisson distribution for the total size n = n 1 + n 2 and a binomial distribution for n 1 between 0 and n (with parameter n 1 /n). Using this, we evaluate the averages The remaining averages indicated by angular brackets only contain the total population size, while the average over the binomial distribution for x 1 has been carried out exactly. With these expressions, the composition isocline equation ∆x 1 = 0 is written as 0 ≈ δα 1 1 − 1/n log S0ϕ n + δϕ 1 1/n , from which we extract the first order approximation reported in the main text, by neglecting the logarithmic dependence on n, and assuming inoculum sizes are large enough that 1−1/n ≈ 1 and setting 1/n ≈ 1/n. Numerical observations suggest that this scaling also contains a coefficient close to 2, due to these neglected terms. However, we are mostly interested in how this position of the composition isocline changes upon variation of the different parameters in the model. A dependence on composition x (τ ) 1 will only appear in higher orders of δα 1 and δϕ 1 , which is consistent with the very straight (orange) lines depicted in Fig. 3 of the main text. We also see that decreasing growth rate differences will push this isocline to larger inoculum sizes. This strong dependence on δα i of the composition isocline is in stark contrast to the position of the population size isocline, which is independent of δα i to first order. The dilution rate d plays no role in the position of the composition isocline, as it cancels already out when calculating the ratio in the cycle mapping, Eq. (S1.6b).

E. Coexistence of multiple strains
This section presents an argument on how multiple strains can coexist in the dynamics. These considerations benefit from a description in absolute numbers n = (n 1 , n 2 , . . . ) for each of the strains, and not use the non-linear transformation of computing fractions x i = n i /(n 1 + n 2 + . . . ).
One way to define coexistence in our model for multiple strains is when all single strain fixed points are unstable with respect to invasion by every other strain. This instability of all single strain fixed points implies that the system will have an additional coexistence fixed point. If only a subset of strains fulfill these pairwise condition, such a subset can coexist over cycles. The example for two strains shown in Fig. 3 depicts such cases.
Formally, this occurs when the Jacobian evaluated at these single strain fixed points features eigenvalues larger than 1. Entries of the Jacobian J at these points n i ≡ (0, . . . , n i , . . . , 0) can be computed by For the following, we assume that the resident strain has index 1. We already computed the first entry, J 11 = ∂ n1 N 1 (n 1 , n 2 . . . ) (n 1 , 0 . . . ) < 1, in Eq. (S1.12), which just indicates that the resident strain is stable on its own. Other entries in the diagonal indicate an invading strain growing when each deme is seeded with exactly one cell of this invading strain, for i ≥ 2, where we used the auxiliary calculations in Eqs. (S1.9) and (S1.10). Apart from entries on the diagonal, the Jacobian is a very sparse matrix: Only the first row contains other non-zero entries, which are final population sizes of the resident strain 1 upon invasion of each other strains (with one cell of the invading strain in each deme). Every other entry indicates the change in final population sizes for the other strains: These are all zero and do not change -the Jacobian is evaluated at the single strain fixed point of the resident strain. Thus, the structure of the Jacobian implies that all entries on the diagonal, Eq. (S1.24), are already the eigenvalues.
Assuming the system has k strains, each of the k single strain fixed points needs to be unstable with respect to invasion to all (k − 1) other strains, leading to k(k − 1) (pairwise) conditions of the form J ii | n j > 1.
The coexistence regions depicted in Fig. 4 are obtained by numerically evaluating these conditions for two strains.
This set of conditions in Eq. (S1.24) is not constructive: it not indicate how to choose parameters to get coexistence. Only when a set of strains is given by their growth parameters α i and ϕ i , one can check if such a set can coexist. Numerical evaluations indicate the the possible range of values for these two parameters decreases with each additional strain. For example in [1], where a related spatio-temporal model was studied, values were chosen with fine-tuned relations along α and ϕ to find coexistence between multiple strains.

S1.2. GROWTH DYNAMICS WITH ANTIBIOTICS
In the presence of antibiotics an additional equation for the dynamics of its concentration is coupled to the microbial strain growth and resource consumption dynamics (see main text). This section presents the derivation of the expansion factor ξ for these dynamical interactions through collective antibiotic reduction by extracellular enzymes. With ξ computed, we can use the general expressions for the changes in the inoculum, Eqs. (S1.17) and (S1.20). The additional steps needed to derive the approximations for isoclines reported in the main text are found below.

A. From dynamical interactions within demes to expansion factors
In the presence of antibiotics, the growth rate of strains is given by Thus, growth rate is a sigmoidal, decreasing curve for increasing antibiotic concentrations B(t). In this expression, the antibiotic concentration is always measured in multiples of the minimal inhibitory concentration (MIC) µ, which is the concentration where the population switches from overall death to overall growth, such that α(B = µ) = 0.
In the within-deme dynamics, the antibiotic concentration B(t) is reduced by the different microbial strains with rates ρ i . Thus, the set of relevant equations is given by For the moment, we ignore the resource dynamics S(t), and ask if the populations in Eq. (S1.26) will survive or go extinct? In a different manuscript, we show that our ensuing approximations are actually very robust against any changes of molecular details, and the results hold for several different dynamics of antibiotic resistance (manuscript in preparation).
In order to solve the combined dynamics of microbial death and antibiotic reduction, the first step involves defining a logarithmic antibiotic concentration K = log(B/µ). Often, only the magnitude of the concentration is important, and not its exact value. Using this logarithmic concentration, we can approximate α i (B) ≈ αiκγ 1+γ log(B/µ) = −α i λK, where we defined the compound parameter λ = κγ/(1 + γ). This allows to write the dynamics in Eqs. (S1.26) as ∂ t N i = −α i λKN i and ∂ t K = − i ρ i N i . We assume that population sizes decay with a rate determined by the initial antibiotic concentration, N i (t) = n i exp − α i λK 0 t , and this death rate does no change significantly during the relevant time. Numerical integration shows that this approximation is sufficiently close to the actual solution. Thus, with this exponential decay in population size we can integrate the dynamics for K, which yields Now we have both trajectories, N (t) and K(t), such that we define a time T 1 implicitly via N (T 1 ) = 1 and thus indicates when the exponential decay of the population reaches a single cell, 1 = n i exp(−α i λK 0 T 1 ). Solving for t 1 and inserting this expression into the dynamics of the logarithmic antibiotic concentration, K(T 1 ), we can check if antibiotics have been reduced sufficiently, K(T 1 ) < 0 (population survived), or still exceeds this value K(T 1 ) > 0 (population will go extinct): For large enough inoculum sizes, 1/n i 1, we can rearrange this inequality to find (S1.27) This condition in Eq. (S1.27) indicates that the microbial population reduces the antibiotic concentration fast enough before going extinct. If the population survives, the growing populations it will degrade any remaining antibiotic fast, and B(t)/µ drops from 1 to 0. We assume that the mixing time T mix is long enough, such that any surviving population will use up all remaining nutrients. The number of dying cells during this initial phase of population death is likely small compared to final population sizes. Thus, these final sizes will be close to the population sizes obtained from the simple resource consumption model, N ≈ S 0 ϕ. Consequently, the expansion factor can be approximated as (S1.28) The Heaviside-Theta function Θ in this expression indicates the survival condition (S1.27) above: the expansion factor evaluates to zero, if all cells die, and is given by the first order term S 0 ϕ/n of Eq. (S1.16) when they survive.

B. Isoclines of the cycle mapping
In order to fully understand the isoclines for antibiotic reduction in Fig. 5, it helps to investigate the dynamics of single strain first. Fig. S1.2 illustrates this map of inoculum sizes, using a single producing strain. The depicted sigmoidal curve emerges from smearing out the sharp cut-off between surviving and extinct populations in the expansion factor ξ, Eq. (S1.28), due to the Poisson distribution in the seeding step. This shape implies that either one or three fixed points occur in the mapping. Similar to the resource consumption dynamics, extinction is always one of these fixed points, n = 0. In contrast to before, however, it is usually stable. The other two fixed points are close to the survival threshold, n 1 ≈ (α 1 λ/ρ) log(B 0 /µ) 2 , and the already known fixed point where all resources are used up, n 1 ≈ dS 0 ϕ 1 . The latter fixed point is always stable, if it exists, and the former is always unstable and acts as barrier between the basins of attraction to either extinction or survival, see Fig. S1.2B. These two fixed points only exists, if their expressions assume values such that the survival threshold is smaller than the fixed point where resources are used up. Otherwise, extinction is the only possible long-term outcome in the dynamics.
For two strains, the two single strain fixed points indicate the intersections of the population size isocline with the boundary at x 1 = 1, see Fig. 5. From the stable fixed point (for resource depletion, with larger inoculum sizes) the isocline extends down to lower x 1 values along the condition derived before for resource consumption, Eq. (S1.19). Without any differences in yield, δϕ 1 = 0, it is a straight line. The other condition, limited by the antibiotic threat, extends along constant producer inoculum sizes, n 1 = const, and is given by the expression computed above for survival of the population. Thus, this part of the isocline will scale with n ∼ x −1 1 , where the proportionality factor is determined by the threshold population size. Both parts of the isocline will meet, when x 1 becomes small enough, and the isoclines connects these two descriptions.
The population composition isocline, 0 = ∆x 1 , is determined by the balance of the local losses in frequency and the coupling to population size changes, 0 = ∆X 1 + Cov X 1 , N/ N . As we assume δα 1 < 0 for the producing strain, the first term is usually negative. Inspecting the expansion in growth rate differences, see Eq. (S1.20), we observe that only a single term does not depend on δα 1 -contrary to the resource consumption before, factors of δα i and δϕ i in ξ are negligible for antibiotic reduction. This term, x 1 (W − 1) , where W = nξ/ nξ , is responsible for large values in the covariance, that lead to positive values in the overall change ∆x 1 , required for coexistence. Inserting ξ for antibiotics the condition for the isocline, 0 = ∆x 1 , can be written as where Θ abbreviates the full expression of the threshold in ξ. The first bracket indicates the difference between the average inoculum fraction of the producing strain in surviving populations and its average fraction in the inoculum. Such an expression is expected to be positive, as surviving populations likely have a larger fraction of the producing strain. However, with increasing inoculum sizes this term will become smaller, as then more populations survive, and the two averages become closer.
When evaluating this first term in Eq. (S1.29) numerically, we find two regimes: For small values of n 1 , it is an exponential decay in n 1 , until n 1 reaches its threshold value in the argument of Θ. For n 1 larger than this threshold, the difference decays algebraically with an exponent close to −1. In contrast, the dependence on the non-producer inoculum size n 2 is weak. Thus, n 1 already specifies the shape of x 1 (W − 1) across the whole phase plane, given by x 1 Θ / Θ − x 1 ∼ 1/n 1 in the relevant region. For the isocline, the decay of this (positive) expression needs to be balanced by all the other terms, which are of order O(δα 1 ); this coefficient δα 1 already determines the magnitude, other dependencies only add minor variations to this value. Hence, for the isocline we seek curves where the first term in Eq. (S1.29) is small and almost constant, which occurs along n 1 ≈ const, which is the main dependency of the first term. Thus, we approximate the scaling of the isocline as (S1.30) which holds as long as enough cells of both strains are in the inoculum. Here, we introduced the small positive constant 1, which we only observe numerically. It is needed, however, in order to explain that the composition isocline is typically steeper than the inoculum size isocline in phase plane, which also scales as ∼ 1/x 1 .

C. Stability of the coexistence fixed point
The two curved isoclines will in general have multiple intersections, which are all fixed points of the cycle mapping. However, we find that at most two of them are stable: A population of only producing strains can survive, as shown for the single strain case above. Such a fixed point appears on the boundary of the phase space, x 1 = 1.
Coexistence can occur, if the two isoclines intersect in a specific way: on the inoculum size isocline, the intersection has to be on the (vertical) part described by resource limitation. On the composition isocline, the intersection needs to occur on the part following the extinction threshold. Only for an intersection on these parts of the isocline, we find a stable fixed point. In Fig. 6 we showed how these isoclines change when changing either the dilution rate d or the growth rate differences δα 1 . Fig. S1.3 extends these results and shows the stability of this potential coexistence fixed point for a much larger range of parameters. The color encodes the result of a linear stability analysis, obtained Stability of fixed points within-deme dynamics involving antibiotic hazards. Displayed are color-coded stability properties of the coexistence fixed point, as long as it exists (gray symbols indicate it does not). Specifically, the colors indicate stable (•), stable with complex eigenvalues (•), unstable with complex eigenvalues (•) and unstable (•) fixed points. In the complex unstable regime, close to the boundary to stability, we also find stable limit cycles, which are not captured in the linear stability analysis of computing eigenvalues of the Jacobian. These limit cycles appear through a Neimark-Sacker bifurcation [5]. The change in shape and position for the two isoclines of inoculum size and inoculum composition are shown in more detail in Fig. 6, where purple and brown boxes in panel (C) correspond to the similar colored borders.
by computing eigenvalues of the cycle map evaluated at the fixed point: for constant dilution rate d and increasing the growth rate difference δα 1 the fixed point goes from stable (green solid), to stable complex (purple solid), to unstable complex (purple circle), and finally to unstable (red circle). Clearly, the actual position in phase plane of these fixed point will change as well, but this information is not displayed in Fig. S1.3. The purple and brown boxes enclose the parameter range shown in more detail already in the main text in Fig. 6. Numerically, we also found stable limit cycles in the unstable complex regime, close to the transition between stable and unstable, from a Neimark-Sacker bifurcation [5]. This transition occurs when the intersection of the two isoclines moves on the inoculum size isocline from the description by resource limitation to the description by antibiotic reduction (at its smallest x 1 values).

D. Additional parameter dependence of the cycle mapping
In the main text, we focused on how isoclines change upon varying the two parameters δα 1 and d. Each of them changes one of the two isoclines (population size or population composition, respectively), while leaving the other isocline almost unaffected. However, the models for antibiotic reduction also features additional parameters, which we explore here.
The dynamics of antibiotic reduction clearly depends on the antibiotic concentration B 0 /µ present at the beginning of the growth phase, and how fast antibiotics is reduced. From the threshold between growth/nogrowth, as stated in Eq. (S1.27), we see that both ρ i /α i λ and B 0 /µ give additional important parameter combinations. Fig. S1.4 illustrates the effects of changing either of those. In general, as long as B 0 /µ > 1, the effect of changing either production rate or initial antibiotic concentration only changes the threshold inoculum size for producing strains that would survive. Changing this threshold only shifts the position of the isoclines, but does not significantly change their shapes. However, when B 0 /µ < 1 the non-producing strain will be able to survive on its own, although it might grow significantly slower in demes where antibiotic concentration is not reduced. Due to this survival of the non-producing strain, the part of the composition isocline that indicated the extinction of the non-producing strain (asymptotic approach towards the x 1 = 1 axis) disappears, and the composition isocline directly intersects with the pure-producer axis of the phase plane, as shown in Fig. S1.4AB.

S1.3. GROWTH DYNAMICS WITH PYOVERDINE
Similar to the public good dynamics with antibiotics, we compute here an expansion factor ξ that takes the effect of pyoverdine production into account. This allows to use the same approximations employed in previous sections to derive the isoclines. Here, we provide the steps in the derivation of these approximations. In contrast to before, however, we mostly focus on the effect of the public good, and neglect several corrections due to the strain differences in δα i and δϕ i .

A. From dynamical interactions within demes to expansion factors
Pyoverdine changes yield ϕ on a population level due to enhanced availability of iron. We assumed its effect is given by, (S1.32) where the latter term in brackets is the time-dependent Y (t) in Eq. (S1.2b). The production of pyoverdine P is given by ∂ t P = i ρ i N i . In order to compute the effects on the expansion factor ξ, we need to evaluate the full resource-use Eq. (S1.15), which contains an integral over these time-dependent changes in yield, dtS(t)(∂ t Y )(t).
To this end, we combine the two ODEs for S(t) and P (t), which are both linear in N (t). First, we rearrange both to bring N (t) to one side of the respective equation, and then equate both expressions, which leads to . (S1.33) In anticipation of the result we wrote the dynamics of the pyoverdine concentration P together with yield ϕY (P ), which is influenced by this concentration. The population composition X i (t) changes only slowly in time, thus we assume that X i (t) ≈ x i for the purpose of the ensuing approximations, and collect all these dependencies in the term R, which we assume to be roughly constant over time. Then, we can integrate this equation from initial conditions S(0) = S 0 and P (0) = 0 up to a time t to get the relation between S(t) and P (t), which is an almost linear expression of the form S(P ) ≈ S 0 − R ϕσ P . As a next step, we take a derivative with respect to time on ϕ(P ), Eq. (S1.32), to obtain ∂ t ϕ(t) = ϕ(σ − 1)(∂ t P ) exp −P (t) . (S1.35) With these two expressions, we transform the integration variable, dt(∂ t P ) = dP , together with adjusting the integration bounds, in the full resource-use dynamics, Eq. (S1.15). Then, we can evaluate this equation to obtain The first line in this solution indicates the expected yield at the end of growth phase. The expression of the second line denotes an correction due to the fact that the concentration of P is increasing during this growth phase, and did not yet start at its final value.
Having this expression, we can simplify further, and neglect the second line to obtain the main scaling of the expansion factor. To this end, we define the abbreviations ξ 0 = S 0 ϕσ/n and H = i xiξ δα i 1+δϕi . The full resource use equation can then be written as which we need to solve for ξ. Note that R also contains a logarithmic dependency on ξ, which does not influence the overall scaling. The term H generates the corrections in δα i and δϕ i , leading to the term G reported in Eq. (S1.16). this generates only multiplicative factors, which is the reason that we solve for ξH first. We obtain For small enough arguments, the negative branch of the Lambert-W function W is linear in this argument, which is the case here: nξ 0 /R is usually large. Thus, we can simply leave out W to obtain a scaling of the expansion factor for pyoverdine production, (S1.36) which assumes H ≈ 1. In order to include larger deviations from unity in H, we will obtain terms similar to the corrections G, as in Eq. (S1.16) For the following, we define A(x) = S 0 ϕσ i x i ρ i /α i , which is the full argument in the exponential.

B. Computing the shape of isoclines
The approximation for the expansion factor ξ in Eq. (S1.36) can be inserted in the general expressions for the population size in Eq. (S1.17) and for the population composition in Eq. (S1.20). At first, we observe that saturation of the public goods dynamics occurs as long as A(x) 1. Then, the exponential term in Eq. (S1.36) is small and negligible and we have ξ ≈ S 0 ϕσ/n. Accordingly, the single strain fixed point of a producer strain is n 1 ≈ dS 0 ϕσ. Since we assume that pyoverdine is not essential for growth, we also find a fixed point for the non-producer strain at n 2 ≈ dS 0 ϕ. For only two strains, the population size isocline then connects these two single strain fixed points: From the producer fixed point it extends to lower producer fractions along a line given by Eq. (S1.19), with the replacement ϕ → ϕσ to indicate saturated public good effects. At a fraction given by A(x 1 ) ≈ O(10 0 ), saturation ceases and final population sizes drop to their base level in absence of pyoverdine. Illustrations in Fig. 7 clearly show this behavior, where the (blue) population size isocline exhibits a sharp bend at small producer fractions.
For the composition isocline, we insert the expression for ξ into Eq. (S1.20). Similar to the case with antibiotics, the main term contributing to a large positive covariance term is x 1 (W − 1) , which evaluates to . (S1.37) Interestingly, when numerically computing these averages over the Poisson distribution, we find that this expression decays exponentially close to exp(−n 1 ), when A(1) 1. The factor σ only plays a role in the coefficient of this scaling. If A(1) 1, the dynamics effectively describes the original resource consumption, as then the pyoverdine production does not influence the dynamics. For intermediate values, A(1) ≈ O(10 0 ), this expression starts to decay much slower. We focus on the simplest case with a significant effect of the public good, where we can use a similar reasoning as with antibiotic reduction: Eq. (S1.37) needs to be balanced by terms of order O(δα 1 ), and thus we again seek curves in phase space where this expression is small and almost constant. Thus, with the exponential decay ∼ exp(−n 1 ), such curves are given by n 1 = const, where the position is strongly determined by the growth rate difference δα 1 . For smaller δα 1 the condition 0 = ∆x 1 is met further out the tail of the exponential decay. Approximately, this balance is given by exp(−n 1 ) ∼ δα 1 . Thus, the scaling of the composition isocline can be approximated by ∆x 1 = 0 ⇔ n ∼ log |δα 1 | + C x 1 , (S1.38) with a logarithmic dependence of δα 1 on its position, and C the coefficient between terms in the balance.  Figure S1.5. Influence of other parameters on the isoclines in the dynamics with pyoverdine. Panels show the isoclines for (A) average inoculum size n and (B) average inoculum composition x1, upon variation of the factor σ that increases yield due to presence of pyoverdine. In panel (A) we additionally illustrate the influence of the production rate ρ1 as different shades of blue. This production rate does not significantly change the isoclines in (B). Parameters dS0ϕ = 10, δα1 = −10 −2 are constant for all plots.

C. Other parameter dependence
The main text only explored the effects of changing growth rate differences δα 1 and the dilution rate d, which are the main parameters for the position of isoclines. Fig. S1.5 explores the effects on isoclines if either of the production rate ρ 1 or the increase in yield σ.

S1.4. NUMERICAL METHODS
Several results in our manuscript have only been obtained numerically, although we tried to explain these results with scaling arguments. Here, we briefly present numerical methods used throughout our work.
In general, the two levels of selection of growth within demes and the mixing cycles are also reflected in our code. We first solve the within-deme dynamics, on a grid of all possible initial conditions (n 1 , n 2 ). These initial conditions are chosen up to values n i,max , such that the probability of finding this combination in the Poisson distribution are tiny and negligible. Often, we use the condition P n i,max n i /P n i n i 10 −5 to determine these maximal inoculum sizes. For most of the results shown in Figs. 5 and 7 this maximal inoculum size is roughly n i,max = 200, which is the value we use in our figures. Solutions to within-deme dynamics are obtained by a standard Runge-Kutta 4th order integration scheme up to the time T mix . After final population sizes N i (T mix ; n 1 , n 2 ) are computed, they are stored as lookup-table. Cycles of growth, mixing and reseeding are realized by computing the sums ni,max n1,n2=0 P n 1 n 1 P n 2 n 2 N i (T mix ; n 1 , n 2 ) .
Fixed points of the mapping (e.g. Fig. S1.3) are obtained via iteration of a multi-dimensional Newton-Raphson scheme, which requires usually only a few steps.