PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of procaThe Royal Society PublishingProceedings AAboutBrowse by SubjectAlertsFree Trial
 
Proc Math Phys Eng Sci. 2016 April; 472(2188): 20150798.
PMCID: PMC4892277

The development of biofilm architecture

Abstract

We extend the one-dimensional polymer solution theory of bacterial biofilm growth described by Winstanley et al. (2011 Proc. R. Soc. A 467, 1449–1467 (doi:10.1098/rspa.2010.0327)) to deal with the problem of the growth of a patch of biofilm in more than one lateral dimension. The extension is non-trivial, as it requires consideration of the rheology of the polymer phase. We use a novel asymptotic technique to reduce the model to a free-boundary problem governed by the equations of Stokes flow with non-standard boundary conditions. We then consider the stability of laterally uniform biofilm growth, and show that the model predicts spatial instability; this is confirmed by a direct numerical solution of the governing equations. The instability results in cusp formation at the biofilm surface and provides an explanation for the common observation of patterned biofilm architectures.

Keywords: biofilm architecture, instability, biofilm structure, biofilm pattern

1. Introduction

Bacterial biofilms are virtually ubiquitous. In a biofilm, bacterial cells attach to one another—and usually to a surface or interface—with a slime matrix which largely immobilizes the cells and creates a distinct micro-environment. Combined with associated phenotypic changes, this provides a measure of protection from physical and chemical attack which allows biofilm to persist in even the most extreme environments [1]. This resistance also makes biofilm a major consideration in many applications of human interest, with both positive effects (biogeochemical cycles, gut flora, wastewater treatment, bioreactors) and negative (medical devices, pathology, industrial fouling and corrosion).

While early mathematical models of biofilm implicitly assumed that it forms a uniform, smooth layer (e.g. [25]), improved microscopic techniques in recent decades have revealed a rich diversity of biofilm structure varying from smooth flat biofilm to forms described as towers, mushrooms, streamers, pores and channels (e.g. [6]). Repeatability is notoriously challenging in biofilm experiments, but results (reviewed by Stoodley et al. [7]) broadly suggest that such architectural forms are favoured under scarce nutrient conditions and relatively low fluid flow. By contrast, high fluid flow and/or high nutrient conditions favour relatively flat biofilm.

The need to understand biofilm architecture stems from the interdependence of the surface architecture, transport of dissolved species (including nutrients, metabolites, disinfectants) across the biofilm/fluid interface, growth rate and fluid flow. For a given surface density of biofilm biomass, the overall chemical exchange, biofilm permeability, flow resistance in confined flow paths, cell attachment/detachment and sloughing may all be presumed to depend on the surface architecture and may affect a range of practical applications. Detachment and sloughing further affect the spread of biofilm and clogging of downstream flow paths.

Many modelling studies of biofilm growth have used discrete approaches to describe the biomass, including cellular automata and individual-based models (reviewed by Laspidou et al. [8], for example). These provide a useful means to perform simulations that readily allow consideration of complex scenarios such as multi-species biofilm, multi-component biochemistry and interactions with fluid flow. However, the physical properties of biofilm are dominated by the slime matrix rather than the cellular component, and in the interests of computational tractability these discrete biomass models take a simplified approach to mechanical deformation. The formation of surface architecture, in particular, depends on the manner in which growth-induced stresses are accommodated. A growing number of experimental studies aim to characterize the mechanical properties of biofilm (reviewed by Böl et al. [9]), which behaves as a viscoelastic solid relaxing to a viscous fluid on times of the order of a minute. For a time scale relevant to biofilm growth (hours to days), a viscous fluid description is appropriate. There is therefore a distinct breed of continuum biofilm models which focus, somewhat more theoretically, on the mechanical description of the biofilm, and in particular the slime matrix. By contrast, Ben Amar & Wu [10] present a model for the formation of biofilm pattern which relies on a characterization of the biofilm as an elastic medium.

Dockery & Klapper [11] used a simple Darcy flow model to describe the growth of biofilm with a single substrate supplied across a fixed-width diffusive boundary layer. Their analysis of a one-dimensional travelling wave solution, corresponding to uniform thick biofilm, revealed a fingering instability mechanism selecting surface structure at a wavelength equal to the length scale of substrate penetration into the biofilm.

More recently, several continuum models have been proposed based on polymer solution theory, to reflect the composition of the slime matrix as a solution largely of extracellular polymeric substances (EPS) in water. Polymer solutions differ from simple solutions in that the conformational entropy of polymer chains even in dilute solution can have significant impact on the rheology. For a solution with volume fraction ϕ of polymer, Flory–Huggins theory [12] provides the simplest theoretical description of a free energy of mixing governing the composition-dependent tendency of a polymer solution to imbibe further solvent. In particular, the fact that biofilm tends not to swell indefinitely suggests the existence of a non-trivial equilibrium composition and is consistent with the ‘poor solvent’ regime of polymer solutions. An additional gradient free energy term (cf. [13]) can become relevant in spatially inhomogeneous solutions and in describing the dynamics of phase separation, which has recently been suggested to have relevance to biofilm [14].

Cogan & Keener [15] used a two-fluid model for the biofilm matrix, with an EPS polymer phase occupying only a small volume fraction ϕ relative to the solvent water phase. They ignored the cellular component of biofilm, taking growth of the polymer phase as a proxy for bacterial growth. Both phases were modelled as Newtonian viscous fluids, and the difference between the phase-averaged pressures was constituted via an osmotic pressure term based on a Flory–Huggins-like free energy of mixing. The model was simplified by neglecting inertial terms and interphase drag (thereby decoupling the phases), and assuming a dominant balance between the EPS phase viscous stress and the osmotic pressure term. Analysing the stability of a one-dimensional travelling wave solution with specified nutrient concentration at the biofim surface (rather than across a boundary layer as in [11]), they identified an instability mechanism with mode selection at a finite wavelength dependent on the strength of surface tension.

