Bacterial social interactions drive the emergence of differential spatial colony structures

Background Social interactions have been increasingly recognized as one of the major factors that contribute to the dynamics and function of bacterial communities. To understand their functional roles and enable the design of robust synthetic consortia, one fundamental step is to determine the relationship between the social interactions of individuals and the spatiotemporal structures of communities. Results We present a systematic computational survey on this relationship for two-species communities by developing and utilizing a hybrid computational framework that combines discrete element techniques with reaction-diffusion equations. We found that deleterious interactions cause an increased variance in relative abundance, a drastic decrease in surviving lineages, and a rough expanding front. In contrast, beneficial interactions contribute to a reduced variance in relative abundance, an enhancement in lineage number, and a smooth expanding front. We also found that mutualism promotes spatial homogeneity and population robustness while competition increases spatial segregation and population fluctuations. To examine the generality of these findings, a large set of initial conditions with varying density and species abundance was tested and analyzed. In addition, a simplified mathematical model was developed to provide an analytical interpretation of the findings. Conclusions This work advances our fundamental understanding of bacterial social interactions and population structures and, simultaneously, benefits synthetic biology for facilitated engineering of artificial microbial consortia. Electronic supplementary material The online version of this article (doi:10.1186/s12918-015-0188-5) contains supplementary material, which is available to authorized users.

1 Computational Framework 1.1 Mechanical Forces Time Evolution Equations: In a two-dimensional space, a bacterial cell (i) can be modeled as a rigid rod with center of mass coordinates ( r i ), center of mass velocities ( v i ), orientation with respect to some specified axis (φ i ), density (ρ i ), and length (l i ). Here we assume that the intracellular density is constant for all cells (ρ i = ρ) so that the mass of a given cell can be calculated using: With this assumption, we can write down the fundamental equations that govern cellular motion, namely Newton's laws, for the force and torque of a rigid body moving in a viscous substance: where the brackets ({}) denote a function of all coordinates and velocities and β denotes the drag from the surrounding viscous medium. The torque is calculated from the cross product of the the lever arm and the corresponding force.
The unique aspects of simulating a bacterial population lie in the facts that cells grow and divide, that the particle number of the system is not conserved, and that l i is a function of time. In order to successfully model the dynamics of such a system, the division time for a given bacterial cell is chosen as the typical time scale of the system and the diameter of the cell is chosen as the basic spatial length scale. Here we introduce two corresponding new S2 variables: where r c and t c have dimensions of length and time respectively. With this change of variables, length in the system can be written as Our equations of motion can thus been written as where the forces and torques are also functions of dimensionless variables. Such nondimensionalization allows to arrange the terms in a way that emphasizes the contribution of the inertia term to the overall motion as To implement simulation, we need to specify the force and torque; however, we notice there are now two clear timescales for the system, namely 1/β which governs that amount of time for a given particle to change velocity and another timescale associated with the force which governs the typical time for an interaction. The detailed forms are presented below.
Hertzian Forces: The following equation for the magnitude of contact force was used to describe the expansion of a bacterial colony on a solid substrate (S1) .
where E is a constant describing the rigidity of the rod, d is the interaction diameter of the rod, and h ij = d − | r ci − r cj |. Here, r ci and r cj denote the closest points between the two rods. In addition, the direction of the force is in the same direction as r ci − r cj . Using this force form, the equations of motion for rods in contact can be expressed as: where the sum over contacts includes cells that have a nonzero contact force (i.e. h ij > 0).
Similarly, the equations for torque can be written as: From the above equations, we can explicitly see the two timescales in the problem, namely: Alternatively, we can use the combined parameter ζ = βρ instead. For bacteria (E. coli)

