Mechanical forces drive a reorientation cascade leading to biofilm self-patterning

In growing active matter systems, a large collection of engineered or living autonomous units metabolize free energy and create order at different length scales as they proliferate and migrate collectively. One such example is bacterial biofilms, surface-attached aggregates of bacterial cells embedded in an extracellular matrix that can exhibit community-scale orientational order. However, how bacterial growth coordinates with cell-surface interactions to create distinctive, long-range order during biofilm development remains elusive. Here we report a collective cell reorientation cascade in growing Vibrio cholerae biofilms that leads to a differentially ordered, spatiotemporally coupled core-rim structure reminiscent of a blooming aster. Cell verticalization in the core leads to a pattern of differential growth that drives radial alignment of the cells in the rim, while the growing rim generates compressive stresses that expand the verticalized core. Such self-patterning disappears in nonadherent mutants but can be restored through opto-manipulation of growth. Agent-based simulations and two-phase active nematic modeling jointly reveal the strong interdependence of the driving forces underlying the differential ordering. Our findings offer insight into the developmental processes that shape bacterial communities and provide ways to engineer phenotypes and functions in living active matter.

Cell division is modeled as the instantaneous replacement of a mother cell when it reaches length by two equal sized daughter cells. The slight volume loss after division is ignored to avoid contact overlapping which may cause unphysical reorientations. It follows that the initial length of a daughter cell is . Given the growth law, the doubling time can be calculated as log .

Cell-cell repulsion
We only consider the repulsive, elastic contact forces between cells, mediated by the soft exopolysaccharides. Linear elastic Hertzian contact theory is applied to quantify the repulsive contact force on cell by cell , written as where is the effective cell stiffness, is the overlapping distance, and is the unit vector parallel to the distance vector . The distance vector is given by connecting two contact points characterizing the smallest distance between two cell centerlines, as shown in Supplementary Fig.  6. The overlapping distance is then calculated by and contact only occurs if 0. Taking the center of mass as the reference point, the moment of contact force is where is the cell director and is the parametric coordinate, along the center-line of the cell, of the contact point.

Cell-to-substrate interactions
The experimental substrate is modeled as an infinite, two-dimensional rigid plane located at 0. We again assume linear elastic Hertzian contact mechanics to quantify the repulsive contact between cells and the substrate. We also add an attractive force which depends on the contact area to account for the matrix protein mediated cell-to-surface adhesion.
To formulate the repulsive contact between the cells and the substrate, we apply a generalized Hertzian contact formula that smoothly accounts for the cell orientation dependent contact energy. In this case, the elastic deformation energy is given by , / / and the equivalent penetration depth is given by where , is the projection of the th cell director onto the -axis. The overlap function denotes the overlapping distance between the cell and the substrate at the local cell-body coordinate -/2 /2. Then, the net force , and moment , from the cell-substrate elastic repulsion can be given by where , and , are the total force and moment vector, and is the moment of inertia. All of the variables are expressed in the body-fixed coordinate system, then transformed into the global coordinate system. We add a small random noise to the net force and moment vectors of the cells at every timestep (10 for forces and 10 for moments).

Choice of parameters for cellular dynamics
The list of physical constants, which were used to successfully capture the verticalization instability in biofilm-dwelling cells in prior work 1 , are used here and given in Table S2.

Geometry of coarse-grained gel system
We model the agarose gel using a coarse-grained particle-based approach. We treat the gel as a collection of spherical particles, each with radius . We treat the interactions of the spheres using a spring network model to recapitulate the elastic behavior of the hydrogel. The pairwise interaction energy between the gel particles is , Σ , where is the distance between particle and , is the equilibrium distance, and is the spring constant. Furthermore, to impart a shear modulus to the system, we also include a three-body interaction energy where , Σ , where is the bond angle formed by particle , , and .

Interactions between gel particles and cells
We again apply linear elastic Hertzian contact theory to describe the repulsive interaction between the coarse-grained gel particles and the cells. Similar to the contact between cells, the elastic contact energy can be written as , / / , where is the contact stiffness between gel and cell, and is the overlapping distance between the th gel particle and th cell. The contact stiffness can be given by the relation , where and are Young's modulus and Poisson's ratio, respectively.

