Steering ecological-evolutionary dynamics to improve artificial selection of microbial communities

Calculating landscape, attractor, and restrictor

In this work, we considered communities with commensal, mutualistic, and exploitative interactions. Below, we describe the differential equations for each type of interaction, and how we calculate the corresponding community function landscape, species-composition attractor, and Newborn restrictor.

Commensal H–M community: The model community for most simulations is the same commensal H–M community used in our previous work¹⁵. The community function landscape plots P(T) as a function of ϕ_M(0) and ({overline{f}}_{P}(0)). Assume that a Newborn community has 100 biomass units, that all cells have the same genotype (all M cells have the same ({f}_{P}={overline{f}}_{P}(0))), that death and birth processes are deterministic, and that there is no mutation. P(T) can then be numerically integrated from the following set of scaled differential equations for any given pair of ϕ_M(0) and ({overline{f}}_{P}(0))¹⁵:

where

and R_M = R/K_MR and B_M = B/K_MB. Unless otherwise specified, landscapes in this paper are obtained by integrating Equations (1–5) from t = 0 to t = 17.

Equation (1) states that Resource R is depleted by biomass growth of M and H, where c_RM and c_RH represent the amount of R consumed per unit of M and H biomass, respectively. Equation (2) states that Byproduct B is released as H grows, and is decreased by biomass growth of M due to consumption (c_BM amount of B per unit of M biomass). Equation (3) states that Product P is produced as f_P fraction of potential M growth. Equation (4) states that H biomass increases at a rate dependent on Resource R in a Monod fashion (Equation (6)) and decreases at the death rate δ_H. Note that Agricultural waste is not a state variable here as it is present in excess. Equation (5) states that M biomass increases at a rate dependent on Resource R and Byproduct B according to the Mankad and Bungay model (Equation (7)⁵¹) discounted by (1 − f_P) due to the fitness cost of making Product, and decreases at the death rate δ_M. In the Monod growth model (Equation (6)), g_Hmax is the maximal growth rate of H and K_HR is the R at which g_Hmax/2 is achieved. In the Mankad and Bungay model (Equation (7)), K_MR is the R at which g_Mmax/2 is achieved when B is in excess; K_MB is the B at which g_Mmax/2 is achieved when R is in excess.

Mutualistic H–M community: If Byproduct is harmful for H, then the community is mutualistic: H and M promote the growth of each other. Such a mutualistic community can still be described by Equations (1–5) and (7), but Equation (6) is replaced with

where larger B₀ indicates lower sensitivity, or higher resistance of H to its Byproduct B.

Exploitative H–M community: If M releases an antagonistic byproduct A that inhibits the growth of H, then the interaction is exploitative: H promotes the growth of M, but M inhibits the growth of H. Besides Eqs (1–5) and (7), we then need to add an equation that describes the dynamics of A

where r_A is the amount of A released when M’s biomass grows by 1 unit. We can then normalize (widetilde{A}) with r_A

so that

We also need to modify the growth rates for H:

where larger A₀ indicates lower sensitivity, or higher resistance of H to M’s Antagonistic by product A.

To calculate the community function landscape, species attractor, and Newborn restrictor, all phenotype parameters, except ({overline{f}}_{P}(0)) take the value from the Bounds column in Table 1. To construct the landscape such as in Fig. 2c, we calculated P(T) for every grid point on a 2D quadrilateral mesh of 10⁻² ≤ ϕ_M(0) ≤ 0.99 and (1{0}^{-2} le {overline{f}}_{P}(0) le 0.99) with a mesh size of Δϕ_M(0) = 10⁻² and ({{Delta }}{overline{f}}_{P}(0)=1{0}^{-2}). To construct the landscapes in Fig. 5b(ii) and b(iii), P(T) was similarly calculated on a 2D grid with a finer mesh of Δϕ_M(0) = 5 × 10⁻³ and ({{Delta }}{overline{f}}_{P}(0)=1{0}^{-4}).

