## Abstract

Biofilms are communities of bacteria adhered to surfaces. Recently, biofilms of rod-shaped bacteria were observed at single-cell resolution and shown to develop from a disordered, two-dimensional layer of founder cells into a three-dimensional structure with a vertically aligned core. Here, we elucidate the physical mechanism underpinning this transition using a combination of agent-based and continuum modelling. We find that verticalization proceeds through a series of localized mechanical instabilities on the cellular scale. For short cells, these instabilities are primarily triggered by cell division, whereas long cells are more likely to be peeled off the surface by nearby vertical cells, creating an ‘inverse domino effect’. The interplay between cell growth and cell verticalization gives rise to an exotic mechanical state in which the effective surface pressure becomes constant throughout the growing core of the biofilm surface layer. This dynamical isobaricity determines the expansion speed of a biofilm cluster and thereby governs how cells access the third dimension. In particular, theory predicts that a longer average cell length yields more rapidly expanding, flatter biofilms. We experimentally show that such changes in biofilm development occur by exploiting chemicals that modulate cell length.

## Main

Biofilms are groups of bacteria adhered to surfaces^{1,2,3}. These bacterial communities are common in nature, and foster the survival and growth of their constituent cells. A deep understanding of biofilm structure and development promises important health and industrial applications^{4,5}. Unfortunately, little is known about the microstructural features of biofilms due to difficulties encountered in imaging individual cells inside large assemblies of densely packed cells. Recently, however, advances in imaging technology have made it possible to observe growing, three-dimensional (3D) biofilms at single-cell resolution^{6,7,8}.

In the case of *Vibrio cholerae*, the rod-shaped bacterium responsible for the pandemic disease cholera^{9,10}, high-resolution imaging revealed a surprisingly complex biofilm developmental program^{7,8}. Over the course of 12–24 h of growth, an individual founder cell gives rise to a dome-shaped biofilm cluster, containing thousands of cells that are strongly vertically ordered, especially at the cluster core. Notably, this ordering is an intrinsically non-equilibrium phenomenon, as it is driven by growth, not by thermal fluctuations^{7,11}.

An important clue to understanding the emergence of vertical order in *V. cholerae* biofilms comes from genetic analyses that established the biological components relevant for biofilm development^{7,8,12}. To facilitate their growth as biofilms, *V. cholerae* cells secrete adhesive matrix components: *Vibrio* polysaccharide (VPS), a polymer that expands to fill gaps between cells, and cell-to-cell and cell-to-surface adhesion proteins. Cell-to-surface interactions enable vertical ordering by breaking overall rotational symmetry. However, despite previous work on the orientational dynamics of bacterial cells and related types of driven active matter^{13,14,15,16,17,18,19,20,21,22,23,24,25}, the nature of this physical process remains unclear.

Here, we establish the biophysical mechanisms controlling *V. cholerae* biofilm development. We show that the observed structural and dynamical features of growing biofilms can be reproduced by a simple, agent-based model. Our model treats individual cells as growing and dividing rods with cell-to-cell and cell-to-surface interactions, and thus serves as a minimal model for a wide range of biofilm-forming bacterial species. By examining individual cell verticalization events, we show that reorientation is driven by localized mechanical instabilities occurring in regions of surface cells subject to high in-plane compression. These threshold instabilities explain the tendency of surface-adhered cells to reorient rapidly following cell division. We incorporate these verticalization instabilities into a continuum theory, which allows us to predict the expansion speed of biofilms as well as overall biofilm morphology as a function of cell-scale properties. We verify these predictions in experiments in which we use chemicals that alter cell length. Our model thus elucidates how the mechanical and geometrical features of individual cells control the emergent features of the biofilm, which are relevant to the survival of the collective.

## Biofilm radius and vertical ordering spread linearly in time

How do cells in *V. cholerae* biofilms become vertical? Biofilms grown from a single, surface-adhered founder cell initially expand along the surface (Fig. 1a and Supplementary Movie 1). This horizontal expansion occurs because cells grow and divide along their long axes, which remain parallel to the surface due to cell-to-surface adhesion^{7}. After about 3 h, progeny near the biofilm centre begin to reorient away from the surface (Fig. 1c). Reorientation events typically involve a sharp change in a cell’s verticality *n*_{
z
}, defined as the component of the cell-orientation vector $\mathbf{n}$ normal to the surface (inset Fig. 1c). At later times, the locations of the reorientation events spread outward, and eventually the biofilm develops a roughly circular region of vertical cells surrounded by an annular region of horizontal cells. Both of these regions subsequently expand outward with approximately equal, fixed velocities. The radial profile of verticality versus time shows that the local transition of cells from horizontal to vertical occurs rapidly, taking 10–30 min for cell-sized regions to develop a vertical majority.

## Agent-based model captures verticalization

To understand how the behaviour of living biofilms arises from local interactions, we developed an agent-based model for biofilm growth (left inset Fig. 1b, Methods and Supplementary Fig. 1). The model extends existing agent-based models^{18,26,27,28} by incorporating the viscoelastic cell-to-surface adhesion^{29,30} that is crucial for *V. cholerae* biofilm formation^{9,12}. Specifically, we treat the cells as soft spherocylinders that grow, divide and adhere to the surface. We simulated biofilms by numerically integrating the equations of motion starting from a single, surface-adhered founder cell (Supplementary Movie 2 and Methods). To make the computations more tractable for systematic studies, we simulated only the surface layer of cells by removing from the simulation cells that become detached from the surface. This quasi-3D model is a reasonable approximation of the full 3D model at early times (Supplementary Fig. 2), and closely matches the dynamics and orientational order observed over the full duration of the experiment (Fig. 1b,d). The emergence of this distinctive orientational–temporal pattern demonstrates that our simple agent-based model is sufficient to capture the physical interactions that underpin the observed ‘verticalization’ transition of experimental biofilms.

