Shape of the growing front of biofilms

The spatial organization of bacteria in dense biofilms is key to their collective behaviour, and understanding it will be important for medical and technological applications. Here we study the morphology of a compact biofilm that undergoes unidirectional growth, and determine the condition for the stability of the growing interface as a function of the nutrient concentration and mechanical tension. Our study suggests that transient behaviour may play an important role in shaping the structure of a biofilm.


Introduction
The stability of a uniform front to small disturbances is a framework for understanding pattern formation in many physical and biological systems [1], with a well-known example in material science being the fingering pattern formed due to supercooling of an alloy, as first characterized by Mullins and Sekerka [2]. In nonequilibrium physics, front instabilities are prevalent [3][4][5] such as classic examples in Hele-Shaw cells [6][7][8][9] and crystal growth [10]. In contrast, the self-organization and collective behaviour of living and synthetic active systems have been intensely studied in recent years [11][12][13][14][15][16], with an area of specific focus being spatial patterns generated by microbial systems [17][18][19][20][21][22][23][24][25]. One particular example in the biological sciences that is receiving much recent attention concerns the growth and spatial structure of biofilms, which are densely packed bacterial communities [26][27][28][29][30].
Bacteria have been experimentally observed to form different patterns in the form of growing colonies when cultured on agar plates at different levels of nutrient concentration [19,[31][32][33][34]. Specifically, the surface of growing colonies form circular (or flat) patterns when the nutrient concentration is high, while the patterns are fractal (or rough) when the nutrient concentration is low. The pattern formation driven by nutrient availability has been theoretically studied [35][36][37] using various models such as the Fisher-Kolmogorov equation [38,39], which combines bacterial diffusion, bacterial growth and nutrient diffusion, all in the dilute limit. Recent studies have also highlighted the importance of the mechanical interactions between the cells [40,41].
Given the wide diversity of microbial systems and their impact on both medical and natural systems, it is important to provide quantitative guidelines for instabilities that may influence the three-dimensional structure of growing biofilms. In this paper we present steps in this direction by analysing the influence of two measurable quantities, namely the nutrient concentration and the effective surface tension that results from active mechanical interactions between bacteria [42], on the instability of a planar growing front of bacteria.
The shape instabilities introduced by nutrient factors have been recently studied using numerical simulations [19,35,40,43,44], which capture the main features of patterning in the experiments. In this paper, we study the stability of the growing front of biofilms in a unidirectional planar growth using a perturbative analysis. We delineate various growth and patterning behaviours as a function of two key control parameters. With the mathematical criterion for the stability analysis, our study illustrates when and how the growing front of biofilms becomes unstable and thus provides insight for understanding microbial growth.

Description of the system
Consider the growth of a biofilm made of a single bacterial species. The scenario for culturing the system is depicted in figure 1(a), where nutrient is supplied from the top of the domain. Denote the nutrient concentration as c x y z t , , , ( ) , where x, y, z are the spatial coordinates and t the time coordinate. The density of bacteria is b 1 3 r~, with b the characteristic length of a single bacterium. Within the biofilm, nutrient is consumed at a rate k c ( ) by each cell, where k c k is a Michaelis-Menten form. This is a nonlinear form that describes crossover from a reaction-limited regime where the nutrient is abundant to a diffusion-limited regime where nutrient is scarce. While at the top layer of the growing biofilm, any of these regimes could be dominant, depletion of nutrient by every layer necessitates that at some depth there will be a crossover to the diffusion-limited regime. For the convenience of analysis, we apply k c k c 0 » ( ) everywhere as an approximation. Because a diffusion process is involved, the nutrient concentration satisfies the equation: where D denotes the diffusion coefficient, x q ( ) is the Heaviside step function, and L t x, H ( )is the biofilm surface at position x y x , = ( )and time t ( figure 1(a)). Denote the velocity of the growing front as V t x, ( ), then Supposing that N u nutrient molecules are consumed on average to make a single bacterium, we can calculate the velocity as where H is the depth of the active growing region within the biofilm ( figure 1(a)). The nutrient concentration and flux should be continuous at the biofilm surface (note that we have chosen the diffusion coefficients equal for simplicity). Nutrient diffusion is subject to boundary conditions where L H + denotes the boundary just above the biofilm surface, L Hthe boundary below the surface, and^the direction perpendicular to the surface. In the initial state, the nutrient concentration is homogenous within the culturing system, i.e. c C t 0 = = ¥ | ( figure 1(a)).

Growth of a flat front
Let us assume that the front grows in a steady state. To faciliate analysis, we approximate the system as semiinfinite in the z-direction, and then we transform to a moving reference frame, x y z t , , , ). In this case, the nutrient concentration satisfies the equation For a uniform stationary moving front, the one-dimensional description is (here k c k c Combined with the boundary conditions, the solution can be determined (see appendix A for details): where the coefficients are In addition, the front velocity has a solution of the form: The velocity in steady state is zero when the nutrient concentration C C c < ¥ ¥ , while V grows monotonically with C ¥ when above the threshold C c ¥ . Equation (9) is plotted in figure 2 for three different values of ζ.

Growth rate for deformation modes
, consider the case that a growing front is slightly perturbed from the flat geometry, with the deformation described by a height profile function h t x, ( ). For convenience, we define C  as c¢ above (+) or below (−) the biofilm surface. Then, for a small perturbation of h t x, ( ) the boundary conditions at the biofilm surface can be approximated as We can construct a general solution for equation (5) of the form , e e , ) and q a  ( ) satisfy the following equations: Combined with equation (10), one finds that A t q,  ( ) is of order h. Approximating these equations to the first order of h, we get the Fourier coefficients as As deformation is involved, there are two sources of contributions to the local front velocity v t x, source) that resists deformation in the growing front, as might be expected owing to cell-cell adhesion or the influence of type IV pili [42]. For simplicity, we use the where ν is the effective surface tension coefficient (see appendix D for details). We note that the stresses that determine the dynamics of a bacterial biofilm in the presence of extracellular matrix, such as exopolysaccharides (EPS), and surface attachment will require a more complex treatment than a simple surface tension description can provide. However, our focus here is on the growing front of a biofilm, which we can reasonably assume will not be influenced by the EPS production and surface attachment. Hence, a simple scalar description of the surface stress coming from a surface tension term should provide a sufficiently accurate starting point in our preliminary study. Furthermore, while at the length scale of single bacteria the granularity of the medium will be important and the stresses will be affected by the anisotropies of the underlying structure, such effects will be sub-dominant as we coarse-grain to the level where the interface can be described by a smooth function, which is the low q limit. Therefore, our work is consistent with this approximation since we are focused on the large length scale behaviour of the surface deformation modes.  (9)).
Developing an approximation to the first order of h in equation (14) yields , .
Combined with equations (14)- (17), in the Fourier space, we find that Combining equation (19) with equations (12) and (13), we find that the unidentified functions q a  ( ) are subject to the following restrictions: Equations (20) and (21) are in a closed form, from which we can obtain q a  ( ). Furthermore, from equation (19) it is clear that the stability of mode q in the growing front is determined by the sign of q q , the deformation mode increases with time, and thus leads to an instability. Consequently, a growing front is stable only under the condition that there is no unstable mode, i.e.
for all q. The dependence of instability on q and C ¥ is shown in figure 3 with different values of ζ and D n . we find that in stable regions (figure 3), q 2 l n peaks at q 0 = . Thus, we can obtain the stability behaviour (and thus the shape) of the growing front via the analysis of small q regimes q 0 , using perturbation analysis, we find the following asymptotic behaviour when q 0  : , there is only a single root, and we denote it ). Then, one finds that the biofilm surface is stable for C C P  ¥ ¥ (Region I, figure 1 figure 1(b)). On the other hand, if D 1 n > , the growing front is always stable since for all q (Region IV, figure 1(b)). The perturbative calculation near q=0 gives us analytical insight into the condition of instability at the largest length scale across the biofilm. However, figures 3(a) and (b) show the instability persists up to a finite threshold in q, so it will be important to examine the fastest growing mode which corresponds to the maximum growth rate in q space. Using a numerical solution of equations (20) and (21), we have calculated the overall growth rate as a function of q, as shown in figure 3(c), with the dependence of q max on the nutrient concentration shown in figure 3(d). The results show that the characteristic length scale of the growing pattern q 2 max p exhibits sensitivity to the nutrient concentration in a narrow range, and disappears when the nutrient concentration is higher than the initial growth threshold by only 30%-50%.

