Bacterial strains

The B. subtilis strains used in this study are listed in Supplementary Table 1. The capBCAE deletion strain was created by amplifying the 5′ arm and 3′ arm fragments of the operon using PCR and assembling the fragments into a pER449 chromosomal integration vector using Gibson assembly. The capBCAE overexpression strain was created by amplifying the 5′ arm and 3′ arm fragments of the operon and assembling the fragments into an ECE174 chromosomal integration vector using Gibson assembly. The vector backbone contains the IPTG-inducible hyperspank promoter. All primers used are indicated in Supplementary Table 2, with overlapping regions for Gibson assembly indicated by lowercase letters. The hag deletion strain was derived from previous work60 with the Phag-mCherry motile cell reporter transformed into the strain. All transformations were performed using a standard protocol61. Where appropriate, growth media were supplemented with antibiotics at the following concentrations: 5 mg l−1 chloramphenicol, 9 mg l−1 neomycin, 300 mg l−1 spectinomycin, 6 mg l−1 tetracycline and 9 mg l−1 kanamycin.

Growth conditions

Biofilms were grown on MSgg media (5 mM potassium phosphate buffer (pH 7.0), 100 mM 3-(N-morpholino)propanesulfonic acid (MOPS) buffer (pH 7.0, adjusted with NaOH), 2 mM MgCl2, 700 μM CaCl2, 50 μM MnCl2, 100 μM FeCl3, 1 μM ZnCl2, 2 μM thiamine HCl, 0.5% (v/v) glycerol and 0.5% (w/v) monosodium glutamate), previously reported to promote biofilm formation24. We used a version of MSgg with reduced buffer (1 mM MOPS) to enable the pH of the overall medium to be lowered to pH 4 (adjusted with HCl). This MSgg with 1 mM MOPS was previously used and shown to have no significant effects on B. subtilis biofilm growth62. The hyperspank promoter was fully induced using 1 mM IPTG.

Growth in microfluidic device

We grew biofilms in commercially available CellASIC B04F microfluidic plates (Millipore-Sigma) following a previously published protocol63. The day before the experiment, we streaked the desired strains from −80 °C glycerol stocks onto Luria-Bertani (LB) agar plates with appropriate antibiotics and incubated the plates overnight at 37 °C. The following day, a single colony was inoculated into 3 ml LB liquid media and incubated with shaking at 37 °C for 5 h. Cells were pelleted at 2,100 × g for 1 min and resuspended in fresh MSgg. We loaded the cells into the microfluidic chamber following the manufacturer’s protocol. The B04F microfluidic plates have two pillars where cells can be trapped, providing an anchor for the biofilm to grow. After cell loading, we purged cells stuck under one of the two pillars in the growth chamber, leaving biofilm growth at only one of the two pillars. Cells were incubated at 30 °C with media being supplied at 1.5 psi.

To induce biofilm dispersion, we followed previous work25,26,27, which removed nutrients from biofilms and increased nutrients. We incubated biofilms in the microfluidic device using standard MSgg. After biofilms reached the desired size, typically around 30 h, we introduced MSgg media without glycerol or glutamate (MS), the sole carbon and nitrogen sources, respectively. After 16 h, we removed MS and flowed in standard MSgg for the remainder of the experiment, up to 72 h.

To aid the growth of the γ-PGA overexpression strain, we supplemented MSgg with 30 mM d-glutamic acid throughout the experiment. Given that γ-PGA is comprised of both l- and d-glutamic acid, adding d-glutamic acid to the medium provides a substrate for increased γ-PGA production.

Time-lapse microscopy

Biofilms were observed using time-lapse microscopy with an Olympus IX83 inverted microscope (Evident Scientific) with a 10× objective (Olympus) and an X-Cite NOVEM light source (Excelitas Technologies). Images were taken at intervals of 30 s (Fig. 4) or 15 min (remaining figures) using an ORCA-Fusion BT camera (Hamamatsu).

γ-PGA production measurement