## Mechanical instabilities cause cell verticalization

Why do the experimental and modelled cells segregate into horizontal and vertical orientations, with transitions from horizontal to vertical proceeding rapidly? To investigate the local mechanics that drive verticalization, we considered the dynamics of a single model cell of cylinder length $\ell $ adhered to the surface. The surface provides a combination of attractive and repulsive forces that, in the absence of external forces, maintain the cell at a stable fixed point with elevation angle *θ* = 0 (that is, horizontal) and penetration into the surface *δ*_{0}. However, when additional forces are applied to the cell, the cell may become unstable to vertical reorientation.

We determined the onset of this instability by performing a linear stability analysis for a cell under constant external forces (inset Fig. 2a). For simplicity, we took the external forces to be applied by a continuum of rigid, spherical pistons that are distributed uniformly around the cell perimeter. The pistons compress the cell in the *x*–*y* plane with an applied surface pressure *p*. For values of *p* larger than a threshold surface pressure *p*_{t}, the cell becomes linearly unstable to spontaneous reorientation (Supplementary Figs. 3 and 4). Our model yields a value of *p*_{t} that increases with $\ell $. In particular, over a broad range of $\ell $, we find a linear increase of *p*_{t} with $\ell $ (Fig. 2a). Intuitively, this increase occurs because the surface adhesion of the model cell scales with its contact area, creating an energy barrier to reorientation that increases with cell length.

To determine whether this simple model can predict verticalization events in the biofilm surface-layer simulations, we examined the forces acting on individual modelled cells throughout the development of a biofilm. Specifically, we computed the reorientation ‘surface pressure’ *p*_{r}, defined as the sum of the magnitudes of the in-plane cell–cell contact forces on a cell, normalized by the perimeter of its footprint, at the instant it begins to reorient (Supplementary Fig. 5). We found that the average reorientation surface pressure $\u27e8{p}_{\mathrm{r}}\u27e9$ increases with $\ell $, as expected from the compressive instability model (Fig. 2a). Furthermore, the predicted value *p*_{t} is in good agreement with the observed $\u27e8{p}_{\mathrm{r}}\u27e9$ for short cells. However, for long cells, $\u27e8{p}_{\mathrm{r}}\u27e9$ saturates more rapidly than *p*_{t}.

## Verticalization depends on cell length

How can longer cells become vertical at surface pressures much lower than the threshold values predicted by the compressive instability model? The large extent of the discrepancy suggests that, for long cells, the in-plane forces alone are insufficient to cause the instabilities (Fig. 2a). Indeed, incorporating the numerically observed distribution of in-plane forces acting on cells into the compressive instability model does not significantly improve the prediction for $\u27e8{p}_{\mathrm{r}}\u27e9$ (Supplementary Fig. 6). Thus, we hypothesized that, in the case of long cells, forces acting in the *z* direction might play an important role in triggering verticalization.

To explore this idea, we returned to the single-cell model and considered the effect of forces in the *z* direction (inset Fig. 2b). Applying small forces in the *z* direction to a cell with fixed centre-of-mass height shifts the equilibrium elevation angle of the cell to a small finite value of *θ*, proportional to the net torque. Under large enough torque, a single end of the cell becomes free of the surface, at which point the cell becomes unstable to further rotation and effectively ‘peels’ off the surface. This nonlinearity, inherent in the geometry of contact, competes with the compressive instability, and under specific conditions can initiate reorientation at much smaller values of surface pressure. Specifically, for a fixed centre-of-mass penetration depth ${\delta}_{0}\ll \ell $, the threshold torque for peeling a cell off the surface due to forces in the *z* direction scales as ${\tau}_{\mathrm{t}}~{\ell}^{2}$ (Supplementary Fig. 4). For the whole-biofilm surface-layer simulation, the average reorientation torque, defined as the total torque due to forces in the *z* direction at the instant of reorientation, closely obeys the predicted ${\ell}^{2}$ scaling for long cells (Fig. 2b). Taken together, our predictions from the compressive and peeling instabilities explain the verticalization of modelled cells over the entire range of cell lengths studied.

## Cell division can trigger verticalization

For both compressive and peeling instabilities, the presence of an energy barrier to reorientation explains the sharp distinction we observed between horizontal and vertical cell orientations. Furthermore, both mechanisms predict larger reorientation thresholds for longer cell lengths. Hence, the model suggests that shorter cells should reorient more readily. To confirm this effect in our simulated biofilms, we compared the distribution of reorientation lengths ${\ell}_{\mathrm{r}}$, defined as the cell cylinder length at the instant of reorientation, to the full distribution of horizontal cell lengths for a series of simulations with different values of the initial cell cylinder length ${\ell}_{0}$, defined as the cell cylinder length immediately following division (Fig. 2c). For all values of ${\ell}_{0}$, we found that the mean reorientation cell length $\u27e8{\ell}_{\mathrm{r}}\u27e9$ is substantially smaller than the mean horizontal cell length. In addition, for simulations with average cell lengths comparable to those in our experiments, most reorientation events occur immediately after cell division. The limited time resolution of the experiments precludes a quantitative analysis of division-induced verticalization; nevertheless, the propensity for cell division to trigger reorientation is clearly observed in the experimental biofilm (Supplementary Movie 1 and Supplementary Fig. 7).

