May 23, 2018, Wednesday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Finger Formation in Biofilm Layers

From MathBio

Jump to: navigation, search


This page details a model governing the growth of a biofilm, as given in the paper Finger Formation in Biofilm Layers by J. Dockery and I. Klapper [1]. The model is developed and is solved in the one dimensional case. A particular example is used to demonstrate how the model would be solved in one dimension.


A biofilm is a substance composed of millions of microorganisms that accumulate on surfaces in flowing aqueous environments [1]. These microorganisms include bacteria, fungi, algae, and protozoa. When these microorganisms form into a biofilm, they become physiologically distinct from the planktonic version of the same cell. Depending on the surrounding conditions, the presence of a biofilm may be detrimental or beneficial.


Biofilm development is characterized in five different stages, as shown in the picture. In the first stage, free-floating microorganism attach to the surface, but the attachment is very weak. Then, if the cell is not immediately separated from the surface, the attachment grows stronger, usually through cell adhesion. The biofilm then begins to grow either through cell division or recruitment of new cells. Once the biofilm has been formed, it may have certain properties such as resistance to antibiotics. The final stage in biofilm development is dispersion. This stage allows the biofilm to spread and form new colonies.

A biofilm model generally contains three principle elements. First, the model must incorporate a transport mechanism for bringing nutrients to the biofilm. Second, the model must require some consumption and growth mechanism. Third, the model must also take into account biofilm necrosis and loss. [1]

The model presented on this page is a simple model of a single substrate limited biofilm growing into a static aqueous environment. [1]


  • Mathematics used: Reaction-diffusion equations, Darcy's Law for modelling flow of a fluid through a porous medium.
  • Type of model: Population of microorganisms in flowing aqueous environments with constant density
  • Biological system studied: Homogeneous biofilms (assumed to be viscous and incompressible fluid of constant density)

Setting Up the Model


In setting up the biofilm model, our three dimensional space is first divided into two regions with a boundary located at z=h(x,y,t). The biofilm region is represented by z<h(x,y,t) while hte aqueous region is given by z>h(x,y,t).

Conditions on the Biofilm Region

The biofilm is treated as a homogeneous, viscous, incompressible fluid of constant density [1] and so satisfies Darcy's law

{\vec  {u}}=-\lambda \nabla p

where {\vec  {u}} is velocity, p is pressure, and \lambda is some constant. Additionally, as stipulated by the second principle element of biofilm models, the biofilm may be growing or decaying. This is represented by

{\vec  {\nabla }}\cdot {\vec  {u}}=g

for some growth function g. Combining these two equations leads to

  -\lambda \nabla ^{2}p=g     (1)

At the interface between the regions, the normal velocity is equal to

{\vec  {n}}\cdot {\vec  {u}}|_{{z=h^{{-}}}}=-{\vec  {n}}\cdot \lambda \nabla p|_{{z=h^{{-}}}}

where {\vec  {n}} is the normal vector and z=h^{{-}} indicates that the interface should be approached from the side of the biofilm. However, from (1), we see that p is actually proportional to \lambda ^{{-1}} and g, which implies that the normal velocity is independent of \lambda and is instead dependent on the growth function g. This conclusion is reasonable since the interface velocity is dependent on conservation of mass and how the biofilm grows or decays. As such, this dependence is represented by g, the growth function, and not by the constant \lambda .

Conditions on the Aqueous Region

Typically, the time for bioflim evolution is on the order of hours. As such, the biofilm is assumed to evolve quasi-statically. For this reason, the aqueous region is assumed to be static, implying that the pressure is a constant for this region. For this model, p is set to 0 in the aqueous region. Because p=0 for z>h and the biofilm evolves quasi-statically, the following boundary conditions for the edges of the biofilm region are obtained:


p_{z}|_{{z=-\infty }}=0

Conditions on the Growth Function

In this model, it is assumed for simplicity that the growth of the biofilm is only dependent on the concentration S(x,y,z,t) of some single, limiting substrate such as oxygen or glucose. Additionally, it is assumed that


where u is a usage function that indicates the rate at which the substrate is used. The usage function u is typically modeled by the Michaelis-Menten equation

u(S)=u_{{max}}{\frac  {S}{S+K_{S}}}

where u_{{max}} represents the maximum usage and K_{S} is the half-saturation. The growth function g is also modeled based on tis form and is often proportional to u.

Model Equations


