Putting theory to the test: An integrated computational/experimental chemostat model of the tragedy of the commons

Cooperation via shared public goods is ubiquitous in nature, however, noncontributing social cheaters can exploit the public goods provided by cooperating individuals to gain a fitness advantage. Theory predicts that this dynamic can cause a Tragedy of the Commons, and in particular, a ‘Collapsing’ Tragedy defined as the extinction of the entire population if the public good is essential. However, there is little empirical evidence of the Collapsing Tragedy in evolutionary biology. Here, we experimentally demonstrate this outcome in a microbial model system, the public good-producing bacterium Pseudomonas aeruginosa grown in a continuous-culture chemostat. In a growth medium that requires extracellular protein digestion, we find that P. aeruginosa populations maintain a high density when entirely composed of cooperating, protease-producing cells but completely collapse when non-producing cheater cells are introduced. We formulate a mechanistic mathematical model that recapitulates experimental observations and suggests key parameters, such as the dilution rate and the cost of public good production, that define the stability of cooperative behavior. We combine model prediction with experimental validation to explain striking differences in the long-term cheater trajectories of replicate cocultures through mutational events that increase cheater fitness. Taken together, our integrated empirical and theoretical approach validates and parametrizes the Collapsing Tragedy in a microbial population, and provides a quantitative, mechanistic framework for generating testable predictions of social behavior.


S1. Model parameterization
Initial concentration of gelatin in the chemostat, S(0): Our M9 gelatin medium contains 1 g of gelatin per 100 ml, or equivalently 10 g per liter.
The strength, or stiffness, of gelatin is indicated by its bloom value.The brand of gelatin used is Alfa Aesar Type A with a bloom of 175, and the average molecular weight of gelatin with a bloom of 175 to 225 is in the range 40 to 50 g/mmol according to product information sheet on www.sigmaaldrich.com.
We can estimate (0) using the average molecular weight and quantity of gelatin.For the case where we estimate the average molecular weight as 50 g/mmol, we calculate this as follows: Since the average molecular weight of gelatin is estimated to be between 40 and 50 g/mmol, we get 0.20 mmol l ≤ (0) ≤ 0.25 mmol l .
Initial concentration of P. aeruginosa strains in the chemostat,   () and   (): Sexton et al. determined that an OD600 =1 is equivalent to 0.48 grams of dry weight per liter for P. aeruginosa grown in a chemostat in minimal medium [1].
By inoculating the bioreactor at an OD600 =0.05 we can determine the initial population density in g DW/l (either X1(0) alone or both X1(0) and X2(0) in coculture).
Given that the molecular weight of the LasB protein is 33 g/mmol [5,6] we can determine the average concentration as The range of enzyme measured in that study was 0-80 µg/ml, with zero representing values that were below the detectable limit.Thus, the actual initial condition could lie anywhere between This range represents the amount of protease enzyme we expect to find in the batch culture prior to inoculation into the chemostat bioreactor.The bioreactor was inoculated from stationaryphase LB cultures to an OD600 of 0.05, resulting in a 100-fold dilution.We consider the same dilution factor for mono and cocultures, neglecting the small reduction in enzyme concentration through inoculation with 90% of the WT strain in cocultures.Thus, the amount of enzyme in the bioreactor will be 0.0 mmol l < (0) ≤ 2.42x10 mmol l .
Initial concentration of product in the chemostat, P(0): We refer to amino acids produced by extracellular protease as "products" collectively.In the experimental design, the chemostat is inoculated from LB stationary phase cultures.We assume that all the usable nutrient and product has been exhausted by the organism in the stationary phase cultures.This means that no product, or possibly a negligible amount of product, is being transferred.
Furthermore, proteolysis-deficient lasR mutants show little to no growth in the 1% gelatin medium (Figure 3 of main text) indicating that there are no breakdown products in fresh gelatin medium.Thus, Chemostat dilution rate, D: The equation to convert peristaltic pump speed in revolutions per minute (RPM) to flow rate was determined experimentally by measuring the time it took to fill a beaker with 10 ml of water at four different pump speeds: 5, 15, 35, and 65 RPM (Table S1).
Table S1.1:Duration to achieve specific volume at various pump speeds.
Using linear regression, we were then able to determine the following equation to convert pump speed (x) to flow rate (y) as y = 0.0731x (Figure S1).
The chemostat experiment is operated at an RPM=2.75 with a bioreactor volume of 100 ml.The resulting flow rate is determined by RPM time (min) to fill 10 ml ml/min The growth rate, F(P), models the growth of P. aeruginosa as a Monod equation dependent on product concentration (P).As indicated above, the products are primarily individual amino acids that are the result of extracellular proteolysis.μ is estimated to be 1.38 1/h when grown in minimal medium containing a 1% mixture of amino acids (casamino-acid, CAA) [7].
Of note, we calculated an exponential-phase growth rate of µ = 0.17 from our batch cultures.That calculation is based on substrate degradation (S).Since our system of equations models growth from the product of degraded substrate (i.e.we use F(P) instead of F(S)), it is more accurate to use the growth rate observed from CAA consumption.