To measure how much γ-PGA is produced, we adapted an established protocol64. First, we grew wild-type, ΔcapBCAE and Phyp-capBCAE strains in liquid LB at 37 °C with shaking. After 4 h, all cultures were normalized to OD600 at 0.1 and resuspended in MSgg. We supplemented the Phyp-capBCAE strain with 30 mM d-glutamic acid as above. After 24 h growth at 37 °C with shaking, we centrifuged the culture at 2,100 × g for 10 min and kept only the supernatant. Next, we added 10 μM methylene blue to stain the γ-PGA. After 30 min incubation, we filtered the solution through a 0.22 μm filter. γ-PGA adsorbed the methylene blue but was too large to pass through the filter. We measured the methylene blue remaining in the filtrate using a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific) at 655 nm.

Particle image velocimetry

To analyse the biofilm movies and track the movement of cells, we first spatially aligned the image stack using the MultiStackReg plug-in in FIJI (RRID:SCR_002285). Four images were taken for each biofilm, for each time point, one image for each quadrant. These four image stacks were stitched together using the BigStitcher plug-in for FIJI. To isolate the biofilm region for analysis, the region outside the biofilm was manually excluded. Next, we used the Particle Image Velocimetry (PIV) plug-in in ImageJ for each pair of successive timeframes to analyse the displacement of cells in the biofilm. To calculate the trajectory of cells across the movie, we analysed the PIV output in Python (RRID:SCR_008394). We generated the trajectory paths of each cell across time points in Python by discretely updating each cell’s position based on the PIV-derived displacement vectors. At intermediate positions on our PIV grid, we used SciPy65 to interpolate the displacement vectors.

Particle tracking

To further analyse the movement of cells in a biofilm, we tracked visible features in the phase-contrast time-lapse images. We used FIJI66 to annotate the X and Y coordinates of these features for each time point, and the X and Y coordinates of the edge of the biofilm. These results were then visualized using Python. This technique was used for Extended Data Figs. 1 and 3.

Fluorescent reporter analysis

To track the dispersing region, we monitored a fluorescent reporter for motile cells, Phag-yfp, using time-lapse microscopy. To visualize the dispersing region, we first binarized the fluorescence image using the Default threshold algorithm in FIJI, which allowed us to compare between different experiments and conditions. We saved the resulting image as a table of grey values using the Transform function in FIJI. We calculated and plotted the mean value for every pixel row of the image. This technique was used for Fig. 3b,d,f,h and Extended Data Fig. 7.

γ-PGA overproduction analysis

To quantify the empty spaces inside the biofilm of the wild-type (WT) and γ-PGA overproduction strains, we first inverted and binarized the image using the Default threshold algorithm in FIJI. Areas where cells are located will have a value of 255, and areas devoid of cells will have a value of 0. We excluded the region outside the biofilm and the region under the cell trap for analysis by selecting regions of interest for these regions and reassigning these pixels NaN values using the changeValues function in FIJI. This exclusion will result in images with different total numbers of pixels, as biofilms are of different sizes. To compare between conditions, we divided the number of white pixels (value 255) or the number of black pixels (value 0) by the total number of pixels in the image.

pH modulation analysis

To quantify the expansion and contraction of the biofilm owing to pH changes, we created a ‘difference movie’ by subtracting the pixel intensities of successive time points for each frame of the phase-contrast time-lapse movie. The pixel shifts shown by this difference movie were tracked using the Multi Measure function in FIJI.

Cluster analysis

To compare dispersion from a WT biofilm versus a γ-PGA overproduction biofilm grown without starvation, we analysed the number of connected clusters for each biofilm. Using Python, we binarized the phase-contrast image to differentiate where cells are present or absent. Using the NetworkX library, we generated a network based on the maximal distance between biomass pixels. We quantified the number of connected clusters in the network for each biofilm using the KDTree function from SciPy65. To eliminate debris from the microfluidic device, we disregarded clusters below a threshold size.

Statistics and data reproducibility

For all experiments, at least n = 3 independent biofilms were used unless otherwise noted. Independent experiments were performed on separate days with new media prepared for each experiment. The exact number of independent biofilms is indicated in the figure legends where the data appear. No tests of statistical significance were necessary and were thus not performed.

Mathematical model

