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.

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.

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⁷. 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.

a,b, Top-down and perspective visualizations of the surface layer of experimental (a) and modelled (b) biofilms, showing the 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. The upper-left panel of a shows a confocal fluorescence microscopy image, and the upper-right panel shows the corresponding reconstructed central cluster using the positions and orientations of surface cells. The upper-left panel of b shows a schematic representation of modelled cell–cell (orange) and cell–surface (yellow) interactions, which depend, respectively, on the cell–cell overlap δ _ij (purple) and cell–surface overlap δ _i (red) (see Methods for details). Scale bars, 5 μm. c,d, 2D growth of biofilm surface layer for experimental biofilm (c; same as shown in a) and modelled biofilms (d). The colour of each spatiotemporal bin indicates the fraction of vertical cells at a given radius r from the biofilm centre, averaged over the angular coordinates of the biofilm (grey regions contain no cells). In d, each spatiotemporal bin is averaged over ten simulated biofilms. In c,d, the horizontal dashed pink lines show the onset of verticalization. The black dashed lines show the edge of the biofilm. Insets show the distribution of cell orientations at time t = 300 min, with colour highlighting horizontal and vertical orientations.

Full size image

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.

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 δ₀. 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.

a,b, Properties of individual cells at the time t_r of reorientation, defined as the time of the peak of total force on the cell prior to it becoming vertical. Analyses are shown for all reorientation events among different biofilms simulated for a range of initial cell lengths {\ell }_0. a, Distributions of reorientation ‘surface pressure’ p_r, defined as the total contact force in the x–y plane acting on a cell at time t_r, normalized by the cell’s perimeter, versus cell cylinder length \ell. The white dashed curve shows the average reorientation surface pressure ⟨p_{\mathrm{r}}⟩ as a function of \ell. The magenta dashed curve shows the threshold surface pressure p_t from linear stability analysis for a modelled cell under uniform pressure, depicted schematically in the inset. b, Distributions of the logarithm of reorientation torque τ_r, defined as the magnitude of the torque on a cell due to cell–cell contact forces in the z direction at time t_r, for different cell cylinder lengths \ell. The white dashed curve shows the average values ⟨{\tau }_{\mathrm{r}}⟩ as a function of \ell. The orange dashed curve shows the scaling {\tau }_{\mathrm{t}}~{\ell }^2 of the threshold torque for peeling from linear stability analysis for a modelled cell, depicted schematically in the inset. c, Mean reorientation length ⟨{\ell }_{\mathrm{r}}⟩ (red), defined as the average value of cell length at t_r, and mean cell cylinder length ⟨\ell ⟩ (grey), defined as the average length of all horizontal cells over all times of biofilm growth, averaged over ten simulated biofilms, each with initial cell cylinder length {\ell }_0, plotted versus {\ell }_0. The inset shows the distribution of reorientation lengths (red) and horizontal surface–cell lengths (grey) for {\ell }_0 = 1 μm. d, Mean avalanche size ⟨N⟩, defined as the average size of a cluster of reorienting cells that are proximal in space and time (Supplementary Figs. 8–10), versus the initial cell length {\ell }_0 for the experimental biofilm (red triangle) and the modelled biofilm (red circles). The open grey triangle and circles indicate the corresponding mean avalanche sizes for a null model. The inset shows a side view of cell configurations in the x–y plane at times t_r for all reorientation events in a simulated biofilm with {\ell }_0 = 2.5 μm. Reorientation events are coloured alike if they belong to the same avalanche. Scale bars, 10 μm and 1 h.