## Verticalization is localized

We next investigated how the surface compression and peeling instabilities influence the propagation of reorientation through a biofilm. First, we generalized our model for the surface compression instability to the multicellular level. A linear stability analysis of the model suggests that reorientation events should be independent and spatially localized for short cell lengths (Supplementary Fig. 3). By contrast, for long cell lengths, the tendency of neighbouring vertical cells to trigger reorientation suggests an (inverse) domino-like effect in which one cell standing up can induce its neighbours to stand up.

To quantify the extent of cooperative verticalization, we computed the size *N* of reorientation ‘avalanches’, defined as groups of verticalization events that are proximal in space and time^{31} (inset Fig. 2d and Supplementary Fig. 8). We found that for long cell lengths, the mean avalanche size $\u27e8N\u27e9$ is significantly larger than the value predicted by a null model that consists of randomizing the angular positions of cells within the biofilm (Fig. 2d and Supplementary Fig. 8). Thus, the dependence of $\u27e8N\u27e9$ on cell length is consistent with the prediction of the inverse domino effect. Interestingly, however, the distribution of avalanche sizes decays roughly exponentially for all values of cell length we studied (Supplementary Fig. 9), with only a modest number of cells (*N* ~ 2–3) involved in typical avalanches (Fig. 2d).

A natural explanation for the small mean avalanche size comes from the reduction in cell footprint that occurs on reorientation, which rapidly alleviates the local surface pressure responsible for verticalization (Supplementary Fig. 5). This effect combines with the disorder of the contact geometries and forces throughout the biofilm, which separates horizontal cells near the verticalization threshold into disconnected groups (Supplementary Movie 3 and Supplementary Fig. 10). Thus, although the inverse domino effect transiently increases verticalization cooperativity, avalanches quickly exhaust the local supply of horizontal cells that are susceptible to becoming vertical. Consequently, verticalization occurs throughout the biofilm in scattered, localized regions.

## Two-fluid model describes propagation of verticalization

To understand how localized cell verticalization gives rise to the global patterning dynamics of the biofilm, we developed a 2D continuum model that treats horizontal and vertical cell densities as two coupled fluids (Fig. 3a and Methods). The local horizontal cell density *ρ*_{h} grows in the plane at a rate *α* and converts to vertical cell density *ρ*_{v} in regions of high surface pressure at a rate *β*. The total 2D cell density is ${\stackrel{\u0303}{\rho}}_{\mathrm{tot}}$ ≡ *ρ*_{h} + *ξρ*_{
v
}, where *ξ* is the ratio of vertical to horizontal cell footprints, and we approximate the surface pressure as linearly proportional to the areal deformation ${\stackrel{\u0303}{\rho}}_{\mathrm{tot}}-{\stackrel{\u0303}{\rho}}_{0}$, where ${\stackrel{\u0303}{\rho}}_{0}$ is the close-packed, but uncompressed, cell density. These interactions yield the following equation for the change in ${\stackrel{\u0303}{\rho}}_{\mathrm{tot}}$ in regions of non-zero surface pressure:

where *γ* is the ratio of the Young’s modulus of the biofilm to the surface drag coefficient, Θ is the Heaviside step function and ${\stackrel{\u0303}{\rho}}_{\mathrm{t}}$ is the threshold surface density for verticalization. We simulated this continuum model, and found that for *α* < *β*, the biofilm generically develops into a circular region containing both horizontal and vertical cells (‘mixed interior’) surrounded by an annular region containing horizontal cells (‘horizontal cell periphery’), closely matching both the experimental biofilm and the agent-based model biofilms (Fig. 3b, Supplementary Movie 4 and Supplementary Note). In this regime, the biofilm front spreads linearly in time at a fixed expansion speed *c*^{*}. Furthermore, the total cell density and the surface pressure are constant in the mixed interior. This constancy is stabilized by the competing effects of cell growth and cell verticalization, and occurs provided that *α* < *β*(1 − *ξ*) (Supplementary Fig. 11). When this condition is satisfied, verticalization can reduce cell density faster than cell density can be replenished by cell growth and cell transport due to gradients in surface pressure. Physically, this results in the cell density rapidly fluctuating around the verticalization threshold. This rapid alternation effectively tunes the verticalization rate down to *α* = *β*(1 − *ξ*), and thereby ensures a constant total cell density and surface pressure in the mixed interior. The resulting ‘dynamical isobaricity’ (constancy of pressure) provides the boundary condition for the horizontal cell periphery that determines the horizontal expansion speed *c*^{*}, independent of *β* and *ξ*. In the limit of slow expansion ${c}^{*}\ll \sqrt{\alpha \gamma}$, *c*^{*} is given by:

where ${c}_{0}^{*}=\sqrt{2\alpha \gamma}$. This expansion speed increases with ${\stackrel{\u0303}{\rho}}_{\mathrm{t}}$ until it saturates at ${c}^{*}=2\sqrt{\alpha \gamma}$ for ${\stackrel{\u0303}{\rho}}_{\mathrm{t}}\gg {\stackrel{\u0303}{\rho}}_{0}$ (Supplementary Note). Intuitively, higher values of the verticalization threshold density ${\stackrel{\u0303}{\rho}}_{\mathrm{t}}$ sustain a wider periphery of horizontal cells, which results in a higher rate of increase in the total number of surface cells. Thus, our continuum model reveals how the geometrical and mechanical properties of individual cells influence the global morphology of the growing biofilm.

## Longer cells yield more rapidly expanding, flatter biofilms

