Figures
Abstract
Chemotactic bacteria form emergent spatial patterns of variable cell density within cultures that are initially spatially uniform. These patterns are the result of chemical gradients that are created from the directed movement and metabolic activity of billions of cells. A recent study on pattern formation in wild bacterial isolates has revealed unique collective behaviors of the bacteria Enterobacter cloacae. As in other bacterial species, Enterobacter cloacae form macroscopic aggregates. Once formed, these bacterial clusters can migrate several millimeters, sometimes resulting in the merging of two or more clusters. To better understand these phenomena, we examine the formation and dynamics of thousands of bacterial clusters that form within a 22 cm square culture dish filled with soft agar over two days. At the macroscale, the aggregates display spatial order at short length scales, and the migration of cell clusters is superdiffusive, with a merging acceleration that is correlated with aggregate size. At the microscale, aggregates are composed of immotile cells surrounded by low density regions of motile cells. The collective movement of the aggregates is the result of an asymmetric flux of bacteria at the boundary. An agent-based model is developed to examine how these phenomena are the result of both chemotactic movement and a change in motility at high cell density. These results identify and characterize a new mechanism for collective bacterial motility driven by a transient, density-dependent change in motility.
Author summary
Bacteria growing and swimming in soft agar often aggregate to form elaborate spatial patterns. Here we examine the patterns formed by the bacteria Enterobacter cloacae. An unusual behavior of these bacteria is the collective movement of cells after the initial aggregation into a tiny spot. Despite the majority of the cells within an aggregate being immotile at any point in the time, the flux of cells entering and leaving the aggregate, as motility is lost and regained in individual cells, led to a net, collective movement of the aggregate. These spots sometimes run into each other and combine. By looking at the cells within these spots under a microscope, we find that cells within each spot stop swimming. The process of switching back and forth between swimming and not swimming causes the movement and fusion of the spots. A numerical simulation shows that the migration and merging of these spots can be expected if the cells swim towards regions of space with high concentrations of attractant molecules and stop swimming in locations crowded with many cells. This work identifies a novel process through which populations of bacteria cooperate and control the movement of large groups of cells.
Citation: Courcoubetis G, Gangan MS, Lim S, Guo X, Haas S, Boedicker JQ (2022) Formation, collective motion, and merging of macroscopic bacterial aggregates. PLoS Comput Biol 18(1): e1009153. https://doi.org/10.1371/journal.pcbi.1009153
Editor: Kiran Raosaheb Patil, University of Cambridge, UNITED KINGDOM
Received: May 27, 2021; Accepted: December 7, 2021; Published: January 4, 2022
Copyright: © 2022 Courcoubetis et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript and its Supporting Information files.
Funding: JB acknowledges support from NSF Directorate for Mathematical and Physical Sciences CAREER Award PHY-1753268. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
In populations of chemotactic bacteria, coupling of the directed movement of individual cells in response to nutrients or chemical stimuli gives rise to spatio-temporal collective phenomena, including swarm bands and aggregates [1–3]. These macroscopic structures are the result of an emergent pattern of bacterial cell density that forms due to the coordinated movement and metabolic activity of billions of bacterial cells in an initially uniform environment. Both swarm band and aggregate formation rely on chemotaxis [1–4]. These collective phenomena can even be predicted with analytical mathematical considerations and recreated with detailed computational models [3,5–10]. Decades of work has resulted in a detailed and predictive understanding of bacterial collective phenomena, based mostly on work with bacterial species Escherichia coli, Salmonella typhimurium and Myxococcus xanthus [1–3,11]. It is unclear whether the collective properties within these three species encompass the full range of collective behaviors observed in all chemotactic bacteria, or if our understanding of these behaviors extends to other species of chemotactic bacteria.
A variety of cellular motility rules have been attributed to the emergence of larger macroscopic properties via the collective organization of high local densities of cells. A main driver of bacterial pattern formation is chemotaxis; individual cells utilize flagella to move in a combination of runs in a straight line, interrupted by tumbles to randomly change the direction of swimming [12,13]. At the molecular level, the switch between the run state and the tumbling state allows the bacteria to navigate towards increasing chemoattractant gradients [14]. For instance, many swarm bands of Enterobacteriaceae species are the result of cells migrating up concentration gradients of nutrients following local depletion [15]. The collective responses can also depend on the cell density and shape. For example, the surface swarming of Bacillus subtilis exhibits different morphologies, depending on the aspect ratios and surface densities of the cells [16]. Moreover, single cell properties such as adhesion, also direct the nature of collective phenomena. During the initial stages, the fruiting bodies of Myxococcus xanthus spontaneously assemble through the adhesion of cells, when two collide with each other [17]. This diversity of mechanisms governing the movement of individual cells has given rise to multiple macroscopic dynamics.
Bacterial cells are internally driven motile agents, and thus belong to the category of active matter. Active matter is a branch of non-equilibrium physics that considers microscopic rules and emergent macroscopic phenomena of energy consuming motile agents [18]. Driven inorganic matter also lies in the realm of active matter and has been used to probe the effect of individual properties of agents on the collective. For instance, self-propelled colloidal particles form aggregates, with size that linearly increases with particle speed [19]. Often living and non-living systems obey similar individual rules, for example swarming and swimming bacteria at high densities and shaken granular materials belong to the same active matter category of self-propelled rods [20]. Local motility interactions, which induce a distance dependent velocity alignment of moving agents, as dictated by the Vicsek model, give rise to emergent collective motion in multicellular organisms [21]. This density driven motility transition has been observed in the schooling of fish, flocking of birds, cells and insects [22]. The concepts of active matter systems can help to understand complex biological and physical processes, and even help to develop micromachines and nanomachines for practical applications [23,24].
Here we report new collective properties observed in the bacterial species Enterobacter cloacae, i.e., the formation, long-distance movement, and merging of macroscale bacterial aggregates. Chemotactic pattern formation was reported in this bacterium as part of a recent study of pattern formation in wild bacterial isolates [6]. Here we quantify a pattern of spots that emerges within an initially well-mixed culture of cells in soft agar and track the movement of each spot over multiple hours. In contrast to bacterial aggregates of Escherichia coli, which have been reported as mostly stationary with slight “jiggling” over time [1], aggregates formed by the bacterium Enterobacter cloacae migrate over distances up to four times their diameter. In addition, this movement results in approximately 36% of the aggregates merging with another aggregate to form a single spot. Our aim is to quantify and explain the spatial characteristics of these aggregates, their motility, and the underlying kinematics of the merging phenomena. High magnification, time-lapse imaging of individual cells within the spots reveals the microscopic mechanism that enables the collective motion of spots, namely a transition of individual cells between the motile and non-motile state. Chemotactic agent-based simulations, which include a novel immotile to motile transition, recreate the spatial order observed in experiments. This transition gives rise to a novel type of collective motility in bacteria.
Results
Aggregate formation by Enterobacter cloacae on soft M9 + glucose agar
Enterobacter cloacae is a bacterial species previously reported to be capable of large-scale pattern formation, including moving bands of high cell density and aggregate formation [6,25]. Here we focus on the aggregate formation, using a large culture dish and uniformly mixing cells into the culture medium at the beginning of the experiment. Cells grown overnight in Luria Bertani media at 37°C, and 180 rpm were uniformly mixed into soft minimal agar media supplemented with glucose at 0.07% inoculum. The media containing cells were poured into 22 cm x 22 cm dishes, with a lid and incubated at room temperature, as depicted in Fig 1A. Pictures of the dish were taken over time using a DSLR camera at 1X magnification. Initially, no macroscopic aggregates could be observed. After 5 hours, aggregates began to form, and at around 20 h a pattern of spots emerged across the plate, as shown in Fig 1B and 1C. Images of aggregates were analyzed using image analysis software as discussed in Materials and Methods.
(A) Sketch of experimental setup. (B) Image of a 100x100 mm2 subregion of the experimental plate before and after aggregate formation. (C) Zoomed in 20x20 mm2 region showing the aggregates at 22.5 h. (D) Number of aggregates on the plate versus time for the duration of the 44 hour experiment. The number of spot aggregates peaks at 22.5 hours, marked by the vertical line.
Aggregates form over about 10 hours, reaching a maximal number of aggregates at 22.5 h, as shown in Fig 1D. A video of spot formation over the entire plate can be found in S1 Video. As shown in the video, many of the aggregates migrate on the plate after formation, and some even merge together. After 36 h, the aggregates begin to dissolve and are no longer visible on the plate. In order to confirm the reproducibility of our results, we used identical culture conditions in two additional, independent experimental set- ups. We observed the emergence of aggregates, the movement of aggregates, and the occurrence of merging events between aggregates (S2 Video). After confirming the observed phenomenon, we proceeded to analyze the formation and the spatial patterns of these aggregates, the migration of aggregates on the plate, as well as the observed aggregate merging process.
Aggregate spatial structure analysis
As shown in Fig 1, thousands of aggregates appear on the plate. We analyzed the distributions of sizes and nearest neighbor distances of the aggregates at the time of the maximum number of aggregates. Replicate experiments showed similar spatial structure (S1 Fig). As shown in S2 Fig, aggregate sizes and the overall spatial distribution did not vary significantly between the time when spots started filling the plate and when the spots began to dissolve. As observed in S3 Fig, there was a large-scale variation in the local density of aggregates, the densest regions containing 25 aggregates per cm2 and the sparsest regions having 5 aggregates per cm2. Despite this variation in the density of aggregates, the trends observed in analysis of the center of the plate, shown in Fig 2, are similar to trends in the larger plate (S2 Fig). The subset of the plate, highlighted in S3 Fig, was used to obtain a higher resolution of the local spatial structure and to avoid edge effects. As shown in Fig 2A, the aggregates were fairly uniform in size, with an average area near 1 mm2. A small fraction of aggregates is larger and does not belong to the same peak in the histogram. Nearest neighbor distances were calculated using Voronoi tessellation, considering all nearest neighbors [26]. As shown in Fig 2B, the distribution of nearest neighbors is a slightly skewed Gaussian, with a characteristic nearest neighbor distance around 2.7 mm. To compare the aggregate point pattern with a random point pattern, 100 random point patterns were generated in a 100x100 region with 1340 points, equal to the number of aggregates found in the 100 x 100mm2 subregion analyzed in Fig 2. Using the Kolmogorov-Smirnov two-sided test, the resulting nearest neighbor probability distribution of the random point patterns (Fig 2B) was found to be non-identical with that of the aggregate pattern with a probability distance of 0.150 and a p-value of 3.69e-155.
A-D) Spatial analysis of the aggregate pattern at 22.5 h for a subregion of the plate. (A) The aggregate area distribution is depicted, with a mean of 1.01 mm2 and standard deviation of 0.25 mm2. (B) Histogram of the nearest neighbor distance distributions of the aggregate pattern and a random point pattern for comparison. For the bacterial aggregates, a pronounced single peak shows the existence of small-scale structure in the pattern. The average is 3.33 mm, with standard deviation of 2.33 mm. The distribution of the random pattern had an average of 3.35 mm, with standard deviation of 3.00 mm and it was found to be statistically different from that of the aggregate pattern, with a two-sided Kolmogorov-Smirnov test distance 0.15 and a p-value of p = 3.69e-155. (C) Structure factor versus wavenumber. The structure factor saturates to 1 at 3 mm-1, suggesting absence of long-range order. The dotted horizontal line represents the structure factor of a random point pattern for large system size. (D) Radial pair correlation function. The pair correlation is zero for distances less than 1.7 average aggregate diameters, as no other spots are detected in that proximity. The pair correlation peaks at 3.37 mm, exhibiting short range order, and saturates to one at 7 mm, consistent with disorder at longer distances. The dotted horizontal line represents the pair correlation factor of a random point pattern for large system size.
To further analyze the overall spatial pattern of aggregates, the pair correlation function and the structure factor of the subregion of the plate were calculated (see Materials and Methods for details). As observed in Fig 2C, the structure factor shows a single peak at 2.5 mm-1, and lacks additional peaks, indicating an absence of long-range order. This type of structure factor resembles that of a liquid, characterized by a single peak at low wavenumber, k, that repeats with a diminishing amplitude for multiples of nearest neighbor distances. However, it is unclear whether the fluctuations for large k are due to noise and finite sample size or long-range order. Fig 2D shows the results of the radial correlation analysis. Here, the ultra-short aggregate-to-aggregate distances appear to be absent, indicating an exclusion zone approximately two times larger than the average aggregate diameter. The short-range structure identified in the pattern matches that of a liquid, which is described by an exclusion region mediated by short range repulsion. In the long-range limit, the pattern of aggregates resembles a gas, as there is no long-range order. In summary, the aggregate pattern can be described as having no Bragg peaks, i.e., no long-range order (Fig 2C), but instead short-range correlations and a restriction on the minimal spacing between aggregates.
Hard sphere models have been widely used to model liquids and successfully capture their quasi-universal spatial structure [27]. In the hard sphere model, each particle is defined as a sphere with a fixed radius, which cannot overlap with the other spheres in the system. Interestingly, the spatial structure of the bacterial aggregates is consistent with a system of closed packed hard spheres whose radii are distributed according to a Gaussian [28]. Specifically, the coefficient of variation for spot size distribution was determined to be η = 0.3, and using the same coefficient of variation for a Gaussian, size distributed hard sphere packing, one retrieves an identical functional form for the radial pair correlation function [28]. This is consistent with the experimental aggregates size distribution. Calculating the pair correlation from a hard sphere model with constant radius of exclusion is enough to recover the qualitative characteristics of the pair correlation function in Figs 2D and S4. However, contrary to the hard-sphere model, the bacterial aggregates are not closely packed, and nearest neighbors are on average separated by approximately three times the average aggregate diameter. The explanation for this difference is that the aggregates form by recruiting bacteria in proximity, thus extending their exclusion region beyond the physical aggregate size.
Aggregate motility analysis
After aggregate formation, a displacement of some of the spot aggregates over time was observed, as shown in Fig 3. In Fig 3A, the green circle shows the location of each aggregate at 14.5 h, the magenta circle shows the location of each aggregate at 37.3 h, and the yellow line indicates the aggregate trajectory. Over twenty hours, the aggregates were displaced by up to a few millimeters, without a significant change in aggregate area or shape. Note that the upper left aggregate dissolves prior to 37.3 h.
(A) Snapshots of a subregion of the plate at T = 14.5 h and T = 37.3 h. The aggregates are depicted along with the trajectory identified by the tracking algorithm. (B) Distribution of the average speed of all aggregates (n = 495). (C) Distribution of path deflection, quantified as the ratio of distance to displacement (n = 495). The path taken by the aggregates slightly deviates from a straight line. D) Log-log plot of mean square displacement versus time. The fit yields a coefficient of 1.54 indicating that the trajectories are, on average, superdiffusive. Brownian motion, corresponding to diffusion, yields a coefficient of 1.0.
Time-lapse movies of the dish were used to quantify the trajectories of 495 aggregates, obtained by tracking the aggregates in six 26 mm x 26 mm subregions, shown in S5 Fig. The complete trajectory and positional information output for each subregion can be found in S1 Appendix. An analysis of the aggregate trajectories reveals that a subset of 35.4% of the aggregates participated in a coalescing or merging event. Specifically, 31.1% aggregates merged with another aggregate, and 4.2% of aggregates merged with two other aggregates resulting in 17.0% reduction in the number of aggregates. In the next sections we discuss the merging of two or more aggregates, while here we focus on the motility of all aggregates. Fig 3B shows the distribution of speeds for all the aggregates. The average speed of aggregate migration is 0.050 mm/hr. The path taken by the aggregates is quantified by dividing the distance over displacement in Fig 3C. Here, a value of 1, the minimal value, indicates a straight-line trajectory. The ratio of distance to displacement shown here, with an average value of 2.0, indicates that aggregates move in a directed manner. Plotting the average mean squared displacement versus time yields that the aggregate trajectories lie in the super-diffusive regime with a power law coefficient of 1.54 (Fig 3D). Thus, the aggregate trajectories are directed. Isolating the non-mergers, a weak but statistically significant negative correlation is found between the average minimum nearest neighbor distance and average speed for each aggregate trajectory, with a Spearman rank-order correlation coefficient correlation -0.12, p-value 0.034 (S6 Fig). The negative correlation indicates that proximity with other aggregates increases the speed of a given aggregate.
Microscopic analysis of aggregates
To shed more light on the microscopic behavior of the cells within as well as around the aggregates, we scaled down the experimental system using the Lab-Tek chamber, so that the aggregate movement could be captured by microscope at regular intervals.
For this experiment, cultures were prepared by tagging 5% of Enterobacter inoculum with RFP [29], to monitor single cell movement. Fluorescent portion of the populations enabled us to track movement of the boundary of an individual aggregate, as high cell density regions of aggregates have higher fluorescence than the rest of the focal plane (Fig 4A and 4B and S3 Video). When measured motility for Enterobacter cells expressing RFP was compared with Enterobacter host cells, we found no significant difference (S8A and S8B Fig), ensuring the reproducibility of the experiment in Lab-Tek chambers with mixed populations.
The Lab-Tek chamber loaded with a growth medium containing Enterobacter cells was incubated at room temperature for 24 hours and then mounted on a microscope after the visible appearance of aggregate (S7 Fig). As observed in liquid cultures in M9 + 0.4% glucose, comprising similar percent of cells in inoculum as that in plate and Lab-Tek chamber experiments, both Ecc1 and Ecc1 + RFP populations not only have similar growth rates (0.48 ± 0.13 h-1 and 0.44 ± 0.04 h-1, respectively) during exponential phase but also reach stationary phase approximately 10 hours after inoculation (S8C and S8D Fig). Reduced cell divisions in the stationary phase implies that the growth likely does not influence the observed pattern formation during the course of microscopy.
High magnification imaging of cells in an aggregate. Aggregates were formed within a coverslip bottom chamber to enable 100X imaging of cells within the aggregate. (A) Motile cells exist in the vicinity of the aggregate and a subset consolidates with the aggregate and becomes immotile. Representative trajectories of individual cells that are motile at t = 0. Yellow lines denote cells that swim towards the spot and become immotile within 0.24 s. Green lines denote cells that remain motile for 0.24 s. (B) Immotile cells within the aggregate. Shown in blue are representative cells within the aggregate that do not change position over 10 seconds. In A and B, trajectories are labeled with respective cellular velocities in μm/sec, and the scale bar is 40 μm. (C) Aggregate movement on the microscale. Progression of the boundary of a typical aggregate over 80 mins. The yellow solid line marks the periphery of the spot in real time, whereas the cyan solid line indicates the position of the spot front at t = 0 mins. Direction of the spot movement is shown with the white arrow. (D) Distribution of aggregate front speeds. By tracking the position of an aggregate boundary over time, the front speed is calculated for n = 40 aggregates. The histogram shows the measured distribution of front speeds with an average value of 0.023 mm/hr.
We pooled data from imaging of 40 individual aggregates recorded over 80 mins in 8 different microscopy sessions. ImageJ was used to analyze aggregate and single cells in aggregate through time. Displacements calculated for single cells reveal the existence of three sub- populations in the milieu (Fig 4A and 4B). As shown in Fig 4A, a typical aggregate was observed to be surrounded by motile cells. A subset of such cells is identified with green circles. Interestingly, some of the motile cells coalesce with aggregate and lose their motility, thus, switching into non- motile cells (yellow circles). Fig 4B, on the other hand, shows cells within the boundary of the aggregate, which were observed to be non- motile for the duration of the video (cyan blue circles).
We then analyzed the progression of the aggregate boundary through time (Fig 4C). We marked and tracked the edge of the aggregate using ImageJ. The difference between the XY coordinates of the centroid of the boundary in the initial frame (solid cyan blue line) and the final frame (solid yellow line) was calculated to obtain the speed of aggregate boundary. An analysis of 40 aggregates revealed that the average speed exhibited by the aggregate boundary is 0.023 mm/hr (Fig 4D), which is comparable to the aggregate speed calculated for a large plate experiment.
Since the microscope was focused on a segment of the aggregate boundary chosen randomly, there is no guarantee that the trailing or receding edge of the aggregate was recorded. Therefore, measuring the boundary propagation speed corresponds to obtaining a component and not the magnitude of the aggregate speed.
Nevertheless, these microscopic observations suggest spot movement must be driven by the flux of single bacterial cells leaving and joining the aggregate. In the direction of movement, individual cells join the aggregate and lose the motility (S3 Video), whereas at the trailing edge motility is regained and cells leave the aggregate (S4 Video).
Aggregate merging analysis
Next, we analyze the merging of two or more aggregates. The subregions studied for quantifying spot motility were also used for analyzing merging events (S5 Fig). As shown in Fig 5A, the movement of aggregates sometimes results in the combination and merging of multiple aggregates into a single aggregate. Fig 5A shows merging of two sets of aggregates. Upon analysis of the trajectories for 495 aggregates, 64.6% of aggregates did not merge, 31.1% aggregates merged with another aggregate, and 4.2% of aggregates merged with two other aggregates. Examples of mergers involving more than two aggregates are shown in S9 Fig. The remainder of the analysis is done for merging events with one other aggregate, which we define as two-spot mergers. Within this subset, only 46 out of the 77 trajectories were included in the analysis, i.e., leaving out merging trajectories with less than 20 positional data points whose dynamics could not be determined with a similarly high degree of accuracy.
(A) Snapshots of a 10 mm x 9mm subregion of the plates at T = 17.5 h and T = 31.5 h. Two instances of two-aggregate merging events are highlighted in the boxed regions. The trajectories are plotted in blue for non-merging aggregates and in red for mergers of two aggregates. (B) Relative distance vs. time for all merging processes, color coded with total aggregate area. (C) The acceleration is positively correlated with the total aggregate area (spearman’s rho 0.30 and p-value 0.040). (D) Distribution of aggregate speed for two-spot mergers and non-merging spots. The average (shown with a dashed line) illustrates that aggregates that eventually merge, on average, move faster. (E) Distribution of trajectory distance divided by displacement for two-spot mergers and non-merging spots. On average, the path taken by aggregates that eventually merge is more direct.
The trajectories taken by merging aggregates are shown in Fig 5B. The relative distance of the two merging aggregates is defined as , where and are the two-dimensional position vectors of aggregate A and B at time t. The trajectory taken by most merging aggregates was similar, exhibiting an increase in aggregate velocity as the distance between aggregates decreased. With this observation in mind, we hypothesized the existence of a distance dependent force law that governs the attraction between aggregates. However, both exponential and power law fits to the acceleration versus relative distance of the trajectories resulted in a wide distribution of exponents and coefficients (S10 Fig). Thus, to quantify the dynamics of the merging spots, we took a simple approach, and the trajectories were treated with a constant acceleration model.
To quantify the merging dynamics, the relative distance versus time for each trajectory was fit to a quadratic function. The quadratic fit was excellent for the majority of trajectories, with an average error of 0.15 pixels. The fit revealed the acceleration of each aggregate, and as shown in Fig 5C, the acceleration was found to be correlated with aggregate size (Spearman rank-order correlation coefficient correlation 0.3, p-value 0.040). Quantifying the relative position of the aggregates relative to point of collision gave similar results, see S11 Fig. Curiously, unlike inertia dominated systems, the smaller spots do not necessarily move more than larger spots (S12 Fig). As shown in Fig 5D and 5E, the aggregates that merged moved fast on average and took a more direct path. The non-merging aggregates had a mean speed of 0.042 mm/hr, whereas the mergers had a mean speed of 0.050 mm/hr, with a statistically significant difference of 0.007 mm/hr (p = 0.0001 for an unpaired t-test). The fact that merging aggregates are on average faster and follow a more direct path suggests the presence of an effective attractive force between aggregates.
Simulation
Based on the experimental observations, a 2D agent-based model was developed to explain the phenomena observed in aggregates of chemotactic bacteria. Our simulations were based on a model proposed in (9), with algorithms discussed in [30–32]. Briefly, 100,000 cells move in space, with movement biased by a chemoattractant molecule produced by the cells. There are two responses to the chemoattractant: 1) at low concentration of attractant the cells move towards regions of higher chemoattractant and 2) at high concentrations of the attractant the cells transition to an immotile state. Cells are revived from the immotile state after a period of time and regain motility. Cells do not grow or divide within the simulation, as the patterns observed in experiments are formed after the cultures enter the stationary phase.
In more detail, 100,000 agents (cells) perform random Brownian motion while experiencing a chemotactic drift force. The chemotactic drift force on a given agent is equal to the gradient of the chemoattractant evaluated at the agent’s current position multiplied by a coefficient of chemotactic sensitivity. Regarding chemotactic sensitivity, receptor law sensitivity was used [9], i.e., the chemotactic force magnitude increases with the slope of the chemoattractant gradient but decreases with the concentration, a formulation motivated by the saturation of bacteria surface receptors. Chemoattractant is produced and expelled by each agent, and the chemoattractants undergo diffusion and decay over time. Within the model, chemoattractant undergoes degradation within the environment. Although not tested here, chemoattracts could be taken up and removed from the environment by cells, which could alter the chemoattractant gradients at short length scales and potentially alter larger scale patterns of cell density. The model does not incorporate cell division, death or food gradients. Further rules were implemented to make the simulation more realistic and capture experimental results: agents were attributed a finite size and they experience a chemoattractant dependent motility transition (See S2 Appendix for detailed model description).
Based on our observation that cells within the aggregates are immotile, the model was generalized to include a mechanism whereby an agent can transition from being motile to immotile and vice versa. To that end, a chemoattractant threshold was introduced, beyond which bacteria transition to an immotile state. Since the chemoattractant concentration is proportional to the local density inside the aggregates, using a chemoattractant threshold corresponds to a density-dependent motility transition. As shown in S13 Fig, introduction of the density-dependent motility transition resulted in a greater fraction of motile cells within the population. Including only the density-dependent transition to the immotile state would lead to the eventual absorption of the vast majority of cells in aggregates. In experiments, however, the aggregates coexist with freely swimming cells throughout the course of the experiment. Furthermore, the analysis of aggregate movement in Fig 4 indicates examples of the departure of cells from the aggregate interface. Therefore, a rule for reactivating motility was needed, even within high cell density regions. There are examples in the literature of a timed motility switch [33], so a rule was incorporated for motility to be regained after a random interval of time.
More specifically, an agent stays in the immotile state for a fixed interval of time, after which there is a constant probability of regaining motility per iteration. To enable newly mobilized cells to leave regions of high cell density after regaining motility, cells ignore the chemoattractant gradient for a short period, performing Brownian motion, after switching from the immotile to motile state. The expected value for the time scale for an agent to regain motility was set to 33 mins (3 minutes of a fixed immotile state and an expected value of 30 minutes for stochastic regain of motility) and the time scale for Brownian motion after regaining motility was set to 15 mins. Lacking experimental data for the details of the motility transition, we had to choose the aforementioned parameters. The duration of Brownian motion was chosen such that an agent can transverse a distance of the order of a millimeter after regaining motility. Finally, the timescale for regaining motility was set to approximately double the duration of Brownian motion such that motile agents make up a significant fraction of the agent population. For a given aggregate, assuming that agents get immediately reabsorbed after regaining motility and performing Brownian motion, this choice approximately sets one third of agents to be in the motile state and two thirds in the immotile state.
The spatial scale of the simulation was set by equalizing the average aggregate radius of the simulation with the experiment. The time scale was set by equalizing the average speed of an agent with that of Enterobacter cloacae, as measured in [6]. To perform this scale calibration, the average speed of every motile agent was calculated from the simulation and equated to the experimental value. The parameters were chosen by numerical testing and are shown along with a laconic description in S1 Table. The goal of the numerical exploration was to retrieve formation of aggregates while qualitatively capturing the relative aggregate size and spacing seen in the experiment. Exploration of systematically varying different parameters and components of the model are shown in S14, S15 and S16 Figs. The model has 13 parameters and is non-linear and stochastic, so it is simply not feasible computationally to systematically test the entire parameter space and analyze the emergent properties. The goal of these simulations was not to extract and use precise numerical values for unknown parameters, but instead to explore whether an established model of chemotactic behavior combined with the experimentally observed motility transition would be able to qualitatively reproduce the experimentally observed phenomena.
Our simulations, using these rules for motility switching, reproduce the formation of bacteria aggregates observed in experiments. Over time, high cell density aggregates composed of mainly immotile cells form. After the formation of aggregates and throughout the simulation, motile and immotile cell populations coexist with approximately 40% of cells retaining their motility at any given time (S17 Fig). As a consequence of setting the spatial scale using the experiment, these aggregates have a diameter of around 1 mm. The size distribution resembles experimental results but exhibits higher variance. As shown in Fig 6B and 6C, the spatial order of aggregates observed in the simulations is similar to the experimental measurements. The pattern retains short range order, described by an exclusion region, resembling a liquid, and disorder at long spatial distances, resembling a gas.
A)- D) Spatial Structure of Aggregates. (A) The simulated spot area distribution has an average of 0.97mm2and a standard deviation of 0.38 mm2. (B) The simulated nearest neighbor distribution has an average of 3.30mm and a standard deviation of 2.49 mm. (C) The spatial correlation function of the simulated aggregates mirrors the saturation of the experimental structure, with saturation to 1 at 3.50 mm. (D) Snapshots before (T = 5.0 h) and after (T = 46.5 h) the merging of two segmented aggregates (see Materials and Methods for segmentation details). E)-H) Merging dynamics in the simulations. (E) Relative distance vs. time for all merging processes, color coded with total aggregate area. (F) The acceleration is not found to be correlated with aggregate size (spearman’s rho 0.01 and a p-value of 0.962). (G) Distribution of aggregate speed for two-spot mergers and non-merging spots. The average (shown by a dashed line) illustrates that mergers, on average, move faster. The speeds of the simulated aggregates are on the same order of magnitude as found in the experiments. (H) Distribution of the ratio of trajectory distance divided by displacement for two-spot mergers and non-merging spots. On average, the paths taken by spots that merge are more direct. Non-merging spots follow a less direct path as compared to the experiments.
In addition to forming a large-scale pattern of aggregates, a collective motility of cells within the aggregate and aggregate merging is also observed. In Fig 6E, 6F, 6G and 6H, we analyze the merging trajectories in the simulation. In this case, the motion is also well described by a constant acceleration model. However, in contrast to the experiment, the aggregate size dependence of the acceleration is not found to be statistically significant, with a Spearman rank-order correlation coefficient of 0.01 and a p-value of 0.962. Nonetheless, the average speed and trajectory distance divided by displacement probability densities, for both merging and non-merging aggregates (Fig 6G and 6H) resemble the experiment. Notably, the trajectories in the simulation for non-merging aggregates are less directed. Finally, the merging frequency in the simulation is lower than in the experiment, with only 6% of spots merging in the simulation, compared with 35% in the experiment.
The model is not able to quantitatively capture all aspects of the experiment. The motility profile of non-merging aggregates reported in the experiments (Fig 5D and 5E), and the respective motility profile in the simulations (Fig 6G and 6H) are quantitatively different. Using distance over displacement as a measure, merging aggregates in the simulation take on average a 47% more direct path while in the experiment they take a 28% more direct path. Furthermore, merging aggregates in the simulation are on average 77% faster than non-mergers while in the experiment the merging aggregates were only 17% faster. This discrepancy illustrates a limitation of the computational model to capture the motility profile of non-merging aggregates. Notably, the average speed of merging aggregates in the simulation, 0.047 mm/hr, was only 6% from the respective value in the experiment (0.050 mm/hr). Finally, the distribution of areas in the simulation is quantitatively different from the experiment. Specifically, the standard deviation of the aggregate areas in the simulation is larger than the experiment by 41%. This difference might be explained by the presence of fewer cells in simulated aggregates as compared to the aggregates observed in the experiments. In both the simulation and experiment, agents can regain motility and exit the aggregates. However, due to the small number of agents per aggregate in the simulation, a small number of agents leaving an aggregate can lead to significant fluctuations in aggregate size. When the ability to transition to motile state or the random walk after regaining motility is removed (S14 Fig), the resulting area distributions become narrow, with respective standard deviations 23% and 31% lower than their experimental counterpart.
Aspects involving the motility switch were investigated systematically (S14, S15 and S16 Figs). Variation or withdrawal of aspects of the motility mechanism do not change the spatial characteristics of the pattern (S14 and S15 Figs) This highlights the fact that the spatial order of the pattern is a direct consequence of chemotaxis. An analysis of the effect of motility on merging was also conducted (S16 Fig). Removing the motility switch altogether results in approximately twice the merging events, but merging and non-merging aggregates does not yield distinct average speed distributions. Thus, the motility switch improves the agreement between simulation and experiment. In addition, the decrease of merging events due to the motility transition prolongs the lifetime of the aggregate pattern. Finally, removing the ability of agents to random walk after regaining motility leads to approximately three times the merging events, but the non-merging aggregates are effectively immotile. This is expected, as cells are less likely to escape an aggregate after regaining motility and thus leading to less motile aggregates.
In the simulation, merging events are caused by an imbalanced flux of agents into and out of the aggregate. An example of a typical merging event is shown in Fig 7A, 7B and 7C. To keep track of the cells in each aggregate in this example, Fig 7D and 7E show the changes in the number of motile and immotile cells over time. When a sufficient number of cells exits the spot and its vicinity, the core of immotile cells collectively dissolves, as the chemoattractant drops below the motility threshold (Fig 7D arrow). This aggregate consists of only motile cells, as shown in Fig 7B, and it moves towards and joins the other aggregate.
A-C) Three frames of the same region at different times in the simulation. The two aggregates, both consisting of motile (blue) and immotile cells (red), eventually merge. (D) The number of cells in a 2 mm x 2 mm square region encompassing the initial position of the bottom aggregate from A) vs. time. The aggregate gradually loses both motile and immotile cells until the chemoattractant concentration drops below the threshold at 11.8 hours (arrow). Then, all cells regain motility and collectively move towards the proximal aggregate. (E) The number of cells in a 2 mm x 2 mm square region positioned on the initial position of the upper aggregate from (A) vs. time. The aggregate experiences fluctuations in the number of cells and eventually receives the motile cells from the proximal aggregate.
In the simulations, prior to fully merging, one of the colonies becomes enriched in motile cells, which then join the other colony. This may account for the accelerated colony velocity during merging events observed in both simulations and experiments, but the experimental data does not directly show this process. Furthermore, it is consistent with the observation that merging events in the experiment lead to the formation of a new radially symmetric aggregate. The merging process depicted in Fig 7 is ubiquitous—more examples of mergers, with visualization of the motility of the agents, are shown in S5 Video. This emergent phenomenon of collective regain of motility is not responsible for the non-merging aggregate speeds reported in Fig 6G. Non-merging aggregates in the simulations consist of immotile cells that move due to the exchange of agents.
Discussion
Many examples of pattern formation and collective motility within bacterial populations have been reported, and here we examine such phenomena within populations of Enterobacter cloacae. This system displays two unique characteristics that have not been reported previously in the context of bacterial emergent behavior: the formation of bacterial aggregates via a transient motility transition and the collective motility of aggregates of swimming bacteria. Such behaviors are reminiscent of other bacterial systems, notably swarming and fruiting body formation in Myxococcus xanthus. Fruiting bodies are complex aggregates of cells that form as a result of bacterial swarming, and work over the years has shown that aggregates form as a result of contact-dependent and diffusive signaling, and that cells within the aggregate display a reduced velocity and changes in cell alignment [34–36]. Enterobacter cloacae accomplishes similar complex behaviors using different molecular and physical mechanisms. There have been several previous reports related to chemotaxis in E. cloacae, including work on phosphate chemotaxis [37,38], although these pathways are not as well characterized as those found in E. coli. Prior work has shown multiple isolates closely related to the strain used in this study form swarm rings and other chemotactic patterns [6]. Genomic sequences have revealed E. cloacae contain genes related to chemotactic response and flagellar synthesis [39], and key chemotactic proteins in Enterobacter cloacae had 87–95% homology with those in E. coli [40,41]. Enterobacter genes were also able to restore swarm ring formation in E. coli mutants [42]. We observed that Enterobacter aggregate motility on soft agar surface is a result of flagella- driven swimming and not because of swarming or gliding of bacterial cells on semi- solid surface [43]. Secondly, aggregate formation and movement are associated with a density-dependent motility transition in which cells temporarily lose motility upon entering regions of high cell density. Because of these important differences, Enterobacter cloacae could serve as an important new model system for studying bacterial collective behavior.
The observed motility transition that occurs in regions of high cell density is potentially a unique mechanism of bacterial aggregate formation. A transition between motile and immotile cells is more typically associated with cells of multicellular organisms, such as the epithelial-mesenchymal transition [44]. The active matter community has examined the cell motility and pattern formation [45–47]. Cells used here demonstrated a complete loss of motility, as opposed to reduced motility in regions of high cell density, suggesting the transition is not due to physical processes such as crowding and jamming. Such a complete loss of motility has been reported in Escherichia coli due to depletion of oxygen [48]. The loss of motility led to the formation of a ring of bacteria. Here a similar motility transition was associated with aggregate formation, and in fact the transition was shown to be transient. As shown in S4 Video, cells that were initially immotile within an aggregate boundary regained motility as the aggregate boundary receded. Transient changes in motility have been reported for Bacillus subtilis, in which molecular “clock” switches immotile cells in a chaining phenotype back to a motile state after several hours [49]. In that system, the timing of the transition was regulated by protein dilution via cell division. In the Enterobacter cloacae, it is unclear if cell division is occurring within dense aggregates, especially given the likelihood of low nutrient conditions within aggregates, so potentially another molecular mechanism would be needed to regulate such a transition.
The results herein constitute a first step towards reporting, quantifying and understanding the novel phenomena of bacterial aggregate motility and merging; future work is warranted for deeper insight into mechanism. The quantitative analysis of macroscopic components such as aggregate motility and merging serve primarily as a quantification of the reported phenomena. For instance, it remains unclear how aggregate motility and merging dynamics can be explained by a set of rules or analytical equations. Perhaps future studies can uncover underlying mechanisms by analysis of the existing data or by perturbing the system in new ways that more clearly reveal the underlying principles of aggregate dynamics. The simple approach of performing a quadratic fit of merging trajectories was able to infer that acceleration during merging correlates with aggregate size and can provide insight in future considerations and model building. In the context of future work, the role and effect of vertical direction of the system should be investigated. For instance, it is unclear if aggregate area is a good description of size as compared to aggregate volume. Furthermore, there could be oxygen gradients that form in the vertical direction that affect the system, as oxygen can play an important role in bacterial motility [48]. There are also multiple additional components that can control the aggregate migration dynamics: agar variations, oxygen, chemoattractant and nutrient gradients that can play a role in the details of migration and merging. Notably, these quantities are hard to measure and new experiments should be developed at this front. Finally, in the microscopic scale, future work is warranted to uncover the details behind the mechanism that controls the observed cellular motility transition that drives motility in the scale of aggregates, potentially through analysis of mutations that modulate aggregate formation and motility.
The potential benefits of aggregate formation, especially aggregates composed of immotile cells, remain unclear. Previous studies have proposed that the clustering of cells into high density aggregates and biofilms is a stress response that enables cells to survive under harsh conditions such as low availability of nutrients or exposure to toxins [47,50]. Would additional benefits be conferred to aggregates of immotile cells? Aggregate formation does not require a motility transition, although as suggested by the agent-based model developed here, the motility transition did produce a more stable pattern of aggregates. Immotile cells may be more energetically efficient. The length scale of the aggregates could also be set by adjustment of the molecular mechanisms that set the threshold density and timing of the motility switch, see S3 Appendix. The transient nature of the observed motility switch did enable migration and merging of the aggregates once formed. The scheme for migration suggested by the experiments is that unbalanced fluxes of cells on different sides of the aggregate resulted in collective movement of the aggregate. Whether such symmetry breaking was random or the result of asymmetry in the chemical and environmental conditions surrounding the cell remains unclear. Potentially, if the direction of the movement of individual aggregates over multiple aggregate lengths is biased by chemical conditions, this may serve to position cells in a more favorable location. The merging process after formation may also be a mechanism to optimize the pattern of aggregates for maximum benefit to the cells.
Materials and methods
Bacterial strains and culturing conditions
Enterobacter cloacae Ecc1 used for the experiment was isolated from the Caltech turtle pond and confirmed using 16s rRNA analysis [6]. Primary cultures of Enterobacter cloacae Ecc1 were grown overnight in Luria- Bertani broth (Difco) at 37°C under constant shaking at 180 rpm. Cultures were then washed thrice with 1X PBS (VWR, life science) and resuspended in 5 ml 1X PBS before inoculating in M9 minimal agar [6].
Aggregate formation in BioAssay plate
To observe aggregate formation, soft agar was prepared by mixing M9 medium (Difco) supplemented with 4% glucose (aMResco) with 0.26% of bacteriological agar (Sigma-aldrich). Primary culture of Enterobacter cloacae Ecc1 was added to sterilized minimal agar to final OD of 0.0002. 95.6 ml of the mix was poured into a 20 X 20 cm Bioassay dish (Thermo Fisher Scientific) with a lid. The BioAssay plate was allowed to set at room temperature for 1 hour before sealing it with parafilm.
The plate was incubated at room temperature for 2 days on a glass table. Reflection of the base of the Bioassay plate using a mirror was captured automatically by camera (Canon, EOS REBEL) every 15 minutes.
Aggregate formation in Lab-Tek chamber
Culture was made visible by mixing Enterobacter cloacae Ecc1 expressing RFP with Enterobacter cloacae Ecc1 to a final concentration of 5%. Culture was mixed with 10 ml soft minimal agar to final OD of 0.0002. 2 ml of mix was allowed to set in the Lab-Tek chamber (Nunc, Thermo Fisher Scientific) for 40 minutes, before sealing the chamber with parafilm. The entire assembly was incubated at room temperature for 24 hours, at that time visible aggregates had formed. The setup was mounted on a microscope for imaging.
Aggregates were imaged in phase contrast and RFP channels using 100X oil objective (1.5 NA CFI plan apochromat) on Nikon Eclipse Ti-E microscope with a sCMOS Camera (Zyla 5.5 sCMOS, Andor) at a 1 minute time interval with an exposure time of 30 and 100 ms respectively.
Measurement of growth rate
To measure the growth rate of E. cloacae Ecc1 transformed with pZE25-RFP and wild-type E. cloacae Ecc1, primary cultures of both the strains were inoculated in 1 ml M9 + (0.4%) glucose media to adjust final OD approximately at 0.0002 similar to cultures used in experiments with bioassay plate and Lab-Tek chamber. 200 μl culture from this stock were distributed per well in 96- well plate (Costar, Corning incorporated) to start three parallel experiments for each culture. OD at 600 nm was measured at room temperature with intermittent shaking using plate reader (TECAN, infinite M200PRO) at the interval of 30 mins for 20 hrs. Readings obtained were fit to logistic equation in MATLAB (R2020a, MathWorks) using cftool to obtain growth rate for each strain.
Estimation of cellular motility
Primary cultures of E. cloacae Ecc1 wild- type as well as transformed with pZE25-RFP were diluted ten times. 5 μl of diluted culture was spotted on a microscope slide, which was then covered with glass coverslip (22 mm X 22 mm). Movement of an individual cell in a drop was recorded using 40X objective on Nikon Eclipse Ti-E microscope with a sCMOS Camera (Zyla 5.5 sCMOS, Andor) with a time interval of 0.1 sec.
Image analysis
The 1X plate frames were processed with background subtraction and enhanced contrast such that only 0.01 percent of pixels are saturated using imageJ [51,52]. The processed frames were segmented using Weka, a machine learning algorithm for microscopic pixel segmentation [53], example segmentation shown from frame T = 22.5h in S18 Fig. Finally, dim segmented aggregates were filtered according to maximum pixel brightness, a process depicted in detail in S6 Video. Plotting the distribution of maximum pixel intensity, the threshold value of 0.45 was chosen to remove detected regions that were not high-density aggregates.
To perform the motility and merging analysis, tracking analysis was performed in imageJ via the Trackmate plugin [54]. The software is capable of producing tracks for all the segmented aggregates. Trajectories that did not lead to a merger were analyzed to determine the motility of non-merging spots. Trajectories that only contained merging events with strictly two aggregates were subject to the merging analysis. Although the aggregate size was also identified by Trackmate, the aggregate size as determined from the segmentation was used instead for greater precision and accuracy.
Images acquired at 100X were analyzed to evaluate the speed of aggregate movement using ImageJ 1.53a [51]. Initial and final time frames were Otsu thresholded to find aggregate boundary, which was then marked using a freehand line tool. XY-coordinates of the centroid of the resultant line were extracted and used to calculate the distance covered by the progressing aggregate front using the Euclidean distance formula.
To measure motility of an individual cell of E. cloacae Ecc1 and E. cloacae Ecc1 + RFP, time-lapse movies of bacterial movement were opened in ImageJ 1.53a. The centroid of a single cell was tracked manually over time and its XY-coordinates were used to calculate temporal displacement of the cell, which in turn was used to measure the velocity of the cell per second. Movements of 25 cells were analyzed manually over time for each strain.
Radial pair correlation and structure factor
The pair correlation function also known as the radial distribution function was calculated by considering the point patterns generated by the aggregate positions. All aggregates were considered in a given square region, except for aggregates that lie within a radius of rmax from the boundary. For each point, the distance between all other points was calculated. Then, the unnormalized probability density distribution was obtained for all the distances obtained. To obtain the probability density, the bins were set according to the desired resolution, defined from 0 up to and including rmax in increments of dr. Furthermore, the probability density was divided with the number density of the pattern, ρ = Ntotal/L2, where Ntotal is the total number of points and L is the length of the region. Finally, the value in the nth bin was divided with the area of the respective annulus that spans from the bin range, with an inner radius of rin = n * dr and an outer radius of rout = (n + 1) * dr. Finally, the values of the resulting probability density was averaged for all points to retrieve the pair correlation function. Furthermore, the structure factor was obtained by taking the fourier transform of the correlation. Thus, by integrating the pair correlation function we can obtain the structure factor: [55].
Simulation
The simulation code was written in Cuda, compiled in gcc version 4.9.4 and run on Cuda version 9.2.88, executed on a single GPU node. System size was set to a 384x384 pixel grid with periodic boundary conditions. At the beginning of the simulation, 100,000 motile agents were placed on random points in the grid. The chemoattractant concentration was initialized to a value of zero. The simulation final time was set to 2000 with a step size of Δt = 0.005. Chemoattractant concentration, agent position and velocity were outputted in text file format for further analysis and plotting. The main parameter set used can be found in S1 Table and is termed the control parameter set and the complete code can be found in S4 Appendix.
The output of the simulation was segmented and tracked. To segment the aggregates from the simulation results, a MATLAB script was written. For each pixel, the number of agents within a radius of 3 pixels was determined. When the number of agents exceeded 40, the pixel was considered to belong to an aggregate, an example segmentation is shown in S19 Fig. In addition, aggregates that consisted of less than 10 pixels were discarded. The parameters chosen to segment aggregates were empirically chosen by comparing segmentations with the respective agent density plot as in S19 Fig. To track the aggregates, the segmented images were entered into imageJ where the Trackmate [54] plugin was used to track the aggregates and output the trajectories.
Supporting information
S2 Video. Experimental replicates.
Video shows spot formation in three different set ups under identical growth conditions at the time interval of 15 mins. A few merging aggregates in three different plates have been indicated with cyan arrows. The middle video shows the plate that was used for all analysis shown in paper. Scale bar- 10 mm. Time stamp- hrs: mins.
https://doi.org/10.1371/journal.pcbi.1009153.s002
(AVI)
S3 Video. Aggregate advancing edge.
Microscopic images of aggregate movement captured in RFP channel with 1 min time interval. Solid yellow line indicates the progression of the aggregate front. Scale bar- 40 μm. Time stamp- hrs: mins.
https://doi.org/10.1371/journal.pcbi.1009153.s003
(AVI)
S4 Video. Aggregate receding edge.
Microscopic images of the aggregate movement recorded in phase contrast channel with frames 1 min apart. Trailing front of the spot has approximately been denoted by the solid yellow line. Scale bar- 40 μm. Time stamp- hrs: mins.
https://doi.org/10.1371/journal.pcbi.1009153.s004
(AVI)
S5 Video. Motility and merging.
In the agent-based simulations, the process that leads to merging in simulations is ubiquitously described by a collective regain of motility of the aggregates. In the video, multiple merging events take place. Cells marked in red are immotile and in blue are motile.
https://doi.org/10.1371/journal.pcbi.1009153.s005
(MOV)
S6 Video. Thresholding dimmer aggregates.
In this video of a zoomed in region of the plate, the processing step of thresholding dimmer aggregates is shown for the course of the experiment. The left panel highlights aggregates that were kept for further analysis with blue and aggregates that were not considered with red. The right panel shows the same image as the left, but without highlighting spot boundaries, for a clearer view of the aggregates.
https://doi.org/10.1371/journal.pcbi.1009153.s006
(AVI)
S1 Appendix. Experimental data for aggregate trajectories.
https://doi.org/10.1371/journal.pcbi.1009153.s007
(ZIP)
S2 Appendix. Mathematical description and simulation components.
https://doi.org/10.1371/journal.pcbi.1009153.s008
(PDF)
S3 Appendix. A model for steady state colony size.
https://doi.org/10.1371/journal.pcbi.1009153.s009
(PDF)
S2 Table. Data used to create main text figures.
https://doi.org/10.1371/journal.pcbi.1009153.s012
(XLSX)
S1 Fig. Distribution of areas, nearest neighbors and pair correlation functions of aggregate positions for replicate experiments.
Spatial quantities are calculated for three different plates at 22.5 hrs. The first two columns are unused replicate experiments, and the third column is obtained from the plate used throughout the study. The quantities are calculated for a 10x10cm subregion in the center of the plate to account for edge effects. A-C) Aggregate area distributions. D-F) Nearest neighbor distance distributions. G-H) Pair correlation functions.
https://doi.org/10.1371/journal.pcbi.1009153.s013
(TIF)
S2 Fig. Distribution of areas, nearest neighbors and pair correlation functions for aggregate positions at different times for both the whole experimental plate and the subregion of interest.
Spatial quantities for three time points 17.5 h, 22.5 h and 27.5 h from left to right. The quantities were plotted with blue for the whole plate and green for the subregion used in the main text, shown in S3 Fig. A-C) Aggregate area distributions. D-F) Nearest neighbor distance distributions. G-H) Pair correlation functions.
https://doi.org/10.1371/journal.pcbi.1009153.s014
(TIF)
S3 Fig. Spatial aggregate density fluctuations.
The aggregate number variations are visualized for the frame used for the spatial analysis, at T = 22.5 h. The highlighted region marks the spatial subset of the plate analyzed in Fig 2.
https://doi.org/10.1371/journal.pcbi.1009153.s015
(TIF)
S4 Fig. Pair correlation for hard sphere model point pattern.
(A) A hard sphere model pattern in a 100x100m2 region generated by assigning 650 points a random position with the condition that they do not lie within a distance of 3.33 mm from any neighboring points. The value 3.33mm is taken from the average nearest neighbor distance in the experimental subregion analyzed in Fig 2. (B) Pair correlation of the respective point patterns. The dotted horizontal line represents the pair correlation factor of a random point pattern for large system size. The pair correlation of the hard sphere model captures the spatial characteristics observed in the pair correlation of the bacterial aggregate pattern. Note that a deterministic radius of exclusion makes the peak at the average nearest neighbor distance more pronounced than in the experiment and leads to a smoother curve.
https://doi.org/10.1371/journal.pcbi.1009153.s016
(TIF)
S5 Fig. Regions of the plate used for analysis of motility and merging.
https://doi.org/10.1371/journal.pcbi.1009153.s017
(TIF)
S6 Fig. Trajectory characterization.
Average speed versus average minimum nearest neighbor distance for all analyzed aggregate trajectories. To obtain the latter quantity for a single trajectory, the nearest neighbor distances were calculated for each frame, the minimum was extracted and averaged over all frames. Spearman rank-order correlation coefficient correlation yielded -0.12, with a p-value of 0.034.
https://doi.org/10.1371/journal.pcbi.1009153.s018
(TIF)
S7 Fig. Depiction of experimental set- up using bioassay plate and Lab-Tek chamber.
The height of agar in both the vessels was 2 mm. Representative images taken by using DLSR camera and 4X objective on Nikon microscope for respective set- up has been shown. Spot formation in the Lab-Tek chamber was imaged in time-lapse microscopy using 100X oil/1.5 NA objective.
https://doi.org/10.1371/journal.pcbi.1009153.s019
(TIF)
S8 Fig. Control measurements for RFP producing populations.
Enterobacter cloacae Ecc1 cells transformed to express RFP protein were compared with Enterobacter cloacae Ecc1 cells for single cell motility (A and B). (C) Their respective growth in M9 + Glucose at room temperature was fit to the logistic growth equation (solid lines) to calculate their respective growth rates as shown in (D). n = 3. Error bars- Standard deviations.
https://doi.org/10.1371/journal.pcbi.1009153.s020
(TIF)
S9 Fig. Three aggregate merging instances.
Merging of three aggregates. 4.2% of aggregates took part in a merger involving three aggregates. The top and bottom row show two instances of a three-aggregate merger over time.
https://doi.org/10.1371/journal.pcbi.1009153.s021
(TIF)
S10 Fig. Merging dynamics and distance dependent force laws.
Probability density distributions of the coefficients obtained by fitting acceleration versus relative distance for each two-spot merging trajectory in Fig 3 with an exponential (A) and a power law (B). Exponents were included for fits that yielded exponents with one standard deviation error that was less than 50% of their respective value. (A) Distribution of the decay exponent, n, obtained by assuming an exponentially dependent force: . (B) Distribution of the power law exponent, n, obtained by assuming a power law dependent force: a(drel) = A0/dreln. It is not possible to distinguish whether a power or exponential law best describes the acceleration profile. Both distributions of the inferred coefficients are wide, a result that does not support the existence of an underlying exponential, or a power law dependent force.
https://doi.org/10.1371/journal.pcbi.1009153.s022
(TIF)
S11 Fig. Aggregate acceleration during merging when analyzed with respect to collision point.
This is an additional analysis performed on two-spot merging trajectories observed in experimental measurements. In this instance, the distance of each spot relative to the collision point versus time is fitted to a quadratic equation. Approximating the trajectory as constant acceleration motion, the acceleration is obtained from the quadratic equation. Finally, it is plotted with respect to the other merging aggregate area. The two quantities have a spearman rank-order correlation coefficient of 0.31 with a p-value of 0.002. The color of each data point corresponds to a merging event taking place in the respectively colored subregion in S5 Fig.
https://doi.org/10.1371/journal.pcbi.1009153.s023
(TIF)
S12 Fig. In this plot of experimentally derived quantities, for each two-spot merging event, the difference in average speed is plotted versus the difference of aggregate size between the two aggregates.
Small aggregates do not necessarily move faster than big aggregates during two-spot merging events. This is the case since the difference of speeds is not always positive, even when the size difference is significant.
https://doi.org/10.1371/journal.pcbi.1009153.s024
(TIF)
S13 Fig. In the simulation, introducing a motility transition dramatically increased the proportion of cells not bound in aggregates over long time scales.
(A) With no motility transition. (B) With the motility transition.
https://doi.org/10.1371/journal.pcbi.1009153.s025
(TIF)
S14 Fig. Effect of varying essential parameters of the simulation on spatial structure.
(A) Chemotactic attraction magnitude (B) Immotility transition threshold (C) Timescale of regaining motility (D) Timescale of random walk after regaining motility.
https://doi.org/10.1371/journal.pcbi.1009153.s026
(TIF)
S15 Fig. Effect of removing aspects of the simulation on spatial structure.
(A) Chemotactic attraction magnitude (B) Immotility transition threshold (C) Timescale of regaining motility (D) Timescale of random walk after regaining motility.
https://doi.org/10.1371/journal.pcbi.1009153.s027
(TIF)
S16 Fig. Effect of removing aspects of the simulation on motility and merging.
(A) Original simulation results for comparison (B) Removing the introduced motility transition (C) Removing random walk process after cells regain motility.
https://doi.org/10.1371/journal.pcbi.1009153.s028
(TIF)
S17 Fig. Evolution of motile and immotile cells in the simulations.
During the simulations, the distribution of motile and immotile cells reaches and retains a mixed steady state.
https://doi.org/10.1371/journal.pcbi.1009153.s029
(TIF)
S18 Fig. Identification of aggregates using Weka segmentation.
A) Experimental image of the entire plate at T = 22.5 h. B) Aggregates identified after using Weka software.
https://doi.org/10.1371/journal.pcbi.1009153.s030
(TIF)
S19 Fig. Segmentation of the simulation.
A) Density heat map of simulation results. For every pixel, the number of agents within a radius of 3 pixels is used to calculate the local cell density. B) Segmentation results after only considering pixels with a density value greater than 40 and aggregates with an area greater than 10 pixels.
https://doi.org/10.1371/journal.pcbi.1009153.s031
(TIF)
Acknowledgments
We would like to thank Marcel Meyer for helpful discussions and sharing the computational framework. Computation for the work described in this paper was supported by the University of Southern California’s Center for High-Performance Computing (carc.usc.edu). Osman Kahraman was also helpful in getting this project off the ground.
References
- 1. Budrene EO, Berg HC. Complex Patterns Formed by Motile Cells of Escherichia Coli. Nature. 1991 Feb 14;349(6310):630–3. pmid:2000137
- 2. Budrene EO, Berg HC. Dynamics of formation of symmetrical patterns by chemotactic bacteria. Nature. 1995 Jul;376(6535):49–53. pmid:7596432
- 3. Woodward DE, Tyson R, Myerscough MR, Murray JD, Budrene EO, Berg HC. Spatio-temporal patterns generated by Salmonella typhimurium. Biophys J. 1995 May;68(5):2181–9. pmid:7612862
- 4. Adler J. Chemotaxis in Bacteria. Science. 1966;153(3737):708–16. pmid:4957395
- 5. Curk T, Marenduzzo D, Dobnikar J. Chemotactic Sensing towards Ambient and Secreted Attractant Drives Collective Behaviour of E. coli. PLOS ONE. 2013 Oct 3;8(10):e74878. pmid:24098352
- 6. Lim S, Guo X, Boedicker JQ. Connecting single-cell properties to collective behavior in multiple wild isolates of the Enterobacter cloacae complex. PLoS ONE. 2019 Apr 4;14(4). pmid:30947254
- 7. Brenner MP, Levitov LS, Budrene EO. Physical Mechanisms for Chemotactic Pattern Formation by Bacteria. Biophys J. 1998 Apr 1;74(4):1677–93. pmid:9545032
- 8. Keller EF, Segel LA. Traveling bands of chemotactic bacteria: A theoretical analysis. J Theor Biol. 1971 Feb 1;30(2):235–48. pmid:4926702
- 9. Meyer M, Schimansky-Geier L, Romanczuk P. Active Brownian agents with concentration-dependent chemotactic sensitivity. Phys Rev E. 2014 Feb 13;89(2):022711. pmid:25353513
- 10. Igoshin OA, Mogilner A, Welch RD, Kaiser D, Oster G. Pattern formation and traveling waves in myxobacteria: Theory and modeling. Proc Natl Acad Sci. 2001 Dec 18;98(26):14913–8. pmid:11752439
- 11. Sager B, Kaiser D. Intercellular C-signaling and the traveling waves of Myxococcus. Genes Dev. 1994 Dec 1;8(23):2793–804. pmid:7995518
- 12. Manson MD. Introduction to bacterial motility and chemotaxis. J Chem Ecol. 1990 Jan 1;16(1):107–18. pmid:24264900
- 13. Berg HC. Chemotaxis in Bacteria. Annu Rev Biophys Bioeng. 1975;4(1):119–36.
- 14. Sourjik V, Wingreen NS. Responding to chemical gradients: bacterial chemotaxis. Curr Opin Cell Biol. 2012 Apr 1;24(2):262–8. pmid:22169400
- 15. Zhang X, Si G, Dong Y, Chen K, Ouyang Q, Luo C, et al. Escape band in Escherichia coli chemotaxis in opposing attractant and nutrient gradients. Proc Natl Acad Sci. 2019 Feb 5;116(6):2253–8. pmid:30674662
- 16. Be’er A, Ilkanaiv B, Gross R, Kearns DB, Heidenreich S, Bär M, et al. A phase diagram for bacterial swarming. Commun Phys. 2020 Apr 3;3(1):1–8.
- 17. Holmes AB, Kalvala S, Whitworth DE. Spatial Simulations of Myxobacterial Development. PLOS Comput Biol. 2010 Feb 26;6(2):e1000686. pmid:20195493
- 18. Ramaswamy S. The Mechanics and Statistics of Active Matter. Annu Rev Condens Matter Phys. 2010;1(1):323–45.
- 19. Buttinoni I, Bialké J, Kümmel F, Löwen H, Bechinger C, Speck T. Dynamical Clustering and Phase Separation in Suspensions of Self-Propelled Colloidal Particles. Phys Rev Lett. 2013 Jun 5;110(23):238301. pmid:25167534
- 20. Bär M, Großmann R, Heidenreich S, Peruani F. Self-Propelled Rods: Insights and Perspectives for Active Matter. Annu Rev Condens Matter Phys. 2020;11(1):441–66.
- 21. Vicsek T, Czirók A, Ben-Jacob E, Cohen I, Shochet O. Novel Type of Phase Transition in a System of Self-Driven Particles. Phys Rev Lett. 1995 Aug 7;75(6):1226–9. pmid:10060237
- 22. Giardina I. Collective behavior in animal groups: Theoretical models and empirical studies. HFSP J. 2008 Aug 1;2(4):205–19. pmid:19404431
- 23. Needleman D, Dogic Z. Active matter at the interface between materials science and cell biology. Nat Rev Mater. 2017 Jul 20;2(9):1–14.
- 24. Bechinger C, Di Leonardo R, Löwen H, Reichhardt C, Volpe G, Volpe G. Active particles in complex and crowded environments. Rev Mod Phys. 2016 Nov 23;88(4):045006.
- 25. Balaban V, Lim S, Gupta G, Boedicker J, Bogdan P. Quantifying emergence and self-organisation of Enterobacter cloacae microbial communities. Sci Rep. 2018 Aug 17;8(1):12416. pmid:30120343
- 26.
Okabe A, Boots B, Sugihara K, Chiu SN. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. John Wiley & Sons; 2009.
- 27.
Hansen J-P, McDonald IR. Theory of Simple Liquids. 4th ed. Elsevier, 2013.
- 28. Lochmann K, Oger L, Stoyan D. Statistical analysis of random sphere packings with variable radius distribution. Solid State Sci. 2006 Dec 1;8(12):1397–413.
- 29. Mittal N, Budrene EO, Brenner MP, Oudenaarden A van. Motility of Escherichia coli cells in clusters formed by chemotactic aggregation. Proc Natl Acad Sci. 2003 Nov 11;100(23):13259–63. pmid:14597724
- 30.
Micikevicius P. 3D finite difference computation on GPUs using CUDA. Proceedings of 2nd Workshop on General Purpose Processing on Graphics Processing Units. Washington, D.C., USA: Association for Computing Machinery; 2009. p. 79–84. https://doi.org/10.1145/1513895.1513905
- 31. Januszewski M, Kostur M. Accelerating numerical solution of stochastic differential equations with CUDA. Comput Phys Commun. 2010 Jan 1;181(1):183–8.
- 32. Green S. Particle Simulation using CUDA. 2010;12.
- 33. Lord ND, Norman TM, Yuan R, Bakshi S, Losick R, Paulsson J. Stochastic antagonism between two proteins governs a bacterial cell fate switch. Science. 2019 Oct 4;366(6461):116–20. pmid:31604312
- 34. Cotter CR, Schüttler H-B, Igoshin OA, Shimkets LJ. Data-driven modeling reveals cell behaviors controlling self-organization during Myxococcus xanthus development. Proc Natl Acad Sci. 2017 Jun 6;114(23):E4592–601. pmid:28533367
- 35. Zhang Z, Igoshin OA, Cotter CR, Shimkets LJ. Agent-Based Modeling Reveals Possible Mechanisms for Observed Aggregation Cell Behaviors. Biophys J. 2018 Dec 18;115(12):2499–511. pmid:30514635
- 36. Sliusarenko O, Zusman DR, Oster G. Aggregation during Fruiting Body Formation in Myxococcus xanthus Is Driven by Reducing Cell Movement. J Bacteriol. 2007 Jan 15;189(2):611–9. pmid:17098901
- 37. Kusaka K, Shibata K, Kuroda A, Kato J, Ohtake H. Isolation and characterization of Enterobacter cloacae mutants which are defective in chemotaxis toward inorganic phosphate. Journal of Bacteriology. 1997 Oct 1;179(19):6192–5. pmid:9324271
- 38. Kato J, Nagata C, Yang L, Kuroda A, Ikeda T, Takiguchi N, et al. Isolation and Characterization of the Enterobacter cloacae cheR Mutant Defective in Phosphate Taxis. Bioscience, Biotechnology, and Biochemistry. 2001 Jan 1;65(2):456–8. pmid:11302189
- 39. Sun E, Liu S, Hancock REW. Surfing Motility: a Conserved yet Diverse Adaptation among Motile Bacteria. Journal of Bacteriology. 2018. 200(23):e00394–18. pmid:30224438
- 40. Liu W-Y, Wong C-F, Chung KM-K, Jiang J-W, Leung FC-C. Comparative Genome Analysis of Enterobacter cloacae. PLOS ONE. 2013 Sep 12;8(9):e74487. pmid:24069314
- 41. Shastry RP, Welch M, Rai VR, Ghate SD, Sandeep K, Rekha PD. The whole-genome sequence analysis of Enterobacter cloacae strain Ghats1: insights into endophytic lifestyle- associated genomic adaptations. Arch Microbiol. 2020 Aug 1;202(6):1571–9. pmid:32166358
- 42. Dahl MK, Boos W, Manson MD. Evolution of chemotactic-signal transducers in enteric bacteria. J Bacteriol. 1989 May 1;171(5):2361–71. pmid:2496104
- 43. Kearns DB. A field guide to bacterial swarming motility. Nat Rev Microbiol. 2010 Sep;8(9):634–44. pmid:20694026
- 44. Kalluri R, Weinberg RA. The basics of epithelial-mesenchymal transition. J Clin Invest. 2009 Jun 1;119(6):1420–8. pmid:19487818
- 45. Meacock OJ, Doostmohammadi A, Foster KR, Yeomans JM, Durham WM. Bacteria solve the problem of crowding by moving slowly. Nat Phys. 2021 Feb;17(2):205–10.
- 46. Bi D, Yang X, Marchetti MC, Manning ML. Motility-Driven Glass and Jamming Transitions in Biological Tissues. Phys Rev X. 2016 Apr 21;6(2):021011. pmid:28966874
- 47. Grobas I, Polin M, Asally M. Swarming bacteria undergo localized dynamic phase transition to form stress-induced biofilms. Giardina I, Walczak AM, editors. eLife. 2021 Mar 16;10:e62632. pmid:33722344
- 48. Douarche C, Buguin A, Salman H, Libchaber A. E. Coli and Oxygen: A Motility Transition. Phys Rev Lett. 2009 May 12;102(19):198101. pmid:19518998
- 49. Norman TM, Lord ND, Paulsson J, Losick R. Memory and modularity in cell-fate decision making. Nature. 2013 Nov;503(7477):481–6. pmid:24256735
- 50. Lyons NA, Kolter R. On the evolution of bacterial multicellularity. Curr Opin Microbiol. 2015 Apr 1;24:21–8. pmid:25597443
- 51. Schindelin J, Rueden CT, Hiner MC, Eliceiri KW. The ImageJ ecosystem: An open platform for biomedical image analysis. Mol Reprod Dev. 2015;82(7–8):518–29. pmid:26153368
- 52. Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al. Fiji: an open-source platform for biological-image analysis. Nat Methods. 2012 Jul;9(7):676–82. pmid:22743772
- 53. Arganda-Carreras I, Kaynig V, Rueden C, Eliceiri KW, Schindelin J, Cardona A, et al. Trainable Weka Segmentation: a machine learning tool for microscopy pixel classification. Bioinformatics. 2017 Aug 1;33(15):2424–6. pmid:28369169
- 54. Tinevez J-Y, Perry N, Schindelin J, Hoopes GM, Reynolds GD, Laplantine E, et al. TrackMate: An open and extensible platform for single-particle tracking. Methods. 2017 Feb 15;115:80–90. pmid:27713081
- 55.
Schmool DS, Kachkachi H. Chapter One—Collective Effects in Assemblies of Magnetic Nanoparticles. In: Camley RE, Stamps RL, editors. Solid State Physics. Academic Press; 2016. https://doi.org/10.1038/srep22872 pmid:26976721