Toward a realistic model of multilayered bacterial colonies

Bacteria are prolific at colonizing diverse surfaces under a widerange of environmental conditions, and exhibit fascinating examples of self-organization across scales. Though it has recently attracted considerable interest, the role of mechanical forces in the collective behavior of bacterial colonies is not yet fully understood. Here, we construct a model of growing rod-like bacteria, such as Escherichia coli based purely on mechanical forces. We perform overdamped molecular dynamics simulations of the colony starting from a few cells in contact with a surface. As the colony grows, microdomains of strongly aligned cells grow and proliferate. Our model captures both the initial growth of a bacterial colony and also shows characteristic signs of capturing the experimentally observed transition to multilayered colonies over longer timescales. We compare our results with experiments on E. coli cells and analyze the statistics of microdomains.


Introduction
Microbial life in general exists not as cells freely floating in an aqueous environment (planktonic state), but rather as sessile, surface associated colonies [1].These sessile, self-organized communities often develop into biofilms, defined as matrix-enclosed colonies adhering to each other and/or to surfaces [1,2].This matrix, called extracellular polymeric substances (EPS), is composed of biopolymers, nucleic acids, and lipids which produce cohesion of the biofilm, adhesion to the substrate, and a number of protective functions [3].Starting from the initial attachment, followed by the growth, proliferation, and ultimately colonization, biofilms represent self-regulated active microbial communities where an interplay of geometry, order and growth-mediated mechanics determine their emerging structure and topology [4][5][6][7][8][9][10].From dental surfaces and open wounds, cancer and gut environments to the biological growth observed on maritime structures (biofouling), biofilms can grow and thrive under diverse settings, thereby regulating infectious diseases and human health, as well as determining the durability of maritime vessels [11] and medical devices such as intravenous catheters, prosthetic heart valves, joint prostheses, peritoneal dialysis catheters, and cardiac pacemakers [12].Importantly, biofilms can also impact the respiratory mucous and the airways of patients with cystic fibrosis, thus regulating patient health and response to medical treatment [13].Because of the specificity of the physical and biological interactions among its constituent cells, biofilms are often considered an emergent form of microbial life, as the properties and function of biofilms may outperform those of the constituent cells [12].
The first step in the development of a biofilm is the growth of a colony of cells multiplying on a substrate, initially at an exponential rate, until nutrient limitations cause a plateauing of the colony size.In this initial phase of rapid expansion, biofilms self-organize into monolayers wherein contact forces between neighboring cells and with the substrate play a prominent role [14,15].These contact forces lead to an orientational transition reminiscent of liquid-crystalline phase transitions.As the degree of confinement decreases, rod-like cells transition from a long-ranged nematic state to a globally disordered state characterized by microdomains with short-range nematic order [14,[16][17][18].
Considerable work has focused on the transition from a two-dimensional monolayer of cells to a three-dimensional (3D) biofilm.A crucial step in the development of a 3D structure is the verticalization of cells [15,19], where the growth of mechanical stresses among the growing cells leads to a buckling instability that turns some cells away from a planar configuration parallel to the substrate [15].
Recent work has shown that, following verticalization, the second layer spreads over the first layer, and thereafter the third layer starts to develop on the second layer and so on [20].Over time, this monolayer stacking process leads to the 3D morphology of the bacterial colony emerging from multiple flat layers of dividing cells parallel to the underlying flat substrate.Developing a model for this monolayer stacking transition is the main goal of this work.

Bacterial model
We model each bacterial cell as a spherocylinder [21] of length  (), depending on time , and diameter  capped at both ends by hemispheres of the same diameter .Each cell is described by the vector to its center-of-mass position  and a unit vector for its orientation ê.In the overdamped limit, the equations of motion read where Γ and Ω are characteristic translational and rotational drag coefficient, respectively,  is the identity matrix [22].In equations (2.1)-(2.2)  = −0.5 ẑ represents a constant adhesion force applied to center of each rod and  = −( + ) ( ê • ẑ) ẑ is a restoring torque that returns our rods to planar alignment, respectively, and ẑ is the unit vector normal to the substrate.Both  and  physically stem from the adhesive forces within the EPS matrix produced by and surrounding the cells.Additionally, to model thermal fluctuations, we include a small fluctuating torque with amplitude ∝ 10 −7 .In equation (2.1), , where    ≡ (  −   )/   , and    = |  −   | is the distance between the center of the cells;  H represents a Hertzian repulsion force between overlapping spherocylinders with cell stiffness  0 = 100 (dependent on Young's modulus) [19,23]; this repulsion is calculated by means of virtual spheres of diameter  centered at the closest points between the axes of spherocylinders so that the cylinders are in contact whenever the distance between the centers of the virtual spheres r  is such that r  −  ≡    < 0. The smallest distance between two spherocylinders is computed using the algorithm by Vega and Lago [24].The substrate also exerts a Hertzian repulsion  Hs =  0  1/2  3 2 s ẑ, where  s is the overlap of the spherocylinder with the solid substrate.In equation (2.2),  H = ( c − ) ×  H , where  c is the vector to the contact point.
We use the following potentials to model the cell-cell interaction  cc acting between each bacterium and their neighboring cells and the cell-substrate potential  cs with the confining surface [22] In equation (2.3),  cc represents a short-range cell-cell attraction with relative strength  cc = −10 −3 , where  = 2.83 is the attraction position,  = 0.16 the width of the potential [22], whereas  cs represents the attraction between cell and substrate modelled through a Yukawa-like potential with  cs = −2.5,where  = 1 is the depth of Yukawa potential, ℎ is the relative distance between cell and substrate, and is ℎ 0 the position of the substrate in our simulations.
Bacterial cell proliferation is the driving force underpinning mechanical interactions in a colony [25][26][27].Away from division events, the cell length growth obeys a linear law (with coefficient   ) modified to include a mechanoresponse () of the cell to the growing pressure exerted by other cells due to proliferation stresses d d

.4)
A cell divides when it reaches a split length   = 6, a small random offset in the daughter cells' length is also used.In equation (2.4),   = 10 −7 is a growth rate constant.It has been observed that when colonies grow, the increasing pressure exerted by each cell on its neighbors leads to a slowdown in the cell length growth, akin to a contact inhibition [23,25]; the term in square brackets is a sigmoid function that models such mechanoresponse to the growing pressure, where  = 60 represents a scale parameter regulating the steepness of this mechanoresponse,  is the total overlap between a cell and its neighbors, and  0 = 0.12 is the inflection point of this dependence.To avoid artificial synchronization effects in cell divisions across the colony, when a cell divides, the two daughter cells differ in length by a random amount with 5% maximum amplitude.
Figure 1 shows a schematic illustration of two spherocylinders; the thin yellow circles represent the virtual spheres used to compute the overlap and the Hertzian repulsion  H between the bacterial cells.We perform molecular dynamics simulations using an Euler integration scheme with a time step Δ = 10 −5 .

Results and discussion
Figure 2 shows a sequence of snapshots from a simulated colony growing from a starting configuration of two randomly oriented cells that, after a number of cycles of growth and division reaches a maturation stage comprising multiple layers [15].In agreement with experimental evidence [28,29], the colony initially grows as a flat, circular monolayer.We note that our simulations reproduce the tangential alignment of cells at the outer edge of the colony.This is due to the fact that a tangentially-oriented rod is mechanically stable to the forces of neighboring cells, whose torques cancel out statistically [28].
As the colony grows, the mechanical stresses mutually exerted on the cells increase, until a threshold value is met and the verticalization transition is crossed [19], resulting in the extrusion of some cells from the initial monolayer.Figure 2(d) shows that some cells have formed the seeds of a flat, second layer, growing directly on top of the initial layer with all cells aligned parallel to the underlying substrate [20].Figure 3 shows the configuration of a colony which has developed a large second layer, visible from the central darker region in figure 3(a) and from the edge-view in figure 3(b).
Close inspection of the colonies reveals that neighboring cells exhibit a large degree of local alignment.Microdomains, best captured in bacterial monolayers, are defined as clusters of cells that display a large local nematic order parameter [14].Figure 4(a) shows a micrograph of bacterial colony captured experimentally, using E. coli (strain C600-WT) bacteria grown on nutrient-rich agar surface.The nonmotile strain of E. coli bacteria were grown at around 30℃, showing a doubling time of approximately 40 minutes.The duration of our experiments were short enough (up to a few hours), ensuring that the expanding monolayers remained nutrient-replete throughout the observation time-scale.The growth of a single bacterium into bacterial monolayers was imaged using time-lapse microscopy (phase contrast mode), acquiring images at 5 minute intervals.Initially, the colony expanded as a monolayer (2D), spanning multiple generations, before extruding into the third dimension to reach a complex 3D structure.The color coding shows that a large degree of local nematic order is present in the monolayer system, and the local orientation can be used to define microdomains.As can be observed in Figure 4(a, b), the intersection of three or more nematic microdomains (i.e., intersection of three or more colored clusters), typically represents topological defect sites [30], the dynamics of which have been discussed in detail elsewhere [6,20].In order to investigate more closely the presence of microdomains, we employ a density-based clustering algorithm.We consider two cells to belong to the same cluster if the angle between their orientations    < 5 • and their center-of-mass distance    < 0.75  .Using these metrics, we use the density-based spatial clustering of applications with noise (DBSCAN) algorithm [31].Figure 4(b) shows the results of the DBSCAN cluster analysis on a representative model colony.Similarly to the experimental results, there is a broad distribution of microdomains, differing in size and detailed morphology.Both simulated and experimental microdomains display a large local nematic order and also large local smectic order [32,33] (though the latter is less pronounced in the experiments, due to the intrinsic biological variability of the cell properties).
Figure 5 shows the probability distribution of the microdomains' area ( ) resulting from the clustering analysis described above and averaged over 39 independent simulations.Our simulations show an exponential distribution ( ) ∼ e − / * in agreement with experimental measurements [14]; furthermore, our estimate for the characteristic domain area  * ≈ 63 2 , which implies that the typical size of these microdomains is  * ≈ 8, which also agrees with the measurement [14].

Figure 1 .
Figure 1.(Colour online) Schematic of two spherocylinders representing rod-like bacteria (e.g., E. coli cells) of length  (excluding the hemispherical caps) and diameter .Thin yellow circles represent the virtual spheres used to compute the Hertzian repulsion  H based on the overlap of the spheres    .

Figure 2 .
Figure 2. (Colour online) Different stages of the colony's life.Starting from an initial condition with two cells randomly oriented (a), the cells undergo different cycles of growth and division (b) and then (c).Eventually, the colony assumes a circular shape with most cells at its boundary aligned tangentially to the colony (d).Each snapshot indicates the number of cells  present.The underlying substrate is not shown.

Figure 3 .
Figure 3. (Colour online) (a) Top view of a large bacterial colony which underwent mono to multiplayer transition.The second layer of cells is visible as the darker region in the center.(b) Edge view of the same colony showing the presence of extruded cells which have divided to form a flat second layer.The number of cells  is indicated in the snapshot.

Figure 4 .
Figure 4. (Colour online) Comparison of microdomains in experimental and simulated bacterial colony.(a) Micrograph of E. coli (strains C600-WT) bacteria grown on nutrient-rich agar plate.The cells are color-coded to visualize local cell orientations.(b) Simulated colony of our rod-like cell model for a single monolayer.The color encodes membership to microdomains, identified via a density-based clustering algorithm (see text for details).Microdomains are characterized by both a large local nematic as well as local smectic order.

Figure 5 .
Figure 5. (Colour online) Probability distribution of the microdomains' area in mature 2D colonies before the first verticalization event.Data are averages over 39 independent simulations.The solid line is an exponential fit to the data.