To calculate the species composition attractor, we integrated Equations (1–5) to obtain ϕ_M(T) − ϕ_M(0) for each grid point on the 2D mesh of ϕ_M(0) and ({overline{f}}_{P}(0)). The contour of ϕ_M(T) − ϕ_M(0) = 0 is then the species attractor (blue dashed curve in Fig. 2b).

The attractor-induced Newborn restrictor at a given ({overline{f}}_{P}(0)) is calculated from its definition: if ϕ_M(0) of a parent Newborn is on the restrictor, then so is the average ϕ_M(0) among its offspring Newborns. Under no spiking, since the average ϕ_M(0) among offspring Newborn is the same as ϕ_M(T) of their parent Adult, the Newborn restrictor coincides with the species attractor (Fig. 3b and Fig. 5b ii). Under x% H spiking, x% of the biomass in Newborns is replaced with H cells. Thus if the parent Adult’s fraction of M biomass is ϕ_M(T), the average ϕ_M(0) among its offspring Newborns is (1 − x%)ϕ_M(T) under x% H spiking. The Newborn restrictor therefore is the contour of (1 − x%)ϕ_M(T) − ϕ_M(0) = 0 (teal curve in Fig. 5a ii and b iii, Fig. 2d ii). Compared with the orange restrictor under no spiking, the teal restrictor is shifted down.

Parameter choices

Details justifying our parameter choices are given in the Methods section of our previous work¹⁵. Briefly, our parameter choices are based on experimental measurements of microorganisms (e.g., S. cerevisiae and E. coli). To ensure the coexistence of H and M, M must grow faster than H for part of the maturation cycle since M has to wait for H’s Byproduct at the beginning of a cycle. Because we have assumed M and H to have similar affinities for Resource (Table 1), the maximal growth rate of M (g_Mmax) must exceed the maximal growth rate of H (g_Hmax), and M’s affinity for Byproduct (1/K_MB) must be sufficiently large. Moreover, metabolite release and consumption need to be balanced to avoid extreme species ratios. We assume that H and M consume the same amount of Resource per new cell (c_RH = c_RM) since the biomass of various microbes shares similar elemental (e.g., carbon or nitrogen) compositions. We set consumption value so that the input Resource can support a maximum of 10⁴ total biomass. The evolutionary bounds are set, such that evolved H and M could coexist for f_p < 0.5, and that Resource was on average not depleted by T to avoid cells entering stationary phase.

In our simulations, we define “mutation rate” as the rate of nonneutral mutations that alter a phenotype. For example in yeast, mutations that increase growth rate by ≥2% occur at a rate of ~10⁻⁴ per genome per generation (calculated from Fig. 3 of Levy et al.⁵²), and mutations that reduce growth rate occur at a rate of 10⁻⁴ ~ 10⁻³ per genome per generation^53,54. Moreover, mutation rate can be elevated by as much as 100-fold in hypermutators. In our simulations, we assume a high, but biologically feasible, rate of 2 × 10⁻³ phenotype-altering mutations per cell per generation per phenotype to speed up computation. At this rate, an average community would sample ~20 new mutations per phenotype during maturation. When we simulated with a 100-fold lower mutation rate, evolutionary dynamics slowed down, but all of our conclusions still held¹⁵. Among phenotype-altering mutations, tens of percent create null mutants, as illustrated by experimental studies on protein, viruses, and yeast^53,55,56. Thus, we assumed that 50% of phenotype-altering mutations were null (i.e., resulting in zero maximal growth rate, zero affinity for metabolite, or zero f_P). Among nonnull mutations, the relative abundances of enhancing versus diminishing mutations are highly variable in different experiments. We based our distribution of mutation effects on experimental studies on S. cerevisiae where the fitness effects of thousands of mutations were quantified under various nutrient limitations in an unbiased fashion⁵⁷. The relative fitness changes caused by beneficial (phenotype-enhancing) and deleterious (phenotype-diminishing) mutations can be approximated by a bilateral exponential distribution with means s₊ = 0.050 ± 0.002 and s₋ = 0.067 ± 0.003 for the positive and negative halves, respectively.

Simulating community selection with small population size