We show our biofilm system using a mathematical model based on partial differential equations for the densities of the two cellular populations (motile and non-motile), a hydrogel field and a nutrient field. The custom Python-based implementation of the simulations is publicly available at https://github.com/dsb-lab/dispersion. The parameter values used in the model are listed in Supplementary Table 3. We consider a two-dimensional continuum model for a growing biofilm composed of two cell phenotypes: a motile cell population with density \({\rho }_{m}\left({\bf{x}},t\right)\) and a non-motile (matrix) cell population with density \({\rho }_{n}\left({\bf{x}},t\right)\). Cell growth depends on the local availability of a nutrient field \(C\left({\bf{x}},t\right)\). The motile population produces an extracellular hydrogel field of γ-PGA, denoted by \(H\left({\bf{x}},t\right)\).

Biomass growth and γ-PGA hydrogel production act as local volumetric sources that generate an effective pressure field \(P\left({\bf{x}},t\right)\), giving rise to a Darcy-type velocity field \({\bf{U}}\left({\bf{x}},t\right)\), which we assume to be unique for all the fields considered below67. The velocity field is therefore not divergence-free; instead, its divergence is determined by local biomass growth and γ-PGA production. In the absence of growth-induced volumetric sources, the velocity field becomes divergence-free.

The transport of \({\rho }_{m}\), \({\rho }_{n}\) and \(H\) is governed by advection by \({\bf{U}}\), supplemented by diffusive terms assuming a small miscibility between fields68. The nutrient field undergoes diffusion, consumption by the biofilm and replenishment from an external reservoir67.

The final model used in this paper was written in a non-dimensional form; the corresponding dimensional model and the non-dimensionalization procedure are detailed in the following sections.

Dimensional model

We are working with a fluid chamber of small height, so we can consider a quasi-two-dimensional system. Then, let \({\bf{x}}\in \Omega \subset {{\mathbb{R}}}^{2}\) and \(t\ge 0\). We define the following dimensional fields:

  • \({\rho }_{m}\left({\bf{x}},t\right)\), \({\rho }_{n}\left({\bf{x}},t\right)\): densities of motile and non-motile (matrix) cells

  • \(H\left({\bf{x}},t\right)\): extracellular hydrogel (γ-PGA) concentration

  • \(C\left({\bf{x}},t\right)\): nutrient concentration

  • \(P\left({\bf{x}},t\right)\): effective pressure

  • \({\bf{U}}\left({\bf{x}},t\right)\): local velocity

We denote the total biomass by \({\rho }_{{\rm{tot}}}\left({\bf{x}},t\right)={\rho }_{m}\left({\bf{x}},t\right)+{\rho }_{n}\left({\bf{x}},t\right)\).

Dimensional model—biomass dynamics

The dynamics of the cell densities are given by

$$\frac{\partial {\rho }_{m}\left({\bf{x}},t\right)}{\partial t}=\alpha \,{\rho }_{m}{\phi }_{\rho }(C)-\nabla \cdot ({\bf{U}}{\rho }_{m})+{D}_{\rho }{\nabla }^{2}{\rho }_{m},$$

$$\frac{\partial {\rho }_{n}\left({\bf{x}},t\right)}{\partial t}=\alpha \,{\rho }_{n}\,{\phi }_{\rho }\left(C\right)-\nabla \cdot \left({{\bf{U}}\rho }_{n}\right)+{D}_{\rho }{\nabla }^{2}{\rho }_{n}.$$

Here \(\alpha\) denotes the maximum growth rate, the advective term \(-\nabla \cdot \left({{\bf{U}}\rho }_{i}\right)\) describes the transport of biomass by the Darcy velocity field, and \({D}_{\rho }\) is a small effective diffusivity that captures a small miscibility. Cell growth is regulated by a saturating Michaelis–Menten-type response to nutrient availability (where Kρ is the growth-rate half-maximal nutrient concentration):

$${\phi }_{\rho }\left(C\right)=\frac{C}{{K}_{\rho }+C}.$$

Dimensional model—hydrogel dynamics

The extracellular polymer field evolves according to

