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][2][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][2][3][4]. These collective phenomena can even be predicted with analytical mathematical considerations and recreated with detailed computational models [3,[5][6][7][8][9][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][2][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.
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 and1C. Images of aggregates were analyzed using image analysis software as discussed in Materials and Methods. 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. 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 100mm 2 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.
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 longrange 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 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. https://doi.org/10.1371/journal.pcbi.1009153.g002 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.
After aggregate formation, a displacement of some of the spot aggregates over time was observed, as shown in 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. 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.
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. 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.
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 and4B). 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).
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 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 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.
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][31][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 densitydependent 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 and6C, 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.
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 and6H, 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 rankorder 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 nonmerging aggregates (Fig 6G and6H) 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 and5E), and the respective motility profile in the simulations (Fig 6G and6H) 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 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.
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][35][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 epithelialmesenchymal transition [44]. The active matter community has examined the cell motility and pattern formation [45][46][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.
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].
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.
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.
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.
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.
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, timelapse 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.
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 r max 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 r max in increments of dr. Furthermore, the probability density was divided with the number density of the pattern, ρ = N total /L 2 , where N total 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 r in = n � dr and an outer radius of r out = (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: sðkÞ ¼ 1 þ 2pr=k Z r max 0 ðgðrÞ 1ÞsinðkrÞdr [55].
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.
PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1009153 January 4, 2022