Individual-based stochastic simulation codes used in this work are largely similar to those in our previous work, except for the modification to simulate species spiking. Below, we briefly recapture the flow of the simulation, which can be found in our previous work¹⁵.

Each simulation begins with n_tot = 100 identical Newborns. The total biomass of each Newborn is BM_target. In most simulations, BM_target = 100 consisting of 60 M cells and 40 H cells of biomass 1. Each Newborn is supplied with abundant agriculture waste and a fixed amount of Resource that supports the growth of 10⁴ total biomass. Unless otherwise specified, maturation time is set to T = 17 (~6 generations) to avoid Resource depletion (i.e., stationary phase in experiments) and cheater takeover.

Each maturation cycle is divided into time steps of length Δτ = 0.05. During each time step, the biomass of each M and H cell grows deterministically according to Equations (4) and (5) (or the corresponding equations for mutualistic and exploitative communities) without the death terms, while the concentration of Resource, Byproduct and Product changes according to Equations (1–3). At the end of each Δτ, each M and H cell dies with a probability of δ_MΔτ and δ_HΔτ, respectively.

Among the survived cells, if a cell’s biomass exceeds the threshold of 2, the cell divides into two identical daughter cells. Each daughter cell then mutates with a probability of P_mut. For a M cell, its f_P, g_Mmax, K_MR and K_MB can mutate independently. For a H cell, its g_Hmax and K_HR can mutate independently. Additionally, H’s resistance to Byproduct B₀ in a mutualistic H–M community and H’s resistance to M’s antagonistic byproduct A₀ in an exploitative H–M community can mutate independently. In our simulations, these biological parameters are inherited upon cell division, and thus we refer to them as phenotypes or genotypes interchangeably. In most simulations of the simple scenario, only f_P of M mutates, while other phenotypes are held at their bounds whose values are shown in the “Bounds” column of Table 1. In simulations where six phenotypes of the commensal H–M community or seven phenotypes of the mutualistic H–M community could be modified by mutations (e.g., Fig. 7 and Supplementary Fig. 12 for the commensal and Supplementary Fig. 20 for the mutualistic community), these phenotypes start from ancestral values shown in the “Ancestral” column of Table 1. In simulations where seven phenotypes of the exploitative H–M community could be modified by mutations (Supplementary Fig. 21), M’s f_P and H’s sensitivity to M’s antagonistic byproduct A₀ start from ancestral values shown in the “Ancestral” column of Table 1. The other five growth phenotypes (2 maximal growth rates and 3 affinities to B and R) start from evolutionary bound shown in the “Bounds” column of Table 1. This choice allows us to speed up the simulation, since if we initiate these five growth phenotypes from ancestral values, there is not enough biomass in Adult communities to perform heritability check until after more than 1000 cycles.

If a mutation occurs, it could be a null mutation with a probability of (frac{1}{2}). A null mutation reduces f_P, g_Mmax, g_Hmax, A₀, or B₀ to zero, while increases K_MR, K_MB and K_HR to infinity (equivalent to reducing affinities to zero). If a mutation is not null, it modifies each phenotype by ~5–6% on average. Specifically, each phenotype is multiplied by (1 + Δs), where Δs is a random variable with a distribution

Here, s₊ = 0.05 and s₋ = 0.067 are the average percentage by which a mutation increases or decreases a phenotype, respectively (for parameter justifications, see our previous work¹⁵).

At the end of a maturation cycle, the amount of Product P accumulated in the Adult, P(T), is the community function. In some simulations, measurement noise is added to the true P(T) to yield the measured community function. For the simple scenario where only f_P is modified by mutations (e.g., Fig. 6(e–h)), measurement noise is a normal random variable with 0 mean and standard deviation of 100, approximately 10% of the community function of Cycle 1. For the complex scenario where 6 or 7 phenotypes of H and M are modified by mutations (e.g., Fig. 7), measurement noise is a normal random variable with 0 mean and standard deviation of 50. The magnitude of noise is comparable to the ancestral commensal and mutualistic community function, and ~1/4 of the ancestral exploitative community function. Top n_chosen Adult communities with the highest measured function are chosen to be reproduced. Sometimes, more than n_chosen Adults might be needed to obtain n_tot = 100 Newborns for the next cycle if there is not enough biomass in n_chosen Adults.

