## Abstract

A general feature of mature biofilms is their highly heterogeneous architecture that partitions the microbial city into sectors with specific micro-environments. To understand how this heterogeneity arises, we have investigated the formation of a microbial community of the model organism *Bacillus subtilis*. We first show that the growth of macroscopic colonies is inhibited by the accumulation of ammoniacal by-products. By constraining biofilms to grow approximately as two-dimensional layers, we then find that the bacteria which differentiate to produce extracellular polymeric substances form tightly packed bacterial chains. In addition to the process of cellular chaining, the biomass stickiness also strongly hinders the reorganization of cells within the biofilm. Based on these observations, we then write a biomechanical model for the growth of the biofilm where the cell density is constant and the physical mechanism responsible for the spreading of the biomass is the pressure generated by the division of the bacteria. Besides reproducing the velocity field of the biomass across the biofilm, the model predicts that, although bacteria divide everywhere in the biofilm, fluctuations in the growth rates of the bacteria lead to a coarsening of the growing bacterial layer. This process of kinetic roughening ultimately leads to the formation of a rough biofilm surface exhibiting self-similar properties. Experimental measurements of the biofilm texture confirm these predictions.

## 1. Introduction

Because they evolve in fluctuating environments and compete for limited resources, wild-type strains of bacteria have developed various cooperative skills such as the ability to swarm or form surface-attached communities known as biofilms [1,2]. These ‘cities of microbes’ [3] are produced in response to various cues and have been shown to provide an advantage in many situations, from nutrient retention to protection against predators through increased antibiotic resistance. The ability to form a biofilm relies on the existence within the community of microorganisms able to secrete extracellular polymeric substances (EPSs) that ‘glue’ cells together and onto the surface [4,5].

The presence of this sticky matrix allows the biofilm to adopt a three-dimensional shape [1,6], as illustrated in figure 1*b*,*d*, similarly to the biological tissues of multicellular organisms [7]. This heterogeneous architecture, which consists of towers and channels, is one of the hallmarks of the biofilm lifestyle. By contrast, bacterial strains which do not produce an extracellular matrix form colonies with a rather homogeneous microstructure as illustrated in figure 1*a*,*c*. Quite surprisingly, this alteration of the microstructure turns out to have a profound impact on the ability of microorganisms to survive in their natural environment. As a result of the structural heterogeneity, the microbial city is partitioned into sectors with specific micro-environments. Microorganisms in turn respond to the local environmental conditions by tuning their metabolic activity and genetic expression. This leads to a dynamic functional regionalization of the biofilm [8] which is sometimes argued to represent a minimal form of multicellularity [9–11]. In a variable environment, the coexistence of multiple phenotypic states within the same community increases the probability that some offsprings are well adapted to the current environment [12]. In such conditions, this bet-hedging strategy ultimately results in an increased fitness of the heterogeneous population over the homogeneous one [13,14]. Besides being strongly coupled to the development of metabolic and phenotypic heterogeneities, the physical structure of the biofilm also serves specific functions. In wild-type strains of *Bacillus subtilis*, the formation of aerial structures is a necessary condition for the production of spores [5] and, thus, species survival.

However widespread, the ability to produce EPS is controlled by species-specific genetic circuitries. In the Gram-positive bacterium *B. subtilis*, a model organism for biofilm formation, the pathway controlling the production of this sticky matrix has been well characterized [15–18]. The matrix typically represents ∼5% of the total (dry) mass of the biofilm [19] and is primarily composed of exopolysaccharides and a protein (TasA). At a somewhat coarse-grained level, it appears clearly that the production of these biomolecules occurs in conditions of low nutrient and is enhanced when a signalling molecule, surfactin, reaches a critical threshold [20]. This lipopeptide is produced by surface colonies of *B. subtilis* and acts both as a surfactant [21] and as a signalling molecule allowing the cells to roughly estimate their density. While not universal, quorum sensing mechanisms are frequently involved in biofilm formation. Although important progress has been made to understand the molecular mechanisms that lead to biofilm formation and development, it is still unclear how the production of EPS results in a complex three-dimensional architecture rather than a flat aggregate as commonly observed on rich medium (as exemplified in figure 1*a*,*c*) or with domesticated strains that have lost the ability to form biofilms [22]. To understand how this structure is built and on which mechanisms it relies to do so, we investigate experimentally and theoretically the biomass production and reorganization in a microbial community of *B. subtilis* cells.

The structure of this paper is as follows. In §2, we analyse the growth kinetics of a macroscopic biofilm to identify the fields that contribute to the spreading of the biofilm. We find that growth saturates in response to the accumulation of metabolic by-products. In §3, we analyse the heterogeneities of the air–biofilm interface and find that (i) it roughens over time and (ii) it is self-similar in the range 1–100 μm. In §4, we draft a theoretical model for the biomass expansion within the framework of continuum mechanics. In order to discriminate between several modelling hypotheses and measure one of the model parameters at the single-cell level, we grow bacterial monolayers in §5. The microscopic observation of these quasi-bidimensionnal biofilms reveals that biomass density is approximately constant and that spreading occurs as a consequence of the mechanical pushing created by the growth and division of the bacteria. Based on these observations, we close our theoretical model in §6. The experimental data on the biomass velocity field are nicely reproduced. We then theoretically analyse in §7 the stability of the air–biofilm interface in the presence of noise in the reproduction rate of the bacteria. We predict a phenomenon of kinetic roughening of the biofilm free surface that leads to the formation of a self-similar heterogeneous architecture. The theoretical critical exponents are in agreement with our experimental data. Our results are then discussed in §8.