Transient growth behaviour
When C C c  ¥ ¥ , in the approximation that L » ¥, the growing front eventually stops (equation (9)), yet we can identify the shape of the front by analysing the transient behaviour before it stops. To study the transient behaviour, we apply numerical studies on the growth process using finite difference methods. Difference equations are converted from equations (1)-(4), and in order to make variables and parameters dimensionless, we assume H/b to be a constant and use H as the unit length in the z-direction, while we define H D 2 t º as the unit time interval. For convenience, define the following dimensionless variables: In our calculations (figure 4), the time step for update is 0.002t, and it makes little differences when we reduce the time steps. We consider the case that a biofilm grows from a thin layer (e.g L t 0 0.1 H * = = ( ) ) towards the nutrient source. As is shown in figure 4(a), the growing front speeds up at first owing to incorporation of more layers of bacterial growth until L H H » (denoted as Phase I). Then the translation speed decreases as the consumption of nutrient by the bacteria overwhelms the supply of diffused nutrient from the source z L = ( ) (denoted as Phase II). Finally, when L H approaches L, the front speed recovers since nutrient supply is increased in the region near the source z L = ( ) (denoted as Phase III). To study the transient behaviour before the growing front stops, we ignore the growth process of Phase I by setting L t 0 1.
). The evolution of the nutrient profile at representative times is shown in figure 4(b), where the quasi-steady state nutrient profile ( figure 4(b), blue dot line) agrees with equation (7). To analyse the behaviour of L  ¥, we simulate the growth process with different system sizes L. When C C c > ¥ ¥ , V min , defined as the minimum value in the velocity profile, approaches the analytical results calculated from equation (9) as system size L increases ( figure 4(d)). When C C c < ¥ ¥ , however, V min decreases inversely with system size (figure 4(c) and (d)), i.e. V L 1 min~( figure 4(d)). We find that , which agrees with the analytical results. In fact, when the velocity is below a threshold, the moving front actually stops (see appendix A for details), and thus there is no Phase III when L is large.
We next focus on the evolution process of Phase II. The stability of the growing front is determined by the local nutrient concentration around the biofilm surface. In the initial state, c x y z C , , = ¥ ( ) and L t 0 1.1 H * = = ( ) , so that the local nutrient concentration around the growing front is much higher than that of the steady state (equation (7)). When t is large, we find the velocity V decreases with time as V t 0.55 especially if L is large (e.g. L 250 * = in figure 4(c)). This result appears to be in good agreement with the experimental studies of biofilm growth reported in [45]. As an approximation, we use the time dependent velocity V t ( ) in this case as a quasi-steady state quantity to measure the local nutrient adequacy around the growing front. From equation (9) (and figure 2), one finds that there is a bijective mapping relation between V and C ¥ when V 0 > , and thus we obtain V V C P P º ¥ ( )and V V C c c º ¥ ( ), the mapping velocity of C P ¥ and C c ¥ in a quasi-steady state. As V t ( ) measures the local nutrient adequacy, so the growing front is stable ) before it stops, thus the biofilm surface is flat (Region V, figure 1(b)). However, if D 1 ) and the velocity decreases with time in Phase II, then we can find times t . Consequently, the growing front is unstable before it stops, resulting in a rough surface (Region III, figure 1(b)).