The substrate diffuses through the aqueous region into the biofilm region. In the biofilm region, the substrate diffuses and is also consumed as as well. For this model, it is assumed that there exists some function H(x,y,t) such that H\geq h and that for z\geq H, the substrate concentration is constant. With the introduction of H(x,y,t), we can now characterize the substrate behavior in a maximum of three regions. In the region z\geq H, the substrate is constantly replenished so that its concentration remains at some constant value S_{\infty }. This is meant to model a bulk flow region away from the biofilm [1]. The region h\leq z\leq H is contained withing hte aqueous region, so here the substrate just diffuses. The region z\leq h exists within the biofilm region. Here the substrate both diffuses and is consumed. The behavior of the substrate in each of these regions is therefore modeled by the following equations:

  S=S_{\infty },z\geq H
  S_{t}-D_{1}\nabla ^{2}S=0,h\leq z\leq H     (2)
  S_{t}-D_{2}\nabla ^{2}S=-u(s),z\leq h

The boundary and matching conditions for these equations are as follows:

S|_{{z=H}}=S_{\infty } S_{z}|_{{z=-\infty }}=0 S|_{{z=h^{{+}}}}=S|_{{z=h^{{-}}}} D_{1}\nabla S|_{{z=h^{+}}}=D_{2}\nabla S|_{{z=h^{-}}}

namely, the substrate completely diffuses through the region z=H and z=-\infty and is continuous at the interface z=h between the aqueous and biofilm regions.

Because of the slow evolution of the system, the substrate concentration can also be assumed to evolve quasi-statically. With this assumption, the time-derivative terms in the equations can be neglected. The equations and associated boundary conditions listed in (2), together with (1), define the biofilm model that will be studied.

Scaling and Nondimensionalization of the Model

Now that the model equations have been developed, the next step is to scale and nondimensionalize them. First, the horizontal system width w is introduced. To find the distance the substrate is able to diffuse in the biofilm region before being consumed, we first note that the substrate consumption rate near the top of the biofilm is given by {\frac  {S_{\infty }}{u(S_{\infty })}}. With this, it follows that the diffusion length is equal to l_{s}={\sqrt  {{\frac  {D_{2}S_{\infty }}{u(S_{\infty })}}}}. The time scale T that measures the time for the biofilm to evolve an amount of O(w) is also introduced. Then, we make the following settings

{\tilde  {x}}={\frac  {x}{w}}, {\tilde  {h}}={\frac  {h}{w}}, {\tilde  {H}}={\frac  {H}{w}}, {\tilde  {t}}={\frac  {t}{T}}, {\tilde  {S}}={\frac  {S}{S_{\infty }}},

{\tilde  {u}}={\frac  {u}{u(S_{\infty })}}, {\tilde  {g}}={\frac  {gS_{\infty }}{u(S_{\infty })}}, {\tilde  {p}}={\frac  {pG^{{-1}}\nu \lambda }{D_{2}}}, \epsilon _{i}={\frac  {w^{2}}{D_{i}T}}

where G is the nondimensional growth number, equal to

G={\frac  {w^{2}}{l_{s}^{2}}}={\frac  {w^{2}u(S_{\infty })}{D_{2}S_{\infty }}}

\nu is the nondimensional ratio of the system growth and substrate evolution times, equal to

\nu ={\frac  {TS_{\infty }}{u(S_{\infty })}}

and \epsilon _{i} are nondimensional numbers comparing system diffusion times to the system evolution time. With these settings and dropping the tildes, the scaled model equations become

  \nabla ^{2}p=-g(u(S)),z\leq h,
  \epsilon _{1}S_{t}-\nabla ^{2}S=0,h\leq z\leq H,
  \epsilon _{2}S_{t}=\nabla ^{2}S=-Gu(S),z\leq h

with the interface velocity

  -{\vec  {n}}\cdot {\vec  {\nabla }}p_{{z=h^{-}}}.

Again, because of the slow evolution of the biofilm, \epsilon _{i} can be considered small and set to zero, thereby removing the time derivatives from the model equations. The boundary and matching conditions for the scaled and nondimensionalized system are


p_{z}|_{{z=-\infty }}=0,


S_{z}|_{{z=-\infty }}=0,


K{\vec  {\nabla }}S|_{{z=h^{+}}}\cdot {\vec  {n}}={\vec  {\nabla }}S|_{{z=h^{-}}}\cdot {\vec  {n}}

where K={\frac  {D_{1}}{D_{2}}}.

The 1D Problem

In one dimension, the model equations reduce to