## 2. The growth kinetics of macroscopic *Bacillus subtilis* biofilm is inhibited by the accumulation of metabolic by-products

In order to develop a macroscopic model describing the biomass production and reorganization, we first inoculated large Petri dishes containing the biofilm-promoting nutritive medium Msgg (minimal salts glycerol glutamate) [23] gelified with 1.5% agarose. In this defined minimal medium, glycerol is a carbon source and glutamate can serve as both nitrogen and carbon source. Using a scanner/incubator, we followed the macroscopic evolution of biofilms for up to 150 h (see Material and methods). Representative pictures of the development are presented in figure 2.

First, we plotted the biofilm area *A* against time *t*. It reveals a linear time dependence *A* ∼ *t*. This observation indicates a slowdown in global growth compared to the classical case where growth is nutrient-limited. Indeed, in this classical case, nutrients penetrate the colony through the boundary and get consumed, giving rise to a ring of active growth at the outer boundary of the colony. The width *w* of this ring is fixed and is approximately given by , where *d* is the diffusion coefficient of the limiting nutrient and *u* the uptake rate of this nutrient by the bacteria. If the width of the ring is small compared with the radius *R* of the colony, then the area of this ring is given by 2*πwR*. From this, it follows that the area *A* = *πR*^{2} of the colony evolves according to the differential equation d*A*/d*t* = *b*(2*πwR*), where *b* is the growth rate of the bacteria. From this equation, we deduce that, when growth is nutrient-limited, *R* evolves linearly in time and thus *A* ∼ *t*^{2} [24]. Our experimental data indicate that the area grows slower than *t*^{2} and thus growth is not limited by the rarefaction of nutrient. To identify which other factor was inhibiting growth, we therefore conducted additional experiments. First, we noted that the pH of the agar plate at the end of the experiments reached a value of 9 ± 0.5 (initial value 7.5), a value at which the growth of *B. subtilis* cells is indeed inhibited [25]. Second, lowering the concentration of the 3-(*N*-morpholino)propanesulfonic acid (MOPS) buffer in the nutritive medium yielded smaller biofilms, again with a pH of 9 after biofilm growth. Together, these results indicate that a chemical is produced by the bacteria which changes the pH and is responsible for the inhibition of growth. This is because, at lower buffer concentration, a lower amount of chemical is needed to reach the same inhibiting pH of 9 and thus smaller biofilms are formed. Because of the high concentration of the MOPS buffer (100 mM) in the medium, only other compounds produced at high enough concentration could give rise to a change in pH from 7 to 9. The two other compounds initially present at high concentration in the nutritive medium are glycerol (approx. 54 mM) and glutamic acid (∼31 mM). Glycerol is transformed into pyruvate before entering the Krebs cycle (with no by-product) while glutamic acid is converted into alpha-ketoglutaric acid, producing ammonia as a by-product in the process. The presence of ammonia in the agar after biofilm growth was indeed confirmed using Nessler's reagent [26]. Since ammonia was indeed detected and was the only compound that could be produced in sufficient amounts to shift the pH, we concluded that ammonia was responsible for the inhibition of growth.

We then measured the velocity field of the biomass during biofilm development. To this end, we took advantage of the microstructure formed at the surface of the biofilms and followed over time the position of morphological features selected outside of the long radial ridges. A representative example of the biomass motion is given in figure 2*f*–*h*, where the position of a single granularity is indicated by a black circle at three consecutive times. In all cases the trajectories of the protrusions were mostly radial and the orthoradial displacements were of the order of the size of the protrusions. Note that, while this method allows us to quantify the macroscopic velocity field, more complex motions occurring on a smaller scale, such as rolling motions of the cells within the protrusions, cannot be detected. Several of those trajectories are represented in figure 2*i*, together with the position of the biofilm edge. As expected, the biomass velocity is higher at the boundary than closer to the centre of the biofilm. The boundary reaches a maximum speed of ∼250 μm h^{−1} at the edge of the biofilm, roughly ∼8 h after inoculation. The velocity then monotonically decreases over time. As can be seen in figure 2*i* however, growth is not restricted to the immediate vicinity of the biofilm edge. For example, at *t* = 18 h, the biomass radial velocity is ∼200 μm h^{−1} at the boundary, ∼100 μm h^{−1} at a distance of 2 mm away from the biofilm boundary and still of the order of ∼20 μm h^{−1} deep inside the biofilm, 5 mm away from the biofilm boundary.

## 3. The air–biofilm interface is self-similar and roughens over time

Next, we quantified the spatio-temporal distribution of the aerial protrusions that are the hallmark of biofilm formation. To this end, we used a mutant strain in which matrix-producing bacteria also synthesize a yellow fluorescent protein (see the Material and methods section). After 24 h of development, these cells make up 80% of the biomass and are homogeneously distributed [8]. Fluorescent images of a 24 h old biofilm that had not developed long radial ridges were collected with a microscope and stitched together to produce a high resolution map of the domains comprising matrix-producing bacteria (figure 3). This cartography reveals that domains of matrix-producing bacteria nucleate at the outer boundary of the biofilm, ripen and merge. This intensity map was then converted into a map of the biofilm height *h*(** r**) (see figure 4

*a–c*and Material and methods). Because the biofilm area increases linearly in time, we can reconstruct the time evolution of the biofilm–air interface

