Global dynamics of microbial communities emerge from local interaction rules

Abstract

Most microbes live in spatially structured communities (e.g., biofilms) in which they interact with their neighbors through the local exchange of diffusible molecules. To understand the functioning of these communities, it is essential to uncover how these local interactions shape community-level properties, such as the community composition, spatial arrangement, and growth rate. Here, we present a mathematical framework to derive community-level properties from the molecular mechanisms underlying the cell-cell interactions for systems consisting of two cell types. Our framework consists of two parts: a biophysical model to derive the local interaction rules (i.e. interaction range and strength) from the molecular parameters underlying the cell-cell interactions and a graph based model to derive the equilibrium properties of the community (i.e. composition, spatial arrangement, and growth rate) from these local interaction rules. Our framework shows that key molecular parameters underlying the cell-cell interactions (e.g., the uptake and leakage rates of molecules) determine community-level properties. We apply our model to mutualistic cross-feeding communities and show that spatial structure can be detrimental for these communities. Moreover, our model can qualitatively recapitulate the properties of an experimental microbial community. Our framework can be extended to a variety of systems of two interacting cell types, within and beyond the microbial world, and contributes to our understanding of how community-level properties emerge from microscopic interactions between cells.

Fig 1. A mathematical framework to scale up from molecular-level properties, to individual-level properties, to community-level properties.

We previously measured the local interaction rules for a cross-feeding community (B) and showed that these can be derived (arrow 1) from the molecular mechanisms of the interaction (A). Here we developed a mathematical model that derives community-level properties (C) either from measured local rules (arrow 2) or directly from the underlying molecular mechanisms (arrow 3). (A) The community consists of two types of Escherichia coli: ΔP is unable to produce the amino acid proline and ΔT is unable to produce the amino acid tryptophan. Cells exchange amino acids with the environment through active uptake (with rate r^u) and passive leakage (with rate r^l). Amino acids are exchanged between cells through diffusion in the environment (with rate D). All rates differ between the two amino acids. (B) Local interaction rules can be fully described by two fundamental quantities: the size of the interaction neighborhood (r_ΔT, r_ΔP, left); and the growth function of a cell (characterized by $\widehat{\mu }$, right). Each dot corresponds to the measured growth rate of a single cell, n = 2610 for ΔP and n = 2162 for ΔT, the line shows the result of a linear regression, data reproduced from [11]. (C) We derive analytical expressions for steady state community-level properties, such as the equilibrium frequency of the two types, their spatial arrangement, and growth rate.

Fig 2. Pair approximation allows for a quantitative description of spatial systems.

(A) Pair approximation assumes that a system of entities living in space and interacting with others close by (left side) can be fully described by tracking the number of all pairwise links between entities (right side). For example, a link counted in N_X←Y indicates that focal cell X interacts with its neighbor Y. (B-D) A system is fully described by: (B) the interaction neighborhoods of both types, characterized by the neighborhood sizes r_A and r_B; (C) the replication neighborhood, assumed to be identical to the smallest interaction neighborhood; and (D) the growth functions of both types, characterized by the maximum growth rates ${\widehat{\mu }}_A$ and ${\widehat{\mu }}_B$. (E) Pair approximation has three dynamical variables that describe the global composition, P(A), and local composition P(B|A, r_A) and P(A|B, r_B) of the system.

Fig 3. Community-level properties depend on local interaction rules.