S4
growing on a solid substrate, the estimated values are (S1) To estimate t 1 and t 2 , we also need the density of a bacterial cell which can be approximated as ρ ≈ 10 −15 kg 3 × 10 −6 m ≈ 3 × 10 −10 kg/m (14) From the above estimates, we can obtain the two time scales as from which we find that the timescale for the velocity to change is much shorter than the timescale for cellular interaction. Therefore, for a practical simulation, it is thus appropriate to neglect the inertial term and choose a longer physical time step, which results in the first order equations for force and torque below: Notably, the timescale for interactions (t 1 ) also provides a point of reference for choosing the integration timestep during simulations because it tells the order of magnitude for time that two rods will contact. Thus, to ensure successful simulation of cellular collisions, we shall S5 choose a timestep that is smaller than t 1 . A reasonable choice would be: where we used 0.5 hr (typical division time for E. coli in exponential growth) for t c . Notice that this corresponds to the interval in real time: It is important to note that this separation of timescales for t 1 and t 2 is a direct consequence of slow movement of cells on a solid substrate. For liquid environments, the timescales may be of similar magnitude during which eliminating the acceleration term in the equations of motion will not be appropriate.
Numerical Integration: The Euler method was used to integrate the time evolution of the force and torque equations. Due to the differences in growth rates of the cells with different interactions, we found it helpful to choose a variable timestep instead of choosing a minimal fixed time interval. The variable timestep was determined by ensuring that the maximal movement of any cell was a small fraction of the cell diameter. This strategy enabled substantially faster simulations for the case of populations with slow overall growth rates due to deleterious interactions (i.e. competition). The force calculations were implemented utilizing the following procedures: (1) cells were classified to a spatial grid; (2) cells in nearest neighbor grids were tested for possible overlaps by finding the closest points for each respective pair of cells; (3) cells that overlapped were assembled into a list; (4) forces and torques at the closest contact points were calculated in parallel and output to a force list; (5) forces and torques were combined for each cell. Some sample code from (S2) was used to aid in development.

S6
where D c is the diffusion constant, α c and f (ρ c , c) determine production (or consumption) by cells, ρ c is the density of cells that interact with the chemical, and β c determines degradation.
Ω is a two-dimensional rectangular region and ∂Ω is the boundary; c 0 is the concentration on the boundary. For the shared nutrition source (n), we assume that there is no autodegradation (i.e. β n = 0). We also re-scale the nutrition concentration so that its value on the boundary equals 1. For the chemicals that mediate cellular interactions, we impose a reactive boundary (c 0 = 0). Similarly, we scale the concentration so that α c = 1. Moreover, in all cases, the two dimensional region simulated is much larger than the colony size, so that S7 the reactive boundary condition for the diffusible chemicals should have negligible effects.
Nutrition: Based on the Monod function (S3), the exact equation for nutrition is given where n is nutrition, ρ is the total density of all species at a grid point in space, defined as the area of all cells within a grid divided by the area of the grid. Parameters were chosen to approximately reproduce the experimentally determined width of actively growing cells of bacterial colonies (30-40µm) (S4), which are in agreement with previous simulations (S1).
The non-dimensional values are D n = 250, α n = 1, and κ = 0.333. The boundary of the simulated area serves as a constant source of nutrient.
Chemicals for Interactions: For chemicals other than nutrition, its production by a given species is assumed constant. Therefore, the corresponding dimensionless reactiondiffusion equation has a form of where c is the chemical, ρ c is the density of the producer species at a grid point in space. In this work, the diffusion (D c ) and degradation (β c ) parameters (D c = 1000 and β c = 10) were chosen so that cellular interaction range is on the order of few spatial grid points. Notably, for a well occupied spatial grid in the absence of diffusion, the steady state value for a given chemical follows S8 which provides an estimate for the level of chemical within the densely packed region of the expanding colony.
Numerical Integration: The Euler method was employed for simulating all of the reaction-diffusion equations. A discretized version of the diffusion operator from a Taylor expansion to second order in grid size was used. A grid size of 5 was used throughout simulations, which is comparable to the max length for a single cell. The simulated region was a square with a side length of 1000, which is much larger than the typical final radius of the colonies (≈ 200) Table S2: Summary of Diffusion Equations. The value for nutrient sensitivity is taken with respect to the maximal concentration at the simulation boundary. Degradation, consumption, and diffusion parameters were chosen to give a desired length scale for active cell growth and interactions.

Description Equation
Nutrients