$$\frac{\partial H\left({\bf{x}},t\right)}{\partial t}=\beta {\phi }_{h}\left(C\right){\rho }_{m}-\nabla \cdot ({\bf{U}}H)+{D}_{h}{\nabla }^{2}H,$$

where \(\beta\) is the γ-PGA production rate, the advective term \(-\nabla \cdot ({\bf{U}}H)\) describes the transport of hydrogel by the Darcy velocity field, and \({D}_{h}\) is the γ-PGA diffusivity. γ-PGA production is activated by nutrient availability (where Kh is the γ-PGA half-maximal nutrient concentration) through

$${\phi }_{h}\left(C\right)=\frac{C}{{K}_{h}+C}.$$

We assume that polymer swelling saturates on a timescale much shorter than that of the biofilm growth and mechanical dynamics considered here. Consequently, we neglect explicit autocatalytic or self-swelling terms in the hydrogel dynamics. Under these assumptions, hydrogel dynamics is captured by a balance between production, advection and diffusion.

Dimensional model—nutrient dynamics

The nutrient field obeys

$$\frac{\partial {\rm{C}}\left({\bf{x}},t\right)}{\partial t}=\eta q({\rho }_{{\rm{tot}}},H)\left({C}_{0}-C\right)-{\gamma }_{\rho }{\phi }_{\rho }\left(C\right){\rho }_{{\rm{tot}}}-{\gamma }_{h}{\phi }_{h}\left(C\right){\rho }_{m}+{D}_{c}{\nabla }^{2}C.$$

The first term in this equation represents nutrient replenishment from an external reservoir with concentration \({C}_{0}\) (the set point of the microfluidic nutrient flow), where η is the nutrient replenishment rate from external flow. This influx is modulated by biomass and hydrogel loading through the factor

$$q\left({\rho }_{{\rm{tot}}},H\right)=\frac{{K}_{c}}{{\rho }_{{\rm{tot}}}+H+{K}_{c}},$$

where Kc is the nutrient penetration coefficient, which reduces nutrient exchange in regions of high biomass or hydrogel density, thereby mimicking the limited penetration of the external nutrient flow into the biofilm.

We assume that nutrient transport occurs in a diffusion-dominated regime and therefore neglect the advective transport of nutrients by the biofilm velocity field. The parameter \({\gamma }_{\rho }\) denotes the nutrient consumption rate associated with cellular growth, whereas \({\gamma }_{h}\) accounts for nutrient consumption linked to γ-PGA production.

Dimensional model—Darcy’s flow and pressure equation

We model biofilm mechanics using a Darcy-type closure,

$${\bf{U}}=-M\left({\rho }_{{\rm{tot}}},H\right)\,\nabla P,$$

where \(P\) is the pressure, calculated as discussed below, and the mobility is \(M\), which depends on the biomass (total cell density) and on the hydrogel according to

$$M\left({\rho }_{{\rm{tot}}},H\right)={M}_{\rho }{\rho }_{{\rm{tot}}}+{M}_{{\rm{h}}}H+{M}_{0}\theta \left(1-{\rho }_{{\rm{tot}}}-H\right).$$

The parameter \({M}_{0}\) in this expression represents the mobility of the surrounding fluid; the \({M}_{\rho }\) and \({M}_{{\rm{h}}}\) terms account for the contributions of biomass and hydrogel to the effective mobility, respectively; and \(\theta\) is the Heaviside function. In particular, we consider \({M}_{h} < {M}^{0}\), reflecting the polymeric, hydrogel-like nature of γ-PGA, which has a lower intrinsic mobility (that is, higher effective viscosity) than the surrounding fluid. Importantly, we have \({M}_{\rho } < {M}_{h}\), meaning that the hydrogel phase is more mobile than the densely packed biomass. Therefore, in regions initially dominated by dense biomass, the production and swelling of γ-PGA lead to an increase in local mobility. Biologically, this process is associated with a redistribution of biomass and a reduction in local cell density, effectively increasing the porosity of the biofilm. This results in an increase in effective large-scale mobility, consistent with our model. Biomass growth and γ-PGA production act as volumetric sources of the biofilm material. We define the total volumetric source as