Interactions between the gel particles and the substrate
There exist two types of interactions between the coarse-grained gel particles and the substrate. First, the gel particles make elastic contact with the substrate, again described by Hertzian contact theory for the contact between a sphere and a flat surface, where the elastic contact energy is given by , / / . Second, we introduce adhesion between the gel particles and substrate which provides an energy barrier to delamination; we take the energy to be of the form , , Σ , , where Σ is the adhesion energy density. The equivalent contact area is given by , π , where is the overlap between the gel particle and the substrate.

Equation of motion
We only consider the three translational degrees of freedom for gel particles, neglecting the rotational degrees of freedom. Therefore, the equations of motion are given by Newton's second law , , where , is the net force and is the acceleration. To prepare the initial amorphous stress-free geometry, we begin with a body-centered cubic crystalline geometry with lattice parameter , where 1.3 . Subsequently, we assigned the system with an initial temperature of 300 and annealed it until it reached a final configuration that is amorphous and stress-free.

Choice of parameters
The spring constant and equilibrium length : The Young's modulus of the coarse-grained gel system is given by , , ≃ 2 , under the condition ≪ . Generally, a smaller leads to a denser gel system and better approximation to a continuum solid.
Here we choose 0.6 μm as a result of a trade-off between simulation quality and computational cost, as the simulation time is proportional to . This value leads to 6 10 Nm , corresponding to the experimentally measured 5 . Radius of coarse-grained gel particle : In order to mimic the continuum constraints posed by the hydrogel in the experiment, the coarse-grained system should not have significant defects larger than the volume of a single cell, which means that should be larger than . On the other hand, cannot be significantly larger than the cell radius , as this will introduce unphysical contacts at the biofilm-gel interface. Taking both requirements into consideration, we choose 1.0 μm, which is nearly double the equilibrium distance 0.6 μm and we keep . Cell-gel contact modulus : Hertzian contact mechanics gives the equivalent contact modulus between two elastic bodies by , where is the Young's modulus and is Possion's ratio of the two contacting elastic bodies. Given the bulk rheology measurement 500 Pa, μ μ 0.49, and 5 kPa 3,4 , it follows 600 Pa. In the simulations, we choose 1500 Pa, somewhat larger than , to avoid the unphysical situation in which cells penetrate into the gel during growth. The three-body interaction parameter and equilibrium angle : To make the amorphous coarsegrained gel system stable, we choose 120° to force the system to deviate from its original cubic configurations. We select such that we attain incompressible solid behavior in the gel /3 , where is the shear modulus and is Young's modulus. These parameters are listed in Table S3.