*h*(

**,**

*r**t*) (figure 4

*d*). The time evolution of the biofilm interface is then analysed using scaling arguments. This scaling approach is usually appealing in physical systems [27] and, more recently, biological systems [28–30]. It has been indeed recognized that the dynamics of many interfaces can be essentially characterized by a few exponents that are largely independent of the physical parameters of the systems but depend rather on the nature of the physical processes driving the growth of the interface. To each general mechanism of growth (such as ballistic deposition, diffusion limited aggregation, etc.) there corresponds a set of exponents known as a universality class that can be used to infer some knowledge about the underlying growth process. First the mean height increases as (figure 4

*e*), where the average is taken over the size

*L*of the system.

In addition to the thickening of the biofilm, its interfacial morphology evolves in time. This feature can be quantified using the typical distance between the valleys and peaks of the surface . This global interface width increases in time as *w*(*t*) ∼ *t ^{β}* with

*β*= 0.5 ± 0.1, indicating a gradual roughening of the interface (figure 4

*f*).

Using a single number to quantify the properties of the interface is somewhat restrictive because it gives only a measurement of the roughness at the global scale and does not indicate whether the peaks and the valleys themselves are rough. In other words, a more refined measure of the roughness should depend on the scale of observation. To investigate this dependance, let us consider the local surface width where the average is taken over all observation windows of size ℓ. The log–log plot of the local width *w*(*ℓ*, *t*) against ℓ for various times reveals essentially two different behaviours separated by a crossover length *ℓ*_{×} (figure 4*g*). For , the curves are linear, thus indicating a power law *w*(*ℓ*, *t*) ∼ *ℓ ^{αlocal}* = 0.6 ± 0.1 with

*α*

_{local}= 0.6 ± 0.1 and the system is rough at all length scales smaller than

*ℓ*

_{×}. Above

*ℓ*

_{×}, the width reaches a saturation value

*w*

_{sat}(

*t*) identical to the global interface width (averaged over the size

*L*of the system) which, as noted above, scales with the age of the biofilm as

*w*

_{sat}(

*t*) ∼

*t*. The crossover between the two regimes essentially occurs when the averaging windows are larger than the size of the largest feature of the interface.

^{β}*ℓ*

_{×}is therefore also a dynamic quantity and scales with time as

*t*

^{1/z}. This suggest the following scaling:

*w*(

*ℓ*,

*t*) ∼

*w*

_{sat}(

*t*)

*f*(

*ℓ*/

*ℓ*(

_{×}*t*)). Upon replacing

*w*

_{sat}(

*t*) and

*ℓ*

_{×}(

*t*) by their respective scaling form, all our data indeed fall on the same master curve

*f*(

*u*) (figure 4

*h*). We can also summarize those scaling relations using the Family–Vicsek ansatz: 3.1with the numerical values

*α*

_{local}= 0.6 ± 0.1,

*β*= 0.5 ± 0.1,

*β*−

*α*

_{local}/

*z*= 0.3 ± 0.2. Do the exponents measured here fall into any of the known university classes and can we learn anything from those measurements? It turns out that similar exponents have indeed been reported in the context of surface growth of metal [31] and polymer [32] films and form a class of scaling behaviour known as

*intrinsic anomalous roughening*[33]. Interestingly, such exponents are thought to stem from non-local effects of the growth process.

As we have shown in §2, growth is not restricted to the immediate vicinity of the boundary in our experimental set-up. Furthermore, local rearrangement of the cells is hindered by the high cohesiveness of the biomass. Consequently, we might expect that bacteria dividing deep inside the biofilm also affect the motion of the free boundary, leading to a non-local growth process.

## 4. Biomechanical model of biofilm growth

In order to explain our experimental observations on biofilm spreading and surface heterogeneities, at least two fields are necessary to build a minimal model. We need to describe how the cell density *c*(** r**,

*t*) and the concentration

*ρ*(

**,**

*r**t*) of the toxic by-product evolve in space

**and time**

*r**t*. The evolution of the by-product concentration can be unambiguously described by a reaction–diffusion equation of the form 4.1where Δ is the Laplace operator. For the sake of simplicity, we assume that the waste is produced by the bacteria at a constant rate

*P*and diffuses away from the biofilm with a diffusion coefficient

_{ρ}*D*. The local microbe concentration

_{ρ}*c*(

**,**

*r**t*) evolves according to two different mechanisms: bacteria reproduce at a rate

*κ*and move with a velocity

**(**

*v***,**

*r**t*), thus creating a flux

*c*. Accordingly,

**v***c*obeys the following continuity equation: 4.2In general, the growth rate

*κ*is a function of many factors such as nutrient concentration, pH, temperature, pressure [34], etc. As mentioned in §3, the limiting factor for the growth process within our experimental set-up is the accumulation of toxic by-products of the metabolism and we will only take into account the dependance of

*κ*on

*ρ*. Various models can be used to fit the data and we choose 4.3where

*κ*

_{0}is the maximum growth rate (i.e. in absence of any waste products),

*ρ*

_{i}is a typical inhibitory value for

*ρ*(i.e. when

*ρ*=

*ρ*

_{i}, the growth rate is half its maximum value) and

*m*characterizes how steeply growth is inhibited. It now remains to identify the physical mechanism(s) underlying the spreading of the biomass, i.e. what is the appropriate equation satisfied by the velocity field

**(**

*v***,**

*r**t*)?