Zhang et al. [16,17] similarly ignored the cellular volume fraction and modelled biofilm as an EPS–water polymer solution. But rather than neglect the water phase velocity, they used a drift flux approach to the two-phase flow, identifying a mixture velocity and an additional polymer network drift flux relative to the mixture. The drift flux was taken as proportional to the gradient of the osmotic pressure, in which they also included a Cahn–Hilliard-like gradient energy term. The model provided a framework for comparing various material models to constitute the EPS–water mixture stress state, including Newtonian viscous and Johnson–Segelman viscoelastic rheology. Results of steady-state analysis for the various models suggest linear instability with a finite selected wavelength, consistent with results of earlier models [11,15,18]. Numerical simulations and extensions of this model have provided the basis for subsequent studies such as those of Lindley et al. [19] and Zhang [20].

Seminara et al. [21] performed experiments on wild-type and EPS- and flagella-deficient strains of Bacillus subtilis. Disc-shaped colonies on agar plates were modelled with a two-fluid polymer solution model assuming dominant momentum balance between osmotic pressure and EPS viscous stress. Further assumptions included swelling quasi-equilibrium (small variation in biofilm composition from swelling equilibrium) and a lubrication theory approximation to reduce the spatial dimension of the model based on the thin disc morphology of the observed biofilm colonies. The model reduced to a degenerate porous medium-type equation for the radial biofilm depth variation, assuming a leading pre-wetted film rather than dealing explicitly with the contact line. The model supported the proposition that biofilm spread is controlled by swelling due to growth-generated osmotic pressure, rather than direct inter-neighbour cell jostling.

Winstanley et al. [22] started from a two-fluid model with Newtonian viscous stresses, and used estimates of model parameters to argue that the viscous stress terms are in fact negligible: at time scales larger than minutes the dominant momentum balance is between interphase drag and osmotic pressure (contrary to the assumptions of [15]). An estimate of the gradient energy parameter suggested that a sharp interface approximation is appropriate for model domains with spatial scale larger than a few micrometres. The study investigated one-dimensional solutions and outlined a potential simplification for the model in higher dimensions.

In this paper, we elaborate on the reduction of the Winstanley et al. [22] model in one lateral dimension in the simplifying regime, where EPS concentration variation through the biofilm is small. We analyse the stability of one-dimensional solutions and show that numerical solutions are consistent with the stability results.

2. Biofilm model

We begin by recalling the polymer/solvent model presented by Winstanley et al. [22]. The EPS has volume fraction ϕ, which is taken to include also the bacteria cells, which themselves produce the EPS. The phase-averaged velocity of the EPS is denoted v, so that the EPS flux is ϕv, and the conservation of EPS is described by the equation

ϕt+.(ϕv)=ϕg(c).
2.1

Here g is a specific growth term which represents the net fractional volumetric rate of extrusion of EPS from the bacterial cells (as well as the growth of the cells themselves). While it may plausibly also depend on the EPS volume fraction, we assume it depends only on the concentration of a rate-limiting nutrient, which for example could be a terminal electron acceptor such as oxygen. Winstanley et al. proposed the commonly adopted Monod form:

g=GcK+c,
2.2

where c is the nutrient concentration, and G and K are constants.

In a similar manner, conservation of the water phase is described by

ϕt+.[(1ϕ)w]=0,
2.3

where w is the phase-averaged water velocity. Uptake and transport of nutrient is given by

[(1ϕ)c]t+.[(1ϕ)cw]=.[(1ϕ)Dc]ϕr(c);
2.4

the equation represents a conservation law for the concentration c of nutrient, allowing also for molecular diffusion, as well as the uptake of nutrient by the bacteria; D is the diffusion coefficient, and we take the uptake function to be

r=RcK+c,
2.5

so that G/R is a yield coefficient for the EPS growth.

The differential equations of the model are completed by a momentum equation for each phase:

0=μEPS.[ϕT]μwk(1ϕ)2(vw)Ψϕpand0=μwk(1ϕ)2(vw)(1ϕ)p,}
2.6

where the tensor T is

T=v+vTθ(.v)δ
2.7

[23],1 where δ is the unit tensor, μEPS and μw are the viscosities of EPS and water, respectively, k is the medium permeability and p is the fluid pressure. These equations are properly derived from models of slow, two-phase viscous flow [24], in which the inertial terms are neglected, and there is an interactive drag between the phases, here chosen as a term proportional to the velocity difference, and such that the water momentum equation is equivalent to Darcy's law. Winstanley et al. additionally included a viscous term for the water, but showed that it was negligible and inconsequential, and so we omit it here.

In two-phase flow theory, the pressures of each phase are generally different and a relationship between them must be prescribed, based on the microscale description of the media. In the present case, consideration of the Flory–Huggins theory of polymer interaction [12] leads to a description of this relation in terms of an osmotic pressure Ψ which can be approximately given for small values of ϕ by

ΨEL[(χ12)ϕ2+13ϕ3],
2.8

where EL is the monomer site energy density and χ is the Flory interaction parameter.2 Evidently, (2.6)2 allows us to define an EPS pressure of the form

ps=p+ϕeqϕΨ(ϕ)dϕϕ.
2.9