$$Q=\alpha {\phi }_{\rho }\left(C\right){\rho }_{{\rm{tot}}}+\beta {\phi }_{h}\left(C\right){\rho }_{m},$$

where the first term on the right-hand side represents the contribution from cellular growth and the second term corresponds to hydrogel production.

These volumetric sources determine the divergence of the velocity field according to

$$\nabla \cdot {\bf{U}}=Q.$$

Substituting the Darcy relation \({\bf{U}}=-M\nabla P\) into this constraint yields a variable-coefficient Poisson equation for the pressure,

$$\nabla \cdot \left(M\nabla P\right)=-Q.$$

Non-dimensionalization of the model

We now define non-dimensional variables as

$$\begin{array}{cccc} c=\frac{C}{{C}_{0}}, & {\tilde{\rho }}_{m}=\frac{{\rho }_{m}}{{\rho }_{0}}, & {\tilde{\rho }}_{n}=\frac{{\rho }_{n}}{{\rho }_{0}}, & h=\frac{H}{{H}_{\rho }}\end{array},$$

$$\begin{array}{cc} \widetilde{t}=\frac{t}{T}=\alpha t, & \widetilde{{\bf{x}}}=\frac{{\bf{x}}}{L}\end{array}.$$

We also define \({\widetilde{\rho }}_{{\rm{tot}}}={\widetilde{\rho }}_{m}+{\widetilde{\rho }}_{n}\). For notation simplicity, we will omit the tilde in cell density in what follows.

The scaling factors arise from characteristic scales for nutrient, biomass, hydrogel, time and length:

  • Nutrient field: \({C}_{0}\) (reservoir concentration)

  • Biomass: \({\rho }_{0}\) (packing density set by the initial conditions)

  • Hydrogel: \({H}_{\rho }\) (hydrogel concentration equivalent in volume to \({\rho }_{0}\), set by the initial conditions)

  • Time: \(T=\frac{1}{\alpha }\) (growth timescale)

  • Length: \(L\) (arbitrary reference scale)

    The scaling relations defined above lead us to the following dimensionless parameters:

  • Diffusion

    $$\begin{array}{rcl} {\widetilde{D}}_{c}=\frac{{D}_{c}}{\alpha {L}^{2}}, & {\widetilde{D}}_{h}=\frac{{D}_{h}}{\alpha {L}^{2}}, & {\widetilde{D}}_{\rho }=\frac{{D}_{\rho }}{\alpha {L}^{2}}\end{array}.$$

  • ECM production

    $$\widetilde{\beta }=\frac{\beta {\rho }_{0}\,}{\alpha {H}_{\rho }}.$$

  • Nutrient consumption

    $$\begin{array}{cccc}{\tilde{\gamma }}_{\rho }=\frac{{\gamma }_{\rho }{\rho }_{0}}{\alpha {C}_{0}}, & {\tilde{\gamma }}_{h}=\frac{{\gamma }_{h}{\rho }_{0}}{\alpha {C}_{0}}, & \tilde{\eta }=\frac{\eta }{\alpha }, & {\tilde{K}}_{c}=\frac{Kc}{{\rho }_{0}}.\end{array}$$

  • Nutrient influx

    $$\widetilde{q}\left({\rho }_{\mathrm{tot}},h\right)=\frac{{\widetilde{K}}_{c}}{{\rho }_{\mathrm{tot}}+\left(\frac{{H}_{\rho }}{{\rho }_{0}}\right)h+{\widetilde{K}}_{c}}=\frac{{\widetilde{K}}_{c}}{{\rho }_{\mathrm{tot}}+h+{\widetilde{K}}_{c}}.$$

    where we express hydrogel in the same normalized units as biomass (that is, \({H}_{\rho }\approx {\rho }_{0}\)).

  • Biomass and hydrogel growth factors

    $$\begin{array}{cc}{\widetilde{\phi }}_{\rho }\left(c\right)=\frac{\left({\chi }_{\rho }c\right)}{1+\left({\chi }_{\rho }c\right)},\, & {\widetilde{\phi }}_{h}\left(c\right)=\frac{\left({\chi }_{h}c\right)}{1+\left({\chi }_{h}c\right)},\end{array}$$

    with \({\chi }_{\rho }=\frac{{C}_{0}}{{K}_{\rho }}\) and \({\chi }_{h}=\frac{{C}_{0}}{{K}_{h}}\).

  • Mobility

    $$\widetilde{M}\left({\rho }_{{\rm{tot}}},h\right)={\rho }_{{\rm{tot}}}+{\widetilde{M}}_{h}h+{\widetilde{M}}_{0}\theta \left(1-{\rho }_{{\rm{tot}}}-h\right),$$

    where \({\widetilde{M}}_{0}={M}_{0}/{M}_{\rho }\) represents the mobility of the surrounding fluid relative to the biomass scale and \({\widetilde{M}}_{h}={M}_{h}/{M}_{\rho }\) quantifies the relative (dimensionless) mobility of the γ-PGA hydrogel.

  • Darcy’s law

    $${\bf{u}}=-\widetilde{M}\nabla p,$$

    where \(\widetilde{\mu }\) is defined above, and the dimensionless pressure is

    $$\begin{array}{cc}\displaystyle p=\frac{P}{\widetilde{P}}, & \mathrm{with}\,\widetilde{P}=\frac{{L}^{2}\alpha }{{M}_{\rho }}\end{array}.$$