p_{{zz}}=-g(u(S)),z\leq h(t) S_{{zz}}=0,h(t)\leq z\leq H(t) S_{{zz}}=Gu(S),z\leq h(t)

and the evolution of the interface is given by

{\dot  {h}}=-p_{z}|_{{z=h^{-}}}

In this case, the quantities S, p, h, and H are assumed to be only functions of z and t and the boundary and matching conditions are the same as before. Furthermore, the function H(t) is rewritten as


for some constant L.

Substrate Behavior

The first step necessary in solving this model is to determine how the substrate behaves in the aqueous and biofilm regions.

Aqueous Region

We begin by solving the substrate equation in the aqueous region because it is simpler. The equation that needs to be solved is


The general solution to this equation is


By applying the boundary equation S(H)=S(h+L)=1, we can solve for the constant C_{2} and find that


So, rewriting C_{1} as A, the substrate equation in the aqueous region is given by

S(z,t)=1+A(z-(h+L)),h\leq z\leq h+L

where A is a constant to be determined.

Biofilm Region

In the biofilm region, the equation governing substrate behavior is


The solution to this equation must satisfy the conditions

S_{z}(-\infty ,t)=0 S(h,t)+LK^{{-1}}S_{z}(h,t)=1.

The second requirement is the result of combining the following matching conditions:

S|_{{z=h^{+}}}=S|_{{z=h^{-}}} KS_{z}|_{{z=h^{+}}}=S_{z}|_{{z=h^{-}}}.

If we use the previously defined substrate equation in the aqueous region for the terms involving h^{+} in the matching conditions, we can solve for an expression for A and derive our second requirement. More specifically, if this is done, the matching condition S|_{{z=h^{+}}}=S|_{{z=h^{-}}} becomes


and the matching condition KS_{z}|_{{z=h^{+}}}=S_{z}|_{{z=h^{-}}} becomes


From this we see that


and the solution must satisfy


Analysis of Substrate Solution in the Biofilm Region

Before making statements about the substrate solution S(z,t), restrictions must be placed on the usage function u(S). First, if no substrate is present, then no substrate is used:


Second, if the substrate concentration is positive, then the usage of the substrate is also positive:

u(S)\geq 0forS\geq 0.

Then, for convenience, the nondegeneracy condition is assumed

u'(0)\geq 0

and u(S) is assumed to be Lipschitz continuous on S\geq 0 as well.

With these restrictions, the solution S(z,t) is unique. To verify this, a phase space analysis at fixed t is performed with the variables S(z) and Q(z)=S_{z}(z). The system to be satisfied is then

{\dot  {S}}=Q

{\dot  {Q}}=Gu(S).

The solution in this system will be represented as a trajectory that starts on the line


and arrives at the S axis at z=-\infty . If the trajectory meets these requirements, it satisfies the conditions placed on the substrate solution.


The only fixed point of this system is located at the origin when S=Q=0. It is a saddle and has two eigenvectors; namely, the stable eigenvector (1,-{\sqrt  {Gu'(0)}}) and the unstable eigenvector (1,{\sqrt  {Gu'(0)}}). The solution to this system is represented by the trajectory that starts at z=h from the intersection of the unstable eigenvector and the line S+LK^{{-1}}Q=1. This trajectory runs along the unstable eigenvector and therefore arrives at the origin at z=-\infty , thereby satisfying the two requirements.

This solution is unique because no other trajectory satisfies the stipulated conditions. Trajectories that start on the line S+LK^{{-1}}Q=1 below the intersection point with the unstable eigenvector go to infinity as z\rightarrow -\infty . Trajectories that start on the line above the intersection point go to S=0 at some finite value of z. Because S_{z}\geq 0 for these trajectories, they cannot be matched to a 0 solution.

So, because there is only one trajectory that satisfies the appropriate conditions, the substrate solution is unique.

General Solution for Pressure and Interface Evolution

Given the substrate solution, a general solution to the pressure and interface evolution equations can be found. Solving the pressure equation,


yields the general solution

p(z,t)=\int _{{z}}^{{h(t)}}\int _{{-\infty }}^{{r}}g(\eta ,t)d\eta dr

where g(u(S)) has been written as g(z, t). Given this solution, the equation governing interface evolution can be easily found

{\dot  {h}}=-p_{z}(h(t),t)=\int _{{-\infty }}^{{h(t)}}g(s,t)ds.

From this equation, we see that the biofilm interface changes depending on the growth of the biofilm.