As the threshold surface density for verticalization ${\stackrel{\u0303}{\rho}}_{\mathrm{t}}$ increases with cell length, we expect biofilms composed of longer cells to maintain a wider periphery of horizontal cells and to thereby expand faster along the surface than biofilms composed of shorter cells. To test this notion in our agent-based model, we computed the expansion speed of the modelled biofilms for a range of initial cell cylinder lengths ${\ell}_{0}$ (Fig. 4b). On fitting the continuum model parameters to those of the agent-based model (Supplementary Figs. 12 and 13), we found that the expansion velocities of the two models were equal to within a few percent. In living, experimental biofilms, we can increase or decrease the average cell length using chemical treatments^{10,32} (Fig. 4a, top row). Similar to the agent-based model biofilms, in experimental biofilms, the surface expansion speed increases with increasing cell length (Fig. 4b and Supplementary Movie 5). Furthermore, when the model biofilm parameters are fitted to experiment (Supplementary Fig. 1), the experimental and model biofilm speeds agree to within 20% or better. Taken together, these observations support the conclusion that self-organized dynamical isobaricity governs the expansion of *V. cholerae* biofilms along surfaces.

How do different surface expansion speeds influence the ensuing biofilm development into the *z* direction? After a few hours, living biofilms grow into roughly semi-ellipsoidal shapes with volume $V=\left(2\mathrm{\u2215}3\right){R}_{\mathrm{B}}^{2}H$, where *R*_{B} is the basal radius and *H* is the height. For equal rates of total volume growth, we expect a biofilm that expands more rapidly along the surface to develop a lower aspect ratio *H*/*R*_{B} than a biofilm that expands less rapidly along the surface. We verified that the chemical treatment does not affect the total volume growth rates of the experimental biofilms (Supplementary Fig. 14), and found that the measured aspect ratio *H*/*R*_{B} does indeed decrease with cell length over a wide range of volumes (Fig. 4a, bottom row; and Fig. 4c). Thus, our results show how the elongated geometries of individual cells govern the global morphology of the collective, an important determinant of biofilm survival in nature.

## Discussion

Our results suggest that bacteria have harnessed the physics of mechanical instabilities to enable the generation of spatially and orientationally patterned architectures. Here, we showed that cell verticalization begins to occur when the local effective surface pressures that arise from cell growth become large enough to overcome the cell-to-surface adhesion that otherwise favours a horizontal orientation. Subsequently, the reduction in cell footprint that occurs on cell verticalization, which acts to reduce the effective surface pressure, provides a mechanical feedback that controls the rate of biofilm expansion. We expect that individual cell parameters have evolved in response to selective pressures on global biofilm morphology (for example, during resource competition)^{6,33,34,35,36}. Since optimal morphology may be condition-dependent, cells may also have evolved adaptive strategies that alter biofilm architecture, which could be investigated experimentally by screening for environmental influences on cell size, shape and surface adhesion^{37}.

For simplicity, we focused on flat surfaces, nutrient-rich conditions and *V. cholerae* strains that have been engineered to have simpler interactions than those in wild-type biofilms (Methods). Moreover, our agent-based model does not explicitly incorporate the VPS matrix secreted by cells^{2,27,38}. Understanding the modifying effects of the VPS matrix, cell and surface curvature (Supplementary Fig. 15), cell-to-cell adhesion (Supplementary Fig. 16) and chemical feedback^{39} will be important directions for future studies. More broadly, we must develop a systematic method to account for the diversity of architectures that can be produced by local mechanical interactions (Supplementary Discussion).

Our study of a two-fluid model for verticalizing biofilms led us to discover a novel type of front propagation. Interestingly, in the biofilm surface layer, the front profile of cell density is precisely uniform starting at some finite distance from the edge, whereas previous models of front propagation saturate asymptotically toward uniformity^{40,41,42,43}. The self-organized nature of this process yields a universal dependence of the expansion speed on the cell geometrical and mechanical parameters that is robust to details of the mechanical feedback. We have focused on the mean-field behaviour of biofilms, but an open question is to understand the role of fluctuations in the ‘pressure’ acting on cells (for example, from a jamming perspective^{44}, a fluctuating hydrodynamical perspective^{45,46} or a combination of approaches).

In summary, we have elucidated the physical mechanism underlying a complex developmental program observed at the cellular scale in bacterial biofilms. The relative biochemical and biophysical simplicity of this prokaryotic system allowed us to quantitatively understand the developmental pathway from the scale of a single cell to the scale of a large community assembly. Going forward, we expect bacterial biofilms will take on increasingly important roles as tractable models that can be used to understand how living systems generate and maintain their structures.

## Methods

### Growing and imaging experimental biofilms

#### Strains and media

The *V. cholerae* strain used in this study is a derivative of wild-type *Vibrio cholerae* O1 biovar El Tor strain C6706, harbouring a missense mutation in the *vpvC* gene (VpvC W240R) that elevates c-di-GMP levels and confers a rugose biofilm phenotype^{47}. Additional mutations were engineered into this strain using *Escherichia coli* S17-*λpir* carrying pKAS32. Specifically, the biofilm gene responsible for cell–cell adhesion, *rbmA*, was deleted. To avoid the effects of cell curvature, we deleted *crvA* encoding the periplasmic protein CrvA responsible for the curvature of *V. cholerae* cells. Biofilm experiments were performed in M9 minimal medium, supplemented with 2 mM MgSO_{4}, 100 μM CaCl_{2} and 0.5% glucose. When indicated, cefalexin (Sigma Aldrich) was added at 4 μg ml^{−1} and A22 (a gift from the Gitai group) was used at 0.4 μg ml^{−1}. These concentrations were experimentally determined to modulate cell morphology without affecting overall mass accumulation.