In response to growth, microorganisms that are not bound together by an extracellular matrix move and collide with each other [35,36]. At a macroscopic scale, the biomass motion is described as diffusive with a velocity field proportional to the local cell concentration gradient [37–40].^{1} In the context of biofilm growth, both the sticky extracellular matrix and the stiffness of the substrate impair the motility of the bacteria on the surface. Furthermore, in the specific case of *B. subtilis* biofilms, the transcription factor Spo0A regulating biofilm formation also represses bacterial motility [8]. Indeed, although motile cells are present at some stages of *B. subtilis* biofilm development [8], experiments conducted on a flagella-null mutant strain show that motility does not significantly affect the expansion rate of *B. subtilis* biofilms [41]. What then is the physical mechanism driving the biomass expansion? This question has been raised in other problems of biological growth (such as tumour growth [42]) and several mechanisms have been suggested. On the one hand, dividing cells create a pressure on their neighbours, therefore moving the microorganisms from high to low pressure regions [43]. On the other hand, some biological tissues, and biofilms in particular [44,45], are cohesive materials. While they are usually incompressible (as most biological tissues comprise a high volume fraction of water), it also requires some energy to separate their constituents. These two effects typically result in the existence of an equilibrium cell density where this cohesive energy is minimal and any departure from this value gives rise to a restoring force [46,47]. As in the case of diffusive cells, this leads to a velocity field proportional to the local cell concentration gradient (albeit with a proportionality coefficient that depends on the cell concentration and can change sign).

Because biofilms are thick (up to hundreds of micrometres) three-dimensional structures, it has been so far difficult to assess the relative influence of those different effects and identify the relevant constitutive equation to be satisfied by ** v**. To better pinpoint the force driving the biomass expansion, we have grown

*B. subtilis*biofilms in an approximately two-dimensional geometry. Growing bacteria in monolayers also allows one to determine the doubling time of the bacteria in the biofilm state which is another parameter of our model.

## 5. The density of cells is constant in two-dimensional micro-colonies

To this end, we inoculated *B. subtilis* bacteria on a thin gel layer as represented in figure 5*a*. In addition, because the ability to form complex three-dimensional domains is tied to the production of an extracellular matrix, we have also monitored EPS production by using a fluorescent reporter protein to follow the activity of the *yqxM-sipW-tasA* operon. This procedure allowed us to identify the bacterial cells producing the protein TasA (the major component of the EPS). Between 1 and 10 motile cells (present either as singleton or doublet) were inoculated on a glass slide coated with a gel layer of biofilm-promoting nutritive medium Msgg [23]. The inter-bacterial distance was a few hundreds of micrometres at the beginning of an experiment (see Material and methods). Cells were then allowed to grow at room temperature under the microscope and both bright-field and fluorescent images were taken at 10 min intervals. Figure 5*b*–*e* shows the typical time course of an experiment.

After a lag phase where bacterial growth rate was slowly increasing (figure 5*b*), the microbes entered an exponential regime (figure 5*c*–*e*) at *t* ∼ 20 h. The first matrix-producing cells then appeared after 33 h and their number increased rapidly (figure 5*f*). Once the cluster of cells was growing exponentially, the division time was then measured on cells at the edge of the microcolonies and, when there was a single layer of bacteria, also at the centre. The doubling time was found to be approximately 85 min, although they were important variations (±15 min) from cell to cell. Note that this value is close to that observed in liquid culture (80 min). We then performed similar experiments with different environmental conditions (the substrate stiffness was varied by changing the agarose concentration of the nutritive gel layer in the interval 0.5–10% and the nutrient was diluted up to ten times). Up to a fourfold dilution, the doubling time of the bacteria did not change appreciably (see figure 5*g*), indicating that nutrient depletion is not the rate-limiting step for the growth process within this range. However, when diluting the nutrient ten times, the division time increased to 100 min. In addition, the length of the bacteria decreased slightly with decreasing nutrient concentration.

Furthermore, we measured bacterial density in windows of 10 × 10 μm^{2} at different locations and times in the two-dimensional microcolonies. We could not find any trend in the temporal or spatial variations of the cell density within the time scale of our experiment (with an average value of cells per 100 μm^{2}) and cells remained tightly packed. Smaller density fluctuations could not be resolved with our experimental set-up. This was true for both matrix-producing and undifferentiated bacteria although their organization was different (figure 5*h*–*j*). Bacteria that did not produce the protein TasA formed tightly packed layers with a polycrystalline two-dimensional order (for small numbers of cells, the alignment between a cell and its lineage is broken by a snapping mechanism (D. Bensimon group, personal communication)). Matrix-producing cells however, tend to form long bacterial chains [15], through a downregulation of the expression of cell wall-degrading enzymes [48,49] associated with the EPS-mediated physical gluing of a mother cell with its lineage. Their number increased and in the late stages (after ∼24 h) of biofilm growth, all the bacteria at the biofilm boundary were producing EPS. At this point, bacteria were mostly aligned in one-dimensional flexible ‘crystals’ or bundles. While there could be systematic variations in the cell density (on a scale larger than the size of the largest observable monolayer, i.e. around 100 μm), these measurements provide an upper bound for the amplitude of the density gradients that could occur inside the bulk of the biofilm. Note that these fluctuations are much weaker than the jump in density occurring at the boundary of the biofilms. From this observation, we concluded that cell density fluctuations in the bulk of the biofilms can be ignored in the description of small two-dimensional microcolonies. In absence of additional data for the density in three-dimensional biofilms, we postulated the same to be true for macroscopic biofilms. As a consequence, the primary mechanism responsible for the biomass spreading in our biomechanical model is the pressure generated by dividing bacteria on the neighbouring cells.

## 6. Closure of the biofilm model and comparison with the macroscopic experiments