A stable equilibrium between the two-phase mixture and the pure solvent occurs when the osmotic pressure is zero and Ψ′>0, and thus

ϕ=ϕeq=3(χ12).
2.10

The boundary conditions for the equations, suitable to describe a growing biofilm in 0<z<s(x,y,t), where z is the distance from the wall at z=0 to which the biofilm is attached, are those of no slip of EPS, no flow of water through the wall and no flux of nutrient at the wall:

v=0,w.n=0,pn=0,cn=0atz=0;
2.11

at the top of the biofilm, suitable conditions are those of a prescribed concentration influx, prescribed ϕ (specifically representing the stable equilibrium Ψ=0 where Ψ′>0, if χ>12, as we assume), continuous water pressure and no traction on the EPS:

dccn=c0c,ϕ=3(χ12),p=0,σ.n=0atz=s,
2.12

where dc is a suitable mass transfer coefficient, representing the thickness of a boundary layer over which the external concentration varies; the stress tensor is

σ=psδ+μEPST;
2.13

in addition, there is a kinematic condition

st+v.(sz)=0,
2.14

which arises from the usual fluid mechanical assumption that particles on the interface remain there, or, in other words, the velocity of the interface is equal to the velocity of particles on it. Winstanley et al. [22] showed that these conditions are sufficient to fully determine the solution for the growth of a laterally uniform biofilm.

Note that this model shares a problem with one of an ordinary viscous fluid, which stems from the combination of the kinematic condition (2.14) and the no-slip boundary condition in (2.11), which together imply that the margin of the biofilm patch is unable to move because the cells adhere to the wall. Rather than address this issue here, in this paper we consider the behaviour of a biofilm patch far from its margin. Fowler & Winstanley [25] have developed a new theory, which they call flux intensity theory, to describe the motion near the margin, and a future piece of work will combine that with the current study.

(a) Non-dimensionalization

The model is non-dimensionalized in exactly the same way as in Winstanley et al. [22]. The dimensionless form of the model can be written as

ϕt+.(ϕv)=ϕg,g=cκ+c,εϕt+.[(1εϕ)w]=0,α(1εϕ)[ct+w.c]=.[(1εϕ)c]ϕg,0=β.[ϕT]ϕ(1εϕ)(vw)Ψεϕp,0=β.[ϕT]Ψp,T=v+vTθ(.v)δandΨ=λϕ2+13ϕ3+O(ε),}
2.15

and the dimensionless parameters are identical to those in the earlier paper. The parameter λ can be defined (taking into account the footnote following (2.8)) as λ=ϕeq/3ϕ0, where ϕ0 is the scale used for the volume fraction of EPS, and is determined by a balance between the osmotic pressure and the interfacial drag, while ϕeq is the equilibrium fraction of EPS, typically 1–5%; Winstanley [26], p. 66 estimated typical values of λ=0.9–4.5 (the values 0.5–2.3 of Winstanley et al. [22] are associated with the incorrect factor 16 which they used in (2.8)).

(i) Boundary conditions

The scaled boundary conditions take the form

v=0,w.n=0,pz=cz=0onz=0
2.16

and

ηcn=1c,p=0,σ.n=0,ϕ=3λonz=s,
2.17

where the scaled stress tensor is

σ=psδ+βTε
2.18

and

ps=p+1ε3λϕΨ(ϕ)dϕϕ.
2.19

The parameter η is defined by

η=dcd,
2.20

where the length scale d is of the order of a typical biofilm thickness, d~100 μm.

Lastly, the kinematic boundary condition on z=s takes the form

st+v.(sz)=0.
2.21

Estimates of the parameters as described by Winstanley et al. [22] are given in table 1. The parameter η was not used by them, as they assumed prescribed concentration c at the free surface, equivalently η=0. Taking η>0 allows for a finite nutrient transfer rate from the medium above.

Table 1.
Typical values of the dimensionless parameters. The values of η and λ are not well constrained.

(b) Reduction of the model

The parameters α, β and ε are all small, and we will eventually consign them all towards zero. Putting α=0 is an assumption of quasi-static reaction kinetics which is commonly made [27], the basis of which is lucidly discussed by Kissel et al. [28]. The asymptotic procedure which we use is novel and non-standard. The neglect of β is a singular perturbation, for which the normal procedure would be to additionally neglect a boundary condition associated with boundary layer behaviour. While we do this, the reduced model obtained is degenerate, and some subtlety is required in elucidating the correct behaviour.

Before beginning our asymptotic simplification, we rescale the equations based on Winstanley et al.'s observation that, even for moderate λ, the assumption of large λ is highly accurate. This suggests that we first rescale the model using

ϕ=3λ(1+νΦ),s,x,v,w13λ,
2.22

where

ν=19λ3.
2.23

This leads to the model equations in the form

νΦt+.[(1+νΦ)v]=(1+νΦ)g,3λενΦt+.[{13λε(1+νΦ)}w]=0,α3λ{13λε(1+νΦ)}[ct+w.c]=.[{13λε(1+νΦ)}c](1+νΦ)g,0=3λβ.[(1+νΦ)T](1+νΦ)[13λε(1+νΦ)](vw)Ψ3λε(1+νΦ)p,0=3λβ.[(1+νΦ)T]Ψp,T=v+vTθ(.v)δandΨ=Φ+O(ν).}
2.24

The boundary conditions take the same form as before, except that now

Φz=0onz=0,Φ=0onz=s,
2.25

and the scaled EPS pressure is

ps=p+13ελ[Φ+O(ν)].
2.26

