Shape of the Growing Front of Biofilms

TL;DR

This paper investigates the morphology and stability of growing biofilm fronts, highlighting how nutrient levels and mechanical tension influence biofilm structure and transient behaviors.

Contribution

It provides a theoretical analysis of the stability conditions of biofilm growth interfaces considering nutrient and tension effects.

Findings

Stability depends on nutrient concentration and mechanical tension.

Transient behaviors significantly influence biofilm structure.

Conditions for stable and unstable growth fronts are identified.

Abstract

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.

Tables

Table 1. Table 1: Behavior of the growing front of biofilm.

Region	Definition	Behavior
I	$C_{\infty }>C_{\infty }^P$,	$V>0$, Flat
II	$C_{\infty }^P>C_{\infty }>C_{\infty }^c$,	$V>0$, Rough
III	$C_{\infty }^c>C_{\infty }$,	$V=0$, Rough
IV	$C_{\infty }>C_{\infty }^c$, $\nu /D>1$	$V>0$, Flat
V	, $\nu /D>1$	$V=0$, Flat

Equations

$${\partial }_{t}c-D{\nabla }^{2}c+\rho k(c)\theta (-z+L_H(\mathbf{x},t))=0,$$

$$L_H(\mathbf{x},t)={\int }_0^{t}d{t}^{\prime }V(\mathbf{x},t^{\prime })+L_H(\mathbf{x},0).$$

$$V(\mathbf{x},t)=\frac{1}{N_u}\underset{(L_H(\mathbf{x},t)-H)\theta (L_H(\mathbf{x},t)-H)}{\overset{L_H(\mathbf{x},t)}{\int }}dzk(c),$$