The non-dimensional pressure can be calculated via

$$\nabla \cdot \left(\widetilde{M}\nabla p\right)=-{Q}_{{\rm{tot}}},$$

where

$${Q}_{\mathrm{tot}}={\rho }_{m}{\widetilde{\phi }}_{\rho }\left(c\right)+{\rho }_{n}{\widetilde{\phi }}_{\rho }\left(c\right)+\beta {\widetilde{\phi }}_{h}\left(c\right){\rho }_{m}.$$

Full non-dimensional model

With these rescaling and dimensionless parameters, and dropping the remaining tildes for readability, the final non-dimensional model is

$${\partial }_{t}{\rho }_{m}={\rho }_{m}{\phi }_{\rho }\left(c\right)-\nabla \cdot \left({\rho }_{m}{\bf{u}}\right)+{D}_{\rho }{\nabla }^{2}{\rho }_{m},$$

$${\partial }_{t}{\rho }_{n}={\rho }_{n}{\phi }_{\rho }\left(c\right)-\nabla \cdot \left({\rho }_{n}{\bf{u}}\right)+{D}_{\rho }{\nabla }^{2}{\rho }_{n},$$

$${\partial }_{t}h=\beta {\phi }_{h}\left(c\right){\rho }_{m}-\nabla \cdot \left(h{\bf{u}}\right)+{D}_{h}{\nabla }^{2}h,$$

$${\partial }_{t}c=\eta q({\rho }_{{\rm{tot}}},h)\left(1-c\right)-{\gamma }_{\rho }{\phi }_{\rho }\left(c\right){\rho }_{{\rm{tot}}}-{\gamma }_{h}{\phi }_{h}\left(c\right){\rho }_{m}+{D}_{c}{\nabla }^{2}c.$$

We note that parameters (for example, α, C0 and Mρ) are used as reference scales in the non-dimensionalization procedure and are therefore absorbed into the definition of the dimensionless variables. Consequently, they do not appear as independent parameters in the final system, which is equivalent to setting them equal to 1.

γ-PGA overproduction in the model

To model γ-PGA overproduction, we introduce an overproduction factor \({\eta }_{{\rm{OE}}}\ge 1\) that multiplies the effective polymer production rate. In the dimensionless formulation, this is implemented by rescaling the γ-PGA production parameter \(\beta\) as

$$\beta \to \,{\eta }_{\mathrm{OE}}\beta .$$

As all cells in the γ-PGA overproduction strain contain the inducible capBCAE operon, we consider that matrix cells are also able to produce γ-PGA in this strain. We assume this production to be heterogeneous. Gene expression in bacterial populations is inherently heterogeneous, even in phenotypically defined subpopulations. In B. subtilis biofilms, this heterogeneity extends to the matrix-producing fraction and persists under inducible expression69,70,71. We represent the fraction of matrix cells that produce γ-PGA by a third cell-density field, which has the same effect on the hydrogel field as the motile cells, but is made up of non-motile cells.

