In silico bacteria evolve robust cooperation via complex quorum-sensing strategies

Many species of bacteria collectively sense and respond to their social and physical environment via ‘quorum sensing’ (QS), a communication system controlling extracellular cooperative traits. Despite detailed understanding of the mechanisms of signal production and response, there remains considerable debate over the functional role(s) of QS: in short, what is it for? Experimental studies have found support for diverse functional roles: density sensing, mass-transfer sensing, genotype sensing, etc. While consistent with theory, these results cannot separate whether these functions were drivers of QS adaption, or simply artifacts or ‘spandrels’ of systems shaped by distinct ecological pressures. The challenge of separating spandrels from drivers of adaptation is particularly hard to address using extant bacterial species with poorly understood current ecologies (let alone their ecological histories). To understand the relationship between defined ecological challenges and trajectories of QS evolution, we used an agent-based simulation modeling approach. Given genetic mixing, our simulations produce behaviors that recapitulate features of diverse microbial QS systems, including coercive (high signal/low response) and generalized reciprocity (signal auto-regulation) strategists — that separately and in combination contribute to QS-dependent resilience of QS-controlled cooperation in the face of diverse cheats. We contrast our in silico results given defined ecological challenges with bacterial QS architectures that have evolved under largely unknown ecological contexts, highlighting the critical role of genetic constraints in shaping the shorter term (experimental evolution) dynamics of QS. More broadly, we see experimental evolution of digital organisms as a complementary tool in the search to understand the emergence of complex QS architectures and functions.

N (0, SD) to the original trait value, where N (0, SD) is the normal distribution with a mean 0 and a standard deviation SD to be SD p , SD S T h or SD r for different traits. 8) Repeat 2) to 7) until Gen max generations were reached.

Model assumptions
For the computational models presented in the paper, we assume that: 1) A single signal type exists.
2) In the propagule pool, the cellular heterogeneity can be obtained via mutation and selection.
3) Each sub-population is founded by a defined number of founding cells G (1 to 10), picked randomly from the propagule pool of the previous generation. 4) For each sub-population, the expected fitness of each founding cell is assessed across a range of testing environments (cellular density varied from 10 1.5 to 10 5 cells per µL) where the founding cells immediately grow until the total sub-population reaches a stationary phase density defined by its current testing environment. 5) In each testing environment, the sub-population forms a closed system, i.e., no mass transfer. 6) For a given testing environment, an individual can be rewarded with a benefit for cooperation only if the number of cooperators is greater than the defined threshold (N T h ).
7) The signal concentration in each testing environment rapidly reaches equilibrium estimated by Eq. (S1) or Eq. (S6), depending on whether invoking the auto-regulation mechanism.

Computational model of quorum sensing without auto-regulation
The computational model of signal dynamics for quorum sensing without auto-regulation is given as below: where S is the local signal concentration, t is time, N is the local cell density, p is the basal signal production rate, and u is the signal decay rate. The equilibrium of Eq. (S1) is given by: where G is the number of mixing genotypes in a sub-population, N j is the local cellular density in the j th sub-population testing environment, and p g and S T hg are the signal production rate and signal response threshold of the genotype g (g = 1, 2, · · · , G) in the sub-population, respectively. Similarly, the function of cooperative benefit of the individual i in the sub-population testing environment j is defined as: where N T h is the cellular density threshold (defined as the median cellular density across all testing environments).

Computational model of quorum sensing with auto-regulation
The computational model of signal dynamics for quorum sensing with auto-regulation is given as below: where S is the local signal concentration, t is time, N is the local cell density, p is the basal signal production rate, r is the ratio of auto-regulation production to basal signal production, K is the half concentration value, and u is the signal decay rate. The equilibrium of Eq. (S6) is given by: the purpose of computational convenience, we used 1 as the Hill function exponent, which can lead to a close form solution, Eq. (S7). When invoking the auto-regulation mechanism, the individual genotype's overall cooperation payoff across all testing environments is assessed as follows: where i (i = 1, 2, · · · , N pop ) represents an individual genotype, j (j = 1, 2, · · · , N env ) represents the index number of a sub-population testing environment, B 0 is the baseline payoff, B coop , C coop , C sig are constants for the benefit of cooperation, cost of cooperation and cost of signaling, respectively, p i and r i are the basal signal production rate and auto-regulation ratio of the genotype i, respectively, and . The function of cooperation cost of the individual i in the sub-population testing environment j is defined as: where G is the number of mixing genotypes in a sub-population, S * gj (calculated by Eq. (S7)) and S T hg are the equilibrium signal concentration and signal response threshold of the genotype g (g = 1, 2, · · · , G) in the sub-population, respectively. The function of cooperative benefit of the individual i in the sub-population testing environment j, and H Bij is defined as the same as in Eq. (S5).

Adding noise to signal
To investigate how clonal populations cope with signal noise to sustain cooperation, we added noise to the equilibrium signal. In the simulations, the noise signal is drawn from a normal distribution with mean S * and standard deviation κ · S * , i.e., N (S * , κ · S * ), where S * is the equilibrium signal calculated from Eq. (S2) or Eq. (S7) and κ is a constant indicating the strength of noise. Note we set all negative values for the noise signal to be 0.