Full size image

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 ⟨p_{\mathrm{r}}⟩ 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 ⟨p_{\mathrm{r}}⟩ for short cells. However, for long cells, ⟨p_{\mathrm{r}}⟩ saturates more rapidly than p_t.

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 ⟨p_{\mathrm{r}}⟩ (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.

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 ⟨{\ell }_{\mathrm{r}}⟩ 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).

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³¹ (inset Fig. 2d and Supplementary Fig. 8). We found that for long cell lengths, the mean avalanche size ⟨N⟩ 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 ⟨N⟩ 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.

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 {\tilde{\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 {\tilde{\rho }}_{\mathrm{tot}}-{\tilde{\rho }}_0, where {\tilde{\rho }}_0 is the close-packed, but uncompressed, cell density. These interactions yield the following equation for the change in {\tilde{\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 {\tilde{\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 {\tilde{\rho }}_{\mathrm{t}} until it saturates at c^{*}=2\sqrt{\alpha \gamma } for {\tilde{\rho }}_{\mathrm{t}}\gg {\tilde{\rho }}_0 (Supplementary Note). Intuitively, higher values of the verticalization threshold density {\tilde{\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.

a, Schematic illustration of the two-component continuum model. Horizontal cells (blue) and vertical cells (red) are modelled, respectively, by densities ρ_h and ρ_v in two spatial dimensions. The total cell density {\tilde{\rho }}_{\mathrm{tot}} is defined as ρ_h + ξρ_v, where ξ is the ratio of vertical to horizontal cell footprints. b, Radial densities ρ of vertical cells (ρ_v, red), horizontal cells (ρ_h, blue) and total density ({\tilde{\rho }}_{\mathrm{tot}}, black), versus shifted radial coordinate \tilde{r}, defined as the radial position relative to the boundary between the mixed interior and the horizontal cell periphery. Results are shown for the continuum model (left; radial cell density in units of μm⁻²), the experimental biofilm (middle; radial cell density in each micrometre-sized bin averaged over an observation window of 50 min) and the agent-based model biofilm (right; radial cell density in each micrometre-sized bin averaged for ten biofilms over an observation window of 6 min). For the continuum model and the agent-based model biofilms, the parameters were chosen to match those obtained from the experiment (Supplementary Figs. 12 and 13). The inset in the left-most panel shows the fraction of vertical cells in the continuum model at a given radius from the biofilm centre (grey regions contain no cells; colour scale is the same as in Fig. 1).

Full size image

As the threshold surface density for verticalization {\tilde{\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.

a, Top-down (upper row) and side views (lower row) of experimental biofilms grown with 0.4 μg ml⁻¹ A22 (magenta), without treatment (yellow) and with 4 μg ml⁻¹ cefalexin (cyan), following overnight growth (upper row) and 7 h after inoculation (lower row). Scale bar, 10 μm. Insets show magnifications of 10-μm²-sized regions of top-down views taken from the peripheries of biofilms. b, Expansion speed c^*, defined as the speed of the biofilm edge along the surface, versus the initial cell cylinder length {\ell }_0 for experimental biofilms (A22, magenta; no treatment, yellow; cefalexin, cyan), agent-based model biofilms (black circles) and the continuum model (dashed black curve). Expansion velocities were determined from a linear fit of the basal radius R_B of the biofilm versus time, where R_B is defined at each time point as the radius of a circle with area equal to that of the biofilm base. For experimental biofilms, the boundary was extracted from the normalized fluorescence data (see Methods for details). For each treatment, the vertical error bars show the standard error of the mean of the expansion speed and the horizontal error bars bound the measured initial cell cylinder length (Supplementary Fig. 1). Inset: model cells with lengths and radii corresponding to the averages for different treatments. c, Biofilm aspect ratio H/R_B for experimental biofilms grown under different treatments, where the biofilm height is defined as H=3V\mathrm{∕}2R_{\mathrm{B}}^2, the height of a semi-ellipsoid with a circular base of radius R_B and volume V equal to that of the biofilm. The error bars show the standard error of the mean. Inset: overlay of biofilm outlines from the bottom row of a. The colour designations and treatments are the same as in a.

Full size image

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{∕}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.

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³⁷.

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³⁹ 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⁴⁴, 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.

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⁴⁷. 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₄, 100 μM CaCl₂ and 0.5% glucose. When indicated, cefalexin (Sigma Aldrich) was added at 4 μg ml⁻¹ and A22 (a gift from the Gitai group) was used at 0.4 μg ml⁻¹. 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⁸. 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⁸. 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⁴⁸. 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⁻⁸ E₀R² to the force acting on the centre-of-mass r and 10⁻⁸ E₀R³ 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²⁹. 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⁴⁹, 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 δ₀.

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 {\tilde{\rho }}_{\mathrm{tot}} defined as the sum ρ_h + ξρ_v, where ξ is the ratio of vertical to horizontal cell footprints. In regions where {\tilde{\rho }}_{\mathrm{tot}} exceeds the close-packing density {\tilde{\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 η₁, where surface pressure is approximated as linearly proportional to the areal deformation {\tilde{\rho }}_{\mathrm{tot}}-{\tilde{\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 {\tilde{\rho }}_{\mathrm{tot}}>{\tilde{\rho }}_{\mathrm{t}}, where {\tilde{\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⁵⁰. 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⁸. 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.

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

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.

Affiliations

1. Joseph Henry Laboratories of Physics, Princeton University, Princeton, NJ, USA

   • Farzan Beroz

2. Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ, USA

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

3. Department of Molecular Biology, Princeton University, Princeton, NJ, USA

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

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

   • Yigal Meir

5. ICS-2/IAS-2, Forschungszentrum Jülich, Jülich, Germany

   • Benedikt Sabass

6. Howard Hughes Medical Institute, Chevy Chase, MD, USA

   • Bonnie L. Bassler

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

   • Ned S. Wingreen

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.

1. Supplementary Information

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

2. Reporting Summary

3. 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² and the total duration is 8 hours with 10 min time steps.

4. 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.

5. 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.

6. 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 ({\tilde{\rho }}_{\mathrm{tot}}= ρ_h + ξρ_v, black), versus radial coordinate r for isobaric regime with {\tilde{\rho }}_0 = 1 m⁻², {\tilde{\rho }}_{\mathrm{t}} = 1.5 m⁻², β = 2.5α, and ξ = 0.5. Dashed grey line shows {\tilde{\rho }}_{\mathrm{t}}.

7. 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.

Publication history

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.