Rate of enzyme-substrate turnover, 𝒌 𝒄𝒂𝒕 :
A lower bound for the rate of enzyme-substrate turnover was estimated to be 1 mass % S / (mass % E • min) [7], where mass % S and mass % E are the biomass concentrations (in w/v) of the substrate and enzyme respectively.1 mass % equals 1 g of E or S per 100 ml of solution, or 10 g per l.
Here, we convert mass % to molar concentrations (in mmol/l).For S we get  = ⋅ 0.02 , where we consider the molecular weight of gelatin to be 50 g/mmol.
The molar concentration of  = ⋅ using the molecular weight of LasB protein to be 33 g/mmol [5,6].
Converting this to kcat we get: 1 mass %  mass %  ⋅ min = 10 g/l ⋅ 0.02 mmol/g 10 g/l ⋅ 0.03 mmol/g ⋅ min ⋅ 60 min h = 40 1/h Given the molecular weight of gelatin was estimated to be between 40 and 50 g/mmol (see above), this yields an estimation of: The authors suggest that the conditions of their experiment likely underestimated kcat, and that the true kcat rate is higher.Indeed, many catalysis rates for the enzymatic activity of proteases lie in the hundreds and thousands per hour [8,9], and other enzymes produced by P. aeruginosa have a catalysis rates in the thousands to the tens or hundreds of thousands per hour [10,11].For these reasons, we expanded the upper bound of the range approx.10-fold to get Amount of product present at half-maximal growth rate,   : Ks is challenging to measure.One reason for this is the difficulty of measuring very low concentrations of growth-limiting product at steady-state [7,12].However, in P. aeruginosa, the Ks for CAA has been estimated as 0.01 mass %, or equivalently .
[7].The brand of gelatin used is Alfa Aesar Type A with a bloom of 175, and the average molecular weight of this kind of gelatin is 40 to 50 g/mmol.
2.9 g l ⋅ mmol 50 g ≤  ≤ 2.9 g l ⋅ mmol 40 g 0.058 mmol l ≤  ≤ 0.0725 mmol l Amount of enzyme produced per bacterial biomass, : The concentration of LasB elastase in an LB stationary-phase culture has a range of 0-80 µg/ml (where 0 represents falling below a detectable limit), with an average of 30.9 µg/ml (equivalently 0.0309 g/l) [4] and a growth yield of  between 4 and 5.
If we consider the range for an OD600 between 4 and 5, we get 0.000390 ≤  ≤ 0.000488.However, while the average concentration was 30.90 µg/ml, the range is 0-80 µg/ml -further increasing the uncertainty to 0.0 < η ≤ 0.00126 Nutrient uptake to biomass conversion factor, : We consider amino acids to be the main products of gelatin proteolysis, which entails LasB and other endo-and exopeptidases [2].In minimal medium containing 0.5% (w/v), or 5 g/l, casamino acids, P. aeruginosa grows to an  of approximately 3.7 [14].We can calculate population density in gDW/ml from the OD600 [1] as follows: 3.7•0.48g/l = 1.776 g/l.
The average molar mass of amino acids is 118.89 g/mol = 0.119 g/mmol [13].Thus, the growth yield is: Number of product molecules produced per substrate molecule in enzyme-substrate catalysis, : The brand of gelatin used (Alfa Aesar Type A) has a bloom of 175 and an average molecular weight of 40 to 50 g/mmol.
The number of amino acids produced, if all of the substrate was converted into usable products, is in the range of , or 336 to 420 respectively.
However, when comparing growth yields in M9 gelatin and in M9 CAA media, it is apparent that only a fraction of the gelatin is converted to usable product.The average growth yield of wild type P. aeruginosa with 1% (w/v) gelatin was OD600 = 0.808 (Figure 3 of main text), whereas the average growth yield with 0.5% casamino acids is OD600 = 3.7 [14].
Thus, the conversion factor of usable product is ⋅ . .= 9.16.
Since the breakdown product molecules can only exist as whole units (i.e.half of an amino acid cannot be produced), we extend this parameter range slightly to ignore fractions yielding: Metabolic burden of enzyme production as a fraction of the total growth rate, q: We can determine q from the ratio of growth rates between strains (also termed relative fitness) during the batch phase of chemostat cultivation.However, we know from previous work that the lasR mutant does not enrich at a constant rate during coculture growth [14] and hence the WT does not experience a constant metabolic burden of enzyme production.In fact, following inoculation and prior to reaching a quorum, the lasR mutant initially grows at the same rates as the WT in coculture, presumably because both strains benefit from the products present in the inoculum.These include secreted proteases and proteolysis products carried over from the preculture.We can attempt to more accurately determine relative fitness by disregarding this initial phase and only considering the phase of lasR mutant enrichment.
Based on this previous work [14], the phases of no enrichment vs. enrichment can be distinguished by the total cell density, with a density of approx.2.5-4x10 8 CFU/ml indicating the initial point of lasR mutant enrichment.This point overlaps with the QS threshold determined below.We can then estimate the initial lasR mutant cell densities at this point according to the mutant fraction quantified upon inoculation of our chemostat cocultures, and we can determine the final cell densities directly from our measurements at the end of the batch phase of our chemostat cocultures.
First, we consider the case where the initial total estimated CFU/ml at point of lasR mutant enrichment is 2.5x10 8 .In this case, the first and second chemostat coculture replicates have a lasR mutant density of 3.43x10 7 and 1.83x10 7 respectively.After 49 hours in batch mode, the first replicate reached a total CFU/ml of 9.67x10 8 and a lasR mutant CFU/ml on antibiotic medium of 2.70x10 8 .The second replicate ended batch mode after 55.5 hours with a total CFU/ml of 1.98x10 9 and a lasR mutant CFU/ml of 7.17x10 8 .
Thus, for the first replicate we can calculate the absolute fitness of the WT and lasR strains as: WT abs fit = ln 9.67x10 − 2.7010 2.510 − 3.4310 = 1.17 lasR abs fit = ln 2.710 3.4310 Hence the relative fitness of the lasR mutant is 2.06/1.17=1.76.
Following a similar calculation as above for the second replicate, we can determine the absolute fitness of the WT and the lasR mutant fitness to be 1.70 and 3.66, respectively.
For the second replicate, the relative fitness of the lasR mutant is 3.66/1.70=2.15.
We calculate the reduced growth of the WT, and q, by setting the lasR mutant growth rate equal to 1, or 100% of the possible growth.For the first replicate we get: Similarly, considering the second replicate where the relative fitness is 2.15, we get  = 0.535.
Repeating the steps above for an initial total CFU/ml of 4x10 8 , the first and second replicates will have a lasR density of 5.48x10 7 and 2.92x10 7 , respectively.Then reproducing the above calculations with these new initial densities yields a relative fitness of 2.29 and 2.60 for Replicate 1 and Replicate 2, respectively.Thus, for Replicate 1 q=0.563 and for Replicate 2 q=0.615.
Hence, we will consider q to be within the range of the minimum and maximum estimated q: 0.432 ≤  ≤ 0.615.
However, to simplify the numerical analysis we consider: 0.425 ≤  ≤ 0.625.
Density of cooperators needed for quorum sensing, QSmin: It is estimated that the QS threshold occurs somewhere in the range of 1.0 x 10 8 ≤ CFU/mL ≤ 5.0 x 10 8 based on QS gene expression data [14].The corresponding value in gDW/mL can be determined by first converting CFU/mL to OD600 according to the correlation in Fig. 1C, and then converting OD600 to gDW/mL using the conversion factor as before.The correlation in Fig. 1C yielded a best fit equation of where the slope is 1.031 to 1.482 and the intercept is 9.272 to 9.320.Considering the average slope and intercept, 1.257 and 9.296 respectively, CFU/mL values of 1 x 10 8 and 5 x 10 8 yield OD600 values of approx.0.0931 and 0.335, respectively.We can then multiply the OD600 by the conversion factor 0.48 to obtain gDW/mL [1].This gives us a range of 0.0447 ≤  ≤ 0.161.