A notional estimate is 3λ≈10, and thus we suppose that the terms in α and ε are negligible in (2.24), but we temporarily retain the terms in β. Additionally, we suppose that ν[double less-than sign]1, which requires only that λ1. For convenience, we define

β¯=3λβ103.
2.27

Ignoring the terms in α, ε and ν, the equations can now be written in the form

.v=g,g=cκ+c,.w=0,2c=g,vw=p,β¯.T=β¯[2v+(1θ)(.v)]=Φ+pandT=v+vTθ(.v)δ;}
2.28

the upper stress boundary condition is, from (2.17), (2.18), (2.26) and (2.27), allowing for the fact that p=Φ=0 at z=s,

T.n=(v+vT).nθ(.v)n=0,
2.29

and the kinematic condition is

st+v.(sz)=0.
2.30

Our earlier comments on the subtlety of the procedure are now evident in (2.28). If we simply put β¯=0, then we can solve for vw, but not for the velocities individually. The problem is that (2.28)5 implies that

×[2v+(1θ)(.v)]=0,
2.31

but this extra constraint is lost if one immediately puts β¯=0. A little sleight of hand is therefore necessary to allow the limit β¯0. First, we define a function H via

Φ+p=β¯[H2Φ+(1θ).v];
2.32

substituting this into (2.28)5 then yields, using also (2.28)4, and without approximation,

H=2w+β¯2[2v+(1θ)(.v)].
2.33

If we now ignore the term in β¯, we have the equations of Stokes flow for w:

.w=0andH=2w,}
2.34

where the first equation is simply (2.28)2.

Equations (2.28)4 and (2.28)5 combine to yield

vwΦ,
2.35

and then (2.28)1, (2.28)2, again neglecting the O(β¯) term, imply

2Φ.v=g.
2.36

Because of the loss of the viscous terms in (2.28)5, it is possible that some of the boundary conditions on v, and consequently w, may not be satisfied. Boundary conditions for w at z=0 follow from those of no slip for v, and are

w.n=0,w.t=Φ.tonz=0,
2.37

where t is the tangent vector to the plane. Boundary conditions at the surface follow from the normal stress condition (2.29) on v, and this leads to

n.(w+wT)=2n.Φ+θgn.
2.38

This seems fine, in the sense that (2.37) and (2.38) provide two conditions each for the fourth-order Stokes flow model (2.34), but they are non-standard in the sense that there is no condition on the pseudo-pressure H, which is suspicious. In fact, the constraint that Φ=p=0 at z=s, together with (2.36), implies that

H=(2θ)g(c)onz=s,
2.39

which together with (2.38) provides a further condition (thus one too many) for (2.34).

The resolution appears to lie in the fact that the neglect of β¯ in (2.28) is indeed singular, and there is a singular layer in which the viscous term is important. We combine (2.28)4 and (2.28)5 in the form

β¯[2v+(1θ)(.v)]=v(wΦ).
2.40

Taking the normal component and denoting v.n=vn, the relevant boundary layer approximation of this is

β¯2vnn2=vn(wnΦn),
2.41

where [partial differential]/[partial differential]n=n.[nabla], and the second term on the left remains small since [nabla]. v=g, and thus its normal component is β¯g(c)(c/n), which is small as c has no singular behaviour; the approximate solution of (2.41) satisfying the normal component of (2.29), which is

vnn=12θg,
2.42

and decaying far from the boundary is

vn=(wnΦn)+12θβ¯gexp[nβ¯],
2.43

where n=|n|, and we take n pointing out of the biofilm patch. Thus, the normal component of (2.29) is satisfied by means of a boundary layer, and it is not necessary that the outer solution satisfy this condition. The boundary conditions for (2.34) at the surface therefore consist of (2.39), together with the tangential component of (2.38), which takes the form, written in terms of components,

tinj(wixj+wjxi)=2tinj2Φxixj.
2.44

(ii) Summary

The variables c and Φ are determined by the equations

2c=2Φ=g(c),g(c)=cκ+c,ηcn=1c,Φ=0atz=sandcz=Φz=0atz=0.}
2.45

The cell velocity is then determined from

v=wΦ,
2.46

and the solvent velocity is determined by solution of the Stokes flow problem

.w=0andH=2w,}
2.47

with the boundary conditions

wn=0,t.w=t.Φatz=0andH=(2θ)g(c),tinj(wixj+wjxi)=2tinj2Φxixjatz=s.}
2.48

Finally, the kinematic condition at the biofilm surface is simply

st+(wΦ).(sz)=0.
2.49

3. Stability

(a) One-dimensional solution

Winstanley et al. [22] gave a one-dimensional solution of the above system in the case that η=0, though with only tangential discussion of the limit ν→0. If solutions depend only on z and t, then we find H is constant, w=0, and the model reduces to

Φ=csc,s˙=csandc=cκ+c,c(0)=0,ηcs=1cs,}
3.1

where cs=c(s,t), cs=cz(s,t), the overdot denotes a time derivative and the prime denotes a spatial derivative with respect to z. The solution for c is a quadrature,

cwcdc[A(c)A(c0)]1/2=z,
3.2

where we can take

A(c)=20cg(c)dc=2[cκln(1+cκ)]
3.3

(thus A(0)=0), and the basal concentration at the wall, cw, is determined by

cwcsdc[A(c)A(c0)]1/2=s,
3.4

while the surface boundary condition for c gives

η{A(cs)A(cw)}1/2=1cs,
3.5

and s then satisfies the ordinary differential equation

s˙={A(cs)A(cw)}1/2.
3.6