III. Single-cell surface adhesion stability analysis
Previous studies have described how cells in 2D bacterial colonies can be peeled from the substrate through a buckling-like instability [5][6][7][8] ; however, the reason why biofilm-dwelling cells self-organize into a verticalized core, that is, remain attached but oriented away from the substrate, remains to be shown. In this section, we consider the stability of surface-adhered cells subject to compression from neighboring cells to determine why surface adhesion facilitates stably surfaceanchored, verticalized cells. Specifically, we consider a minimal model examining the instability of a single spherocylindrical cell that is compressed by its neighbors, following the same physical model described in the agent-based simulation section above; however, in our theoretical analysis we neglect any external gel confinement, and only consider the effects of bulk viscosity for analytical tractability. We assume that the central cell is not growing and instead we mimic the effects of growth-induced compression by bringing neighboring cells closer together. The central cell initially starts out parallel to the substrate in close contact with its neighbors, and the peripheral neighboring cells are brought closer to each other to squeeze the central cell (see the schematic in Figure 2f). We find that the central cell evolves through four phases upon increasing compression: I. the cell remains horizontal and the elastic potential energy in the cell increases; II. the central cell rotates from horizontal to vertical and all potential energy is consequently released; III. the central cell remains vertical, and the elastic potential energy increases again when the outer cells touch the center cell again; and IV. the central cell becomes unstable and is ejected from the substrate. We confirm this sequence with a simple 2D agent-based model where cells are confined to the -plane. The transition from I to II corresponds to a "verticalization" instability 1 , from II to III corresponds to a trivial geometric transition, and from III to IV corresponds to a "pinch-off" instability. Note that the pinch-off instability is locally (marginally) stable, but under sufficiently large -perturbations, the central cell gets ejected from the substrate. Experimentally, we expect that many factors could contribute to the perturbations including but not limited to: deviation of the cell shape from a perfect spherocylinder, fluctuating protein bonds, and variations in the adhesion protein concentration, etc. In what follows below, we consider these two instabilities in more detail paying specific attention to the role of the cell-to-surface adhesion.
The equations of motion for the position , , , and director , , ∥ of a cell written using Lagrangian mechanics, are: where is a Lagrange multiplier corresponding to the constraint | | 1 0. Here ∑ is the total potential energy and includes cell-to-cell repulsion, as well as cell-to-surface adhesion and repulsion terms. These surface interaction terms are given by In phase I, in the absence of any compression from neighboring cells, the balance of surface adhesion and repulsion leads to an equilibrium cell configuration 0 at an apparent height where Σ / is the dimensionless adhesion coefficient. In the presence of compression, we consider perturbations away from this stable point.
To simplify the analysis, we assume that all neighboring cells compress and lead equally to a buildup of elastic energy or pressure in the central cell of interest ( / / where is the overlap distance between two cells). We assume hexagonal close packing around the caps of the cell and square close packing around the cylindrical portion of the cell leading to a maximum of ∼ 6 / neighbors compressing the central cell (see schematic in Supplementary Fig 10a). Finally, we make the simplifying assumption that the dominant compression that leads to verticalization is due to the two neighbors pushing on the two ends of the central cell. Note that in practice the exact configuration around the central cell will vary, leading to a stochastically varying stress field; our analysis is meant to provide an estimate of the relative energy barriers between the verticalization and pinch-off instabilities.
Combining the cell-surface adhesion, cell-surface repulsion, and cell-cell repulsion terms, we find the evolution equation for small is given by where the sign of Ω determines the stability of the fixed-point. Supplementary Figure 10c shows plots of Ω for different dimensionless adhesion coefficients. For small values of , which are most physically relevant, the zeros of Ω are only weakly dependent on . This analysis suggests that the verticalization instability only weakly depends on adhesion, that is, the energy barrier to verticalization is nearly constant. This counterintuitive result arises because for ≪ 1, adhesion does not impart a significant torque on the cell. Instead, while compression from peripheral cells is destabilizing, cell-surface repulsion acts to provide a stabilizing torque. In phase III, again assuming no compression from neighboring cells, the balance of surface adhesion and repulsion leads to an equilibrium verticalized configuration 1 with height Note that the effective penetration depth (the degree to which the soft cell is flattened against the substrate) is larger for verticalized cells than for horizontal cells. This difference in penetration depth occurs due to the difference in geometry between a hemispherical cap and a cylindrical body and how the cell-surface adhesion and cell-surface repulsion energies scale with the penetration depths. Again, we consider compression due to neighboring cells (at a height , ). In this case, the footprint of the verticalized cell is a circle and so we assume 6 neighboring cells that are hexagonally closed packed around it (Supplementary Fig 10b). In this configuration all neighbors contribute equally to the pinch-off of the central cell. Here we look for finite-sized perturbations away from the fixed point; the evolution equation of the height of the cell of interest is given by , .
We look for points , such that , 0. We only consider points 0 , as points , results in no net vertical force such that the cell of interest remains stably attached to the substrate. Specifically, we focus on points very close to /2 to determine the critical energy required to produce marginally unstable cells. Supplementary Figure  10d shows plots of , /2 0.999 for different adhesion coefficients. In contrast to the verticalization transition, in this case we find a strong dependence of the zeros of on ; therefore, the energy barrier to pinch-off depends strongly on cell-to-surface adhesion.
We combine these two stability analyses into a single plot of the energy landscape experienced by the central cell for different adhesion energies (Supplementary Figure 10e). We see that for very small adhesion energies, pinch-off is energetically favored over verticalization, whereas for larger adhesion energies, pinch-off is energetically more costly than verticalization. Importantly, we find that for the parameters listed in Table S2, the critical energy corresponding to the verticalization instability is much smaller than that of the pinch-off instability, leading to preferential verticalization of horizontal cells rather than pinch-off of verticalized cells in a population of cells in the basal layer (Fig. 2f). This analysis explains the observed positive correlation between the adhesion energy and the fraction of stably verticalized cells ( Supplementary Fig. 5, 7, 8): adhesion energy affects the stability of cells with respect to the pinchoff instability but much less so to the verticalization instability.