To model supplementation with d-glutamic acid, we introduce a rescaling factor \({\eta }_{C}\) for the nutrient field in the dimensional set-point concentration \({C}_{0}\). Guided by the experimental design, we assume that the added d-glutamic acid only increases the ability of cells to produce γ-PGA, without modifying the growth terms directly. After the non-dimensionalization, this supplementation leads to the following parameter rescaling:

$${\chi }_{h}\to {\eta }_{C}{\chi }_{h},$$

$${\gamma }_{h}\to {\eta }_{{OE}}{\gamma }_{h}/{\eta }_{C}.$$

We note that the overproduction factor \({\eta }_{{\rm{OE}}}\) also affects the γ-PGA-associated consumption term parametrized by \({\gamma }_{h}\), as polymer secretion requires metabolic resources.

Initial conditions, boundary conditions and numerical integration parameters

Initial biomass is prescribed as a compact colony centred within the simulation domain, with geometric perturbations at the periphery to break radial symmetry. The simulation domain is a two-dimensional square \(\left[0,{L}_{x}]\times [0,{L}_{y}\right]\) with \({L}_{x}={L}_{y}=20\) (dimensionless units), discretized on a uniform grid of \({N}_{x}={N}_{y}=1,024\) cells. The pressure field is solved using Dirichlet boundary conditions, \(P=0\), at the domain boundary. For the remaining fields, Dirichlet boundary conditions are imposed consistently with the numerical implementation: \({\rho }_{m}={\rho }_{n}=h=0\) and \(c=1\) at the boundary. Advection terms are discretized using a first-order upwind scheme, whereas diffusive and pressure operators are approximated using finite differences on the non-periodic grid.

Initial biomass geometry and phenotype partition

The two phenotypes are initialized in distinct subregions: a core colony region with a protruding plume-like region and a periphery.

  • A baseline circular colony of radius \({R}_{{\rm{bio}}}=1.5\) is placed at the domain centre \(\left({x}_{0},{y}_{0})=({L}_{x}/2,{L}_{y}/2\right)\).

  • The biofilm interface is perturbed by an angular Fourier noise constructed from a finite number of modes (\(k=1,\ldots ,6\)) with random phases. This produces a perturbed radius

    $${R}_{{\rm{bio}}}(\theta )={R}_{{\rm{bio}}}\left[1+{\varepsilon }_{{\rm{core}}}\eta (\theta )\right],$$

    with \({\varepsilon }_{{\rm{core}}}=0.08\), where \(\eta (\theta )\) is a normalized sum of cosine modes. This perturbation breaks the circular symmetry and seeds morphological heterogeneity at \(t=0\).

  • The matrix phenotype \({\rho }_{n}\) is first assigned to the biofilm region and bounded in [0, 1].

  • A plume region is then constructed and assigned to the motile phenotype \({\rho }_{m}\). The plume is defined as the union of the following:

  1. 1.

    a perturbed circular cap of radius \({R}_{{\rm{plume}}}=0.5\,\)(with perturbation amplitude \(0.10\)), and

  2. 2.

    a tapered, cone-like extension of length \({L}_{{\rm{plume}}}=3/2\,{R}_{{\rm{plume}}}\) and base half-width \({R}_{{\rm{core}}}=0.2\), whose width decreases with distance from the base.

Inside the plume mask, \({\rho }_{n}\) is set to 0 and \({\rho }_{m}\) is set to \(1\). In that way, we are able to mimic one of the characteristic branches that are formed by motile cells during biofilm maturation.

In the overproduction case, add \(N=50\) randomly located ‘spots’ of γ-PGA-producing matrix cells (see ‘γ-PGA overproduction in the model’) in the matrix region.

To avoid sharp interfaces, the initial fields are regularized using a smooth hyperbolic tangent profile with characteristic width \(\Delta x\).

Nutrient and hydrogel initialization

The nutrient is initialized uniformly at 0, \(c(x,0)=0\), and the γ-PGA hydrogel is initialized as \(h(x,0)=0\) everywhere.

Reporting summary

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