A Particular Example

To demonstrate how the biofim model may be solved in one dimension, we choose our usage and growth functions to be


g(u)=\alpha u=\alpha S.

Then, because \alpha can be scaled into the pressure, it is set to \alpha =1. With these particular functions, the model equations become greatly simplified.

The Biofilm Region

In the biofilm region, the equation that needs to be solved is then


This has a general solution of

S(z,t)=C_{1}e^{{{\sqrt  {G}}z}}+C_{2}e^{{-{\sqrt  {G}}z}}.

But, after applying the condition S_{z}(-\infty ,t)=0, this becomes

S(z,t)=C_{1}e^{{{\sqrt  {G}}z}}.

At this point, the second condition S(h,t)+LK^{{-1}}S_{z}(h,t)=1 can be applied. This lets us solve for the constant C_{1}, finding

C_{1}={\frac  {1}{1+K^{{-1}}L{\sqrt  {G}}}}e^{{-{\sqrt  {G}}h}}.

With this, we see that the substrate equation in the biofilm region is equal to

S(z,t)={\frac  {1}{1+K^{{-1}}L{\sqrt  {G}}}}e^{{-{\sqrt  {G}}(h-z)}},z\leq h.

The Aqueous Region

In the aqueous region, we know the substrate equation has the form of


We also know that KS_{z}(h^{+},t)=S_{z}(h^{-},t). So, using the derivative of the substrate equation in the biofilm region for S_{z}(h^{-},t), i. e.

S_{z}(h^{-},t)={\frac  {{\sqrt  {G}}}{1+K^{{-1}}L{\sqrt  {G}}}}

we can find an expression for S_{z}(h,t) and determine that the substrate equation in the aqueous region is

S(z,t)=1-{\frac  {{\sqrt  {G}}}{K}}({\frac  {1}{1+K^{{-1}}L{\sqrt  {G}}}})(h+L-z),h\leq z\leq h+L

Pressure and Interface Evolution Equations

The pressure equation now has the form

p_{{zz}}=-S=-{\frac  {1}{1+K^{{-1}}L{\sqrt  {G}}}}e^{{-{\sqrt  {G}}(h-z)}}

and therefore has the general solution

p(z,t)=-{\frac  {1}{G}}({\frac  {1}{1+K^{{-1}}L{\sqrt  {G}}}})e^{{-{\sqrt  {G}}(h-z)}}+C_{1}z+C_{2}.

Applying the condition p_{z}(-\infty ,t)=0 shows that C_{1} must be equal to 0. Applying the condition p(h,t)=0 shows that C_{2} must be equal to

C_{2}={\frac  {1}{G}}({\frac  {1}{1+K^{{-1}}L{\sqrt  {G}}}}).

So, the solution to the pressure equation is

p(z,t)={\frac  {1}{G}}({\frac  {1}{1+K^{{-1}}L{\sqrt  {G}}}})(1-e^{{-{\sqrt  {G}}(h-z)}}),z\leq h.

The equation governing interface evolution easily follows since

{\dot  {h}}=-p_{z}(h,t)

leading to

{\dot  {h}}={\frac  {1}{{\sqrt  {G}}}}({\frac  {1}{1+K^{{-1}}L{\sqrt  {G}}}})t.

Recent Extension

In their paper [2], Alpkvist and Klapper, cite above discussed paper. They propose a multidimensional continuum model for heterogeneous growth of biofilm systems with multiple species and multiple substrates. The fundamental assumptions for their model are a combination of the ones presented in Wanner and Gujer (1986) [3] and above discussed paper. Their model is notable improvement of above discussed model. Authors of the paper have taken a number of steps to provide an adequate numerical approximation of the coupled nonlinear hyperbolic and elliptic systems describing the model. The numerical procedures have been tested and implemented for two possible biofilm systems. Representative simulations for these two systems have been presented in two and three spatial dimensions. Thus they have developed a model capable of describing and capturing heterogeneity in multiple spatial dimensions for multiple substrates and species.

External Links

  1. Dockery, J., Klapper, I., Finger Formation in Biofiom Layers, SIAM J. Appl. Math, 2001, Volume 62, 853-869
  2. E. Alpkvist, I.Klapper, "A multidimensional multispecies continuum model for heterogenous biofilm development," Bull. Math. Biol. 69, 765-789 (2007)
  3. Wanner, O., Gujer, W., 1986. A multispecies biofilm model. Biotech. Bioeng. 28(33), 314–328