This solution settles, as t for finite κ, to a travelling wave solution in which cw→0, s˙V=A(cs)=22κln(1+1/κ), and c is given implicitly by

zs=zVt=ccsdcA(c),
3.7

and cs is given by

ηA(cs)=1cs.
3.8

(b) Linear stability analysis

We study the stability of the travelling wave solution given above. We assume a two-dimensional flow in which the coordinates are x,z, and we can then define a stream function ψ such that

w=(ψz,ψx).
3.9

The normal and tangent vectors are

n=(sx,1)(1+sx2)1/2andt=(1,sx)(1+sx2)1/2,
3.10

and then the boundary conditions at the surface take the form

(1sx2)(ψzzψxx2Φxz)=2sx(ΦzzΦxx+2ψxz),H=(2θ)g(c),Φ=0,η(czsxcx)(1+sx2)1/2=1candst+ψx+Φz+(ψzΦx)sx=0,}
3.11

while the boundary conditions at the base become

c,ψ0,Φcsasz,
3.12

where cs is the constant solution of (3.8).

We change to a moving frame in which

s=Vt+Σandz=Vt+ζ;
3.13

other than changing z to ζ and s to Σ in the equations and boundary conditions, the only change occurs in the kinematic condition, which becomes

V+Φζ+Σt+ψx+(ψζΦx)Σx=0onζ=Σ.
3.14

The basic state for the system is

ψ=0,Σ=0,Φ0=csc,H=(2θ)g(cs),A(cs)=1csη=Vandc=c0(ζ),ζ=ccsdcA(c).}
3.15

Denoting perturbations to ψ,c,Φ,H by Ψ,C,ϕ,h, respectively, we find the linearized equations are

hx=2Ψζ,hζ=2Ψxand2C=2ϕ=g(c0)C,}
3.16

with linearized boundary conditions

C,Ψ,ϕ0asζ
3.17

and

ϕζ+Φ0Σ+Σt+Ψx=0,ΨζζΨxx2ϕxζ=2Φ0Σx,h=(2θ)g(cs)(C+c0Σ),ϕ+Φ0Σ=0and(C+c0Σ)=η(Cζ+c0Σ)onζ=0.}
3.18

We take

Σ=eikx+σt,k>0,
3.19

without loss of generality, and then

Ψ=f(ζ)Σ,C=b(ζ)Σ,ϕ=a(ζ)Σ,h=d(ζ)Σ,
3.20

where we find

f=(A+Bζ)ekζ,d=2iBkekζ,a=Dekζb(ζ)andb=[k2+g{c0(ζ)}]b,}
3.21

and the boundary conditions (3.18) then imply

σ=a0Φ0ikA,2Bk+2Ak22ik(kDb0)=2ikΦ0,ikB=(112θ)g0(b0+c0),Db0+Φ0=0andb0+c0=η(b0+g0),}
3.22

where a0=a′(0), b0=b(0), b0=b′(0), and the other quantities with suffix zero are steady-state solutions evaluated at ζ=0.

Simplification of this yields

σ=(1(1/2)θ)(b0+c0)g0k,
3.23

and the value of b0 is determined from the boundary value problem

b=[k2+g{c0(ζ)}]b,ηb+b=ηg0c0atζ=0andb0asζ.}
3.24

To assess the stability criterion (3.23), we take as a non-restrictive example the function

g(c)=μc;
3.25

then the equation for b (and c0) can be solved, giving

c0=eμζ1+ημandb=μe(k2+μ)1/2ζ1+η(k2+μ)1/2,
3.26

from which we find

σ=(1(1/2)θ)μ3/2η[(k2+μ)1/2μ][1+η(k2+μ)1/2][1+ημ]k,
3.27

which is positive, since we may take θ23. Note that when η=0, σ=0, and the steady state is neutrally stable. For η≠0, the growth rate tends to zero as k→0 and k, and apparently the model is missing some dissipative term which would enable stabilization at high wavenumber.

In more detail, we may illustrate the situation with the toy system

uxx=uvandvt=u+Dvxx,}
3.28

for which perturbations to the steady state u=v=0 have solutions exp(σt+ikx), where

σ=11+k2Dk2.
3.29

The solution is linearly unstable at long wavelengths (small k), but the diffusion term causes σ to become negative at large k, and the solution is well posed. However, if D=0, the second equation becomes hyperbolic, and σ→0 as k.

We suggest that the same sort of behaviour may be occurring here. In fact, it seems the remaining dissipative term β¯ in (2.28) may provide the required stabilizing term. A partial argument for this can be given by re-doing the stability analysis above by replacing (2.47)2 with (2.33). Following the same analysis through, we find that (3.27) is replaced by

σ=μ3/2{1+η(k2+μ)1/2}[(1(1/2)θ)η{(k2+μ)1/2μ}(1+ημ)k12β¯{1+(1θ)μ}k],
3.30

so that the weak viscous term provides a stabilizing effect at large wavenumber. Illustration of the growth rate as a function of k is given in figure 1, both with and without the stabilizing term. Note how small the growth rate is.

Figure 1.
The growth rate σ as a function of wavenumber k, using (3.30) with parameters μ=1, η=0.1, θ=0. The solid line gives the rate for β¯=0, the dashed one for β¯=0.005. (Online version in colour.) ...

4. Numerical solutions

