A multi-scale eco-evolutionary model of cooperation reveals how microbial adaptation influences soil decomposition
Construction of the five-compartment model
Here we explain the construction of the five-compartment model (Fig. 1a). This is step 1 among the four steps described in “Results” section (subsection “Ecosystem dynamics at microsite scale”). We use upper bars in our initial notations to indicate parameters prior to rescaling.
The five-compartment model captures the stochastic processes acting at the level of C, D, M, Z, X entities (molecules, cells) (Fig. 1a) within a microsite. Dynamics of C, D, M, Z, X occur in continuous time. Mt is the number of cells at time t. Ct, Dt, Zt are the numbers of SOC molecules, DOC molecules, and enzyme molecules respectively. Xt is the number of complexes formed by an enzyme molecule binding a SOC molecule. There are constant external sources of SOC and DOC. When a cell dies, a fraction p of the molecules released are recycled into SOC, while the rest is recycled into DOC. A fraction l of dead microbes and deactivated enzymes may be lost due to leaching.
We denote by α the structural cost of a cell, which is the equivalent in number of DOC molecules of one cell (without storage). We denote by (alpha ^{prime}) the energetic cost of a cell, which is the number of DOC molecules consumed to produce the energy needed for the synthesis reactions involved in the production of one cell. We denote by β the equivalent in number of DOC molecules of one SOC molecule, and the structural and energetic cost of producing one molecule of enzyme by ρ and (rho ^{prime}), respectively. We assume that the energetic costs are carbon released by cells as CO2 (cell respiration) that diffuses out of the system instantly. We define the biomass production fraction and enzyme allocation fraction as
$${bar{gamma }}_{M}:=frac{1}{alpha +alpha ^{prime} },quad {bar{gamma }}_{Z}:=frac{1}{rho +rho ^{prime} }.$$
(1)
The event times are given by independent exponential random variables whose parameters are defined by event rates (Supplementary Tables 2–4). These event rates give an approximation of the average frequency of each event. The rates of cell growth and enzyme production depend on the trait φ. Once a cell has doubled its initial size, reproduction occurs by releasing the mother cell at its initial size, and the daughter cell at its same size. The cell must therefore take up and store its structural and energetic cost, ((alpha +alpha ^{prime} )), in DOC molecules in order to reproduce. We denote N the number of uptake events before reproduction. The number of DOC molecules taken up at each uptake event is then ((alpha +alpha ^{prime} )/N), hence the notation ({{bf{1}}}_{{Dge (alpha +alpha ^{prime} )/N}}) which equals 1 if (Dge (alpha +alpha ^{prime} )/N), and 0 otherwise. The same notation, ({{bf{1}}}_{{Dge rho +rho ^{prime} }}), is used for the production event of an enzyme molecule. Uptake is stochastic, but reproduction is deterministic, which means that when a cell has performed N uptake events, it reproduces with probability 1. A larger N means a larger number of uptake events between 2 reproduction events, which also means less DOC molecules taken up at each uptake event. The model tracks the dynamics of the number of cells, SOC, DOC, enzyme molecules, and also of the DOC stored in each cell.
Enzyme–substrate complexes form at rate ({bar{lambda }}^{k}) as one enzyme molecule (e.g. cellulase) bind one SOC molecule (e.g. cellulose). A complex may either dissociate (with no decomposition) at rate ({bar{lambda }}_{-1}^{varepsilon }), or react at rate ({bar{mu }}^{varepsilon }) and convert the molecule of SOC into β molecule of DOC while the enzyme is released and free again to react with new molecules of SOC (Supplementary Table 3).
System size k does not appear in this system of equations, yet it enters the volume-dependent parameters of the model, IC, ID, θ, and Km. We denote V the unit volume of soil that contains on average one microbial cell and the corresponding equilibrium of carbon mass of SOC, DOC and enzymes. The system size k is the number of well-mixed unit volumes in one microsite, which determines the number of cells sharing SOC, DOC, and enzymes. The volume of one microsite is therefore k × V. Increasing k amounts to increasing microsite volume, the number of cells sharing resources in one microsite, the amount of resources per microsite, and the volume-dependent parameters, such as the amount of SOC entering a microsite per unit of time, IC. In our analysis, the unit volume V is fixed, and we vary k to investigate the effect of microsite volume on the system’s eco-evolutionary dynamics. With very large k, the hybrid model can be approximated by a fully deterministic model which takes the form of a system of four ordinary differential equations (see Supplementary Information 3.3 and Supplementary Fig. 1), similar to the microbial decomposition model first introduced by Allison et al.53. However, empirical data suggest that k is of the order of 10–10054. When k = 1, there is only one cell in the microsite, which volume is V defined as the unit soil volume expected to contain a single cell. A value of k greater than 1 means that each microsite contains k cells and k times the amount of SOC, DOC and enzyme molecules of 1 cell; thus, the microsite volume is k × V, and volume-dependent parameters are rescaled by k. Specifically, there are four volume-dependent parameters: the external input of C, ({bar{I}}_{C}^{k}), the external input of D, ({bar{I}}_{D}^{k}), the half-saturation constant of DOC uptake, ({bar{K}}_{m}^{k}), and the encounter intensity of two given SOC and enzyme molecules, ({bar{lambda }}^{k}). The external inputs increase proportionally with the volume of the microsite, while the encounter intensity of two given molecules in a microsite decreases as its volume increases. The half-saturation is inversely proportional to the affinity between a given cell M and a given DOC molecule d, which decreases with increasing microsite volume. We thus obtain the following scaling relationships:
$${bar{I}}_{C}^{k}=k{bar{I}}_{C}, {bar{I}}_{D}^{k}=k{bar{I}}_{D}, {bar{lambda }}^{k}=frac{lambda }{k} {rm{and}} {bar{K}}_{m}^{k}=k{bar{K}}_{m}.$$
(2)
In our simulations, we generally assume that k is equal to 10, to match the empirical observation that (cells) in soil habitat tend to interact with 10 to 100 other cells at all time54.
Derivation of the hybrid model
In Supplementary Information 3, we present the next two steps (2 and 3 in “Ecosystem dynamics at microsite scale” of “Results” section) to derive the hybrid model on which all our results are based. In Supplementary Information 3.1, we explain how the dynamics generated by the five-compartment model can be captured with a reduced model with four-state variables (step 2). In Supplementary Information 3.2, we explain how the stochastic-deterministic PDMP model can be derived from the stochastic four-state variable model (step 3). In the hybrid model (step 4), only cell death remains stochastic, and cell dynamics is measured in unit of number of individuals (M), while other entities are now in carbon mass unit. The rescaled SOC, DOC and enzyme abundances are denoted with lower case letters c, d, and z.
Simulation algorithm for the hybrid spatial model
One technical benefit of the hybrid model is its much greater computational tractability. Here we describe the algorithm used to perform simulations of the hybrid model. We ran the model on single microsite to produce the simulations reported in Fig. 2. We ran the model on a 10 × 10 lattice of microsites for the subsequent figures. The algorithm is based on the Gillepsie algorithm55 as used in Champagnat et al.37, Fournier and Méléard56, which straightforwardly extends to the simulation of PDMPs.
To couple PDMP models across microsites, we account for the DOC and dispersal of cells between adjacent microsites. The DOC diffusion between microsites is modeled by approximating a continuous diffusion with a Euler scheme in which time is discretized with a fixed time step interval, τdiff. τdiff is chosen sufficiently small to provide a fine enough discretization of the DOC diffusion.
A simulation starts with a given amount of M, z, c, and d in each microsite at time t = 0, while the initial amount of DOC stored within each cell is determined uniformly at random. Two stochastic events (death of a cell) may not occur at the same time. Assume that the process has been computed until time ti; to continue the computation to time ti+1, we proceed as follows.
First, we simulate T, an exponential random variable with parameter r(ti) = dMM(ti), which corresponds to the death rate of the total cell population at time ti (M(ti) being the total number of cells on the entire lattice). We then compute
$${t}_{i+1}:={t}_{i}+min left(T,{tau }_{{rm{diff}}}right).$$
To obtain c(ti+1), d(ti+1), and z(ti+1) in each microsite at time ti+1, and the variation in amount of DOC stored within a cell in the corresponding microsite, we use a Euler scheme that solves the dynamical system
$$left{begin{array}{l}hskip -136pt dot{c}(t)={I}_{C}-{l}_{C}c-theta zc,\ dot{d}(t)={I}_{D}-{l}_{D}d+theta zc+(1-l){d}_{Z}z-{V}_{max }frac{d}{{K}_{m}{,}+{,}d}{omega }_{M}M,\ hskip -84pt dot{z}(t)=varphi {gamma }_{Z}{V}_{max }frac{d}{{K}_{m}{,}+{,}d}{omega }_{M}M-{d}_{Z}z,\ hskip -93pt dot{Delta }(t)=(1-varphi ){gamma }_{M}{V}_{max }frac{d}{{K}_{m}{,}+{,}d}{omega }_{M},end{array}right.$$
in each microsite between ti and ti+1, where M is the number of cells in the microsite at time ti, Δ gives the amount variation of DOC stored within a cell, Δ(ti) = 0 and the other initial conditions are the biomass of c, d, z in the corresponding microsite at time ti.
Note that, within a microsite, the variation of stored DOC is the same for all cells and corresponds to Δ(ti+1). Hence, this amount is added to the amount of DOC stored within each cell living in the corresponding microsite. If, for a cell j, the resulting amount Sj(Ti) + Δ(ti+1) is over ωM, a new cell appears. The amount of stored DOC within the new cell and the mother cell is then updated: Sj(ti+1) = (Sj(Ti) + Δ(ti+1) − ωM)/2. Otherwise, Sj(ti+1) = Sj(Ti) + Δ(ti+1). To determine the position of the new cell, the following steps are taken:
A uniform random variable ϑ1 in [0, 1] is simulated.
If ϑ1 More
