Simulating biofilm growth

1
Simulating biofilm growth
Ravindra Duddu Brian Moran David Chopp
Stephane Bordas
IUTAM 2006 modeling of moving interfaces, INSA
Lyon, France
2
University of Glasgow, Scotland (UK)
3
UNIVERSITY OF GLASGOWDepartment of Civil
EngineeringGroup Mechanics and MaterialsHead
of group Regius Professor Nenad Bicanic
4
Department ofCivil EngineeringGroup Mechanics
and MaterialsHead of group Regius Professor
Nenad Bicanic
Rankine
5
BiofilmsAggregates of bacteria in aqueous
environment
6
Biofilms
7
BiofilmsAggregates of bacteria in aqueous
environment
8
Mathematical representation
9
Biofilm Model Equations
Biomass Balance
where
dividing by
(1)
now the biomass flux is defined as
using the above in (1) and applying the product
rule
10
Biofilm Model Equations
For single-species model
where is the net specific rate of biomass
production
But it can be shown that biomass volume-fraction
is a constant given by
Now we define a velocity potential as given below
So for a single-species biofilm with constant
biomass we can the biomass conservation equation
in the form of a second order differential
equation
Similarly writing a mass balance for substrate
gives
11
Biofilm Equations in 2D - Strong Form
Substrate
Euler
Velocity Potential
Poisson
12
Biofilm Equations in 2D - Weak Form
The trial functions are chosen such that
The test functions are defined as
The weak form of the substrate and velocity
potential equations in the biofilm domain
The substrate equation could be treated as
quasi-steady state and so we drop out the
advection terms in the weak form.
13
Biofilm Equations in 2D Linearization
The resulting linearized equations are given by
14
Biofilm Equations in 2D Discretized equations
The discretized equations after using the XFEM
approximation with discontinuous derivative
enrichment and the penalty method for imposing
Dirichlet boundary conditions along arbitrary
iso-level curves are given by
where
15
Biofilm Equations in 2D Newton-Raphson Solver
The tangent and residual operators are defined by
The nodal values of increment of substrate and
velocity potential could be found iteratively as
till the norm of the residuals is less than the
specified tolerance. Now the nodal values of
velocity potential are used to calculated the
interface speed.
16
The eXtended Finite Element Method (XFEM)
Substrate and velocity potential fields are only
C0 continuous in O and there is a jump in their
derivatives across the interface. The interface
is implicitly defined by a level set function,
. Discontinuous enrichment is a piecewise
continuous, C0, enrichment function D, whose
first derivatives are discontinuous across the
interface
The unknown field, say u(x, t) can be
approximated as
17
The Level Set Method
The level set method is based upon representing
an interface as the zero level set of some higher
dimensional function, F. From an Eulerian
viewpoint we use a fixed mesh of grid points to
implicitly represent the biofilm-fluid interface
as a level contour of a function. For stability
purposes, it is preferable to keep F as a signed
distance function. To update the interface we
advance ? in time using the evolution equation
for the level set method as given below
where F is the interface speed given by the
physics of the phenomenon. Since F is defined
only on the interface and we need F to defined at
the nodes for computation, we define a compatible
Fext using velocity extensions.
18
Objectives of velocity extension

• Construct a velocity field on the whole domain
• Based on the field on the interface
• So that the signed-distance property is conserved
• Without changing the speed of the interface
• Interface speed F propagated outward from the
  interface in the direction of the normal (also
  known as causality)

19
Velocity extensions and Fast Marching Methods
The scalar velocity, F, is defined on the
interface as
But the level set evolution equation requires F
to be known everywhere in the domain. So we
define a suitable Fext which is known at all the
grid nodes
Fast marching method to construct the velocity
field as given below. Say we want to find Fi, j
, the value of Fext at the node (i, j), we have
known
Solving for the unknown, Fi, j , we get
gt0 by causality
20
A simple weighted average
Fext is constant along lines normal to the
interface (use norm grad phi 1 and level set
evolution eqn)
21
Velocity extensions and Fast Marching Methods

• To start the FMM, we must directly compute the
  value of Fext on grid points in a neighbourhood
  of the interface. We use bicubic interpolation
  and a Newton-Raphson iterative scheme in the
  neighbourhood of the interface to find the
  closest point on the interface to the grid point
  to find the value of Fext.
• Fext at the grid point is set to the value of F
  at the closest point
• Why bicubic?
• Linear and bilinear causes shrinking of the
  interface upon multiple initializations hence
  preserves mass upon multiple reinitialization
• First order accurate FMM in first layer even if
  globally second order gt need second order
  accuracy in the interface definition
• Easily scales up to higher dimensions as opposed
  to linear interpolation (N-cubic interpolants)
• Note this was used for fatigue crack growth
  (Sukumar and Chopp).