Hill coefficient for QS function, n:
There is evidence that QS in P. aeruginosa functions like an 'on-off' switch, with positive autoregulation and receptor dimerization as the underlying mechanisms [15][16][17].To obtain this type of behavior in our model with the QS function Q(X1), we require an ultrasensitive response with a Hill coefficient n > 1 [18,19].A Hill coefficient of 2 is commonly used to model the rheostat like behavior of bacterial growth, thus, we use n=2 in this study.

Fig S1. 1 .
Fig S1.1.Correlation between pump speed and flow rate.Linear regression was used to establish an equation which converts pump speed (RPM) (x) to flow rate (ml/min) (y).

Fig S1. 2 .Fig S1. 3 .
Fig S1.2.Line graphs illustrating 10,000 simulations with the estimated parameter and initial condition values randomly chosen from the estimated parameter ranges for three experimental data types: (A) WT-only population density, (B) coculture cell density, and (C) the cheater frequency.Experimental data is shown in green and simulations are shown in an opaque black line, thus, darker regions occur where more simulations overlap with one another.This figure shows that growth conditions mirroring our results are possible within our established estimated parameters.

Fig S1. 4 .
Fig S1.4.Comparison of the best fit as determined by the lowest weighted normalized RMSE values (WNRMSE) (green) and the final best fit (black) across three data sets: WT only (left), coculture (center), and cheater frequency (right).Solid light grey lines show the chemostat data for comparison.A vertical dashed grey line indicates the transition from batch mode to chemostat mode.Horizontal dashed grey lines indicate the range (coculture) or minimum (cheater frequency) cutoffs used to find the simulation with the overall best fit.
Maximum growth rate constant,   :