Cell Growth and Division
Cellular growth follows the following equation: S9 where l i is the length of the cell (i), A i is the cell area, g is the maximal growth rate, κ is the parameter associated with nutrient consumption (0.333 as in section 1.2), ξ is an interaction parameter, and T is the concentration of diffusible chemical produced by the other species.
This specific form for elongation is similar to that found in (S1).
As discussed previously, we expect the maximal concentration of the diffusible chemicals to be approximately 1/β c (or 1/10 in the simulations) for a well-packed spatial grid of producing cells. We thus chose the interaction parameter ξ so that the exposure of the cell to the maximal concentration resulted in either cessation of cell growth (toxin) or doubling of cell growth (public good) ( ξ = +10 or −10 respectively). To ensure that the cells in the population never decrease in length, the growth rate was truncated to zero for any possible negative rates that arose.
The growth constant g was determined by considering the case of nutrient saturation (n = 1) at which the doubling time is approximately one time unit in the absence of other diffusible chemicals, i.e., Calculating the time for a cell to grow from a length of 1.5 to 4.0 for division gives: S10 Finally, by setting the above time as 1 and plugging in the parameter value of κ, we obtained For cellular division, we assumed that each cell has a length that increases over time until it reaches a specified division length. The division length for a cell was generated randomly from a Gaussian distribution with a mean of 4 and a standard deviation of 0.3. Any values for division length outside of [3.1, 4.9] were truncated to the extreme values. Notably, the total of cell length plus the mechanical interaction range is conserved during division so that newly divided cells did not overlap. Each daughter cell received a fraction of the parent length between 0.4 and 0.6 (adding up to one in total) randomly.

S11
To interpret the results obtained from the simulation of our community modeling framework, we considered a simplified case of the system-a well-mixed model. As shown below, the model consists of three variables, including nutrition (n) and the populations of two cellular species (u and v).
The parameters are nutrient supply (n 0 ), degradation or flow rate (D), nutrient consumption (m), consumption to growth ratio (γ), nutrient sensititivey (κ), and interspecies interaction strength (ξ 1 and ξ 2 ). Notice that the first term on the right hand side of the equations for both u and v consists of nutrient consumption and interaction terms; the second term results form death or flow out of the system. As we are primarily interested in the system's steady states and their stabilities, we performed the change of variables using S12 which results in the following equations after dropping the bars and rescaling the parameters: The steady states of the system can then be derived by setting the above equations to zero. When considering the solutions to the equations, there are three possible types: u and v both are zeros; one of them is zero; and neither of them is zero. For the case where both are zero, the steady state solution is: For the case where v is extinct and u is nonzero, the steady state solution is: Similarly, the steady state solution for the case of u is extinct and v is nonzero can also be expressed. Furthermore, for the case where both are nonzero, the detail expression can also be derived (It is not listed here as the expression is quite cumbersome).
Next we studied the stability of the system at the steady states for given types of interaction. As the extinction of the both species is trivial, we primarily focused on the cases S13 where at least one species survives, which require the following condition: for positive κ.
We found that, as long as the above condition is satisfied, the state with both u and v extinct is linearly unstable as desired. Evaluating the Jacobian for the other two single species fixed points shows that the u nonzero, v extinct state is stable as long as ξ 2 > 0.
Similarly, the u extinct, v nonzero state is stable as long as ξ 1 > 0.
A numerical exploration of the number of stable steady states in Mathematica confirmed the above analytical results. For parameter values of α = 2 and κ = 1/3, we explored the impacts of social interaction on the system's steady state properties by varying the values of ξ 1 and ξ 2 from −0.2 to 0.2. Here, the values for ξ 1 and ξ 2 were kept relatively small to reflect the realistic case where cellular interactions from toxins or public goods are a secondary effect compared with that from nutrition. As shown in Figure 6 in the main text, the different values for ξ 1 and ξ 2 lead to different outcomes (extinction, bistable extinction, and coexistence). The axes in the plot are marginal cases (with zero eigenvalues), however, the test cases through simulations show extinction (SI Figure S12). For the figures, Figure   6C in the main text and Supporting Figure S12, interaction parameters (ξ 1 and ξ 2 ) were assigned 0, −0.1, or 0.1 depending on the interaction type being tested. S14 Figure S1: The fraction of green cell lineages as a function of the total cells of the communities that are neutral, mutualistic, and competing. It is a supplement to Figure 4C in the main text. Figure S2: Statistical analysis of bacterial communities that are compared according to time evolution. Six simulation runs were performed for each interaction type in order to obtain the statistics. S16 Figure S3: Statistical analysis of the role of social interactions in determining community structure for the case of a medium initial density (64 cells, 1:1 ratio). S17 Figure S4: Statistical analysis of the role of social interactions in determining community structure for the case of a medium initial density (16 cells, 1:1 ratio). Figure S5: Changing the initial density has little effect on colony roughness for both amensalism and competition. S19 Figure S6: Representative structures emerged from amensal communities with three different initial densities: high (256 cells), medium (64 cells) and low (16 cells). As the initial cell density drops, the number of total green cells (victim species) increases, eventually resulting in coexistence.