When they divide, bacteria exert a pressure *p*(** r**,

*t*) on their neighbours. Because bacteria are essentially incompressible and since we did not observe significant cellular density fluctuations in our experimental set-up, bacteria are not compacted nor compressed by this pressure, but rather displaced from regions of high pressure towards regions of lower pressure, with a velocity

**. The simplest linear relationship corresponding to this mechanism is known as Darcy's law [50] 6.1where**

*v**λ*is a material parameter. This relation can be interpreted as the condition of mechanical equilibrium in an overdamped regime, in which a constant force is needed to maintain a steady velocity. Outside the biofilm, there is no pressure generated and thus

*p*= 0 and

**=**

*v***0**. Furthermore, the incompressibility condition implies that the local volume variation is entirely due to cell proliferation: 6.2In addition, the assumption of constant cell density (denoted

*c*

_{0}) in the biofilm reduces (4.2) to (6.2). Plugging (6.1) into (6.2) and introducing the dimensionless variables , with

*L*being a typical length, we obtain the following equations inside and outside of the biofilm (tildes are dropped for clarity):

*Inside*:
6.3*Outside*:
6.4where *C*_{1} = *D _{ρ}*/

*L*

^{2}

*κ*

_{0}is the ratio of the characteristic time scale of growth over the time scale of diffusion of the waste-product and

*C*

_{2}=

*c*

_{0}

*P*/(

_{ρ}*ρ*

_{i}

*κ*

_{0}) is the time scale of growth over the time scale of end-product production.

The set of equations (6.3) and (6.4) form our model for the biomass production and spreading, in conjunction with the diffusion of the by-product. Given an initial biofilm patch, equations (6.3) and (6.4) are solved numerically in two dimensions and the velocity ** v** is used to evolve the patch radius. Although

**depends only on the local pressure gradient, this is a non-local model. Indeed, because the pressure field must satisfy the Poisson equation, the interface velocity depends on the overall amount of mass created throughout the biofilm and not just on the local cell concentration gradient. The maximum growth rate value is**

*v**κ*

_{0}= 1/85 min

^{−1}, as measured in the previous section. The by-product diffusion coefficient

*D*is chosen to be that of ammonia in water, i.e.

_{ρ}*D*= 7.10

_{ρ}^{−9}m

^{2}s

^{−1}. Choosing

*L*= 1 mm (the typical initial radius of the biofilm), we get

*C*

_{1}= 36. Fitting the numerical results to the experimental data for the area gives

*C*

_{2}= 7.5 and

*m*= 1.2 for the stiffness coefficient. As seen in figure 2

*e*, this model correctly describes the (almost) linear time dependence for the area found experimentally except in the first 8 h as bacteria need to build a confluent layer.

Once the parameters of the model are found, we can then predict how a material point in the biofilm moves over time. As can be seen in figure 2*i*, this prediction is in excellent agreement with the data obtained on the displacement of the morphological features at the surface of the biofilm.

It is however important to realize that some growth also occurs in the thickness direction. While the biofilm thickness increases from ∼0.5 μm at the edge to ∼40 μm at 3 mm from the edge, the thickness saturates and becomes approximately constant further inside the biofilm. Instead of growing everywhere in thickness, the central layer grows mostly in two dimensions but long radial ridges are also formed. The growth in the thickness direction is then essentially localized along those ridges. Assuming six radial ridges (see figure 2*c*) with twice the thickness of the biofilm and a height of 1 mm, the ridges account for 10% of the overall biomass for a biofilm diameter of 2 cm.^{2} Because the standard deviation for the area measurement is large (see figure 2*e*), we have ignored those ridges in our model and treated the biofilm as a flat disc. Of course, by neglecting 10% of the biomass, we slightly underestimate the growth rate. Because the maximum growth rate *κ*_{0} has been measured independently, the values of *C*_{2} and *m* are probably slightly overestimated compared to what would be obtained by taking into account the full three-dimensional structure of the biofilm.

## 7. Kinetic roughening of the biofilm interface

Next, we investigated the formation of heterogeneities at the surface of the biofilm by analytically studying the behaviour of a three-dimensional film growing homogeneously in space. Since we focus on the formation of heterogeneities at a mesoscopic scale, we shall neglect the (macroscopic) variations of the by-product field. In that case, the growth function *g* = 1/(1 + *ρ*^{m}) is a function of time only that we write *g* = *g*^{(0)}(*t*). Because of the stochasticity in gene expression, bacterial growth rates fluctuate in space and time and the interface will not remain perfectly flat. Let us consider a bacterial film of thickness *h* in the *z-*direction and of infinite extent in the {*x*, *y*}-plane, evolving according to the set of equations (6.3) and (6.4):

*Inside*:
7.1*Outside*:
7.2The motion of the free boundary of the biofilm evolves according to the kinematic condition
7.3where ** n** is the unit vector normal to the interface. If the thickness is also independent of

*x*and

*y*, i.e. if the surface stays flat, then

*h*,

*p*and

**are functions of**

*v**z*only. Solving (7.1), (7.2) and (7.3) subject to the conditions that at

*t*= 0 the thickness is 1 and that the biofilm remains attached to the agar (i.e.

*v*(0) = 0), we arrive at the following planar solution: 7.4 7.5 7.6This solution corresponds to a flat biofilm growing perfectly homogeneously. Partly not only because of the stochasticity in gene expression, but also because of the heterogeneities of the environment, such a situation is unlikely and we may expect bacteria to divide at slightly different growth rates. In mathematical terms, this implies that the growth function

_{z}*g*is the sum of the space-independent baseline

*g*