22
Key points

• The FMM is about one order of magnitude faster,
  when applicable (monotonic and F is only a
  function of position, not of
• The information about location of the interface
  at a given time must come from an earlier time.
  Information flows from the zero iso-contour to
  positive contours.
• Upwind finite difference
• See also General Upwind Methods by Vladimirsky
  and Sethian and recent work by Chopp to loosen
  requirement on monotonicity

23
Level Set Update
The basic level set evolution equation is a type
of Hamilton-Jacobi equation. We use Upwind
finite differences for approximating
as given below
Update of the level set function is done by
forward Euler in time as given by
24
Upwind Finite Differences
Our finite difference methods must also respect
the directions of the characteristics, then it
follows that at any given grid point, the
difference operator must look in the direction
from which the information is coming
Interface
The direction of the finite differences at points
A, B is indicated in the above figure. This is
the scheme used in computing .
25
Mesh and connectivity
Finite difference grid
The XFEM Engine
The Level Set Engine
Velocity extensions using fast marching
Substrate and velocity potential fields solved
using the XFEM
Level set evolution
The level set function update
Coupling of the XFEM with the level set methods
26
Coupling of the XFEM with the level set methods
27
Computation of the gradient of the velocity
potential
28
Numerical examples

• 1D results interface velocity computation,
  comparison with FD
• 2D results
• One semi-circular biofilm colony
• One semi-circular biofilm colony with
  semi-permeable boundary
• Three semi-circular biofilm colonies
• Wavy slab biofilm

29
1D example
30
XFEM with discontinuous-derivative enrichment vs
Immersed Interface method using finite
differences
1D results for the biofilm problem
Interface speed F, calculated when the height of
the 1D biofilm is 0.2 mm
31
One semi-circular biofilm colony
Water
bacteria
We simulate the growth of a semi-circular
biofilm colony. At t0 the biofilm colony is a
semi-circle of radius 0.01 mm.
32
The interface motion from t 0 to t 44.5 days
33
Substrate Concentration at t 44.5 days, darker
areas are rich in substrate
34
Velocity potential at t 44.5 days, finger
tip-splitting
35
One semi-circular biofilm colonysemi-permeable
substratum
We simulate the growth of a semi-circular biofilm
colony. At t0 the biofilm colony is a
semi-circle of radius 0.01 mm.
36
The interface motion from t 0 to t 7.3 days
37
Substrate Concentration at t 7.3 days, darker
areas are rich in substrate
38
Velocity potential at t 7.3 days, we get a flat
biofilm in this case.
39
Three semi-circular biofilm colonies
We simulate the growth of three semi-circular
biofilm colonies. At t0 the middle colony is a
semi-circle of radius 0.02 mm and other two
colonies are semi-circles of 0.01 mm radius.
40
The interface motion from t 0 to t 28.6 days
41
Substrate Concentration at t 28.6 days, darker
areas are rich in substrate
42
Velocity potential at t 28.6 days, only the
middle biofilm is growing at a much faster rate.
shielding
43
Wavy slab biofilm
We simulate the growth of a slab biofilm of
height 0.01 mm and add a noise term to it as
given by
44
The interface motion from t 0 to t 31 days
45
Substrate Concentration at t 31 days, darker
areas are rich in substrate
46
Velocity potential at t 31 days, we can see
only one of the peaks is growing.
shielding
47
Conclusions

• Nonlinear interior coupled Poissons problem
• Discontinuous derivative enrichment
• Level set and fast marching treatment of
  interface update
• Hybrid finite difference level sets and extended
  finite element meshes
• Quantitatively and qualitatively reproduces
  natural behaviours in biofilms
• Promising for similar classes of problems

48
FUTURE

• 3D
• Fluid-structure interaction
• Complex substratum geometries (e.g. lungs)
• Influence of stress in the biofilm on growth
• Model detachment and reattachment downstream
• Space-time

• Multi-species
• Quorum sensing
• Antibiotics

49
the end
50
Scotland
51
Scotland
52
Scotland
53
(No Transcript)
54
Scotland
55
(No Transcript)
56
(No Transcript)
57
Research projects extended FEM (August 2006)
Crack initiation by loss of ellipticity/hyperboli
city
Void growth and cracks in biomorphic
materials With Svetozara Petrova, University of
Augsburg, Germany.
With Timon Rabczuk,Germany
Dense microfissuration of ageing concrete
With Amor Guidoum and Cyrille Dunant,
EPFL-Material science
Ductile fracture Enrich with HRR fields? Crack
growth criterion?
A-posteriori error estimates and local remeshing
for X-FEM With Le Phong, Ho Chi Minh Ville
Shearography and holography, non-destructive
evaluation With Xavier Schwab, Stuttgart, Germany