Chosen Adults are reproduced into Newborns with different methods. If “cell sorting” is used, then the deviation of a Newborn’s total biomass from the target BM_target = 100 is within 2, and the deviation of a Newborn’s species ratio from that of the parent Adult is within 2%. If “pipetting an inoculum of a fixed total biomass” is used, then a Newborn’s total biomass is within a deviation of two from the target BM_target, while its species composition fluctuates stochastically. If “pipetting” is used, then Newborn’s total biomass and species composition both fluctuate stochastically. The dilution fold of each Adult is adjusted, so that the average Newborn community’s total biomass is BM_target over all selection cycles. If a fraction φ_S of a Newborn’s biomass is to be replaced by M or H cells, each Newborn gets on average a biomass of (B{M}_{{{{{{{{rm{target}}}}}}}}}left(1-{varphi }_{S}right)) from its parent Adult community and on average a biomass of BM_targetφ_S from M or H-spiking mix. Specifically, suppose that the biomass of an Adult is BM(T) = M(T) + H(T) where M(T) and H(T) are the biomass of M and H at time T, respectively. If a fraction φ_S of each Newborn’s biomass is to be replaced by a spiking mix, this Adult is then reproduced into n_D Newborns, where

and (leftlfloor xrightrfloor) is the floor (round-down) function. If n_D is larger than n_tot/n_chosen, only n_tot/n_chosen Newborns are kept. Otherwise, all n_D Newborns are kept and as many additional Adults with the next highest functions are reproduced to obtain n_tot Newborns for the next cycle. These Newborns are then topped off with either M or H spiking mixes so that their total biomass is on average BM_target = 100, as described in the next subsection. Note that the fold of dilution of an Adult is calculated based on biomass, a continuous variable. However, the biomass is composed of individual biomass of discrete cells. During reproduction, integer number of cells is distributed into each Newborn community.

Simulating species spiking when only M’s cost f
_P mutates

In the simple scenario where only f_P of M is modified by mutations, phenotypes of all H cells are the same. Within an Adult community, all H cells also have identical individual biomass L_H, because simulations start with H cells of biomass 1 and because growth is synchronous. To mimic reproducing a chosen Adult through pipetting an inoculum of a fixed total biomass into each Newborn with a φ_S-H-spiking strategy, H and M cells from the chosen Adult are randomly assigned to a Newborn community, until its total biomass comes closest to (B{M}_{{{{{{{{rm{target}}}}}}}}}left(1-{varphi }_{S}right)). If φ_S > 0, the number of H cells supplemented to the Newborn community is the nearest integer to (B{M}_{{{{{{{{rm{target}}}}}}}}}{varphi }_{S}{L}_{H}^{-1}). Because integer number of cells is assigned to each Newborn, the total biomass might not be exactly BM_target but within a small deviation of ~2 biomass units.

To mimic reproducing through pipetting, each M and H cell in an Adult community is assigned a random integer between 1 and dilution factor n_D (Equation (12)). All cells assigned with the same random integer are then dealt to the same Newborn, generating n_D Newborn communities. If φ_S > 0, the number of H cells supplemented into each Newborn is a random number drawn from a Poisson distribution of a mean of (B{M}_{{{{{{{{rm{target}}}}}}}}}{varphi }_{S}{L}_{H}^{-1}).

To mimic reproducing through cell sorting, each Newborn receives a biomass of (B{M}_{{{{{{{{rm{target}}}}}}}}}left(1-{varphi }_{S}right)) from its parent Adult. Suppose that the fraction of M biomass in the parent Adult is ϕ_M(T), then M cells from the parent Adult are randomly assigned to the Newborn, until the total biomass of M comes closest to (B{M}_{{{{{{{{rm{target}}}}}}}}}{phi }_{M}(T)left(1-{varphi }_{S}right)) without exceeding it. H cells with a total biomass of (B{M}_{{{{{{{{rm{target}}}}}}}}}left(1-{phi }_{M}(T)right)left(1-{varphi }_{S}right)) are assigned similarly. If φ_S > 0, the number of H cells supplemented to the Newborn community is the nearest integer to (B{M}_{{{{{{{{rm{target}}}}}}}}}{varphi }_{S}{L}_{H}^{-1}) where L_H is the biomass of individual H cell in the parent Adult. Because each of M and H cells had a length between 1 and 2, the actual biomass of M and H assigned to a Newborn could vary from the target by up to 2 biomass units. Consequently, deviations of BM(0) from BM_target and of ϕ_M(0) from parent Adult’s ϕ_M(T) are only a few percent.