#### Biofilm growth

*V. cholerae* strains were grown overnight at 37 °C in liquid LB medium with shaking, back-diluted 30-fold and grown for an additional 2 h with shaking in M9 medium until early exponential phase (OD_{600 nm} = 0.1–0.2). These re-grown cultures were diluted to OD_{600 nm} = 0.001 and 100 μl of the diluted cultures was added to wells of 96-well plates with No. 1.5 coverslip bottoms (MatTek). The cells were allowed to attach for 10 min, after which the wells were washed twice with fresh M9 medium, and, subsequently, 100 μl of fresh M9 medium was added, with or without drugs. The low initial inoculation density enabled isolated biofilm clusters to form. The locations of the founder cells were identified, and 1 h after inoculation, imaging was begun on the microscope stage at 25 °C.

#### Microscopy

Details of the imaging system have been described elsewhere^{8}. Briefly, images were acquired with a spinning-disc confocal microscope (Yokogawa), a 543 nm laser (OEM DPSS) and an Andor iXon 897 EMCCD camera. For single-cell-resolution imaging, a 60× water objective with a numerical aperture of 1.2 plus a 1.5x post-magnification lens was used. To avoid evaporation, immersion oil with a refractive index of 1.3300 ± 0.0002 (Cargille) was used instead of water. The time difference between each image acquisition was 10 min, and the total imaging time was 8 h. Only the bottom 5 μm of the biofilm was imaged (with a *z* step size of 0.2 μm) to avoid excessive photobleaching and phototoxicity. For coarse-grained imaging, a 20x multi-immersion objective was used without post-magnification. In this case, the time difference between each image acquisition was 30 min, and entire biofilms were imaged (with a *z* step size of 1 μm). At the end of the coarse-grained time course, the biofilm clusters were imaged again with high magnification to determine cell lengths. All image acquisitions were automated using Nikon Element software. All cells harboured mKO fluorescent proteins expressed from the chromosome. Experimental images in Fig. 4 were false-coloured to differentiate between different growth conditions.

#### Image processing

The cell segmentation protocol and Matlab codes have been described elsewhere in detail^{8}. Cell position, cell length and cell orientation were used as input for 3D rendering in Fig. 1a and for further analysis. For biofilm clusters grown in the presence of A22 or cefalexin, cell lengths were manually measured in the bottom cell layers of the biofilms using the Nikon Element software. To define the biofilm shape parameters in the coarse-grained images, the bottom cell layers of the biofilms were first identified by finding the brightest *z*-cross-section, according to the total fluorescence intensity. The contour of the individual biofilm cluster was next identified using the 3D Canny edge detection method implemented in Mathematica. To correct for the inevitable optical stretching in the *z*-direction, we compared the heights obtained from the same cluster in the coarse-grained and the fine resolution images. By imaging a series of biofilm clusters of different sizes, we obtained a curve of a biofilm cluster’s actual height versus its apparent height in the coarse-grained images, which we used to calibrate the heights measured in the coarse-grained images.

### Modelling biofilms

#### Agent-based model

We model the volume occupied by a cell as a cylinder of length $\ell $ with two hemispherical end caps each of radius *R*. We treat cell growth as elongation that increases cell volume at a fixed rate *α* (chosen randomly at birth from a narrow Gaussian distribution to desynchronize cell divisions). Cells are born with an initial cell cylinder length ${\ell}_{0}$ and grow to twice their total length. When a cell reaches the division length, the cell is instantaneously replaced by two identical daughter cells. In our model, cell-to-cell and cell-to-surface overlaps exert repulsive forces according to Hertzian contact mechanics for elastic materials^{48}. In the case of the cell-to-surface overlap, we also include an attractive interaction with an energy proportional to the cell-to-surface contact area (that is, the Derjaguin approximation^{29,30}). To account for symmetry-breaking microscopic irregularities, we add a small amount of random noise to each component of the generalized force at every time step (10^{−8} *E*_{0}*R*^{2} to the force acting on the centre-of-mass **r** and 10^{−8} *E*_{0}*R*^{3} to the generalized force acting on the orientation vector $\mathbf{n}$). Finally, we include two sources of viscous drag: a modest damping of 3D motion through the surrounding fluid and biopolymer matrix, and a much larger damping of sliding motion along the adhesive surface^{29}. We determine the parameters in our model by fitting to the experimental data (Supplementary Fig. 1).

#### Simulation of agent-based model