^{(0)}(

*t*) plus a noisy quantity that depends on space and time

*ε*

*η*(

*x*,

*y*,

*t*), where

*ε*is the amplitude of the noise: 7.7where, for simplicity, we have omitted the

*z*-dependance of the noise. How does this noise affect the biofilm interface? Will it remain flat or will it lose its symmetry? To answer this question, we look for a solution of (7.1)–(7.3) that also departs from the flat solution (7.4)–(7.6) by a quantity of order

*ε*: 7.8 7.9 7.10Plugging (7.7)– (7.10) into the system (7.1)–(7.3) and developing to first order in

*ε*, we obtain the following set of linear differential equations for the fields

*p*

^{(1)},

*v*^{(1)}and

*h*

^{(1)}:

*Inside*:
7.11*Outside*:
7.12together with:
7.13We then solve (7.11) in Fourier space, again with the condition that *p*^{(1)} = 0 at the free boundary and at the base of the biofilm. Plugging the result in (7.13) yields
where is the amplitude of the ** q** component of the perturbation

*h*

^{(1)}(

*t*). Note that the operators

*ω*(

**,**

*q**t*) and

*μ*(

**,**

*q**t*) are both non-local (in real space) and time-dependent. From the deterministic version (i.e.

*η*

**= 0) of equation (7.14), we learn that any perturbation with a wavelength larger than the thickness of the biofilm (|**

_{q}**|**

*q**h*

^{(0)}(

*t*) < 1) will grow over time (

*ω*(

**,**

*q**t*) > 0) while smaller wavelengths are damped (

*ω*(

**,**

*q**t*) < 0). Thus, any finite size system would eventually reach a flat state if the growth process was allowed to continue indefinitely. This relaxational dynamics is unlikely because perturbations are continuously generated in the presence of noise (

*η*

**≠ 0) and an initially flat interface may roughen over time. Assuming a white Gaussian noise such that and , equation (7.14) has the following solution: 7.15The non-locality of the noise is critical for the development of a rough interface. Indeed, if**

_{q}*μ*(

**,**

*q**t*) = 1 and in the limiting case of a growing layer propagating at constant speed

*V*, it can be shown [51] that the surface is logarithmically rough in two dimensions and smooth in three dimensions or more. In the present case, the non-local noise does roughen the Laplacian front in both two dimensions and three dimensions. In the limit , the

*q*-component of the perturbation decays as

*ω*(

**,**

*q**t*) ∼−

*g*

^{(0)}(

*t*)

*h*

^{(0)}(

*t*)|

**| and**

*q**μ*(

**,**

*q**t*) ∼ 1/|

**| which leads to . If now we assume a general form for the time-dependent thickness**

*q**h*

^{(0)}(

*t*) ∼ 1 +

*Vt*, we obtain while in the limit . Inverting this to real space, we find the local interface width in three dimensions: 7.16