$$\left\{\begin{array}[]{l}c\left|{{}_{z=L}}\right.={C_{\infty}},\\ c\left|{{}_{z={L_{H}}^{+}}}\right.=c\left|{{}_{z={L_{H}}^{-}}},\right.\\ {\partial_{z}}c\left|{{}_{z={L_{H}}^{+}}}\right.={\partial_{z}}c\left|{{}_{z={L_{H}}^{-}}},\right.\end{array}\right.$$

$${\partial }_{t^{\prime }}{c}^{\prime }-V{\partial }_{z^{\prime }}{c}^{\prime }-D{\nabla }^2{c}^{\prime }+\rho k(c^{\prime })\theta (-z^{\prime })=0.$$

$$-V\frac{d{c}^{\prime }}{d{z}^{\prime }}-D\frac{d^2{c}^{\prime }}{dz{}^{\prime 2}}+\rho k_0{c}^{\prime }\theta (-z^{\prime })=0.$$

$$\begin{array}[]{lcr}c^{\prime}(z^{\prime})=\left\{\begin{array}[]{ll}{C_{\infty}}\left(1+B\,{e^{{{-Vz^{\prime}}\mathord{\left/{\vphantom{{-Vz^{\prime}}D}}\right.\kern-1.2pt}D}}}\right)&\;(z^{\prime}\geq 0)\\ {C_{\infty}}\left(1+B\right)\,{e^{{\lambda_{1}}z^{\prime}}}&\;(z^{\prime}<0)\end{array}\right.\;;&{C_{\infty}}\geq C_{\infty}^{c},\\ \texttt{{and}}\\ c^{\prime}(z^{\prime})=\left\{\begin{array}[]{ll}{C_{\infty}}\,\frac{{z^{\prime}}}{L}&\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(z^{\prime}\geq 0)\\ 0&\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(z^{\prime}<0)\\ \end{array}\right.\;;&{C_{\infty}}\leq C_{\infty}^{c}.\\ \end{array}$$

$$C_{\infty }^c=\frac{N_u\rho }{1-\exp (-\sqrt{\zeta })}$$

$$B=-\frac{1}{\zeta }{\ln }^2(1-N_u\rho /C_{\infty })$$

$${\lambda }_1=-\frac{1}{H}\ln (1-N_u\rho /C_{\infty })$$

$$\zeta =H^2\rho k_0/D$$

$$\displaystyle V=\left\{\begin{array}[]{lr}\frac{D}{H}\left[{\ln\left({1-\frac{{{N_{u}}\rho}}{{{C_{\infty}}}}}\right)-\frac{\varsigma}{{\ln\left({1-\frac{{{N_{u}}\rho}}{{{C_{\infty}}}}}\right)}}}\right]&;\;{C_{\infty}}\geq C_{\infty}^{c},\\ \\ 0&;\;{C_{\infty}}\leq C_{\infty}^{c}.\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}{C^{+}}\left|{{}_{z^{\prime}=h({\mathbf{x},t})}}\right.={C^{-}}\left|{{}_{z^{\prime}=h({\mathbf{x},t})}},\right.\\ {\partial_{z^{\prime}}}{C^{+}}\left|{{}_{z^{\prime}=h({\mathbf{x},t})}}\right.={\partial_{z^{\prime}}}{C^{-}}\left|{{}_{z^{\prime}=h({\mathbf{x},t})}}.\right.\end{array}\right.$$

$$\left\{\begin{array}[]{l}{C^{\rm{+}}}=C_{0}^{\rm{+}}+\int_{\mathbf{q}}{{A_{+}}({\mathbf{q},t})\;{{\mathop{\rm e}\nolimits}^{i\mathbf{q}\cdot\mathbf{x}}}}\;{{\mathop{\rm e}\nolimits}^{-{\alpha_{+}}(\mathbf{q})z}},\\ {C^{-}}=C_{0}^{-}+\int_{\mathbf{q}}{{A_{-}}({\mathbf{q},t})\;{{\mathop{\rm e}\nolimits}^{i\mathbf{q}\cdot\mathbf{x}}}}\;{{\mathop{\rm e}\nolimits}^{{\alpha_{-}}(\mathbf{q})z}},\end{array}\right.$$

$$\left\{\begin{array}[]{l}{\partial_{t}}{A_{+}}+({V{\alpha_{+}}-D{\alpha_{+}}^{2}+D{q^{2}}})\;{A_{+}}=0,\\ {\partial_{t}}{A_{-}}-({V{\alpha_{-}}+D{\alpha_{-}}^{2}-D{q^{2}}-\rho{k_{0}}})\;{A_{-}}=0.\end{array}\right.$$

$$A_{+}=A_{-}=-\frac{C_{\infty }\left(1+B\right)\rho k_0/D}{{\alpha }_{+}+{\alpha }_{-}}h(\mathbf{q},t).$$

$$v(\mathbf{x},t)=v_J(\mathbf{x},t)+v_b(\mathbf{x},t),$$

$$v_J(\mathbf{x},t)=\frac{1}{N_u\rho }\underset{h(\mathbf{x},t)-H}{\overset{h(\mathbf{x},t)}{\int }}d{z}^{\prime }(V{\partial }_{z^{\prime }}{C}^{-}+D{\partial }_{z^{\prime }}^2{C}^{-}),$$

$$v_b(\mathbf{x},t)=\nu {\nabla }^{2}h(\mathbf{x},t),$$

$$\begin{array}[]{l}{v_{J}}\simeq{V}\left[{1+\lambda_{1}h({\mathbf{x},t})}\right]+\frac{V}{{{N_{u}}\rho}}\int_{\mathbf{q}}{{A_{-}}\;{{\mathop{\rm e}\nolimits}^{i\mathbf{q}\cdot\mathbf{x}}}}({1-{{\mathop{\rm e}\nolimits}^{-{\alpha_{-}}H}}})\\ \;\;\;\;\;\;\;\;\;+\frac{D}{{{N_{u}}\rho}}\int_{\mathbf{q}}{{A_{-}}\;{{\mathop{\rm e}\nolimits}^{i\mathbf{q}\cdot\mathbf{x}}}}({1-{{\mathop{\rm e}\nolimits}^{-{\alpha_{-}}H}}})\;{\alpha_{-}}.\end{array}$$

$$v(\mathbf{x},t)=V+{\partial }_{t}h(\mathbf{x},t).$$

$${\partial }_{t}h(\mathbf{q},t)=\left[\lambda (\mathbf{q})-\nu q^2\right]h(\mathbf{q},t),$$

$$\lambda (\mathbf{q})=V{\lambda }_1-\frac{V{\lambda }_1\left(V+D{\alpha }_{-}\right)}{D({\alpha }_{+}+{\alpha }_{-})}\frac{(1-e^{-{\alpha }_{-}H})}{(1-e^{-{\lambda }_{1}H})}.$$

$$\lambda -\nu q^2=D{{\alpha }_{+}}^2-V{\alpha }_{+}-D{q}^2,$$

$$\lambda -\nu q^2=D{{\alpha }_{-}}^2+V{\alpha }_{-}-D{q}^2-\rho k_0.$$

$$\lambda (\mathbf{q})-\nu q^2\approx \left[\frac{1-\nu /D}{1/f(C_{\infty })-1}-\nu /D\right]D{q}^2,$$

$$f(C_{\infty })\equiv {{\lambda }_1}^2\zeta /H^2-H{\lambda }_1\left(\frac{2}{1+{{\lambda }_1}^2\zeta /H^2}-1\right)\left(\frac{C_{\infty }}{N_u\rho }-1\right).$$

$$V=\frac{1}{N_u}{\int }_{L_H-H}^{L_H}dz\left[k_{0}c(z)-\epsilon \right]\theta \left(k_{0}c(z)-\epsilon \right),$$

$$V=\theta \left(\frac{1}{H}{\int }_{L_H-H}^{L_H}dz{k}_{0}c(z)-\epsilon \right)\frac{1}{N_u}{\int }_{L_H-H}^{L_H}dz{k}_{0}c(z),$$

$$-D{\partial }_{z^{\prime }}^2{c}^{\prime }\left(z^{\prime }\right)+\rho k_0{c}^{\prime }\left(z^{\prime }\right)\theta \left(-z^{\prime }\right)=0.$$

$$\left\{\begin{array}[]{l}c^{\prime}(z^{\prime}=L^{\prime})={C_{\infty}}\\ {C^{+}}(z^{\prime}=0)={C^{-}}(z^{\prime}=0)\\ {\partial_{z^{\prime}}}{C^{+}}(z^{\prime}=0)={\partial_{z^{\prime}}}{C^{-}}(z^{\prime}=0)\end{array}\right..$$

$$\frac{1}{H}{\int }_{-H}^{0}d{z}^{\prime }{k}_0{c}^{\prime }\left(z^{\prime }\right)<\epsilon .$$

Full text

Keywords: biofilm, instability, pattern formation

1 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 contrast, the self-organization and collective behaviour of living and synthetic active systems have been intensely studied in recent years [3, 4, 5, 6, 7, 8], with an area of specific focus being spatial patterns generated by microbial systems [9, 10, 11, 12, 13, 14, 15, 16, 17]. 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 [18, 19, 20, 21, 22].

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 [11, 23, 24, 25, 26]. 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 [27, 28, 29, 30, 31] using various models such as the Fisher-Kolmogorov equation [32, 33], 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 [34, 35].

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 analyzing the influence of two measurable quantities, namely the nutrient concentration and the effective surface tension that results from active mechanical interactions between bacteria [36], on the instability of a planar growing front of bacteria.

The shape instabilities introduced by nutrient factors have been recently studied using numerical simulations [11, 29, 34, 37, 38], 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 criteria for the stability analysis, our study illustrates when and how the growing front of biofilms becomes unstable. Our analysis agrees 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 [39] and may shed light on our understanding of division of proto-cells in early forms of life [39, 40].

2 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 $\rho=1/{b^{3}}$, 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^{\prime}\frac{c}{{c+{K_{m}}}}$ 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)\approx{k_{0}}c$ everywhere as an approximation. Because a diffusion process is involved, the nutrient concentration satisfies the equation:

[TABLE]

where $D$ denotes the diffusion coefficient, $\theta(x)$ is the Heaviside step function, and ${L_{H}}({\mathbf{x},t})$ is the biofilm surface at position $\mathbf{x}=(x,y)$ and time $t$ (Figure 1(a)).

Denote the velocity of the growing front as ${V}({\mathbf{x},t})$, then

[TABLE]

Supposing that ${N_{u}}$ nutrient molecules are consumed on average to make a single bacterium, we can calculate the velocity as

[TABLE]

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. Nutrient diffusion is subject to boundary conditions

[TABLE]

where ${L_{H}}^{+}$ denotes the boundary just above the biofilm surface and ${L_{H}}^{-}$ the boundary below the surface. In the initial state, the nutrient concentration is homogenous within the culturing system, i.e. $c\left|{{}_{t=0}}\right.={C_{\infty}}$ (Figure 1(a)).

To faciliate analysis, we approximate the system as semi-infinite in the $z$-direction, and then we transform to a moving reference frame, $(x^{\prime},y^{\prime},z^{\prime},t^{\prime})$, with $(x^{\prime},y^{\prime},z^{\prime},t^{\prime})=(x,y,z-{L_{H}},t)$ and $c^{\prime}(x^{\prime},y^{\prime},z^{\prime},t^{\prime})=c(x,y,z,t)$. In this case, the nutrient concentration satisfies the equation

[TABLE]

3 Growth of a Flat Front

Let us assume that the front grows in a steady state. For a uniform stationary moving front, the one-dimensional description is (here $k(c)\approx{k_{0}}c$)

[TABLE]

Combined with the boundary conditions, the solution can be determined (see appendix A for details):

[TABLE]

where the coefficients are

[TABLE]

Meanwhile, the front velocity has a solution of the form:

[TABLE]

The velocity in steady state is zero when the nutrient concentration $C_{\infty}<C_{\infty}^{c}$, while $V$ and $C_{\infty}$ is bijective when above the threshold $C_{\infty}^{c}$. Equation (9) is plotted in Figure 2 for three different values of $\zeta$.

4 Growth Rate for Deformation Modes

When $C_{\infty}>C_{\infty}^{c}$, and thus $V(t\to\infty)>0$ , consider the case that a growing front is slightly perturbed from the flat geometry, with the deformation described by a height profile function $h(\mathbf{x},t)$. For convenience, we define ${C^{\rm{\pm}}}$ as $c^{\prime}$ above (+) or below (-) the biofilm surface. Then, the boundary conditions at the biofilm surface are

[TABLE]

We can construct a general solution for equation (5) of the form

[TABLE]

where $C_{0}^{+}\equiv{C_{\infty}}\left(1+B\,{e^{-Vz^{\prime}/D}}\right)$, $C_{0}^{-}\equiv{C_{\infty}}\left(1+B\right)\,{e^{{\lambda_{1}}z^{\prime}}}$, and $\int_{\mathbf{q}}{\equiv\int{\frac{{{d^{2}}{\mathbf{q}}}}{{{{({2\pi})}^{2}}}}}}$. By substituting equation (11) into equation (5), we find ${A_{\pm}}({\mathbf{q},t})$ and ${\alpha_{\pm}}(\mathbf{q})$ satisfy the following equations:

[TABLE]

Combined with equation (10), one finds that ${A_{\pm}}({\mathbf{q},t})$ is of order $h$. Approximating these equations to the first order of $h$, we get the Fourier coefficients as

[TABLE]

As deformation is involved, there are two sources of contributions to the local front velocity $v({\mathbf{x},t})$,

[TABLE]

where $v_{J}({\mathbf{x},t})$ represents the biofilm growth caused by nutrient flux,

[TABLE]

while $v_{b}({\mathbf{x},t})$ is a surface related contribution (a new source) that resists deformation in the growing front, as might be expected owing to cell-cell adhesion or the influence of type IV pili [36]. For simplicity, we use the following generic form

[TABLE]

where $\nu$ is the effective surface tension coefficient. Developing an approximation to the first order of $h$ in equation (14) yields

[TABLE]

Meanwhile ${L_{H}}({\mathbf{x},t})=\overline{{L_{H}}(t)}+h({\mathbf{x},t})$, where $\overline{{L_{H}}(t)}$ denotes the average over $x$ and $y$ axis in ${L_{H}}({\mathbf{x},t})$. By applying ${\partial_{t}}$ on both sides:

[TABLE]

Combined with equations (14)-(17), in the Fourier space, we find that

[TABLE]

where

[TABLE]

Combining equation (19) with equations (12) and (13), we find that the unidentified functions ${\alpha_{\pm}}(\mathbf{q})$ are subject to the following restrictions:

[TABLE]

Equations (20) and (21) are in a closed form, from which we can obtain ${\alpha_{\pm}}(\mathbf{q})$. Furthermore, from equation (19) it is clear that the stability of mode $q$ in the growing front is determined by the sign of ${\lambda(\mathbf{q})-\nu{q^{2}}}$: when ${\lambda(\mathbf{q})-\nu{q^{2}}}>0$, 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. ${\lambda(\mathbf{q})-\nu{q^{2}}}\leq 0$ for all $q$. The dependence of instability on $q$ and $C_{\infty}$ is shown in Figure 3 with different values of $\zeta$ and ${\nu\mathord{\left/{\vphantom{\nu D}}\right.\kern-1.2pt}D}$. One can find that in stable regions (Figure 3), $\lambda-\nu{q^{2}}$ peaks at ${{q}\approx 0}$. Thus, we can obtain the stability behavior (and thus the shape) of the growing front via the analysis of small $q$ regimes $(q\approx 0)$.

Since ${\alpha_{+}(q=0)}=V/D$, ${\alpha_{-}(q=0)}=\lambda_{1}$ and $\lambda(q=0)=0$, using perturbation analysis, we find the following asymptotic behavior when $q\approx 0$:

[TABLE]

where

[TABLE]

When $C_{\infty}\geq C_{\infty}^{c}$, $f({{C_{\infty}}})$ is a monotonically decreasing function of $C_{\infty}$, with $f({C_{\infty}}=C_{\infty}^{c})=1$ and $f({C_{\infty}}\to\infty)=-1$. So if $\nu/{D}<1$, for equation $f({C_{\infty}}){\rm{=}}\nu/{D}$, there is only a single root, and we can denote it as $C_{\infty}^{P}$. Then, one finds that the biofilm surface is stable for $C_{\infty}\geq C_{\infty}^{P}$ (Region I, Figure 1(b)), while it is unstable for $C_{\infty}^{c}<{C_{\infty}}<C_{\infty}^{P}$ (Region II, Figure 1(b)). On the other hand, if $\nu/{D}>1$, the growing front is always stable since $\lambda(\mathbf{q})-\nu{q^{2}}<0$ for all $q$ (Region IV, Figure 1(b)).

The perturbative calculation near $q=0$ gives us an 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_{\rm max}$ on the nutrient concentration shown in Figure 3(d). It is intriguing to find out that the characteristic length scale of the growing pattern $2\pi/q_{\rm max}$ exhibits such extreme sensitivity to the nutrient concentration in a narrow range, and disappears when the nutrient concentration is higher than the initial growth threshold by only only 30%-50% in our typical examples.

5 Transient growth behaviour

When $C_{\infty}\leq C_{\infty}^{c}$, in the approximation that $L\approx\infty$, the growing front eventually stops (equation (9)), yet we can identify the shape of the front by analyzing the transient behavior before it stops. To study the transient behavior, we apply numerical studies on the growth process using difference equations converted from equations (1) and (4). In the simulation (Figure 4), we assume $H/b$ to be a constant and use ${H}$ as the unit length in the $z$-direction, while we define $\tau\equiv{H^{2}}/D$ as the unit time interval. For convenience, define the following dimensionless variables: ${L^{\rm{*}}}\equiv L/H$, ${L_{H}}^{*}\equiv L_{H}/H$, ${t^{*}}\equiv t/\tau$, ${V^{*}}\equiv V\cdot\tau/H$, ${c}^{*}\equiv\frac{{{c}}}{{N_{u}}{\rho}}$ and ${C_{\infty}}^{*}\equiv\frac{{{C_{\infty}}}}{{N_{u}}{\rho}}$.

We consider the case that a biofilm grows from a thin layer (e.g ${{{L_{H}}^{*}(t=0)}=0.1}$) 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}\approx 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 behavior before the growing front stops, we ignore the growth process of Phase I by setting ${{{L_{H}}^{*}(t=0)}=1.1}$ (so that $L_{H}(t=0)>H$). The evolution of the nutrient profile at representative moments is shown in Figure 4(b), where the quasi-steady state nutrient profile (Figure 4(b), blue dot line) agrees with equation (7). To analyze the behavior of $L\to\infty$, we simulate the growth process with different system sizes $L$. ${V_{\min}}$, defined as the minimum value in the velocity profile, decreases inversely with system size (Figure 4(c-d)), i.e. ${V_{\min}}\sim 1/L$ (Figure 4(d)). 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_{\infty}$ and ${{{L_{H}}^{*}(t=0)}=1.1}$, the local nutrient concentration around the growing front is much higher than that of the steady state (equation (7)). 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_{\infty}$ when $V>0$, and thus we obtain ${V^{P}}\equiv V\left({C_{\infty}^{P}}\right)$ and ${V^{c}}\equiv V\left({C_{\infty}^{c}}\right)$, the mapping velocity of ${C_{\infty}^{P}}$ and ${C_{\infty}^{c}}$ in a quasi-steady state. As $V(t)$ measures the local nutrient adequacy, so the growing front is stable when $V(t)>V^{P}$, while unstable when $V^{P}>V(t)>V^{c}$. If $\nu/D>1$, the growing front is always stable ($\forall q$, $\lambda(\mathbf{q})-\nu{q^{2}}<0$) before it stops, thus the biofilm surface is flat (Region V, Figure 1(b)). However, if $\nu/D<1$, ${V^{P}}=V\left({C_{\infty}^{P}}\right)>0$, since $V(t\to\infty)=0$ and velocity decreases with time in Phase II, then we can find times $t^{P}>0$ satisfying $V(t)<V^{P}$ when $t>t^{P}$. Meanwhile ${V^{c}}=V\left({C_{\infty}^{c}}\right)=0$, thus when $t>t^{P}$, $V^{P}>V(t)>V^{c}$. Consequently, the growing front is unstable before it stops, resulting in a rough surface (Region III, Figure 1(b)).

6 Discussion

The growth and patterns formed by a biofilm is summarized in Table 1. According to the mathematical model developed here, in theory, there are five distinct regions. If $\nu/D>1$, the biofilm surface is always flat (Region IV and V), yet the growth is transient when ${C_{\infty}}<C_{\infty}^{c}$ (Region V), while sustainable when ${C_{\infty}}>C_{\infty}^{c}$ (Region IV). In reality, the value of $\nu/D$ is usually significantly smaller than one, so these regions (Region IV and V) may be difficult to observe in experiments. If $\nu/D<1$, the sustainable growth threshold is determined by the nutrient threshold $C_{\infty}^{c}$, while the pattern formation is governed by the nutrient threshold $C_{\infty}^{P}$ ($C_{\infty}^{P}>C_{\infty}^{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 [34]. Furthermore, the patterning in Regions I-III agree well with those of microbial colonies in the experimental studies [11, 23, 24, 25, 26].

To summarize, our study provides a mathematical framework to differentiate growth and patterning behaviors of biofilms and so provides insight for understanding microbial growth.

This work was supported by the Human Frontier Science Program RGP0061/2013. We thank Berenike Maier and Tom Cronenberg for helpful discussions.

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) [41, 42, 43]. Definitely, there is a minimum threshold of nutrient flux to initiate cell growth, and we can define it as $\varepsilon$. To consider this effect, the formula of front velocity changes as:

[TABLE]

where we have applied $L_{H}(t=0)>H$. $\varepsilon$ is small and can be neglected when there are enough nutrients to sustain bacterial growth, while it needs to be taken into consideration when the moving speed approaches to zero. For convenience, we can approximate (A.1) as

[TABLE]

Using (A.2), we find that there is a threshold of velocity for biofilms growth: when $V<\frac{H}{{{N_{u}}}}\varepsilon$, $V=0$. Thus, when ${C_{\infty}}<C_{\infty}^{c}$, as we find in Figure 4(c-d), ${V_{\min}}$ depends inversely on the system size $L$. When $L$ is sufficiently large, so that ${V_{\min}}<\frac{H}{{{N_{u}}}}\varepsilon$, 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\to\infty$) when ${C_{\infty}}<C_{\infty}^{c}$ (note that $V(t\to\infty)=0$). Supposing that the growing front finally stops at $L_{H}(t\to\infty)=L-L^{\prime}$, and using coordinates $(x^{\prime},y^{\prime},z^{\prime},t^{\prime})=(x,y,z-L+L^{\prime},t)$, we obtain the following equation as $t^{\prime}\to\infty$:

[TABLE]

For $z^{\prime}>0$, ${C^{+}}(z^{\prime})={C_{\infty}}\left({A_{c}}+{B_{c}}{\lambda_{1}}z\right)^{\prime}$; whereas $z^{\prime}<0$, ${C^{-}}(z^{\prime})={C_{\infty}}{C_{1}}\,\exp\left(z^{\prime}\sqrt{\frac{{\rho{k_{0}}}}{D}}\right)$. The boundary conditions are:

[TABLE]

In steady state $V=0$, which means [see (A.2)]

[TABLE]

Combined with (A.4), we obtain ${C_{1}}=A_{c}=B_{c}=1/(1+{\lambda_{1}}L^{\prime})$. When $L$ is large, and $\frac{{L-L^{\prime}}}{{L^{\prime}}}\ll 1$, we can approximate $L^{\prime}$ as $L$ and thus obtain the nutrient profile in equation (7).

When is $L$ large enough to be approximated as $\infty$? The criterion lies in (A.5), which requires a threshold value of the system size for the approximation of $L\to\infty$. Specifically, if we take $L_{0}$ as the threshold, then

[TABLE]

where $L_{0}$ exhibits an inverse dependence on $\varepsilon$. When $L>L_{0}$, we can approximate $L$ as $\infty$.

Appendix B Upper bound for the growth speed

We assumed that the uptake rate of nutrient by bacteria takes the Michaelis-Menten form, i.e.

[TABLE]

with the maximum value of $k(c)$ to be $k^{\prime}$. From equation (3), we find that there is an upper bound to the translation speed, i.e.

[TABLE]

Despite the existence of an upper bound to velocity, it is easy to find that the bijective mapping relation between $C_{\infty}$ and $V(t\to\infty)$ still holds for $C_{\infty}\geq C_{\infty}^{c}$. Using the same analysis, we can get similar growth and patterning behavior of a biofilm when we apply the correction that arises due to Equation (B.2).