Supporting Figures
S20 Figure S7: Statistical analysis of the role of social interactions in determining community structure for the case of a high initial density (256 cells) and a 1:7 relative abundance.
S21 Figure S8: Statistical analysis of the role of social interactions in determining community structure for the case of a high initial density (256 cells) and a 1:3 relative abundance.
S22 Figure S9: Statistical analysis of the role of social interactions in determining community structure for the case of a high initial density (256 cells) and a 3:1 relative abundance.
S23 Figure S10: Statistical analysis of the role of social interactions in determining community structure for the case of a high initial density (256 cells) and a 7:1 relative abundance. Figure S11: Fraction of green cells as a function of total cells in mutualistic communities with different initial relative abundance. Mutualism tends to decrease deviations in the relative abundance away from 1:1. Figure S12: Simulations of well-mixed two-species communities with different social interactions. For control, the green cell fraction is determined purely by the initial conditions. For commensalism, the red cells always go extinct. For amensalism, the green cells always go extinct. For mutualism, equal species abundance is achieved. For competition, one of the species goes extinct depending on initial conditions. For parasitism, the red cells always go extinct. This figure supplements Figure 6 in the main text.
S26 Figure S13: Time course images for cell density using the well-mixed initial conditions with high cell density and Neutralism. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point. Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S27 Figure S14: Time course images for cell density using the well-mixed initial conditions with medium cell density and Neutralism. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point. Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S28 Figure S15: Time course images for cell density using the well-mixed initial conditions with low cell density and Neutralism. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point. Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S29 Figure S16: Time course images for cell density using the well-mixed initial conditions with high cell density and Mutualism. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point. Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S30 Figure S17: Time course images for cell density using the well-mixed initial conditions with medium cell density and Mutualism. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point. Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S31 Figure S18: Time course images for cell density using the well-mixed initial conditions with low cell density and Mutualism. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point. Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S32 Figure S19: Time course images for cell density using the well-mixed initial conditions with high cell density and Competition. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point. Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S33 Figure S20: Time course images for cell density using the well-mixed initial conditions with medium cell density and Competition. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point.
Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S34 Figure S21: Time course images for cell density using the well-mixed initial conditions with low cell density and Competition. The total cell density (top row) shows a densely packed quasi-circular colony with a maximum height corresponding to approximately 10 cells packed in one spacial grid (5 µm by 5 µm). The values for cell density were calculated by counting the number of cells with a center of mass at a specified spatial grid point. Since cells with a center of mass at one grid point may extend into a neighboring grid, the values were averaged over nearest neighbors. The cell densities for individual species (i.e. green or red) are shown in rows 2 and 3. The mesh lines shown in the image correspond to approximately 6 spatial grids in the simulation. The boundary of the grid was chosen for display purposes; the actual simulation boundary is more distant from the expanding colony to minimize effects from the boundary conditions. S35 Figure S22: Comparison of the cell number density plots and the single cell image for a top down view of the final state for a simulation with well-mixed initial conditions and Neutralism. The two-dimensional nature of the simulations allows the local density to be inferred by considering the color of the single cells in the images of the far right column.
S36 Figure S23: Comparison of the cell number density plots and the single cell image for a top down view of the final state for a simulation with well-mixed initial conditions and Mutualism. The two-dimensional nature of the simulations allows the local density to be inferred by considering the color of the single cells in the images of the far right column.
S37 Figure S24: Comparison of the cell number density plots and the single cell image for a top down view of the final state for a simulation with well-mixed initial conditions and Competition. The two-dimensional nature of the simulations allows the local density to be inferred by considering the color of the single cells in the images of the far right column.     Table S9: Run statistics for communities simulated using the well-mixed model. The mean and standard deviation of the fraction of green cells were obtained using 1000 runs of simulations with random initial conditions (each of the species concentrations is between 0 and 1 and nutrition is set to be 1).