Discussion
The growth and patterns formed by a biofilm are summarized in table 1. According to the mathematical model developed here, in theory, there are five distinct regions. If D 1 n > , the biofilm surface is always flat (Region IV and V), yet the growth is transient when C C c < ¥ ¥ (Region V), while sustainable when C C c > ¥ ¥ (Region IV). In reality, the value of D n is usually significantly smaller than one (see appendix E for details), so these regions (Region IV and V) may be difficult to observe in experiments. If D 1 n < , the sustainable growth threshold is determined by the nutrient threshold C c ¥ , while the pattern formation is governed by the nutrient threshold C P ¥ (C C P c > ¥ ¥ ). These two nutrient thresholds, obtained naturally from our analytical analysis, can illustrate the origin of thresholds for roughness and branching in the colony patterns of a recent simulation study [40]. Furthermore, the patterning in Regions I-III agree well with those of microbial colonies in the experimental studies [19,[31][32][33][34]. To summarize, our analysis agrees qualitatively with experimental studies concerning patterning of microbial colonies and can illustrate puzzles that are not fully understood in the previous simulation studies. Our study will be relevant to a wide variety of practical questions such as the behaviour of multispecies biofilms, the impact of digestive enzymes that may free nutrients, and the influence of cooperation among cells or the presence of cheater cells on the evolution of a biofilm. Moreover, it is inherently related to the recent stability analysis that has been used to study the chemically driven growth and division of droplets [46] and may shed light on our understanding of division of proto-cells in early forms of life [46,47].

Acknowledgments
This work was supported by the Human Frontier Science Program RGP0061/2013. HAS thanks NSF-MCB-1344191. We thank Berenike Maier and Tom Cronenberg for insightful and helpful discussions, and Carey Nadell, whose seminar at Princeton on the subject during a visit by RG initiated our discussions and inspired our work.

Appendix A. Nutrient threshold for growth
Bacteria are living systems out of equilibrium. When the living environment is harsh, for instance, when nutrients are insufficient, some species of bacteria switch to a protective state named a spore, which is quasi inanimate with significantly lowered energy dissipation (so that the bacteria may survive for even hundreds of years in the harsh environment) [48][49][50]. There exists a minimum threshold of nutrient flux to initiate cell growth, and we can define it as ε. To consider this effect, the formula of front velocity changes as: ) . Due to the consumption by microbes, the nutrient is depleted very quickly when the concentration is below certain levels [51]. Therefore, we can consider ε to be small, and neglect it when there are enough nutrients to sustain bacterial growth. However, ε needs to be taken into consideration when the moving speed approaches to zero. For convenience, we can approximate (A.1) as Using (A.2), we find that there is a threshold of velocity for biofilms growth: when V H N u e < , V=0. Thus, when C C c < ¥ ¥ , as we find in figures 4(c) and (d), V min depends inversely on the system size L. When L is sufficiently large, so that V , then there is no Phase III in the growth process. Furthermore, we can derive a solution to the nutrient profile in the steady state (t  ¥) when C C c < ¥ ¥ (note that V t 0  ¥ = ( ) ). Supposing that the growing front finally stops at L t L L H  ¥ = -¢ ( ) , and using coordinates x y z t x y z L L t , , , , , , ¢ ¢ ¢ ¢ = -+ ¢ ( ) ( ), we obtain the following equation as t¢  ¥: . The boundary conditions are: Table 1. Behaviour of the growing front of biofilm.