Simulating species spiking when both H and M cells evolve

In the more complex scenario, both H and M evolve. We thus need to spike with evolved H and M clones. Additionally, Newborns are spiked with H or M clones from their own lineage as demonstrated in Supplementary Fig. 11a. Below, we describe the simulation code for the experimental procedure (Supplementary Fig. 11a) we simulated.

In all simulations where 6 or 7 phenotypes are modified by mutations, chosen Adults are reproduced through pipetting in a similar fashion as described above. After Newborns are reproduced from a chosen Adult in Cycle C − 1, a preset number of H or M cells are randomly picked from the remaining of this Adult to form H or M-spiking mix for Cycle C. At the end of Cycle C, we choose 10 Adults with the highest functions. Assuming that each chosen Adult is reproduced through pipetting with φ_S-H-spiking strategy, a Newborn receives on average a biomass of (B{M}_{{{{{{{{rm{target}}}}}}}}}left(1-{varphi }_{S}right)) from its parent Adult community and on average a biomass of BM_targetφ_S from H spiking mix generated at the end of Cycle C − 1. Since each chosen Adult usually gives rise to 10 Newborns, the number of cells distributed from the chosen Adult to each Newborn is drawn from a multinomial distribution. Specifically, denote the integer random numbers of cells that would be assigned to 10 Newborns to be {x₁, x₂,…, x₁₀}. If the chosen Adult has a total biomass of BM(T) composed of I_M M cells and I_H H cells (both I_M and I_H are integers), the probability that {x₁, x₂,…, x₁₀} cells are assigned to 10 Newborns, respectively, and x₁₁ cells remain, is

Here, ({p}_{0}=B{M}_{{{{{{{{rm{target}}}}}}}}}left(1-{varphi }_{S}right)/BM(T)) is the probability that a cell is assigned to one of 10 Newborns, p₁₁ = 1 − 10p₀ is the probability that a cell is not assigned to Newborns. Thus, ({x}_{11}={I}_{H}+{I}_{M}-mathop{sum }nolimits_{i = 1}^{10}{x}_{i}) is the number of cells remaining after reproduction, from which H and M cells are randomly picked to generate the spiking mix for Cycle C + 1.

Suppose that the current spiking strategy is φ_S-H, then these 10 Newborns are spiked with H-spiking mix generated in Cycle C − 1. An average of BM_targetφ_S of H biomass is spiked into each Newborn so that the total biomass of Newborns is on average BM_target. Suppose that five H cells from the parent Adult’s lineage are randomly picked at the end of Cycle C − 1, and that they have biomass {L_H1, L_H2, L_H3, L_H4, L_H5}, respectively. The total number of H cells assigned to each Newborn, x_H, is then randomly drawn from a Poisson distribution with a mean of (B{M}_{{{{{{{{rm{target}}}}}}}}}{varphi }_{S}/{overline{L}}_{H}), where ({overline{L}}_{H}=frac{1}{5}mathop{sum }nolimits_{j = 1}^{5}{L}_{Hj}) is the average biomass of the five H cells. Each spiked H cell has an equal chance of being one of the five cells.

Updating spiking percentage based on heritability checks

When the community function landscape is unknown, we can estimate heritability of community function under different spiking percentages through parent–offspring regression. In most simulations (e.g., Fig. 7), heritability evaluation is carried out about every 100 cycles (“periodic heritability check”). In the simulations demonstrated in Supplementary Fig. 17, the average improvement rate in community function is estimated from the chosen Adults over the last 50 cycles. Heritability evaluation is carried out when this average improvement rate becomes negative (“adaptive heritability check”). For both periodic and adaptive checks, heritability evaluation can be postponed until within-community selection improves cell growth sufficiently to provide sufficient biomass for heritability check.