^{n}Although the noise generates perturbations at all possible wavelengths with equal probabilities, we have shown using a stability analysis that each mode of the perturbation evolves at a different growth rate. Each wavelength therefore contributes differently to the interface morphology. For a given realization of the noise, we summed these contributions to find the profile of the interface at all time. Averaging over several realizations of the noise then yielded the scaling exponents: *α*_{local} = *β* = 0.5 in agreement with the experiments (*α*_{local} = 0.6 ± 0.1 and *β* = 0.5 ± 0.1). While this result is a direct consequence of the high cohesiveness of the biomass, it is, on the other hand, weakly dependent on the exact process limiting growth. Therefore, this should hold even if growth is limited by nutrient depletion instead of by-product accumulation (albeit with different scaling exponents, because the mean height might evolve differently in time). Recalling that *n* = 0.6 ± 0.1, our model predicts *α*/*z* = *β* = 0.5 and *β* − *α*_{local}/*z* = 0.20 ± 0.05 in good agreement with the experiments.

## 8. Discussion

As a final note on the morphogenesis of *B. subtilis* biofilm, note that there is one more level of structure at the surface of the biofilm. Long radial ridges such as those seen in figure 2*d* do not result from a kinetic roughening mechanism because they appear at the surface of biofilms over time scales of minutes, much faster than the characteristic time scale of growth. Furthermore, observation of the ridges at the biofilm boundary reveals that they are hollow, at least near the outer boundary of the biofilm. This suggests that they form by a delamination process, similar to that seen in swelling gel layers [52], occurring when mechanical stresses accumulated during growth overcome the bonding strength of the bacteria with the substrate. Recent results indeed show that areas of dying cells appear to focus the surrounding growing tissues. This focusing relieves part of the compressive stresses and further promotes the nucleation of wrinkles at the surface of the biofilm [53].

The importance of mechanical forces has also been highlighted in another recent experiment performed on a slightly different system [54]. Biofilm pellicles formed at the air–water interface are less constrained by the underlying fluid and can develop a complex structure due to a buckling instability driven by the compressive stresses. Interestingly, this also suggests that the adhesion between the biofilm and its substrate (which is not included in our model) is an important parameter controlling the formation of ripples at the colony surface and further work is needed to understand this coupling. Although the development of a complex structure has been tied to the activation of specific cell fate in *B. subtilis* biofilm (such as sporulation), several studies indicate that the biofilm architecture itself may also endow the bacterial biofilm with specific abilities such as increased resistance to wetting [55] as well as enhanced liquid transport [56].

In this work, we have highlighted that matrix-producing cells remain attached by their poles after division and form long flexible bundles. Preliminary experiments also tend to indicate that those bundles can align within an external gradient of nutrients and further work is needed to characterize and understand this coupling. Because mass reorganization is strongly hindered in this case, we have developed a biomechanical model for an incompressible biomass where the motion of the cells is overdamped. Note that the biomechanical model formed by (6.3) and (6.4) corresponds to the sharp interface limit of the more general two-phase mixture theory also used in modelling tumour growth [57,58]. A model built upon this theory has also been recently developed to describe the early stages of biofilm growth (when there could be a radial gradient of ECM concentration) [41] but requires a bacterial doubling time of 150 min to fit the experimental data. However, this value is quite large and not in agreement with our independent measurement of 85 min (see figure 5*g*). While working in the constant density limit is indeed supported mostly by the observation of two-dimensional microcolonies, note that the opposite limit where the density is allowed to fluctuate is not enough to explain the time evolution of the biofilm area. Differential growth must be taken into account to reproduce the experimental observations. This effect was introduced by coupling the biomechanical model to an equation for the production and diffusion of the by-product. The resulting mathematical model correctly describes the velocity field of the biomass throughout the biofilm. While modelling biofilm growth using the theory of mixtures is an interesting alternative, it increases the number of parameters while our minimal mechanical model is parameter-free (the material parameter is lumped inside the pressure term). As one of our goals was to develop a minimal model, we did not pursue this avenue of reflexion within this paper.

When the stochasticity of the growth process is taken into account, the model predicts that the biofilm interface roughens over time. Experimental data confirm that the surface of mature biofilm is indeed rough, with self-similarity exponents close to the theoretical prediction. While our approach allows the derivation of analytical results, alternative descriptions of the biomass are possible. In particular, cellular automata [59] as well as individual-based models [60,61] have been used to model the spreading of biofilms in various experimental conditions. Such models indeed reproduce the cohesive nature of the biomass and, interestingly, also predict an increase in the global roughness of the biofilm interface [62], under conditions of nutrient limitation.

In the case where bacteria are free to move relative to each other, it is necessary to add an additional diffusive term in the equation of motion for the interface. This diffusive term typically acts to smooth the interface and leads to the formation of either weakly (logarithmically) rough or smooth air–colony interface. This might explain why colonies of microorganisms in which the microbes do not produce EPS fail to form a complex architecture.

## 9. Material and methods

### 9.1. Strain and culture conditions

Strains used in this study are the wild-type strain NCIB3610 and its derivative *amyE:: P _{yqxM}*

_{−yfp}(

*spec*). This mutant was generated by fusing the promoter

*P*of the

_{yqxM}*yqxM-sipW-tasA*operon (encoding the TasA protein) to a gene coding for a yellow fluorescent protein and inserting the construct in the

*amyE*locus of the

*B. subtilis*chromosome [8]. For routine growth,

*B. subtilis*cells were streaked from frozen stocks onto LB agar plates and inoculated overnight at 37°C. For biofilm growth, bacteria from a liquid LB culture were collected at the indicated optical density and transferred to Msgg medium [23], either liquid or fortified with agarose. When appropriate, 100 μg ml

^{−1}of spectinomycin was added to the medium.

### 9.2. Nucleation of the biofilm phenotype

Msgg medium supplemented with agarose (0.5–10%) was brought to 75°C after autoclaving. 1 ml of this mixture was filtered through a 220 nm filter and deposited on a clean glass slide within a ∼250 μm thick square frame which was immediately closed with a coverslip and allowed to solidify for 1 h. Cells from a liquid culture were collected at low optical density (OD_{600} ∼ 0.25), checked for fluorescence to ensure that bacteria had not triggered collective behaviour prior to inoculation and diluted such that 0.2 μl contained at most 10 bacteria. The coverslip was then gently removed and the agarose pad was trimmed with a scalpel to leave a ∼6 × 6 mm and ∼250 μm thick square of agar at the centre of the chamber to prevent oxygen shortage. Immediately following this procedure, the agar surface was inoculated with 0.2 μl of the diluted *B. subtilis* culture. After evaporation of the droplet (∼30 s in a fume hood), the chamber was closed with a clean coverslip and sealed to prevent dehydration. The slide was first incubated for 2 h at 30°C, then transferred under a microscope at room temperature for further observation.

### 9.3. Microscopy and image analysis

Microscopic images were acquired using a Zeiss 135 inverted microscope equipped with a 63× dry objective with phase contrast optics and appropriate filter set for fluorescence imaging. Pictures were acquired using a cooled CCD camera (Cooke) controlled with *μ* Manager and analysed with ImageJ.

### 9.4. Critical surfactin concentration

Plates with 48 wells containing liquid Msgg nutritive medium at different concentrations (1×, 0.5×, 0.1×), supplemented with increasing concentrations of surfactin (0, 0.01, 0.1, 0.5, 1, 2, 5, 10 μM) were inoculated with 1 μl of a suspension of *B. subtilis* cells grown to exponential phase (OD_{600} ∼ 0.5) in LB medium. The fluorescence in each of the wells was then recorded using a Fluoroskan plate reader.

### 9.5. Biofilm growth kinetics and inhibition

Msgg fortified with 1.5% agarose plates was dried under the fume hood for 30 min and then inoculated with 1 μl of a suspension of exponentially growing *B. subtilis* cells. The plate was then inverted and placed on a scanner. Water-filled Petri dishes and two fans helped maintain a saturated humidity in the incubator. One image was taken every 20 min for up to 5 days. pH measurements were performed using pH paper.

### 9.6. Kinetic roughening of the biofilm interface

A biofilm was grown for 24 h on Msgg at 30 ± 2°C before being imaged with a 10× objective and fluorescence optics. An array of 4 by 20 images was collected and stitched together using ImageJ software. This intensity map was then converted into a height map. Assuming a homogeneous spatial distribution of fluorescent bacteria, the emitted light by unit thickness *I*_{0} is constant. The light intensity produced by a layer of thickness d*z* at a depth *z* below the surface is *I*_{0}d*z*. The fraction of this light reaching the detector is given by the Beer–Lambert law and reads *I*_{0}e^{(−z/λ)}d*z* where *λ* is the penetration depth. Integrating this quantity between *z* = 0 and *h* where *h* is the thickness of the biofilm yields the total intensity *I* produced by a biofilm of thickness *h* reaching the detector, i.e. *I* = *I*_{0}*λ*(1 − e^{(−h/λ)}). Inverting this relation yields *h* = *λ* log{1/(1 − *I*/(*I*_{0}*λ*))}. The emitted light by unit thickness *I*_{0} was found by assuming the thickness at the biofilm edge to be 1 μm and assuming a penetration depth of λ = 235 μm [41].