(A) The global composition of the community (global frequency of A, given by Eq 1) primarily depends on the ratio of maximum growth rates of the two types. If the neighborhood sizes of both types are large (r_A = r_B = 10, purple line), the equilibrium frequency approaches that of a well-mixed system (black line, given by Eq 2). If the neighborhood sizes are small (r_A = r_B = 3, green line), the type with the higher growth rate attains a higher frequency than in a well-mixed system. When either type becomes fixed (i.e., P(A) = 0 or P(A) = 1) cells can no longer reproduce. In natural communities, where populations sizes can change, this would lead to the collapse of the community; however, in our model we assume a constant population size, so in this case the community persists even though there is no longer any turnover of cells (i.e., community productivity drops to 0). (B) The neighborhood size affects the global composition of the community, when the two types have different maximum growth rates. The solid lines show the frequency of type A for spatial systems where r_A is held constant at 3, while r_B is increased. The equilibrium frequency varies as r_B increases: for larger r_B, the frequency in a spatial system (solid lines) moves closer to the frequency in a well-mixed system (dotted lines). This result holds when the types have different maximum growth rates (red and green curves) but not when they have equal maximum growth rates (black line). (C) The neighborhood size strongly affects the local composition of the community. Here the two types have the same maximum growth rate and neighborhood size. For both types, the local frequency of the partner cells is much lower than the global frequency, when the neighborhood size is small. (D) The neighborhood size strongly affects community productivity. The community productivity (i.e. the gross production of new biomass) is smaller when cells have smaller interaction neighborhoods, because the local frequency of the partner around each cell is lower. This effect is stronger in communities where the types have different maximum growth (the red and green curves are below the black curve).

Fig 4. Small interaction neighborhoods can lead to community collapse.

(A) When the neighborhood size is small (r_A = r_B = 3, purple line), asymmetric communities, where the cell types have different maximum growth rates, have low productivity. When the asymmetry is too large, communities cannot grow (i.e. community productivity is 0) in a spatially structured environment even though they could grow in well-mixed environments. When the neighborhood size is large (r_A = r_B = 30, orange line), the productivity of the spatial community is close to that of the well-mixed community even when there is an asymmetry. (B) Increasing the neighborhood size of just one the types is enough to increase the productivity of spatial communities and prevent collapse. r_A is held constant at 3, while r_B is increased.

Fig 5. The neighborhood size strongly affects the equilibrium frequency when the growth rate of a cell depends on the absolute number of partner cells.

(A) When the growth rate depends on the frequency of the partner type within the interaction neighborhood (frequency dependence), the equilibrium frequency of type A is almost completely determined by the ratio of the maximum growth rates (Eq 1). (B) When the growth rate depends on the number of cells of the partner type within the interaction neighborhood (density dependence), the equilibrium frequency of type A depends both on the ratio of the maximum growth rates and on the ratio of the neighborhood sizes (see Equation 32 in S3 Text). The slow growing type can still dominate the community when it has a much larger interaction range (blue region, top right). In both panels r_A = 10.

Fig 6. Our mathematical framework can explain experimentally observed community properties.

The model parameters for the experimental community were calculated from the biophysical model (S2 Table). (A) Experimental communities approach a stable frequency of ΔT. Individual communities (thin lines, n = 22) and their average value (thick line) are shown. The initial frequency of ΔT in the communities is determined by the initial number of cells that enter the microfluidic growth chambers and the subsequent growth before the start of image acquisition; as a result, it cannot be controlled experimentally. (B) Pair approximation predicts a unique stable equilibrium. The dynamical equations (Equations 12–14 in S2 Text) were solved numerically starting from different initial frequencies. (C) Cellular automaton simulations also reach a unique stable equilibrium. (D) The observed equilibrium frequency of ΔT is consistent with the model predictions. Data: P(A) = 0.19 (95% confidence interval (CI): 0.17–0.20), Pair-approximation: P(A) = 0.20, and cellular automaton: P(A) = 0.18 (CI: 0.18–0.19, evaluated after 100’000 time steps). The difference between the model prediction and data is less than 8%. (E) The frequency of the partner cell type within the interaction neighborhood (local frequency) is lower that the global frequency because cells are surrounded by their own offspring. Dots show measurements for 21 replicate communities, bar shows mean value. Pair approximation predicts (Eq 3) a decrease in frequency of 0.99 for ΔP and 0.89 for ΔT, the experimental values are 0.92 (CI: 0.84–0.99) and 0.85 (CI: 0.82–0.88), respectively. The difference between the model prediction and data is less than 13%. (F) The average growth rate of the community is reduced due to cell clustering. An in-silico analysis (see Methods) shows that the growth rate in clustered communities, with experimentally observed spatial arrangements, is reduced by a factor of 0.87 (CI: 0.85–0.90) compared to randomized communities, where cell clusters were disrupted (Data reproduced from [11]). Pair approximation predicts a decrease by a factor of 0.92 (S3 Table); the difference between the model prediction and data is less than 6%.