During one round of heritability evaluation, heritability of community function is estimated through parent–offspring community function regression under all candidate spiking strategies (Supplementary Fig. 11b). The current spiking strategy is updated if an alternative spiking strategy confers significantly higher community function heritability.

To evaluate heritability under one spiking strategy, up to 100 Newborn communities are generated under this spiking strategy. After these mature into Adults, their functions are the parent functions. Each Adult parent then gives rise to six Newborn offspring under the same spiking strategy. When the six Newborn offspring mature into Adults, the median of their functions is the average offspring function. When offspring functions are plotted against their parent functions, the slope of the least-squares linear regression (green dashed line in Supplementary Fig. 11b) quantifies the heritability of community function. Heritability of a community function is thus similar to heritability of an individual trait, except that we use median instead of mean of offspring functions, because median is less sensitive to outliers. The 95% confidence interval of heritability is then estimated by nonparametric bootstrap^58,59. More specifically, first, 100 pairs of parent–offspring community functions are resampled with replacement. Second, heritability is calculated with the resampled data. Third, 1000 heritabilities are calculated from 1000 independent resamplings, from which the 95% confidence interval is estimated from the 5th and 95th percentile.

An alternative spiking strategy is considered significantly more advantageous than the current spiking strategy if heritability of the alternative spiking strategy is higher than the right endpoint of the 95% confidence interval of the heritability of the current spiking strategy. If more than one alternative spiking strategies are more advantageous, the one with the highest heritability is implemented to replace the current strategy. Similarly, an alternative spiking strategy is considered more disadvantageous if heritability of the alternative spiking strategy is lower than the left endpoint of the 95% confidence interval of the heritability of the current spiking strategy. When implementing random spiking strategy, the current spiking strategy is updated with a strategy randomly picked from candidate spiking strategies.

Simulating community selection with large population size

When the population size of each community is scaled up by 10 or 100 times (Supplementary Figs. 2 and 18b), the simulation codes described above become inefficient. Instead of tracking the biomass and phenotype of each cell in a large population, we divide the cells into categories and track the number of cells from different categories, where a category is defined by a unique combination of cell biomass and phenotype ranges. In our simulations, the biomass of each cell ranges between 1 and 2, f_P of each M cell ranges between 0 and 1. Since H cells do not mutate, H cells are divided into 100 categories. H cells that belong to category i have a biomass between [1 + (i − 1) × ΔL, 1 + i × ΔL] where ΔL = 10⁻². Since only f_P of M cells are modified by mutations, M cells are divided into 100 × 10⁵ categories. M cells that belong to category (i, j) have a biomass between [1 + (i − 1) × ΔL, 1 + i × ΔL] and f_P between [(j − 1) × Δf_P, j × Δf_P] where Δf_P = 10⁻⁵. Every time f_P of a M cell is modified by mutations, this cell jumps from the current category to a new category determined by its new f_P value.

Similar to simulations with small population sizes, each selection cycle starts with n_tot = 100 Newborn communities. Maturation time T is divided into time steps of length Δτ = 0.05. Over each time step, the growth in cell biomass and the changes in metabolites are simulated in a similar fashion as described above. At the end of each time step, the number of cells to die or to mutate in each category is drawn from a bionomial distribution. If f_P of a M cell is modified by mutation, the mutation effect is drawn from the same distribution as described above: (frac{1}{2}) of mutations reduce f_P to 0 and the other (frac{1}{2}) is randomly drawn from the distribution in Equation (11).

At the end of a maturation cycle, top 10 Adults with the highest functions are chosen. Each then reproduces 10 Newborns via pipetting for the next cycle. The fold of dilution is similarly adjusted, so that the average of Newborn total biomass is BM_target over all selection cycles. From each category of a chosen Adult, the number of cells assigned to a Newborn community is randomly drawn from a multinomial distribution.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com