Generating genetic mixing
In the evolution simulations where we varied the genetic relatedness, different numbers of mixing genotypes were introduced to form sub-populations to evaluate the overall cooperation payoff for each genotype in every generation. Unless specified otherwise, the actual number of mixing genotypes in each sub-population in every generation was drawn from a zero-truncated Poisson distribution with the average beingḠ where (λ G ∈ [0, 10]; step size: 0.1). We define the clonal case (G = 1) when λ G = 0 where exact one genotype will be selected to form the sub-population, i.e., no genetic mixing. Note that the number of sub-populations may be different in every generation due to the variation of mixing genotypes in each subpopulation.

Constructing constitutive cooperators
To investigate how decision making interact with social behaviors of cooperation, we compared the overall payoff of cooperation of individuals mediated by QS with those in the absence of collective control. Specifically, we constructed constitutive cooperators which do not have the ability to make social informed choices.
In the clonal case, wild-type individuals will always cooperate. This will incur a penalty to each of such individuals for cooperating in wrong environments 1 . In the genetic mixing scenarios, the cooperative benefits of wild-type individuals will be shared evenly with all group members. In the simulations, all individuals were subject to mutation, switching from a wild-type to mutant, or mutant to wild-type depending on their own initial type. The actual number of replacement individuals was drawn from a Poisson distribution with the mean being 0.01. Note that mutant individuals will always reap the benefits of cooperation without paying for any cost. Formally, the overall cooperation payoff in the constitutive cooperation scenarios can be defined as: where P W T is the proportion of wild-type individuals in the sub-population with a group size G.

Introducing constitutive cheats
To challenge the quorum sensing system, in the simulations of invasion by a cheat phenotype, we replaced a certain number of individuals 2 chosen at random in the population pool with constitutive cheats in every generation at a certain rate. The constitutive cheat is defined as a genotype with a zero basal production rate and a maximum possible signal threshold. In other words, a constitutive cheat does not produce or respond to signal. The actual number of cheats introduced into each generation is drawn from a Poisson distribution with the mean being λ Cheat . Note that the constitutive cheats are both immutable, which means they cannot be eliminated through mutation, and inheritable, which means their offspring are still cheats.

Measuring phenotypic assortment of cooperative investment
To test if the auto-regulation mechanism could be explained by the generalized reciprocity theory, we recoded the mean value of cooperative investment 3 within each sub-population in the genetic mixing scenario where individuals are grouped into small collectives. We then plotted the group mean cooperative investment against individual cooperative investment. Finally, the regression line was fitted using the generalized linear model with a normal distribution. The slope of the regression line indicates the phenotypic assortment of cooperative investment. When the slope is high, the behaviors among individuals shifts closer to each other, investing more in cooperation. Otherwise, the behaviors of investment for cooperation vary among individuals.