Region
Definition Behaviour In steady state V=0, which means (see (A.2)) ). When L is large, and 1 L L L -¢ ¢  , we can approximate L¢ as L and thus obtain the nutrient profile in equation (7). When is L large enough to be approximated as ¥? The criterion lies in (A.5), which requires a threshold value of the system size for the approximation of L  ¥. Specifically, if we take L 0 as the threshold, then where L 0 exhibits an inverse dependence on ε. When L L 0 > , we can approximate L as ¥.
demands that the difference in the normal stress (or pressure) across the interface is balanced by the local mean curvature of the surface times the surface tension; a combination that is called the Laplace pressure. In the nonequilibrium case of a growing biofilm front, while we do not have access to the energetic definition of surface tension, we can use a generalization of the force balance argument to define a quantity that plays the same role as surface tension. Microscopically, we can use the Laplace pressure equilibration argument to define surface tension as the tendency to eliminate any curvature in the absence of normal stress difference. In the case of the biofilm, this mechanisms can be provided by the pulling of the pili: since they are attached to the neighbouring bacteria, any curvature in the profile necessitates an asymmetry in the distribution of pulling pili that will relax the bacterial configuration to a flat surface where there is minimal asymmetry. This crude generalization can be used to define an interfacial normal stress difference g that is proportional to the local curvature h

Appendix E. Parameters estimation and discussions
Two parameters are most important in our model; one is C c ¥ , the other is D n . From equation (8a), C c ¥ is of order N ; u r for common carbon sources (e.g. glucose), N u is of the order of 10 9 [52,53]. While ρ is typically b 3 -. Considering the spaces between bacteria in a biofilm, ρ can be estimated as