Controlling evolutionary dynamics to optimize microbial bioremediation

Abstract Some microbes have a fascinating ability to degrade compounds that are toxic for humans in a process called bioremediation. Although these traits help microbes survive the toxins, carrying them can be costly if the benefit of detoxification is shared by all surrounding microbes, whether they detoxify or not. Detoxification can thereby be seen as a public goods game, where nondegrading mutants can sweep through the population and collapse bioremediation. Here, we constructed an evolutionary game theoretical model to optimize bioremediation in a chemostat initially containing “cooperating” (detoxifying) microbes. We consider two types of mutants: “cheaters” that do not detoxify, and mutants that become resistant to the toxin through private mechanisms that do not benefit others. By manipulating the concentration and flow rate of a toxin into the chemostat, we identified conditions where cooperators can exclude cheaters that differ in their private resistance. However, eventually, cheaters are bound to invade. To overcome this inevitable outcome and maximize detoxification efficiency, cooperators can be periodically reinoculated into the population. Our study investigates the outcome of an evolutionary game combining both public and private goods and demonstrates how environmental parameters can be used to control evolutionary dynamics in practical applications.


SHIBASAKI And MITRI
Despite these exciting advances, these approaches remain only partially successful. In medical microbiology, the emergence of resistant pathogens has led to a life-threatening global health crisis (Tacconelli et al., 2018). On the engineering side, we lack sufficient understanding to maximize the benefits gained from microbes, in terms of production rates or the efficiency of degradation of toxic compounds (Giri et al., 2020).
In all the examples above, there are two obstacles for controlling microbes. First, microbes' response to changes in their environment and to the behavior of neighboring cells may cause the extinction of species that contribute to community function or that protect against pathogens (ecological instability). Second, their large population size and short generation times mean that microbes quickly acquire mutations that can lead to evolutionary instability. Species can thereby gain resistance to antibiotics or lose their ability to perform a desired community function (Akita & Kamo, 2015;Bull & Barrick, 2017;Kumar, Maschke, Friehs, & Schügerl, 1991;Rugbjerg, Myling-Petersen, Porse, Sarup-Lytzen, & Sommer, 2018;Rugbjerg, Sarup-Lytzen, Nagy, & Sommer, 2018). In the bioremediation of heavy metals, for example, mutants that do not degrade these harmful compounds may invade the original population (Ellis, Lilley, Lacey, Murrell, & Godfray, 2007;O'Brien, Hodgson, & Buckling, 2014) and exclude it. To optimize the functional efficiency of a microbial system, therefore, we need to consider both ecological and evolutionary dynamics of microbial populations (Schuster et al., 2010).
In this study, we focus in on a bioremediation problem and develop a mathematical model to investigate the simple case of a single species degrading a toxic compound in a chemostat. We first explore its ecological and evolutionary stability and, second, use this knowledge to optimize the efficiency of toxin degradation over time. In our model, we assume that toxins are harmful to the microbes (e.g., heavy metals), who degrade them by secreting a costly product (e.g., extracellular enzymes). The benefit of degrading toxins is shared by the whole population as, for simplicity, our model does not include spatial structure. Toxin degradation can therefore be regarded as a public goods game (Broom, Pattni, & Rychtář, 2018;Samuelson, 1954), as defined in microbiology (Hummert et al., 2014;Smith & Schuster, 2019;West, Diggle, Buckling, Gardner, & Griffin, 2007). It is important to note that not all bioremediation systems correspond to public goods games (e.g., Röling et al., 2002;Smith, Graham, & Cleland, 1998), but here we focus on a subset of these systems where the compound to be degraded is toxic to the microbes and its degradation is costly. In such a system, the evolutionary instability of bioremediation is expected (Ellis et al., 2007;O'Brien et al., 2014).
In our model, microbes can adopt one of four strategies: They can be product secretors that pay a cost to contribute to the public good (cooperators) or nonsecretors that do not (cheaters) (see O'Brien et al. (2014) for an empirical example). In addition, microbes can be sensitive to the toxins or can acquire resistance, for example, by activating efflux pumps to expel toxins from within the cell (Blair, Webber, Baylay, Ogbolu, & Piddock, 2015;Bottery, Wood, & Brockhurst, 2016;Rojo-Molinero, Macià, & Oliver, 2019), or thickening the cell wall. In essence, public good secretion can also be seen as a form of extracellular resistance to the toxin. In other words, here we consider toxin resistance through private or public means, whereby a cell benefits only itself or also the remaining population, respectively.
The population and evolutionary dynamics are then analyzed using evolutionary game theory, where a strategy is considered to be evolutionarily stable if it is not invaded by mutants with another strategy (Maynard Smith & Price, 1973). Evolutionary game theory typically considers the frequencies of strategies (i.e., frequencies sum to one), for example, in the replicator dynamics (Cressman & Tao, 2014). Here, however, since toxin concentration decreases with the absolute number or density of degrader microbes (cooperators), our model describes the dynamics of the densities of strategies (as in Hauert, Holmes, and Doebeli (2006);Hauert, Wakano, and Doebeli (2008); Gokhale & Hauert, 2016). And since the microbes' death rate depends on toxin concentration in the environment and their resistance level, our model also includes environmental feedback (as in Gong, Gao, and Cao ( 2018); Tilman, Plotkin, and Akçay (2020); Weitz, Eksin, Paarporn, Brown, and Ratcliff (2016)), where each strategy affects the environment differently, and the changing environment affects the fitness of each strategy differently. We use this model to derive a protocol for optimal toxin degradation. We first show that cooperators that secrete toxin-degrading enzymes are excluded by cheaters that do not, if they have the same level of resistance to the toxin. This recapitulates a well-known result that can be explained by the tragedy of the commons (Hardin, 1968).
However, we then show that cooperators can invade a population of cheaters if their level of toxin resistance is different. Since we assume that cheaters are unlikely to quickly acquire double mutations leading to cooperators with a different resistance level, maintaining degradation is only possible if we periodically inoculate these cooperators back into the chemostat. The success of this approach relies on the ability of cooperators to invade cheaters of different resistance. We then calculate the values of the experimentally controllable parameters (inoculation probabilities of cooperators, chemostat dilution rates, and inflowing toxin concentrations) that maximize the cumulative efficiency of detoxification.
In sum, our model combines population dynamics, evolutionary dynamics, and environmental feedback to optimize a population-level function. Integrating ecology and evolution into microbial public goods games is increasingly appreciated in microbial applications (Moreno-Fenoll, Cavaliere, Martínez-García, & Poyatos, 2017;Sanchez & Gore, 2013). And while optimizing bioremediation is the case study we are considering here, our approach of controlling evolutionary dynamics by changing environmental parameters can be applied to many other microbial functions.

| MODEL
In our scenario (Figure 1a), cells can take on one of four strategies depending on whether they produce the enzymes that degrade the toxin (cooperate) or not (cheat), and whether the cells are resistant to the toxin (resistant) or not (sensitive): sensitive cooperator (sCo), sensitive cheater (sCh), resistant cooperator (rCo), and resistant cheater (rCh).
We begin by defining the bacterial population dynamics in our system. The dynamics of the density x of each strategy i ∈ {sCo, sCh, rCo, rCh} in a chemostat are defined by growth, death, and dilution out of the chemostat: [Correction added on 13 August 2020, after first online publication: The equation 1 has been corrected.] where r i is the intrinsic growth rate of strategy i due to nutrients that are not explicitly defined in the model, δ i (T) is the death rate of strategy i given toxin concentration T, and α is the dilution rate. The total densities of the four strategies ∑ i x i should be lower or equal to one in Equation (1); that is, the carrying capacity is equal to one in the absence of death or dilution. In this formulation, a useful proxy for fitness is the ratio between intrinsic growth and death at a given toxin concentration, whether death is by toxin or by dilution: should be satisfied for any strategy i that exists in the chemostat (x i > 0). In addition, when (1) where r is the maximum intrinsic growth rate.
Cellular death rate δ i (T) increases with toxin concentration T, and is represented by a Hill equation as is common in models of death by drugs (Chou, 2006): where d max is the maximum death rate, K i is the half maximal toxin concentration of strategy i, and n is the Hill coefficient, which determines the steepness of the function. Resistance can be modeled either by increasing K i or decreasing n (Sampah, Shen, Jilek, & Siliciano, 2011). Here, we assume that resistant cells have a larger K than sensitive cells K r > K s , such that they reach d max at a higher toxin concentration than the sensitive cells. Note that the toxin concentration T changes over time, as described below.
Due to the dilution in the chemostat, a proportion of cells of each strategy i flows out of the chemostat. The dilution rate into and out of the system is denoted by α To describe the population dynamics of each strategy, however, it is necessary to also formulate the dynamics of the toxin concentration because it affects the microbes' death rate, and because the toxin concentration changes over time as cooperators detoxify it.
The dynamics of the toxin concentration T in the chemostat are defined by the concentration flowing into and out of the chemostat, and detoxification by cooperators: where T in is the toxin concentration flowing into the chemostat, and f(x Co ) is the degradation rate, which is assumed to follow a

Michaelis-Menten function:
where x Co = x sCo + x rCo , that is, the sum of sensitive and resistant cooperators. Whether cooperators are resistant or sensitive has no impact on toxin degradation. K d represents the density of cooperators x Co that gives half the maximum of f(x Co ). All parameters of the model are listed in Table 1 (3b)

TA B L E 1 List of parameters
As the goal of this study is to maximize the efficiency of detoxification, we define detoxification efficiencyϕ as the difference between the toxin concentration flowing into and out of the chemostat multiplied by the dilution rate: With this definition, ϕ is proportional to the amount of detoxified liquid and is composed of the degree of detoxification and the amount of liquid flowing out of the chemostat. Although this equation gives the efficiency at any time t, we mainly focus on the efficiency at an equilibrium.

| All strategies can persist in mono-cultures with no mutation
We first analyze whether cooperators and cheaters can persist in mono-culture. Remember that only cooperators produce public goods that degrade the toxin and thereby increase the survival of all cells in the chemostat. When a few cooperators (either sensitive (a, f) Diagrams of pairwise invasion analysis when a single mutation (a) or double mutations (f) occur. A → B represents that A is invaded by B, and the color of each arrow shows the condition for successful invasion. Black arrows represent successful invasion regardless of the toxin concentration, while pink and blue arrows represent toxin concentration-dependent invasion (that invasion succeeds when the toxin concentration is low or high, or when the toxin concentrations are intermediate, respectively). (b-e, g-j) Examples of an invasion state space for each pair of strategies. A → B represents that A is a resident strategy and B is an invader. Pairs of sCo and sCh, and rCo and rCh are omitted since cheaters are fitter than cooperators in these pairs, regardless of the parameter values. In the yellow areas, the resident strategies are not invaded by the invaders, while the invasion succeeds in the orange areas. Each solid line is a boundary under which the resident strategy persists in mono-culture. The dashed line in each panel represents where the fitness proxy W of resident and invader strategies are equal at an equilibrium reached in a mono-culture of the resident strategy. Note that residents can coexist with invaders under certain conditions (see Appendix S3). Parameter values in (b-e, g-j) are f max = 0.5, K d = 0.2, r = 1, c d = 0.15, c r = 0.2, d max = 1, K s = 0.3, K r = 0.6, and n = 3. Note that panels (b-e, g-j) are examples of the state space given the parameter values; different parameter values will show different invasion landscapes where i is the focal strategy (see Appendix S1 for derivation). By solving dT∕dt = 0, dx i ∕dt = 0, one can find a trivial equilibrium (x i = 0) and one or more nontrivial equilibria T * , x * i > 0 that should satisfy: which we can calculate numerically using Newton's method.
In the absence of cooperators, we assume that the toxin concentration is equal to the incoming toxin because cooperators are the only degrader cells. In the case of a mono-culture of cheaters then, T = T in regardless of cell density, and their density at a stable equilibrium x * i , i ∈ {sCh, rCh} is positive and given by Equation (9b) if inequality (8) holds (see Appendix S1 for details). From now on, we focus on conditions where cell density converges to a positive value in mono-culture (i.e., where inequality (8) holds).

| Cheaters can invade a population of cooperators
Next, we ask what happens when a cheater mutant invades a population of cooperators at its nontrivial equilibrium or vice versa.
Cheater mutants can invade and exclude a population of cooperators at any toxin concentration and independently of whether they are both sensitive or both resistant, as long as the resident and mutant have the same resistance levels (e.g., sCo and sCh) (Figures 1c, and 2a). In contrast, cooperators are unable to invade a population of cheaters of the same resistance level at any toxin concentration. These findings recapitulate the classical result that cheaters will always dominate in a well-mixed environment  because cooperators pay a cost for producing degrading enzymes, but are as sensitive to the toxin as cheaters. In other words, the tragedy of commons (Hardin, 1968) occurs in this case.

| Toxin concentration determines invasion of sensitive and resistant cells
For mutants that differ in their private resistance level (e.g., sCo and rCo), the toxin concentration determines whether invasion succeeds or not (Figure 2b-e). Intuitively, this is because the benefit of being resistant to the toxin is quite low when its concentration is very low. Similarly, when toxin concentration is very high, the death rate of the resistant strain is close to that of the sensitive strain and too high to compensate for the cost of private resistance. Under these conditions, sensitive cells can invade a population of resistant ones ( Figure 1D). Instead, resistant cells can invade a population of sensitive cells when the toxin concentration in the chemostat is intermediate. Two strategies that differ only in their resistance level never stably coexist (see Appendix S3 for derivation).

| Cooperators can invade and coexist with cheaters of different resistance levels
Thus far, we have considered whether mutants can invade a resident population that differs in only one trait, their private or their public resistance (i.e., cooperative toxin degradation). While we assume that the time to reach an equilibrium following a single mutant invasion is shorter than the time for a second mutation to occur, we nevertheless explore invasions by such double mutants here (Figure 2f). Once a mutant has invaded, whether it will coexist with the resident population is unclear because, as cooperators increase, toxin concentration decreases, which changes the fitness landscape. In other words, increasing cooperator density can decrease the fitness difference between cooperators and cheaters. In Appendix S3, we show that cooperators and cheaters of different resistance levels (e.g., rCo and sCh) can indeed stably coexist at certain parameter ranges. Nevertheless, these two coexisting strains can then be invaded by cheaters with different resistance (e.g., rCh), which excludes the other two strategies (see Appendix S4). In other words, cooperators are never evolutionarily stable because they can be invaded and excluded by cheaters of the same resistance level, as we show next.

| In the long-term, cooperators are unlikely to be maintained
Having analyzed the outcomes of the invasion of all mutants in the short-term (Figure 2a,f), we can now predict how the population in the chemostat will change in the long-term. Regardless of which type we start with, as cells mutate, the population will transition between different genotypic "states," which can be represented by the state transition diagram in Figure 3. The probability of cooperators mutating into cheaters or vice versa is given by μ 1 , while μ 2 is the probability to change the level of resistance. For simplicity, we assume that these mutation probabilities are very small. Once a mutation occurs, we use Equation (1) as outlined above to take us to the following equilibrium state. We assume that no further mutations will occur before the equilibrium is reached, but relaxing this assumption does not alter the overall dynamics (Appendix S7). For example, sCh will appear in the population of sCo with probability μ 1 (1-μ 2 ) and exclude it. Then, rCh can appear in the population of sCh with probability (1-μ 1 ) μ 2 , but may invade or not, depending on the toxin concentration ( Figure 3).
Although in principle cooperators can invade a population of cheaters that differ in the level of resistance (e.g., rCo can invade sCh), (a) invasion success depends on α and T in , (b) double mutations are expected to be rare (μ 1 μ 2 is close to 0), and (c) if a cheater mutant of the same resistance as the cooperator invades (e.g., rCh), it will dominate the population and replace the cooperators. Accordingly, it is very difficult to maintain cooperators in the chemostat due to natural selection. This brings us to one of the main findings of the study: Even though cooperators and cheaters can coexist under some conditions, to maintain costly microbial detoxification, it is necessary to inoculate cooperators manually and to change α and T in to favor their survival. Crucially, though, because cooperators are able to invade cheaters of opposite sensitivity, these inoculations can maintain cooperators-and thereby detoxification-in the short-term. In the following sections, we show how to control the values of α and T in and inoculation probabilities to maximize the efficiency of detoxification.

| Culture conditions can be controlled to optimize detoxification efficiency
Ultimately, our goal is to maximize the efficiency of detoxification ϕ, which depends on the absolute abundance of the two types of cooperators in the chemostat. In turn, these abundances can be controlled by changing the culture conditions through two parameters: the dilution rate α and the toxin concentration flowing into the chemostat T in .
To maximize the objective function in Equation (7) (c) T = T * i when only one type of cooperators i is present. In the latter two cases, one can calculate the equilibria (analytically or numerically) and their corresponding detoxification efficiency ϕ for each culture condition (values of α and T in ). We can then find the optimal culture conditions that maximize this efficiency (Figure 4), although the equilibrium can be ecologically unstable for some parameter values.
Intuitively, the maximum efficiency is larger in a mono-culture of cooperators than in a co-culture of cooperators and cheaters of different levels of resistance (see Appendix S6). If cheaters can be excluded from the population by changing α and T in , the optimal strategy for cultivation is (a) to exclude the cheaters by adjusting the culture conditions and then (b) to change the culture conditions to maximize the productivity of a mono-culture of cooperators. F I G U R E 3 Schematic illustration of the state transitions for longterm evolutionary dynamics. Arrows represent state transitions resulting from natural selection. Solid arrows show transitions that are independent of toxin concentration, and dashed arrows transitions that depend on toxin concentration, here depicted for a given T in and α. The transition from sCo to rCo and vice versa does not occur if cooperators coexist with cheaters. Values along the arrows represent state transition probabilities

| Inoculating cooperators to optimize detoxification efficiency
Above, we showed that even though they are unlikely to appear by double mutation, cooperators can invade a population of cheaters if their level of resistance is different (Figure 3). Instead of waiting for these mutants to arise naturally, it would be more efficient to manually inoculate cooperators into the population, and to change α and T in to allow them to invade successfully and to exclude the cheaters.
Assuming that we cannot observe the prevalence of each strategy at will, the problem is how often to inoculate sensitive or resistant cooperators to maximize detoxification efficiency over time. If cooperator inoculation probabilities are too small, cheaters will dominate the population, leading to a detoxification efficiency of zero.
If they are too large, we may inoculate cooperators unnecessarily (e.g., sCo into a mono-culture of sCo) or when they cannot invade (e.g., sCo into a mono-culture of sCh). Such unfavorable inoculations can be costly because they can require a higher inflowing toxin F I G U R E 4 Optimal detoxification at each equilibrium. The efficiency of detoxification at an equilibrium state ϕ(α, T in , T) given α and T in is represented by color in each panel: (a) only sCo, (b) coexistence of sCo with rCh, (c) only rCo, and (d) coexistence of rCo with sCh. The red stars represent the maximum efficiency of detoxification in each. In the areas above the dashed gray lines in panels (a) and (c), the cooperator cannot persist in mono-culture (inequality 8).
We used different values of n in the top and bottom rows to allow cooperators to coexist with cheaters with a different level of resistance Even though the transitions are probabilistic, we assume that the establishment of strategies following mutation or inoculation (i.e., short-term dynamics) is deterministic. Relaxing this assumption by introducing an establishment probability ε into the transition matrix P (m; ) does not change the ergodicity of the Markov chain, and we arrive at the stationary distribution in the same manner. In Appendix S7, we show how to calculate the expected cumulative efficiency of detoxification Φ from the beginning of the cultivation to time step s (defined in Equation (A.67)). This calculation is somewhat cumbersome, but for large s, Φ is approximately proportional (see Appendix S7) to the expected efficiency of detoxification at the stationary distribution * : To maximize the expected cumulative detoxification efficiency Φ, therefore, we can calculate m that maximizes Equation (11). In Box 1, we show an example of how to use this approach in practice.

| D ISCUSS I ON
In this study, we have shown how to control the ecological and evo-

Box 1 An example of optimizing detoxification efficiency.
Imagine that we have set up the experimental system described, and would like to compute the optimal inoculation probabilities.
We first need to define the Markov chain to make predictions, and second, we need to experimentally measure the parameters of our bacterial strains, in particular their degradation efficiency  , 2008Sanchez & Gore, 2013). Rather than increasing the growth rate of others, cooperators in our model degrade toxic compounds, which decreases the death rate of surrounding cells (O'Brien et al., 2014). This scenario enables us to introduce the evolution of resistance to the toxin, for example, through efflux pumps, as a private good, which we base on studies of drug-dose effect and resistance to it (Chou, 2006;Sampah et al., 2011). Unsurprisingly, cheaters always exclude cooperators with the same private resistance level because detoxification is costly . We show, however, that because the benefit of private resistance depends on the toxin concentration in the chemostat, cooperators can invade a population of cheaters that differ in their private resistance. The co-occurrence of two strains that differ both in their degradation ability as well as their resistance level is unlikely to suddenly arise by mutation, especially if we assume that mutations are rare. To maintain the degradation of toxins, therefore, it is necessary to periodically inoculate cooperators into the chemostat while changing the dilution rate and inflowing toxin concentration to guarantee invasion success.
Optimal values for these parameters (cooperator inoculation probabilities, dilution rate, and inflowing toxin concentration) that maximize the detoxification efficiency of the system can be calculated using our model. As input, the model requires experimental measurements of growth and death rates of the chosen microbe and its mutants (i.e., intrinsic growth rate of each strategy r i , maximum death rate d max , Hill coefficient n, median-effect toxin concentrations K s , K r , and degradation efficiencies ϕ of cooperators).
Our model and its results can also apply to problems other than bioremediation that involve survival in toxic environments. In essence, we are studying the evolutionary dynamics of public resistance (which is cooperative) and private resistance (which is not). Consider, analogously, two types of antibiotic resistance mechanisms: Public mechanisms are costly and benefit the producing cell as well as its neighbors such as extracellular secretion of antibiotic-degrading enzymes (e.g., β-lactamases (Yurtsev, Chao, Datta, Artemova, & Gore, 2013)), and private resistance mechanisms only benefit the producing cells, such as efflux pumps. The evolutionary dynamics in a scenario whereby cells can switch between these different resistance mechanisms and being sensitive to the antibiotic correspond to Figure 2a,f. In this case, however, an objective function would aim to minimize rather than maximize the densities of the most resistant strains. Another interesting aspect is that the benefits of resistance depend on toxin concentration in the chemostat, which is affected by the density of cooperators and the toxin concentration flowing into the chemostat.
In other words, the public goods game affects the benefit of the private goods. This is why cooperators can invade a population of cheaters when they differ in their resistance level (Figure 2f).
Of course, our model relies on a number of assumptions and focuses only on a subset of possible bioremediation systems. First, we assume that only cooperators can detoxify. In reality, sensitive cells (cooperators or cheaters) may decrease toxin concentrations by passively absorbing them (Bottery et al., 2016). Including toxin absorption would not change our findings because (a) cooperators decrease toxin concentration more effectively than cheaters with the same resistance levels, (b) the invasion analysis is still valid, and (c) we can still calculate the optimal inoculation probabilities as shown in Section 3.7.
If toxin absorption by sensitive cells is significant, sCh could be better detoxifiers than rCo. In this case, it would be unnecessary to inoculate rCo to optimize the detoxification efficiency. Second, we only consider extracellular toxin degradation (e.g., by enzyme secretion), while toxins can also be degraded inside cells (O'Brien & Buckling, 2015).
For intracellular degradation, a different functional form of detoxification f(x Co ) would be necessary, but we expect similar results as long as this function increases monotonically with the density of cooperators. Indeed, the invasion analysis is independent of the form of f(x Co ).
Similarly, we assume that toxins kill the microbes and that their degradation does not contribute to growth. In reality, many compounds that are undesirable for humans are instead used as substrates by microbes (Atashgahi et al., 2018). This latter case is simpler than the one we consider here, since detoxification is no longer cooperative and there is no risk of cheaters arising and collapsing the system. Finally, detoxification may carry a negligible cost, for example, if it the toxic compound is neutralized by a change in pH, which occurs naturally due to a microbe's metabolism.
Another issue is how to define detoxification efficiency ϕ. Rather than Equation (7), one could, for example, define ϕ as the time needed for the toxins to decrease to a negligible concentration. This would change the optimal culture conditions α and T in , but not the procedure to find the optimal introduction probabilities of cooperators m, which are independent of the formulation of ϕ. Our model also fixes some parameters, such as the Hill coefficient n, which can evolve in reality (Sampah et al., 2011). Similarly, the cost of resistance c r can decrease over time due to compensatory evolution (Andersson & Hughes, 2010;San Millan et al., 2014). Allowing these parameters to evolve would make it more difficult for sCo to invade rCh because the relative fitness of rCh will increase.
We also assume that our system is well mixed and that there are no spatial gradients within the chemostat. Spatial structure, for example, whereby detoxifying enzymes diffuse slowly through the chemostat and have a patchy distribution can favor the coexistence of cooperators and cheaters (Allison, 2005). Indeed, previous empirical bioremediation studies have reported coexistence of cooperators with cheaters (Ellis et al., 2007;O'Brien et al., 2014). Theoretically, this may be due to a difference of resistance levels between cooperators and cheaters as we show here, but a simpler explanation would be the presence of spatial gradients. Relaxing the assumption of a perfectly well-mixed chemostat would make the persistence of cooperators easier. It may also increase the public benefit of toxin resistance, which we have considered to be private here (Rojo-Molinero et al., 2019).
Finally, there may be other ways of periodically introducing cooperators. Experimentally, our constant inoculation probabilities represent a situation where stock strains of cooperators would be manually added into the chemostat. If instead, multiple chemostats are running in parallel, another way of introducing cooperators would be to exchange certain amounts of fluids between chemostats.
Ecologically, this would correspond to migration among patches, and the optimal migration probabilities would depend on the probability distribution of the different strategies in each chemostat.
Comparing the optimal introduction probabilities and the cumulative efficiency of detoxification between the model presented here and a multi-chemostat system is left for future work.
In summary, we have combined an ecological model with evolutionary game theory to develop a protocol for the control and optimization of a bioremediation system by microbes, and guard it against collapse through the emergence of cheaters. More broadly speaking, our scenario motivates the integration of important elements from ecological models, such as population densities and environmental feedback, into evolutionary game theory. In essence, our model can be adapted to any practical applications involving costly microbial traits, where manipulating environmental conditions can be used to control evolutionary dynamics. Achieving this will allow us to better anticipate evolutionary change in microbial systems that we strive to control, whether this involves increasing toxin degradation as we have shown here, the production of public goods such as useful chemicals, or eliminating antibiotic resistant pathogens.

ACK N OWLED G EM ENTS
We thank Masato Yamamichi and Xiang-Yi Li for useful feedback on the earlier versions of the manuscript. SS is funded by the Nakajima

Foundation and the University of Lausanne and SM by European
Research Council Starting Grant 715097.

CO N FLI C T O F I NTE R E S T
There is no conflict of interest.

AUTH O R CO NTR I B UTI O N S
SS and SM conceived of the project, SS built and analyzed the model, and SS and SM wrote the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Python code and csv files used for simulations and analytic calculations are available at https://github.com/Shota SHIBA SAKI/Contr ollin gEvol ution adyDy namic sGitHub.