IV. Continuum model for growth-induced macroscopic cell ordering
In this section, we present a minimal coarse-grained description to explain the growth and self-organization observed in the basal layer of V. cholerae biofilms. We first consider the simpler case with one population of cells that either generate in-plane growth or not, before considering the complete two-phase model which accounts for pressure-dependent cell verticalization and outof-plane growth.

Theoretical framework
We limit our consideration to the basal plane of the biofilm and assume it is a growing, quasi-2D system with mesoscopic nematic order tensor 2〈 ⨂ /2〉 9 . Based on its symmetries, in 2D, can also be written as cos 2 sin 2 sin 2 cos 2 where is the angle of the average head-less director and q is the scalar order parameter, which quantifies the degree to which the cells are locally ordered. In the analysis below we will assume that this scalar order parameter is constant. We assume the bacterial cells grow with a uniform growth rate but that the growth can be out-of-plane, resulting in a spatially varying in-plane growth rate . Taking the cell density as , and the coarse-grained growth-induced velocity field as , mass conservation requires The density is nearly uniform in the experiment, that is, , where ≪ 1. To leading order, Eq. (22) becomes The growth-induced velocity field is found by integrating the spatially dependent growth rate. In this analysis, we will neglect density fluctuations , but it can be related to the fluctuation in local pressure due to varying local configurations. The evolution of the nematic order tensor is described by the Beris-Edward equation 10 where is the flow-alignment parameter and 0 for rod-shaped objects that tend to align in shear. Here ∇ ∇ ∇ ⋅ and ∇ ∇ are the traceless strain-rate and vorticity tensors and is the molecular field. Γ relaxes towards a bulk state with minimal spatial variation in the director. However, we do not expect Γ to play an important role in biofilm ordering because biofilm-dwelling cells secrete exopolysaccharides which are soft and deformable and act as a "cushion" between cells. This can be seen in the biofilm from the ∆BC mutant (Fig.  1e) in which minimal local alignment was observed. Contrast this with the ∆vpsL mutant and other bacterial colonies [11][12][13] where there is significant local alignment (Fig. 1f) when no exopolysaccharides are present. Finally, force balance requires where is the stress tensor in the biofilm, whose divergence is balanced by surface friction. Here we have assumed that energy dissipation is dominated by surface friction rather than viscous dissipation inside the biofilm, corresponding to the so-called "dry" limit 14 . We take the frictional force density to be linearly proportional to the velocity and the friction coefficient to be isotropic and constant. In many active nematic systems, the active stress is anisotropic and dependent on . For instance, in motile cell colonies, cells generate anisotropic active stresses [15][16][17] while in growing bacterial colonies a number of different constitutive relations have been proposed 11,12,18,19 .
Although at the microscopic scale, these forces may be anisotropic, on the mesoscopic scale pertaining to our theory, they are expected to be isotropic on average. This leads to the assumption of an isotropic stress field which is balanced by surface friction, Pressure in this system arises due to compression from neighboring cells, mediated by the exopolysaccharides, which arises from cell proliferation. To close the model, we impose the following set of boundary and initial conditions. By symmetry, we expect that the velocity goes to zero at the center of the biofilm, that is, the center of the biofilm does not move and , .
We also assume that the pressure accumulates towards the center of the biofilm and that at the outer edge of the biofilm the pressure is zero, where denotes the edge of the biofilm. The outer radius evolves through the kinematic condition .

(30)
Finally, we assume that the orientation of the biofilm is described by some initial order parameter tensor which can also be written as , , 0 , . Equations (23), (24), (26) and boundary conditions and initial conditions Eq. (27)-(31) form the basis of the theoretical model. In the following subsections, we consider a few idealized cases to determine under which conditions biofilms self-organize into an aster pattern. Throughout, we will also make the simplifying assumption that all variables except for are axisymmetric and that there is no azimuthal flow, .