Despite the comments in the preceding section, we have not attempted to solve the model where the weak viscous terms are included. Instead we solve (2.45)–(2.47) using a front-fixing transformation (e.g. [29]), but to compare with the stability analysis we have used both the function g(c)=c/(c+κ) and the simpler g(c)=μc; the results are similar, but we use the latter in our figures below. For our computations, we ignore the divergence term by putting θ=0. It can be shown that non-zero values of θ lead to the same results, adjusted by a rescaling of the variables. This might also be inferred from (3.27). The problem is mapped from a domain with a time-dependent boundary to one with a fixed, time-independent boundary. This is done by a transformation of the vertical coordinate system from z to ξ, via ξ=z/s(x,t). There are additional terms in the differential equations as a result; however, we are able to implement a standard finite-difference scheme on a fixed grid after the transformation has been made.

After the coordinate transformation has been performed, we discretize the governing equations and boundary conditions using a finite-difference scheme with semi-implicit time-stepping. We use a staggered, Cartesian, two-dimensional grid with ni×nj discrete points. There are periodic boundary conditions in x and boundary conditions (2.48) and (2.49) at the base and free surface. The nonlinear system of equations is solved using a Newton–Krylov solver, provided by the Portable Extensible Toolkit for Scientific Computation [3032]. At each time step, we solve for c and ϕ via (2.45), before using these solutions in the boundary conditions for the coupled u,w and H solver. The solver finally uses explicit time stepping to solve for the new surface elevation. The solution at each step is accepted when the absolute size of the nonlinear system residual is less than 10−8.

We solve the system in domains of varying widths, all of sufficient depth that the travelling wave limit is an accurate description with c1 at the bed. To test stability of the system, we introduce a perturbation to the initially flat surface. We do this with a regular perturbation in the form of a sine wave. The results are discussed in the following section.

As we shall see, while our numerical method is successful, the front-fixing method will break down when the interface develops corners, which it does, due possibly to the neglect of the small viscous term in (2.33). One could of course simply include the viscous term, but realistic values of β¯ would then yield an impossibly stiff problem. The question then arises whether other methods might avoid this difficulty. Different numerical methods for this problem have been discussed by Du et al. [33]. Specifically, they develop an interface-capturing method where the interface separates a viscous two-phase gel network and a viscous pure solvent. Their method uses an artificial small additional network volume fraction ε in the equations, so that they solve the two-fluid viscous model everywhere, and the interface is ‘captured’ by the condition that the volume fraction reaches zero. With this addition, their numerical method then uses finite differences.

It seems that all such methods using finite differences will come to grief when the interface develops corners, and that a much more subtle approach would be necessary in order to track a non-smooth interface.

5. Discussion

Examples of the solutions are shown in figures 25. There is nothing too interesting in the concentration or velocity distributions, and we focus on the interface shape. Figure 2 shows a typical such solution, which shows the slow growth of spatial instability from an initial sinusoidally oscillating interface. For comparison, figure 3 shows the growth from a more random initial state, which yields a more spotty growth field. As discussed below, one cannot take a completely random initial state, since this model lacks a high wavenumber smoothing property, and the numerical solution would be unable to proceed in that case.