To simulate the agent-based model, we numerically integrate the equations of motion using an explicit embedded Runge–Kutta–Fehlberg method. Our implementation of this method in C++ is adapted from the GNU Scientific Library^{49}, and available freely online at GitHub (https://github.com/farzanb/Agent-based-biofilm-model). The initial conditions consist of a single cell of length ${\ell}_{0}$ at elevation angle *θ* = 0 and penetration depth *δ*_{0}.

#### Continuum model

To describe the radial expansion and orientational dynamics of the biofilm, we treated the horizontal cells and vertical cells as continuous fields *ρ*_{h} and *ρ*_{v}, with the total density ${\stackrel{\u0303}{\rho}}_{\mathrm{tot}}$ defined as the sum *ρ*_{h} + *ξρ*_{v}, where *ξ* is the ratio of vertical to horizontal cell footprints. In regions where ${\stackrel{\u0303}{\rho}}_{\mathrm{tot}}$ exceeds the close-packing density ${\stackrel{\u0303}{\rho}}_{0}$, the rate of change of the horizontal cell density is proportional to the divergence of the flux of cells due to cell transport plus terms due to cell growth and cell verticalization. We assume that cell transport is given by the gradient of surface pressure divided by a surface viscosity coefficient *η*_{1}, where surface pressure is approximated as linearly proportional to the areal deformation ${\stackrel{\u0303}{\rho}}_{\mathrm{tot}}-{\stackrel{\u0303}{\rho}}_{0}$ (that is, following ‘Hooke’s law’). Cell growth in the plane is proportional to the local horizontal cell density and occurs at a fixed rate *α*. We assume that cell verticalization locally converts horizontal cell density to vertical cell density at a fixed rate *β* in regions where ${\stackrel{\u0303}{\rho}}_{\mathrm{tot}}>{\stackrel{\u0303}{\rho}}_{\mathrm{t}}$, where ${\stackrel{\u0303}{\rho}}_{\mathrm{t}}$ is the threshold surface density for verticalization.

#### Simulation of continuum model

We simulated the dynamics of the radial cell densities in Python using FiPy, a finite volume PDE solver^{50}. We performed our simulation on a cylindrical grid containing 30,000 radial points with a spacing of Δ*x* = 1 nm and a temporal step size of Δ*t* = 1 ms.

### Reporting Summary

Further information on experimental design is available in the Nature Research Reporting Summary linked to this article.

### Code availability

The cell segmentation protocol and Matlab codes for image processing have been described in a previous publication^{8}. The simulation used for the agent-based model is available on GitHub (https://github.com/farzanb/Agent-based-biofilm-model). The code used for the analysis of the agent-based model simulations is available from the corresponding author upon request.

### Data availability

The data are available from the corresponding author on request.

## Additional information

**Publisher’s note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## References

- 1.
Ghannoum, M., Parsek, M., Whiteley, M. & Mukherjee, P. K.

*Microbial Biofilms*2nd edn (ASM, Washington DC, 2015). - 2.
Flemming, H. C. & Wingender, J. The biofilm matrix.

*Nat. Rev. Microbiol.***8**, 623–633 (2010). - 3.
Gordon, V. D., Davis-Fields, M., Kovach, K. & Rodesney, A. Biofilms and mechanics: a review of experimental techniques and findings.

*J. Phys. D***50**, 223002 (2017). - 4.
Hall-Stoodley, L., Costerton, J. W. & Stoodley, P. Bacterial biofilms: from the natural environment to infectious diseases.

*Nat. Rev. Microbiol.***2**, 95–108 (2004). - 5.
Costerton, J. W., Stewart, P. S. & Greenberg, E. P. Bacterial biofilms: a common cause of persistent infections.

*Science***284**, 1318–1322 (1999). - 6.
Stewart, E. J., Satorius, A. E., Younger, J. G. & Solomon, M. J. Role of environmental and antibiotic stress on

*Staphylococcus epidermidis*biofilm microstructure.*Langmuir***29**, 7017–7024 (2013). - 7.
Drescher, K. et al. Architectural transitions in

*Vibrio cholerae*biofilms at single-cell resolution.*Proc. Natl Acad. Sci. USA***113**, E2066–E2072 (2016). - 8.
Yan, J., Sharo, A. G., Stone, H. A., Wingreen, N. S. & Bassler, B. L.

*Vibrio cholerae*biofilm growth program and architecture revealed by single-cell live imaging.*Proc. Natl Acad. Sci. USA***113**, E5337–E5343 (2016). - 9.
Teschler, J. K. et al. Living in the matrix: assembly and control of

*Vibrio cholerae*biofilms.*Nat. Rev. Microbiol.***13**, 255–268 (2015). - 10.
Bartlett, T. M. et al. A periplasmic polymer curves

*Vibrio cholerae*and promotes pathogenesis.*Cell***168**, 172–185 (2017). - 11.
Onsager, L. The effects of shape on the interaction of colloidal particles.

*Ann. NY Acad. Sci.***51**, 627–659 (1949). - 12.
Berk, V. et al. Molecular architecture and assembly principles of

*Vibrio cholerae*biofilms.*Science***337**, 236–239 (2012). - 13.
Doostmohammadi, A., Thampi, S. P. & Yeomans, J. M. Defect-mediated morphologies in growing cell colonies.

*Phys. Rev. Lett.***117**, 048102 (2016). - 14.
Duvernoy, M.-C. et al. Asymmetric adhesion of rod-shaped bacteria controls microcolony morphogenesis.

*Nat. Commun.***9**, 1120 (2018). - 15.
Enos-Berlage, J. L. & McCarter, L. L. Relation of capsular polysaccharide production and colonial cell organization to colony morphology in

*Vibrio parahaemolyticus*.*J. Bacteriol.***182**, 5513–5520 (2000). - 16.
Boyer, D. et al. Buckling instability in ordered bacterial colonies.

*Phys. Biol.***8**, 026008 (2011). - 17.
Smith, W. P. J. et al. Cell morphology drives spatial patterning in microbial communities.

*Proc. Natl Acad. Sci. USA***114**, E280–E286 (2016). - 18.
You Z., Pearce, D. J. G., Sengupta, A. & Giomi, L. Geometry and mechanics of micro-domains in growing bacterial colonies. Preprint at https://arxiv.org/abs/1703.04504 (2017).

- 19.
Mamou, G., Mohan, G. B. M., Rouvinski, A., Rosenberg, A. & Ben-Yehuda, S. Early developmental program shapes colony morphology in bacteria.

*Cell Rep.***14**, 1850–1857 (2016). - 20.
Grant, M. A. A., Waclaw, B., Allen, R. J. & Cicuta, P. The role of mechanical forces in the planar-to-bulk transition in growing

*Escherichia coli*microcolonies.*J. R. Soc. Interface*https://doi.org/10.1098/rsif.2014.0400 (2014). - 21.
Su, P.-T. et al. Bacterial colony from two-dimensional division to three-dimensional development.

*PLoS One***7**, e48098 (2012). - 22.
Asally, M. et al. Localized cell death focuses mechanical forces during 3D patterning in a biofilm.

*Proc. Natl Acad. Sci. USA***109**, 18891–18896 (2012). - 23.
Kumar, N., Soni, H., Ramaswamy, S. & Sood, A. K. Flocking at a distance in active granular matter.

*Nat. Commun.***5**, 4688 (2014). - 24.
Narayan, V., Ramaswamy, S. & Menon, N. Long-lived giant number fluctuations in a swarming granular nematic.

*Science***317**, 105–108 (2007). - 25.
Yadav, V., Chastaing, J.-Y. & Kudrolli, A. Effect of aspect ratio on the development of order in vibrated granular rods.

*Phys. Rev. E***88**, 052203 (2013). - 26.
Farrell, F. D. C., Hallatschek, O., Marenduzzo, D. & Waclaw, B. Mechanically driven growth of quasi-two-dimensional microbial colonies.

*Phys. Rev. Lett.***111**, 168101 (2013). - 27.
Ghosh, P., Mondal, J., Ben-Jacob, E. & Levine, H. Mechanically-driven phase separation in a growing bacterial colony.

*Proc. Natl Acad. Sci. USA***112**, E2166–E2173 (2015). - 28.
Lardon, L. A. et al. iDynoMiCS: next-generation individual-based modelling of biofilms.

*Env. Microbiol.***13**, 2416–2434 (2011). - 29.
Chen, Y., van der Mei, H. C., Busscher, H. J. & Norde, W. Viscous nature of the bond between adhering bacteria and substratum surfaces probed by atomic force microscopy.

*Langmuir***30**, 3165–3169 (2014). - 30.
Barthel, E. Adhesive elastic contacts - JKR and more.

*J. Phys. D***41**, 163001 (2008). - 31.
Candelier, R., Dauchot, O. & Biroli, G. Dynamical facilitation decreases when approaching the granular glass transition.

*Eur. Phys. Lett.***92**, 24003 (2010). - 32.
White, C. L., Kitich, A. & Gober, J. W. Positioning cell wall synthetic complexes by the bacterial morphogenetic proteins MreB and MreD.

*Mol. Microbiol.***76**, 616–633 (2010). - 33.
Yang, D. C., Blair, K. M. & Salama, N. R. Staying in shape: the impact of cell shape on bacterial survival in diverse environments.

*Microbiol. Mol. Biol. Rev.***80**, 187–203 (2016). - 34.
Young, K. D. The selective value of bacterial shape.

*Microbiol. Mol. Biol. Rev.***70**, 660–703 (2006). - 35.
Stocker, R. & Seymour, J. R. Ecology and physics of bacterial chemotaxis in the ocean.

*Microbiol. Mol. Biol. Rev.***76**, 792–812 (2012). - 36.
Chang, F. & Huang, K. C. How and why cells grow as rods.

*BMC Biol.***12**, 54 (2014). - 37.
Jubair, M., Morris, J. G.Jr. & Ali, A. Survival of

*Vibrio cholerae*in nutrient-poor environments is associated with a novel “persister” phenotype.*PLoS One***7**, e45187 (2012). - 38.
Yan, J., Nadell, C. D., Stone, H. A., Wingreen, N. S. & Bassler, B. L. Extracellular-matrix-mediated osmotic pressure drives

*Vibrio cholerae*biofilm expansion and cheater exclusion.*Nat. Commun.***8**, 327 (2017). - 39.
Vicsek, T.

*Fluctuations and Scaling in Biology*(Oxford Univ. Press, Oxford, 2001). - 40.
Fisher, R. A. The wave of advance of advantageous genes.

*Ann. Eugen.***7**, 353–369 (1937). - 41.
von Saarloos, W. Front propagation into unstable states.

*Phys. Rep.***386**, 29–222 (2003). - 42.
Fort, J. & Pujol, T. Progress in front propagation research.

*Rep. Prog. Phys.***71**, 086001 (2008). - 43.
Puliafito, A. et al. Collective and single cell behavior in epithelial contact inhibition.

*Proc. Natl Acad. Sci. USA***109**, 739–744 (2012). - 44.
Delarue, M. et al. Self-driven jamming in growing microbial populations.

*Nat. Phys.***12**, 762–766 (2016). - 45.
Bi, D., Lopez, J. H., Schwarz, J. M. & Manning, M. L. A density-independent rigidity transition in biological tissues.

*Nat. Phys.***11**, 1074–1078 (2015). - 46.
Marchetti, M. C. et al. Hydrodynamics of soft active matter.

*Rev. Mod. Phys.***85**, 1143 (2013). - 47.
Beyhan, S. & Yildiz, F. H. Smooth to rugose phase variation in

*Vibrio cholerae*can be mediated by a single nucleotide change that targets c-di-GMP signaling pathway.*Mol. Microbiol.***63**, 995–1007 (2007). - 48.
Landau, L. D., Lifshitz, E. M., Kosevich, A. M. & Pitaevskii, L. P.

*Theory of Elasticity*(Pergamon Press, Oxford, 1986). - 49.
Galassi, M. et al.

*GNU Scientific Library Reference Manual*3rd edn (Network Theory Ltd, Godalming, 2017). - 50.
Guyer, J. E., Wheeler, D. & Warren, J. A. FiPy: Partial differential equations with Python.

*Comp. Sci. Eng.***11**, 6–15 (2009).

## Acknowledgements

We thank B. Bratton, X. Chen, S. Jena, N. Pappireddi and J. Shaevitz for insightful discussions and encouragement, G. Laevsky for assistance with imaging and T. Bartlett and Z. Gitai for sharing the Δ*crvA* strain and reagents. This work was supported by the NIH under grant no. R01 GM082938 (F.B., Y.M. and N.S.W.), the Howard Hughes Medical Institute (B.L.B.), NSF grant MCB-1713731 (B.L.B.), NSF grant MCB-1344191 (H.A.S., B.L.B. and N.S.W.), NIH grant 2R37GM065859 (F.B. and B.L.B.), the Max Planck Society-Alexander von Humboldt Foundation (B.L.B.) and the Eric and Wendy Schmidt Transformative Technology Fund from Princeton University (H.A.S., B.L.B. and N.S.W.). J.Y. holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund.

## Author information

### Affiliations

#### Joseph Henry Laboratories of Physics, Princeton University, Princeton, NJ, USA

- Farzan Beroz

#### Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ, USA

- Jing Yan
- , Benedikt Sabass
- & Howard A. Stone

#### Department of Molecular Biology, Princeton University, Princeton, NJ, USA

- Jing Yan
- , Bonnie L. Bassler
- & Ned S. Wingreen

#### Department of Physics, Ben-Gurion University of the Negev, Beer-Sheva, Israel

- Yigal Meir

#### ICS-2/IAS-2, Forschungszentrum Jülich, Jülich, Germany

- Benedikt Sabass

#### Howard Hughes Medical Institute, Chevy Chase, MD, USA

- Bonnie L. Bassler

#### Lewis-Sigler Institute for Integrative Genomics, Princeton University, Princeton, NJ, USA

- Ned S. Wingreen

### Authors

### Search for Farzan Beroz in:

### Search for Jing Yan in:

### Search for Yigal Meir in:

### Search for Benedikt Sabass in:

### Search for Howard A. Stone in:

### Search for Bonnie L. Bassler in:

### Search for Ned S. Wingreen in:

### Contributions

F.B., J.Y., and N.S.W. conceived the research. F.B. performed all simulations and analysis. J.Y. carried out all experiments, and J.Y. and F.B. contributed to the experimental design. All authors contributed to the manuscript preparation.

### Competing interests

The authors declare no competing interests.

### Corresponding author

Correspondence to Ned S. Wingreen.

## Supplementary information

### Supplementary Information

Supplementary Information, Supplementary Figures 1–16, Supplementary Discussion, Supplementary References 1–20

### Reporting Summary

### Supplementary Video 1

Growth of a

*V. cholerae*biofilm cluster, showing cross-sectional images of the bottom cell layer. The strain constitutively expresses*mKO*. The viewing window is 45 by 45 µm^{2}and the total duration is 8 hours with 10 min time steps.### Supplementary Video 2

Visualization of the surface layer of a modelled biofilm with initial cell cylinder length ${\ell}_{0}$ = 1 µm, showing positions and orientations of horizontal (blue) and vertical (red) surface-adhered cells as spherocylinders of radius

*R*= 0.8 µm, with the surface shown at height*z*= 0 µm (brown). Cells with*n*_{ z }< 0.5 (>0.5) are considered horizontal (vertical), where $\mathbf{n}$ is the orientation vector. Scale bar, 3 µm. The total duration is 10 hours.### Supplementary Video 3

Visualization of the surface layer of a modelled biofilm with initial cell cylinder length ${\ell}_{0}$ = 2 µm, showing horizontal (blue) and vertical (red) cells as spherocylinders, the surface (brown), and cell-to-cell contact forces (yellow lines connecting the centres of cells, with thicknesses proportional to the force). Cells with

*n*_{ z }< 0.5 (>0.5) are considered horizontal (vertical), where $\mathbf{n}$ is the orientation vector. The length of the scale bar is 3 µm, and its thickness corresponds to 300 pN.### Supplementary Video 4

Numerical simulation of the continuum model assuming no vertical cell transport, showing radial densities of horizontal cells (

*ρ*_{h}, blue), vertical cells (*ρ*_{v}, red), and total density (${\stackrel{\u0303}{\rho}}_{\mathrm{tot}}$=*ρ*_{h}+*ξρ*_{v}, black), versus radial coordinate*r*for isobaric regime with ${\stackrel{\u0303}{\rho}}_{0}$ = 1 m^{−2}, ${\stackrel{\u0303}{\rho}}_{\mathrm{t}}$ = 1.5 m^{−2},*β*= 2.5*α*, and*ξ*= 0.5. Dashed grey line shows ${\stackrel{\u0303}{\rho}}_{\mathrm{t}}.$### Supplementary Video 5

Expansion of

*V. cholerae*biofilm clusters grown with the drug A22 at a concentration of 0.4 µg mL^{–1}(left) and the drug Cefalexin at a concentration of 4 µg mL^{–1}(right). Cross-sectional images of the bottom cell layers are shown. The strain constitutively expresses*mKO*. Scale bar, 30 µm. The total duration is 10 hours with 30 min time steps.

## About this article

### Publication history

#### Received

#### Accepted

#### Published

### DOI

https://doi.org/10.1038/s41567-018-0170-4

### Rights and permissions

To obtain permission to re-use content from this article visit RightsLink.