Uniform growth limit
If all growth is in-plane, then is uniform and constant. Integrating Eq. (23) we find the velocity field is In this case, is simply stretched by growth, and there is no tendency for to rotate or align in any given direction. An initially isotropic system will remain isotropic as it grows ( Supplementary  Fig. 14a-c) 12,20 . This is what happens in the 2D expansion of bacterial colonies, which remain macroscopically isotropic as they grow. Although microscopic forces can align cells locally, they cannot introduce any long-range order. Note that in two dimensions axisymmetric growth is a special case where there is no net strain in the system: if instead the system were confined into a rectangular channel, where growth is unidirectional, the strain rate tensor would not be zero and therefore there would be net alignment along the direction of the channel 18,19,21,22 .

Alignment due to a growth void
Now suppose there is a radius inside which the cells have verticalized and are angled out of plane. These cells do not contribute to the growth of the basal plane and so Substituting into Eq. (23), the growth-induced velocity is The corresponding strain rate and vorticity tensors in polar coordinates ( , ) are Under uniform growth / and / are equal to /2 , leading to no net deviatoric strain-rate. In contrast, in the presence of a growth void, these two terms become unequal, leading to a net deviatoric strain-rate. This is because the radial gradient in velocity is no longer balanced by azimuthal stretching due to outward expansion. For a velocity of the form given in Eq. (35), the strain rate and vorticity tensors in polar coordinates are , if Next, the nematic order parameter can be transformed to a polar coordinate system, given by where is the azimuthal angle and cos cos is the coordinate transformation matrix from cartesian to polar coordinates. Finally, given , and , the matrix form of Eq. (24) in polar coordinates is ∂ ∂ ∂ cos 2 θ sin 2 θ sin 2 θ cos 2 θ 2 for , which yields two scalar equations cos 2 θ ∂θ ∂ ∂θ ∂ 0.

(40A,B)
Defining Θ and combining the two scalar equations above, we find that for , By first ignoring the advective term, we see that Θ ∼ sin 2Θ which has stable fixed points at Θ or , characteristic of an aster pattern. In other words, the rightmost term, which arises due to gradients in the growth-induced velocity field, drives cells to reorient towards an aster pattern. The same observation was made during inward growth in bacterial colonies where deviation of the velocity field from that of isotropic expansion resulted in radial alignment of cells 20 .
More generally, if we suppose that instead of a growth void, there is a general differential growth rate: for and Δ for then the evolution equation for Θ for is given by If the biofilm has a faster growing outer rim Δ 0, the bacteria will tend to organize into an aster pattern. If, on the other hand, the biofilm has a faster growing core, i.e. Δ 0, Θ will have stable fixed points at Θ /2, characteristic of a vortex pattern. This means that biofilms that exhibit excess growth in the core will be driven to form a vortex.
Returning to Eq. (41), we solve it by the method of characteristics, assuming some initial maximum biofilm radius , which yields where / and Θ corresponds to the angular representation . Plots of Eq. (43) are given in Supplementary Figures 14d-f. Note that, although cells tend to align in an aster pattern, the aligning torque becomes weaker as cells are advected outwards. Therefore, does not necessarily reach a value of 1 ( Supplementary Fig. 14g). Depending on the initial size of the growth void as well as the flow alignment parameter, different degrees of radial alignment are reached in the final state, with a larger / leading to overall more radial alignment.
To be more general, now instead of a finite void, suppose there is an arbitrary velocity field , corresponding to an arbitrary growth rate , . Following the same procedure as above, we find the evolution equation for Θ is where /4 / quantifies the aligning torque due to gradients in the flow field. Equation (44) thus describes a generalized time-evolution equation for the orientation field due to gradients in the flow velocity.

Two-phase active nematic model for cell verticalization and alignment
We expand the active nematic model in the previous subsection to include the verticalization dynamics in a two-phase model consisting of vertical and horizontal cells. Consider a quasi-2D layer of cells on a substrate where the cells are either parallel (horizontal) or tilted (vertical) with respect to the substrate and let denote the fraction (0 1) of vertical cells and 1 therefore denote the fraction of horizontal cells. Horizontal cells proliferate and generate their progeny in the basal layer of interest, while vertical cells do not contribute to the expansion to the basal layer. The growth-induced stress generated by the horizontal cells also cause them to verticalize with a pressure dependent rate . The evolution equation for the fraction of these two populations of cells is therefore given by: Given that cell verticalization requires a minimum threshold pressure , we make the simple assumption that the conversion rate is a linear function of the excess pressure We initiate the biofilm with some initial radius , and assume all cells are horizontal to begin with, , 0 0. Finally, we assume symmetry at the origin which requires 0 and that is axisymmetric. Equations (26), (44)-(47) with the above mentioned initial and boundary conditions make up the full two-phase continuum model for cellular ordering driven by verticalization.