Figure 2.
Evolution of the interface s at increasing times t=0,0.5,…,3 on a periodic domain of length 2π, where the initial condition for s is s=2+0.05sin5x, and we have used g=μc, with the parameter values being μ=1, η ...
Figure 3.
(a,b) Evolution of the interface s at increasing times t=0,0.2,…,1.2 on a periodic domain of length 2π, where the initial condition for s is s=2+0.01{cos2x+sin3x+0.7cos4x+0.4sin5x+0.5sin6x+ ...
Figure 5.
Growth rate determined from the numerical solution as portrayed in figure 4, compared with the theoretical result in (3.27). The parameters used were μ=1 and η=0.1. The same wavenumbers are indicated as in figure 4, and ...

In seeking to compare our numerical results with the stability analysis, we encounter a problem which in retrospect is not surprising. It is caused by the form of the kinematic equation (2.30) for the interface s. So long as the normal interface velocity vn=v.n is given, then (2.30) can be solved by the method of characteristics (cf. [25]), and if vn>0 (which it certainly is in the one-dimensional solution, see (3.1), and therefore initially), any initially oscillatory interface will have convergent characteristics and thus form cusps; and, indeed, this is what happens, and, when it does, the solution is singular and our numerical method breaks down.

In more detail, the singularity is manifested as the formation of corners in the biofilm interface. There is nothing intrinsically wrong with the model in this respect: it is quite feasible that the surface z=s(x,t) should form such kinks. However, our numerical front-fixing method automatically introduces second derivatives of the free surface into the differential equations, and the method consequently breaks down when the first derivatives of s approach the singularity. A quite different method would be necessary to prolong the numerical solution past this point.

Such an aspiration is beyond the scope of the present study, since we were not even sure that the model would be sufficiently well posed to have a solution in two dimensions. A simpler alternative is to include an artificial diffusion term δsxx in the kinematic condition, and we have implemented this; the problem then is that this provides a stabilizing effect which can remove the instability altogether, unless we choose η to be fairly large. However, such results are difficult to follow, because the method then requires such small time steps that it becomes prohibitive to run long enough.

One simple apparent remedy would be to restrict the initial amplitude A of the spatial variation of s to be very small, since we are only really seeking to validate the linear stability result. The problem with this is evident in figure 4. Since the ordinate measures log10A, the value of A˙/A is the slope of the amplitude curves, and is evidently not constant, particularly for the small wavenumber perturbations. The reason for this is evident in figure 2. Even though the initial mean value of s is sufficiently large for the one-dimensional solution to be accurate (if s were constant), the spatial variation already causes a nonlinear adjustment to the growth rate: what we see in figure 4, most obviously at the lowest wavenumber, is a rapid transient where the growth rate of the mean interface position adapts to its initial distortion, before proceeding to grow pseudo-exponentially. This is why we estimate the slope of the curves some way beyond the initial transient. On the other hand, we cannot go too far, because cusp formation indicates the breakdown of the linear theory. With all these provisos in mind, we consider that the estimated growth rates determined from figure 4 are reasonably in agreement with the theoretical estimate, as shown in figure 5.

Figure 4.
Evolution of the amplitude A of the spatial oscillation of s with time, for varying initial oscillation periods (and with periodic boundary conditions on the corresponding periodic domain). The function g is taken to be g=μc, with μ=1, ...

6. Conclusion

In this paper, we have developed the one-dimensional theory for biofilm growth presented by Winstanley et al. [22] to consider solutions which are not laterally uniform. The development and simplification of the model is non-trivial, and it leads to a highly non-standard type of Stokes flow model, whose boundary conditions preclude confidence. Nevertheless, we also develop a direct numerical solution, which shows that a planar interface is unstable, and the results are in satisfactory agreement with a linear stability analysis.

Our analysis also shows, both numerically and analytically, that the interface will develop cusps in finite time. As explained near the beginning of §5, the formation of cusps is a result of the hyperbolic nature of the kinematic condition (2.49), and we consider this to be a real feature of our model, which unfortunately means that our numerical method which assumes a smooth interface is doomed to fail when the cusps form. Such singularities can be alleviated by incorporating extra physical effects, most notably a surface tension, as was done by Cogan & Keener [15], but we do not consider that there is any physical justification for such a term, as it would essentially be a surface tension for the polymer matrix. Rather, we tentatively suggest that the formation of cusps in the interface may represent the first steps in our theory of the formation of more exotic, non-smooth architectures, such as are seen experimentally.

We mentioned in the introduction that, broadly, biofilm architectures appear to be favoured in conditions of poor nutrient supply (e.g. [34]). In order to compare such a statement with our theoretical results, we need to tease apart the dependence of the parameters μ and η in (3.27) on the nutrient level c0 and the transport layer thickness dc in (2.12). Low nutrient supply means low c0 and/or high dc. The latter implies high η (see (2.20)), while, since in (3.25) g=μc replaces g=c/(c+κ) in (2.28), we essentially have μ~1/κ, and, since κc01 [34], eqn (2.11), low c0 means low μ.

In figure 6, we plot σ given by (3.27) for three sets of combinations of μ and η. Ideally, we would see larger σ for low μ and high η, but the opposite conclusion seems to be the case, at least for fixed wavenumber. However, more important should be the maximum growth rate, and this does not vary very much. Possibly a more robust conclusion from this figure would be that poor nutrient supply favours growth at longer wavelengths, but it may be rash to try and infer too much from such a simple stability theory.

Figure 6.
The growth rate σ(k,μ,η) given by (3.27) for values μ=5, η=0.2; μ=η=1; and μ=0.2, η=5 as indicated. (Online version in colour.)

Footnotes

1The usual value of θ=23. The divergence term in (2.7) was omitted by Winstanley et al. [22], i.e. θ=0. Inclusion of a bulk viscosity μB[23], p. 154 would give a value θ=23μB/μEPS. We will proceed with the general term θ, although, as we shall see, the precise value has little effect.

2The factor 13 was incorrectly given as 16 in equation (2.7) of Winstanley et al. [22], with largely cosmetic consequences.

Authors' contributions

A.C.F. devised and solved the model, and wrote the paper. T.M.K.-S. constructed and solved the numerical model. H.F.W. performed the literature review, and assisted with the analytic and numerical work. All authors worked on the final text.

Competing interests

We have no competing interests.

Funding

This publication has emanated from research conducted with the financial support of Science Foundation Ireland under grant nos. SFI/09/IN.1/I2645 and SFI/13/IA/1923.

References

1. Whitman WB, Coleman DC, Wiebe WJ 1998. Prokaryotes: the unseen majority. Proc. Natl Acad. Sci. USA 95, 6578–6583. (doi:10.1073/pnas.95.12.6578) [PubMed]
2. Wanner O, Gujer W 1986. A multispecies biofilm model. Biotechnol. Bioeng. 28, 314–328. (doi:10.1002/bit.260280304) [PubMed]
3. Wanner O, Reichert P 1996. Mathematical modeling of mixed-culture biofilms. Biotechnol. Bioeng. 49, 172–184. (doi:10.1002/(SICI)1097-0290(19960120)49:2<172::AID-BIT6>3.0.CO;2-N) [PubMed]
4. Tiwari SK, Bowers KL 2001. Modeling biofilm growth for porous media applications. Math. Comp. Model. 33, 299–319. (doi:10.1016/S0895-7177(00)00246-6)
5. Lee MW, Park JM 2007. One-dimensional mixed-culture biofilm model considering different space occupancies of particulate components. Water Res. 41, 4317–4328. (doi:10.1016/j.watres.2007.06.026) [PubMed]
6. Costerton JW. 2007. The biofilm primer, vol. 1. Berlin, Germany: Springer Science & Business Media.
7. Stoodley P, Cargo R, Rupp CJ, Wilson S, Klapper I 2002. Biofilm material properties as related to shear-induced deformation and detachment phenomena. J. Ind. Microbiol. Biotechnol. 29, 361–367. (doi:10.1038/sj.jim.7000282) [PubMed]
8. Laspidou CS, Kungolos A, Samaras P 2010. Cellular-automata and individual-based approaches for the modeling of biofilm structures: pros and cons. Desalination 250, 390–394. (doi:10.1016/j.desal.2009.09.062)
9. Böl M, Ehret AE, Bolea Albero A, Hellriegel J, Krull R 2013. Recent advances in mechanical characterisation of biofilm and their significance for material modelling. Crit. Rev. Biotechnol. 33, 145–171. (doi:10.3109/07388551.2012.679250) [PubMed]
10. Ben Amar M, Wu M 2014. Patterns in biofilms: from contour undulations to fold focussing. Europhys. Lett. 108, 38003 (doi:10.1209/0295-5075/108/38003)
11. Dockery J, Klapper I 2002. Finger formation in biofilm layers. SIAM J. Appl. Math. 62, 853–869. (doi:10.1137/S0036139900371709)
12. Flory PJ. 1953. Principles of polymer chemistry. Ithaca, NY: Cornell University Press.
13. Cahn JW, Hilliard JE 1958. Free energy of a nonuniform system. I. Interfacial free energy. J. Chem. Phys. 28, 258–267. (doi:10.1063/1.1744102)
14. Ghosh P, Mondal J, Ben-Jacob E, Levine H 2015. Mechanically-driven phase separation in a growing bacterial colony. Proc. Natl Acad. Sci. USA 112, E2166–E2173. (doi:10.1073/pnas.1504948112) [PubMed]
15. Cogan NG, Keener JP 2004. The role of the biofilm matrix in structural development. Math. Med. Biol. 21, 147–166. (doi:10.1093/imammb/21.2.147) [PubMed]
16. Zhang T, Cogan NG, Wang Q 2008. Phase field models for biofilms. I. Theory and one-dimensional simulations. SIAM J. Appl. Math. 69, 641–669. (doi:10.1137/070691966)
17. Zhang T, Cogan NG, Wang Q 2008. Phase field models for biofilms. II. 2-D numerical simulations of biofilm-flow interaction. Commun. Comput. Phys. 4, 72–101.
18. Picioreanu C, van Loosdrecht MCM, Heijnen JJ 1998. Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach. Biotechnol. Bioeng. 58, 101–116. (doi:10.1002/(SICI)1097-0290(19980405)58:1<101::AID-BIT11>3.0.CO;2-M) [PubMed]
19. Lindley B, Wang Q, Zhang T 2012. Multicomponent hydrodynamic model for heterogeneous biofilms: two-dimensional numerical simulations of growth and interaction with flows. Phys. Rev. E 85, 031908 (doi:10.1103/PhysRevE.85.031908) [PubMed]
20. Zhang T. 2012. Modeling of biocide action against biofilm. Bull. Math. Biol. 74, 1427–1447. (doi:10.1007/s11538-012-9719-z) [PubMed]
21. Seminara A, Angelini TE, Wilking JN, Vlamakis H, Ebrahim S, Kolter R, Weitz DA, Brenner MP 2012. Osmotic spreading of Bacillus subtilis biofilms driven by an extracellular matrix. Proc. Natl Acad. Sci. USA 109, 1116–1121. (doi:10.1073/pnas.1109261108) [PubMed]
22. Winstanley HF, Chapwanya M, Fowler AC, McGuinness MJ 2011. A polymer-solvent model of biofilm growth. Proc. R. Soc. A 467, 1449–1467. (doi:10.1098/rspa.2010.0327)
23. Batchelor GK. 1967. An introduction to fluid dynamics. Cambridge, UK: Cambridge University Press.
24. Drew DA, Passman SL 1999. Theory of multicomponent fluids. New York, NY: Springer.
25. Fowler AC, Winstanley HF 2012. Movement of a sessile cell colony. Math. Proc. R. Ir. Acad. 112A, 79–91. (doi:10.3318/PRIA.2011.112.08)
26. Winstanley HF. 2011. Mathematical modelling of biofilm growth and bioavailability. DPhil. thesis, University of Oxford, Oxford, UK.
27. Klapper I. 2012. Productivity and equilibrium in simple biofilm models. Bull. Math. Biol. 74, 2917–2934. (doi:10.1007/s11538-012-9791-4) [PubMed]
28. Kissel JC, McCarty PL, Street RL 1984. Numerical simulation of mixed-culture biofilm. J. Environ. Eng. 110, 393–411. (doi:10.1061/(ASCE)0733-9372(1984)110:2(393))
29. Crank J. 1984. Free and moving boundary problems. Oxford, UK: Clarendon Press.
30. Balay S, Gropp WD, McInnes LC, Smith BF 1997. Efficient management of parallelism in object oriented numerical software libraries. In Modern software tools in scientific computing (eds E Arge, AM Bruaset, HP Langtangen), pp. 163–202. Basle, Switzerland: Birkhäuser Press.
31. Balay S. et al. 2013. PETSc users manual. Technical report ANL-95/11 - Revision 3.4, Argonne National Laboratory, Lemont, IL, USA.
32. Balay S. et al. 2014. Portable, Extensible Toolkit for Scientific Computation. See http://www.mcs.anl.gov/petsc.
33. Du J, Guy RD, Fogelson AL, Wright GB, Keener JP 2013. An interface-capturing regularization method for solving the equations for two-fluid mixtures. Commun. Comp. Phys. 14, 1322–1346. (doi:10.4208/cicp.180512.210313a)
34. Heydorn A, Nielsen AT, Hentzer M, Sternberg C, Givskov M, Ersbøll BK, Molin S 2000. Quantification of biofilm structures by the novel computer program COMSTAT. Microbiology 146, 2395–2407. (doi:10.1099/00221287-146-10-2395) [PubMed]

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of The Royal Society