### 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^{15}. 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{dR}{dt}=-{c}_{{RM}}{g}_{M}M-{c}_{{RH}}{g}_{H}H$$

(1)

$$frac{dB}{dt}={g}_{H}H-{c}_{{BM}}{g}_{M}M$$

(2)

$$frac{dP}{dt}={f}_{P}{g}_{M}M$$

(3)

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

(4)

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

(5)

where

$${g}_{H}(R)={g}_{{Hmax}}frac{R}{R+{K}_{{HR}}}$$

(6)

$${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)$$

(7)

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)^{51}) 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

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

(8)

where larger *B*_{0} 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

$$frac{dwidetilde{A}}{dt}={r}_{A}{g}_{M}left(1-{f}_{P}right)M$$

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}

$$A=widetilde{A}/{r}_{A}$$

so that

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

(9)

We also need to modify the growth rates for H:

$${g}_{H}={g}_{H}(R)={g}_{{Hmax}}frac{R}{R+{K}_{{HR}}}frac{{A}_{0}}{A+{A}_{0}}$$

(10)

where larger *A*_{0} 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 (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^{15}. 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^{4} 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^{−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 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^{−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 held^{15}. 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