Non-dimensionalization
We non-dimensionalize Eq. (26), (44)(47) by the following dimensional scales: time 1/ , pressure , length / / and velocity , , which yields with boundary conditions ̃ 0, 0 0, ̃ ̃ , ̃Θ | ̃ 0, and ̃| ̃ 0 and initial conditions Θ , , 0 Θ , and , 0 0. At the boundary, when solving for Θ, we impose a sharp velocity gradient Δ /Δ̃ ̃ / /2 where is the average length of a cell, to account for the fact that across the boundary cells there is a sharp velocity gradient which tends to align cells tangentially to the boundary 7,12 . There are two key parameters that control the evolution of the system: the dimensionless flow alignment parameter /4 and the dimensionless verticalization rate / . The flow alignment parameter dictates how quickly the cells will align to a straining flow, while the verticalization rate dictates how quickly horizontal cells are converted to vertical cells. For simplicity, we assume that the only distinguishing feature of the nonadherent mutant is that 0; however, we would also expect that in the absence of adhesion, would also be smaller and lead to an underestimation of the characteristic length-scale.

Choice of parameters
We choose the following set of parameters when solving for the evolution of the system. Growth rate : The number of bacteria in a three-dimensional biofilm was measured over time and the growth rate was found to be 0.6 hr.

Threshold verticalization pressure
: From the theoretical verticalization analysis, using experimentally calibrated parameters, it was found that the critical cell overlap distance at which point verticalization happens is ~0.125 μm. From this value, the critical cellcell contact pressure can be calculated as 4 / 3π * / where is the contact modulus of the cells which we estimate to be 500 Pa using bulk rheological measurements (shear modulus of a bulk WT* biofilm was measured to be ~200 Pa) and * /2 0.4 μm. This yields a threshold pressure of ∼ 130 Pa. Friction coefficient : Using the surface viscosity defined in Table S2, we estimate the friction coefficient as / where is the surface footprint of the cell which we estimate to be ∼ 3 μm . The friction coefficient is therefore ∼ 7 10 Pa s/μm . Characteristic length scale : The characteristic length-scale / / in the system corresponds to the characteristic size of biofilm when verticalization begins. For the parameters defined above, this length scale is ∼ 3.4 μm.

Mean cell length
: Taking the average of the largest and smallest cylindrical sections of the spherocylindrical cells (corresponding to right before and right after division) gives a length of 1.8 μm. Initial radius of biofilm , : We initialize the biofilm to have an area equivalent to the footprint of a single cell. This gives a value , ∼ 1 μm.
We solve the dimensionless Eq. (48)-(51) using a finite difference scheme with explicit time-stepping, and up-winding of advective terms. Example solutions are given in Supplementary  Figures 19 and 20. We fit the two unknown dimensionless parameters to the experimental results and find that 2.5 and 1.5. These parameters compare favorably to what has been measured before. Namely, Beroz et al. 1 measured in ABS that in a biofilm growing without a confining agarose gel, ∼ 2.5 hr. . Here we find a somewhat smaller value, ∼ 1.5 hr. , likely owing to the fact that confinement due to the agarose gel tends to impart normal stresses that suppress verticalization. Nonetheless, we find that the results are only weakly dependent on for ( Fig. 4a-d; Supplementary Fig. 19d-k). Finally, we estimate by measuring the cell-cell orientation correlation for neighboring cells in the basal layer of the experimental biofilms and find ∼ 0.5 and so ∼ 3. Measuring the intrinsic flow alignment parameter is difficult as it is dependent on many factors including the size, shape, aspect ratio and local nematic order 23 , however, a value of 3 is in reasonable agreement with values of 0.3-2 reported for other growing and non-growing bacterial systems 11,18,24 .