Partitioning selection on cooperative investment
To further uncover the influence of the auto-regulation mechanism on cooperative behaviors in our evolu- where z represents the trait cooperative investment, w i is the number of offspring (fitness) produced by the individual i,w is the mean number of offspring produced, ∆z i represents the difference between the average z value among the individual i's offspring and i's own z value, Cov i (·) and E i (·) denote the expectation and covariance over all individuals i in the population respectively.
By introducing the genetic mixing in the simulations, individuals in the population have been assigned into small groups. We are able to further partition that selection based on cooperative investment to account for individuals that are nested within collectives. Specifically, we can expand Eq. (S11) by substituting it into the expectation term in its right hand side. Note that the groups that form each subpopulation g are the individuals ig. We can re-write the two-level Price equation as follows: wherew g = E(w ig ) andz g = E(z ig ). The first covariance term on the right hand side of the equation indicates the selection on cooperative investment at level of subpopulations (between-group selection), whereas the second expectation term captures the selection at individual level (within-group selection). Figure S1. Evolution trajectories of basal production rate and signal threshold for high cost of signaling. We evolved 5, 000 initially identical genotypes for 5, 000 generations in a patchy, variable density environment (see Main Text Figure 6). In all simulations, the cost of cooperation was fixed and there was no auto-regulation. We used a fixed cost of signaling (C sig = 10 10 ), and recorded the evolution trajectories of basal production rate and signal threshold. All reported results were averaged over 30 replications (shaded area indicates the standard deviation). The remaining parameters used in the simulations can be found in Table S1. Although the evolved traits were at equilibrium, the predicted transients are consistent with our previous mathematical model [4]. It should be noted that when p = 0, S T h will stop evolving due to the cooperation collapse.  Table S1. for QS with no auto-regulation and auto-regulation, respectively. All reported results were averaged over 30 replications. Note that surfaces were smoothed using the spline interpolation method. The remaining parameters used in the simulations can be found in Table S1.  Figure S4. Evolved traits against genetic relatedness. We evolved 5, 000 initially identical genotypes for 5, 000 generations with no auto-regulation (A) and auto-regulation (B), respectively. In all simulations, the cost of cooperation and the cost of signaling were fixed. Each dot represents the evolved mean results (averaged over the last 50 generations) for different average number of genotypes per groupḠ (λ G ∈ [0, 10]; step size: 0.1). The remaining parameters used in the simulations can be found in Table S1. Figure S5. Evolution trajectories of traits for QS-controlled cooperation with no autoregulation. For a population with an intermediate genetic mixing (λ G = 2), we evolved 5, 000 initially identical genotypes for 5, 000 generations with no auto-regulation. We used a fixed cost of cooperation and a fixed cost of signaling, and recorded the evolution trajectories of basal production rate and signal threshold. All reported results were averaged over 30 replications (shaded area indicates the standard deviation). The remaining parameters used in the simulations can be found in Table S1.
Intermediate Relatedness (λG = 2) Intermediate Relatedness (λG = 2) Low Relatedness λG = 10 Low Relatedness λG = 10  Table S1. The full evolution timelapse can be found in Supplemental Material, Videos S1-S6.  Figure S7. Overall cooperation payoff for fixed trait evolution. We used a fixed signal production rate of 0.5 × 10 −8 and a fixed response threshold of 3 µM for QS-controlled cooperation without autoregulation, respectively. In all cases, we evolved 5, 000 initially identical genotypes for 5, 000 generations. Each dot represents the evolved mean results (averaged over the last 50 generations) for different average number of genotypes per groupḠ (λ G ∈ [0, 10]; step size: 0.1). The vertical error bars represent the standard deviation of overall cooperation payoff over 30 replications. The remaining parameters used in the simulations can be found in Table S1. Figure S8. Evolution trajectories of traits for QS-controlled cooperation with auto-regulation. For a population with an intermediate genetic mixing (λ G = 2), we evolved 5, 000 initially identical genotypes for 5, 000 generations with no auto-regulation. We used a fixed cost of cooperation and a fixed cost of signaling, and recorded the evolution trajectories of basal production rate, observed production rate and signal threshold. All reported results were averaged over 30 replications (shaded area indicates the standard deviation). Note that the basal production rate is slightly increasing which implies it is approaching the observed production rate over a longer evolution timescale. However, the observed production rate is largely converged. The remaining parameters used in the simulations can be found in Table S1. Auto-regulation Ratio Figure S9. Evolved auto-regulation ratio against genetic relatedness. We evolved 5, 000 initially identical genotypes for 5, 000 generations with auto-regulation. In the simulations, the cost of cooperation and the cost of signaling were fixed. Each dot represents the evolved mean results (averaged over the last 50 generations) of auto-regulation ratio (r as in Eq. (S6)) for different average number of genotypes per group G (λ G ∈ [0, 10]; step size: 0.1). The vertical error bars represent the standard deviation over 30 replications. The remaining parameters used in the simulations can be found in Table S1.  Table S1.  Table S1. Figure S12. Regression analysis for individual and group mean investment for cooperation (G = 2). For fixed costs of cooperation and signaling with the number of mixing genotypes G = 2, we collected 5, 000 same initial genotypes and evolved them for 5, 000 generations with no auto-regulation (A) and auto-regulation (B). We recorded the individual and group mean investment for cooperation at the last generation over 100 replications. Each blue dot represents an individual's investment against its group mean investment. The red lines are the regression lines fitted using the generalized linear model with a normal distribution. The analysis of covariance shows there is a significant difference between the slope of no autoregulation in (A) and the slope of auto-regulation in (B) (F -test, p = 0.000). The remaining parameters used in the simulations can be found in Table S1.  Figure S13. Selection on cooperative investment within and between groups (G = 2). We evolved 5, 000 initially identical genotypes for 5, 000 generations with no auto-regulation (A) and autoregulation (B), respectively. In all simulations, the cost of cooperation and the cost of signaling were fixed, and the number of mixing genotypes was fixed G = 2. We recorded the two-level Price equation components in every generation. The reported results were the average value over 100 replications. The remaining parameters used in the simulations can be found in Table S1.  Figure S14. Selection on cooperative investment within and between groups (G = 5). We evolved 5, 000 initially identical genotypes for 5, 000 generations with no auto-regulation (A) and autoregulation (B), respectively. In all simulations, the cost of cooperation and the cost of signaling were fixed, and the number of mixing genotypes was fixed G = 5. We recorded the two-level Price equation components in every generation. The reported results were the average value over 100 replications. The remaining parameters used in the simulations can be found in Table S1.  p min /p max the minimum/maximum basal production rate 0/2 × 10 −8 µMs −1 per cell p init initial value for basal production rate 0.5 × 10 −8 µMs −1 per cell r min /r max the minimum/maximum auto-regulation ratio 0.01/50 r init initial value for auto-regulation ratio 10 u signal decay rate 10 −4 µLs −1 m mass transfer rate 0 to 5 × 10 −5 µLs −1 1 FU: fitness unit 2 CTE: cooperative testing environment 3 CB: cooperative behaviour 4 The local cellular densities are evenly spaced within the range 10 1.5 to 10 5 (cells per µL).