Steering ecological-evolutionary dynamics to improve artificial selection of microbial


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 work15. 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))15:







$$frac{dH}{dt}={g}_{H}H-{delta }_{H}H$$


$$frac{dM}{dt}={g}_{M}left(1-{f}_{P}right)M-{delta }_{M}M$$





$${g}_{M}(R, B)={g}_{{Mmax}}frac{{R}_{M}{B}_{M}}{{R}_{M}+{B}_{M}}left(frac{1}{{R}_{M}+1}+frac{1}{{B}_{M}+1}right)$$


and RM = R/KMR and BM = B/KMB. Unless otherwise specified, landscapes in this paper are obtained by integrating Equations (15) from t = 0 to t = 17.

Equation (1) states that Resource R is depleted by biomass growth of M and H, where cRM and cRH 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 (cBM amount of B per unit of M biomass). Equation (3) states that Product P is produced as fP 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)51) discounted by (1 − fP) due to the fitness cost of making Product, and decreases at the death rate δM. In the Monod growth model (Equation (6)), gHmax is the maximal growth rate of H and KHR is the R at which gHmax/2 is achieved. In the Mankad and Bungay model (Equation (7)), KMR is the R at which gMmax/2 is achieved when B is in excess; KMB is the B at which gMmax/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 (15) and (7), but Equation (6) is replaced with

$${g}_{H}(R)={g}_{{Hmax}}frac{R}{R+{K}_{{HR}}}exp left(-frac{B}{{B}_{0}}right)$$


where larger B0 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 (15) and (7), we then need to add an equation that describes the dynamics of A


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


so that



We also need to modify the growth rates for H:



where larger A0 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−2 ≤ ϕM(0) ≤ 0.99 and (1{0}^{-2} le {overline{f}}_{P}(0) le 0.99) with a mesh size of ΔϕM(0) = 10−2 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−3 and ({{Delta }}{overline{f}}_{P}(0)=1{0}^{-4}).

To calculate the species composition attractor, we integrated Equations (15) 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 work15. 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 (gMmax) must exceed the maximal growth rate of H (gHmax), and M’s affinity for Byproduct (1/KMB) 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 (cRH = cRM) 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 104 total biomass. The evolutionary bounds are set, such that evolved H and M could coexist for fp < 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−4 per genome per generation (calculated from Fig. 3 of Levy et al.52), and mutations that reduce growth rate occur at a rate of 10−4 ~ 10−3 per genome per generation53,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−3 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 held15. Among phenotype-altering mutations, tens of percent create null mutants, as illustrated by…

Read More:Steering ecological-evolutionary dynamics to improve artificial selection of microbial

Notify of
Inline Feedbacks
View all comments

Get a Copy of Free Newsletter

Subscribe to our mailing list and get interesting stuff and updates to your email inbox.

Thank you for subscribing.

Something went wrong.