## Acknowledgement

We thank the R. Kolter laboratory for the gift of the *P _{yqxM}* −

*yfp*strain and J. Merrin for sharing the scanner/incubator apparatus. We thank P. Kumar, Y. Maeda, J. Merrin and A. Petroff for discussions.

## Appendix A. Bacterial clusters need to reach a critical area

At *t* ∼ 33 h, the first event of phenotypic switching occurs within the cluster (figure 5*e*), when the population reaches approximately 3000 cells. Owing to both cellular division and additional phenotypic switching from other bacteria, the number of matrix-producing cells rose quickly from 2 to 37 cells in 3 h (figure 5*f*). Switching was not restricted to areas close to already switched bacteria neither at preferential locations in the cluster. Because matrix secretion is under the indirect control of the lipopeptide surfactin, this indicates that the bio-surfactant was homogeneously distributed within the biofilm at the onset of differentiation.

Within the time scale of our experiment, the density of cells does not significantly evolve and cells remain tightly packed. From this observation, we conclude that cellular density may not be the relevant parameter triggering collective behaviour for clustered growth. To better pinpoint the relevant factor initiating biofilm formation, we performed similar experiments with different environmental conditions (the substrate stiffness was varied by changing the agarose concentration of the nutritive gel layer in the interval 0.5–10% and the nutrient was diluted up to 10 times). Surprisingly, those experiments revealed that the emergence of the biofilm phenotype does not occur at a critical number of cells (which varies by up to a factor of 10), neither at a critical density, but at a critical *size* of the cluster of 656 ± 208 μm^{2}.

This critical size requirement can be explained by the following model. A gram of (dry) *B. subtilis* bacteria typically produces 10 mg of surfactin per hour [63]. Because bacteria are tightly packed within a cluster, this corresponds to a volumetric production rate of surfactin of *P _{s}* ∼ 0.01 mol (l)

^{−1}h

^{−1}. Assuming the surfactin molecules degrade or diffuse fast outside of the cluster (which has radius

*R*and thickness

*H*), its concentration

*s*is zero outside and there is a radial gradient of surfactin the magnitude of which is of the order of

*s*/

*R*. The resulting outward flux

*J*of surfactin molecules is then

*J*∼

*D*/

_{s}s*R*where

*D*∼ 10

_{s}^{−10}m

^{2}s

^{−1}is the diffusion coefficient of surfactin. Thus, the variation of the number of surfactin molecules in an interval d

*t*is

*πR*

^{2}

*Hs*(

*t*+ d

*t*) −

*πR*

^{2}

*Hs*(

*t*) =

*πR*

^{2}

*HP*d

_{s}*t*− 2

*πRHJ*d

*t*or in differential form A1from which we find

*s*(

*t*) =

*P*

_{s}R^{2}/2

*D*(1 − exp(−

_{s}*D*/2

_{s}t*R*

^{2})) under the assumption that

*s*(0) = 0, i.e. there is no surfactin initially. Because growth is much slower than diffusion on submillimetre length scales, the concentration

*s*quickly saturates at a value

*P*/(2

_{s}A*πD*) proportional to the cluster area

_{s}*A*. Within this simple model, a critical concentration therefore implies a critical area: A2

Since the critical surfactin concentration *s*_{crit} necessary to trigger biofilm formation was found to be ∼2 μM (see Material and methods), the corresponding critical cluster area is *A*_{crit} ∼ 400 μm^{2}, of the order of the experimental data. This finding implies that there must be at least one surfactin-producing cell present prior to the commitment of the cluster to a sedentary life. Considering that the reported proportion of surfactin-producing cells in mature biofilms is small (approximately 1 in every 1000–3000 cells [64]), it is surprising that, given such a low switching probability, we have observed a well-defined threshold for the nucleation of the biofilm phenotype. It will therefore be interesting to investigate the emergence of surfactin producers to assess whether they are overproduced in the early stages of biofilm formation and, if so, why and how.

## Footnotes

One contribution of 13 to a Theme Issue ‘Biophysics of active systems: a themed issue dedicated to the memory of Tom Duke’.

↵1 In addition, chemotactic microorganisms can sense local chemical concentration and bias their motion accordingly. This leads to an additional term in the velocity field, proportional to the local chemical gradient.

↵2 Note also that because bacteria in those ridges do not have access to nutrients, the ridges can only grow through mass transfer from the flat parts of the biofilm in contact with the solid surface.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.