Fifth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 10, 2003, pp 123–142. http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

A transient density-dependent diffusion-reaction model for the limitation of antibiotic penetration

in biofilms ^{∗}

Hermann J. Eberl & Messoud A. Efendiev

Abstract A mathematical model is presented that describes the disinfection of microbial biofilms by antibiotics. It is the first multi-species/multi-substrate generalization of a continuous prototype biofilm model comprising degenerating as well as fast diffusion. The boundedness of the model solution is established. The dynamic model behaviour in dependence of model parameters is studied by numerical simulations. A characteristic dimensionless parameter, the disinfection number, is derived, allowing an a priori statement whether all active biomass will be removed or not.

1 Introduction

Biofilm Modeling

Bacterial biofilms are accumulations of microorganisms in aqueous systems that grow attached to interfaces and surfaces. The bacteria are embedded in a polymeric matrix which offers them protection against harmful impact from the environment. Therefore, in a biofilm bacteria live in a protected mode of growth which enhances their ability to survive in hostile environments.

124 A transient density-dependent diffusion-reaction model EJDE/Conf/10

For a long time, dynamic biofilm models have been based on the seminal work of Wanner and Gujer [19]. In their approach, the biofilm will grow only in one direction, perpendicular to the substratum. This describes the formation of a biofilm layer with uniform thickness. The accumulation of biomass in this model is described by a convective mechanism, with the convection velocity directly related to the production of new biomass. However, as is observed in laboratory experiments, using advanced techniques such as confocal laser scanning microscopy (CLSM) or scanning electron micrograph (SEM), in realiter biofilms can grow in complicated spatial architectures with voids and channels, rather than as homogeneous layers. Therefore, in recent years several models have been proposed that also allow for a description of the evolution of spatially heterogeneous biofilm structures. The crucial task in multi-dimensional biofilm modeling is to formulate a spatial biomass spreading mechanism that allows for spatial heterogeneities and which is sensitive to the shape determining factors. In most approaches the biomass is considered a continuum. One can distinguish between models based on discrete, stochastic local sets of rules, e.g. [10, 13, 14, 15, 20], and the more recent fully continuous deterministic models based on partial differential equations [6, 8]. These models either treat biomass density as a constant [6, 10, 13] or as a dependent variable [8, 14, 15]. As an alternative to the continuum ansatz, an individual based model was proposed in [12], in which single bacteria are considered instead of a biomass continuum.

∗Mathematics Subject Classifications: 92D25, 92C50, 35K65.

Key words: Biofilms, antibiotic disinfection, simulation, nonlinear diffusion-reaction.

\begin{matrix}c \\ ◯ \end{matrix}2003 Southwest Texas State University.

Published February 28, 2003.

In this study we will generalize the prototype density-dependent diffusionreaction model [8] to a binary biofilm system with two biomass fractions. It comprises degeneracy (as in the porous medium equation) as well as a singularity of the diffusion coefficient. In [8], only a single-species/single-substrate biofilm was considered. Real biofilms, however, consist of more than one population, and often several dissolved substrates must be taken into account. Since the overall mechanisms leading to mixing or separation of several species in a biofilm are not yet fully understood, no general multi-species biofilm model can be formulated at present. Instead, particular biofilm systems must be modelled separately. As a first example we will consider a binary biofilm with active and inert biomass in the presence of antibiotics and nutrients. A thorough mathematical analysis of a complex model like this is currently possible only in a restricted way, in particular due to occurrence of both non-standard diffusion effects. Therefore, mainly computer simulations are deployed to study the qualitative behaviour of the system. In the current study our main interest lies not in the formation of spatially heterogeneous biofilm architectures but in the interaction of the system’s components, mainly in the spreading mechanism in the presence of more than one population. Therefore, we will restrict ourselves to the one-dimensional case. However, the model can be applied to study spatial effects as well.

Antibiotic Disinfection of Biofilms

Many bacterial infections in the human body are caused by biofilms (see [4] for a list of examples), the treatment of which is rather complicated because the sessile bacteria of a biofilm community are less susceptible to antibiotics than their suspended planktonic conspecifics.

Biologists claim several different mechanisms counteracting antibiotic therapy [4]: The bacteria in a biofilm communicate with each other through a cellto-cell signaling mechanism. This enables them to react to changing situations as a population community instead of individually. A mathematical model for this effect called quorum sensing was suggested and recently discussed in [5]. A second mechanism responsible for reduced susceptibility of biofilms to antibiotics is the failure of antimicrobial agents to fully penetrate the biofilm. In the present study we turn to this aspect. Diffusion of antibiotics in the biofilm is slower than in pure water, and it is also slower than diffusion of many feeding nutrients or oxygen. Furthermore, at the biofilm surface the antibiotics are quickly degraded in biochemical reactions. This hampers the full penetration of the biofilm. The bacteria at the interface will be deactivated, but due to the lack of antibiotics the microorganisms in the deeper layers of the biofilm are not affected. Thus, they can survive and multiply. This was experimentally investigated in [1], e.g. Modelling studies on this effect using inherently one-dimensional biofilm models were reported in [7, 16, 17, 18].

The mathematical description of antibiotic disinfection of biofilms is more complicated than a prototype single-species/single-substrate biofilm system, the consideration of which is sufficient for the modelling of spatially heterogeneous biofilm architectures, even if many biofilms in the human body are formed by one single bacterial species: Not only the living organisms must be taken into account but also inert biomass as the product of the disinfection process. Moreover, of course, besides the nutrients (including oxygen), the antibiotic itself must be considered as a dissolved component. Hence, with active and inert biomass, we have a binary biofilm. Previous experimental findings and results obtained with classical modelling techniques allow for a comparison of our model’s behaviour with established results. Therefore, the antibiotic disinfection regime is an ideal candidate to test the applicability of the new biomass spreading mechanism that was originally formulated for a mono-species model of a biofilm system with more than one component.

2 The Basic Biofilm Model

A density-dependent diffusion-reaction model for the formation of spatially heterogeneous single-species/single-substrate biofilms was introduced in [8] for the dependent variables substrate concentration C and biomass density M. The key model features are: (i) there is a sharp boundary between biofilm and bulk liquid, (ii) biomass spreading only takes place if the biomass density approaches a prescribed maximum value which establishes an upper bound, and (iii) new biomass is produced locally as long as there are nutrients available to feed bacteria. The system of evolution equations reads

\begin{matrix}C_{t} = ∇_{x} · (D_{C}∇_{x}C) − \frac{γCM}{K + C} \\ M_{t} = ∇_{x} · (D_{M}(M)∇_{x}M) + \frac{ξ_{1}CM}{K + C} − ξ_{3}M \end{matrix}  (2.1)

where the density dependent biomass diffusion coefficient is given by

D_{M}(M) = d_{M} \frac{M^{β}}{(M_{max} − M)^{α}} , α, β > 1, d_{M} > 0  (2.2)

The parameters in (2.1) are the nutrient diffusion constant D_{C}, the specific maximum consumption rate γ, the Monod half-saturation constant K and the specific maximum biomass growth rate ξ_{1}, which all are positive. The biomass decay rate ξ_{3} is non-negative. The model is considered in a bounded domain Ω ⊂ ℝ^{d}, d = 1, 2, 3. The standard boundary conditions considered are mixed Neumann/Dirichlet conditions.

Figure 1: Formation of a biofilm colony in time according to (2.1): First a wavy layer develops. After nutrients get limited inside the biofilm the bigger colonies start to dominate over the smaller ones and mushroom type structures develop (from [8]).

In this model, the distinction between the actual biofilm and the surrounding liquid phase is given by the biomass density. The liquid region is described by Ω_{1}(t) = {x ∈ Ω | M(x, t) = 0}. The biofilm itself is the region Ω_{2}(t) = {x ∈ Ω | M(x, t) > 0}. Both regions vary in time due to the evolution of the biomass density. The degeneracy D_{M}(0) = 0 is responsible for model properties (i) and the first part of (ii) as in the porous medium equation. The singularity D_{M}(M) → ∞ as M → M_{max} was introduced in order to guarantee the second part of property (ii). Property (iii) is satisfied due to the reaction kinetics applied. In [8] it was demonstrated by computer simulations that model (2.1) is able to describe the formation of spatially heterogeneous biofilm structures, such as experimentally observed mushroom-type colonies, see Figure 1. The heterogeneity of the predicted biofilm architectures is sensitive to the initial and boundary conditions as well as to a dimensionless parameter G that relates biomass production and nutrient availability. In [9], existence, uniqueness, and boundedness of the model solution were proven, in particular it was shown that M indeed obeys the upper bound M_{max}.

3 A Model For Antibiotic Disinfection

Governing Equations

The prototype biofilm model (2.1) only describes mono-species biofilms. In order to model spatio-temporal processes of biofilms with more than one biomass fraction, it must be adapted and generalized in an appropriate manner. A simple example is antibiotic disinfection of biofilms, for which inert biomass must be taken into account in addition to viable biomass. From the modelling point of view, this is treated as a second population and a binary biofilm is obtained. The basic setup for the disinfection process follows the model developed in [16]:

The dissolved substrates B and C are transported due to diffusion in bulk and biofilm, and consumed by active biomass in the biofilm. The disinfection process is considered as a one-to-one conversion of viable biomass into inert biomass, depending on the availability of antibiotics. Viable biomass X grows, as long as nutrients are available, by consumption as in (2.1) and it decays due to disinfection with antibiotics and due to natural death. The inert biomass Y grows due to disinfection. The reaction terms are first order in B and X, and they depend on C according to Monod kinetics. As in [16], we assume that live and dead cells are moved together and we put this into our model framework. Therefore, in our model the nonlinear biomass diffusion coefficient is the same for X and Y and depends on the total density of both fractions. The upper bound M_{max} must be obeyed simultaneously by both viable and inert biomass. Similarly, spatial spreading occurs only when this maximum value is approached by the total biomass density. The resulting biofilm model reads

\begin{matrix}B_{t} = & ∇_{x} · (D_{B}∇_{x}B) − βBX \\ C_{t} = & ∇_{x} · (D_{C}∇_{x}C) − \frac{γXC}{k + C} \\ X_{t} = & ∇_{x} · (D_{M}(X + Y)∇_{x}X) + \frac{ξ_{1}XC}{k + C} − ξ_{2}BX − ξ_{3}X \\ Y_{t} = & ∇_{x} · (D_{M}(X + Y)∇_{x}Y) + ξ_{2}BX \end{matrix}  (3.1)

where oxygen is considered to be the limiting substrate C. The biomass diffusion coefficient D_{M}(M) in (3.1) is given by (2.2). Implicit in this model formulation is the assumption that viable and dead biomass have the same maximum cell density M_{max}. Model equations (3.1) are valid in the domain Ω ⊂ ℝ^{d}, d = 1, 2, 3. All newly introduced parameters are positive. In (3.1), β and ξ_{2} are first order reaction constants describing the interaction of antibiotics and biomass. D_{B} is the diffusion coefficient of the antibiotics. The remaining parameters of (3.1) are as in (2.1).

In order to realistically model antibiotic penetration into the biofilm one must take into account that the diffusion process in the solid biofilm Ω_{2} is slower than in the surrounding liquid phase Ω_{1}. This induces the further model generalization

D_{B}(x) = \begin{cases}D_{B,1} & for x ∈ Ω_{1} \\ D_{B,2} & for x ∈ Ω_{2} \end{cases}  (3.2)

where D_{B,1} and D_{B,2} are constants with τ_{B} := D_{B,2}/D_{B,1} ≤ 1. The same applies for the diffusion coefficient of substrate concentration C. The ratios τ of biofilm and liquid diffusivity depend on the size and molecular weight of the molecules (cf. [2, 18]). For various antibiotics, [18] gives typical values around τ_{B} ≈ 0.8. For small molecules like oxygen one has τ_{C} in a good approximation in the range 0.95 ≤ τ_{C} ≤ 1. Moreover in the case of oxygen and antibiotics, it holds that D_{C} ≥ D_{B}.

The evolution equations given above are completed by appropriate initial and boundary conditions, connecting the system Ω with the external world.

Properties of the Disinfection Model

Boundedness of the Solution

In [9] the boundedness of the biomass density M was proven for the monospecies model (2.1) with parameters satisfying \frac{ξ_{1}C_{0}}{K+C_{0}} − ξ_{3} > 0. This latter condition is no severe restriction for practical considerations since otherwise the production of new active biomass would be impossible. For mixed homogeneous Dirichlet/Neumann boundary conditions for the biomass density, one obtains M < M_{max} almost everywhere.

The boundedness of the total biomass density X + Y < M_{max} (almost everywhere) of the antibiotic model (3.1) can be obtained as a corollary of this theorem as follows.

Proposition 3.1 In a bounded domain Ω, the solution of the disinfection model (3.1) with non-negative parameters satisfying \frac{ξ_{1}C_{0}}{K+C_{0}} − ξ_{3} > 0, boundary conditions

B(t, x) = B_{0} for x ∈ ∂Ω_{DB}, \frac{∂B}{∂n} = 0 for x ∈ ∂Ω_{NB}
C(t, x) = C_{0} for x ∈ ∂Ω_{DC}, \frac{∂C}{∂n} = 0 for x ∈ ∂Ω_{NC}
X(t, x) = 0 for x ∈ ∂Ω_{DX}, \frac{∂X}{∂n} = 0 for x ∈ ∂Ω_{NX}
Y (t, x) = 0 for x ∈ ∂Ω_{DX}, \frac{∂Y}{∂n} = 0 for x ∈ ∂Ω_{NX}

with

∂Ω_{NB} ∪ ∂Ω_{DB} = ∂Ω, ∂Ω_{NB} ∩ ∂Ω_{DB} = ∅, Ω_{DB} ≠ ∅
∂Ω_{NC} ∪ ∂Ω_{DC} = ∂Ω, ∂Ω_{NC} ∩ ∂Ω_{DC} = ∅, Ω_{DC} ≠ ∅
∂Ω_{NX} ∪ ∂Ω_{DX} = ∂Ω, ∂Ω_{NX} ∩ ∂Ω_{DX} = ∅, Ω_{DX} ≠ ∅

and initial data with

C(0, x) ∈ L^{∞}(Ω) ∩ H^{1}(Ω), 0 ≤ C(0, x) ≤ C_{0}
B(0, x) ∈ L^{∞}(Ω) ∩ H^{1}(Ω), 0 ≤ B(0, x) ≤ B_{0}

and

X(0, x), Y (0, x) ∈ L^{1}(Ω), 0 ≤ X(0, x), 0 ≤ Y (0, x),
X(0, x) + Y (0, x) < M_{max}

satisfies

0 ≤ C(t, x) ≤ C_{0}, 0 ≤ B(t, x) ≤ B_{0},

and

X + Y < M_{max} almost everywhere

Outline of proof The following steps are taken:

i) For the solutes B and C the assertion follows from the maximum principle.

ii) X and Y are non-negative by comparison and invariance theorems.

iii) Adding the evolution equations for X and Y in (3.1) yields

(X + Y)_{t} = ∇_{x} · (D_{M}(X + Y)∇_{x}(X + Y)) + X (\frac{γC}{k + C} − ξ_{3})

Since

0 ≤ X (\frac{γC}{k + C} − ξ_{3}) ≤ (X + Y) (\frac{γC}{k + C} − ξ_{3})

it follows that the solution M of the corresponding boundary value problem (2.1) is an upper bound for X+Y by monotonicity and comparison arguments. Hence, ∥X + Y ∥_{L^{∞}(ℝ_{+}×Ω)} ≤ ∥M∥_{L^{∞}(ℝ_{+}×Ω)} < M_{max} follows by the main theorem proven in [9].

Note that in many applications Y (0, x) ≡ 0 everywhere in Ω. Then, the initial conditions for X are the same as the ones for the total biomass M in [9].

Similarly, the finite speed of the propagation of the interface between Ω_{1} and Ω_{2} is established. That is, for small X + Y the diffusion operator in (3.1) behaves like the degenerating density-dependent diffusion operator of the porous medium equation.

Initial Penetration of Solutes

Since the transport and reaction processes governing the evolution of substrates are much faster than the evolution and formation of biomass (cf. [11]), a quasisteady state scenario is investigated under the assumption of frozen biomass X and Y. The substrates C and B relax to an equilibrium very quickly. This difference in characteristic time-scales is the reason why the actual initial conditions for B and C in (3.1) are not of particular importance.

At initial time t = 0 of the disinfection process, we have constant X ≡ X_{0} in Ω_{2} and X ≡ 0 in Ω_{1}, and Y ≡ 0 everywhere. In this case, the equations for C and B decouple. In a one-dimensional setup the steady states of B satisfy a linear second order ordinary differential equation. After rescaling the diffusionreaction equation, the spatially one-dimensional equation reads in dimensionless form

\begin{matrix}0 = b_{xx} − θ_{b}^{2}b for x ≤ λ := L_{f}/L_{z} \\ 0 = b_{xx} for λ ≤ x ≤ 1 \\ b_{x}(0) = 0, b(1) = b_{0} \end{matrix}  (3.3)

with θ_{b}^{2} = βX_{0}L_{f}^{2}/D_{B,2}. In the chemical engineering literature, the dimensionless variable θ_{b} is called the Thiele number. At the interface of Ω_{1} and Ω_{2}, the antibiotic concentration b and its diffusive flux are continuous. These two internal boundary conditions allow for a unique solution of (3.3) that reads

b(x) = \begin{cases}b(λ)\frac{cosh(θ_{b}x)}{cosh(λθ_{b})} & for x ≤ λ \\ b(λ) + \frac{x−λ}{1−λ} (b_{0} − b(λ)) & for λ ≤ x ≤ 1 \end{cases}  (3.4)

The concentration b(λ) at the interface of Ω_{1} and Ω_{2} can be eliminated using the internal boundary conditions described above. This yields

b(λ) = \frac{b_{0}}{1 + (1 − λ)τ_{b}θ_{b} tanh(λθ_{b})}  (3.5)

Solution (3.4) shows how the penetration of antibiotics in the biofilm is hampered by the diffusion-reaction mechanism quantified by the Thiele number. A detailed investigation of this effect was carried out in [17, 18].

In the sequel, we will derive a condition for the net biomass production preventing antibiotic disinfection, based on the analytical solution (3.4). An analytical solution of the corresponding two-point boundary value problem for C is not available. However, in cases where the concentration at the biofilm interface is small, i.e. c(λ) ≪ k (this is the case in some interesting applications, see the parameters used in [16]), the Monod term can be approximated by a first order reaction and one obtains a similar equation as for B, which reads

\begin{matrix}0 = c_{xx} − θ_{c}^{2}c for x ≤ λ, \\ 0 = c_{xx}for λ ≤ x ≤ 1, \\ c_{x}(0) = 0, c(1) = c_{0} \end{matrix}  (3.6)

with θ_{c}^{2} = γX_{0}L_{f}^{2}/(kD_{C,2}). The solution of (3.6) is analogous to (3.4).

If the antibiotic concentration is sufficiently high to dampen the formation of new viable biomass, the biomass production term in the evolution equation (3.1) for X must be dominated by the loss term, i.e.

∫_{0}^{λ} \frac{ξ_{1}}{k} cdz ≤ ∫_{0}^{λ} (ξ_{2}b + ξ_{3}) dz  (3.7)

Upon substituting the analytical solutions into (3.7), one obtains the following condition for antibiotic disinfection that prevents the net production of active biomass

1 ≤ D := \frac{kθ_{c}}{θ_{b}} · \frac{b_{0}ξ_{2} tanh(λθ_{b}) + ξ_{3}λθ_{b}}{c_{0}ξ_{1} tanh(λθ_{c})} · \frac{1 + (1 − λ)τ_{c}θ_{c} tanh(λθ_{c})}{1 + (1 − λ)τ_{b}θ_{b} tanh(λθ_{b})}  (3.8)

If D < 1, the production of new viable biomass is faster than the disinfection process.

One-dimensional Dynamic Simulations

Boundary Conditions and Numerical Scheme

As pointed out earlier, a complete analytical treatment of the model equation (3.1) is difficult. Therefore, we carry out computer simulations to study the model behaviour, restricting ourselves to the one-dimensional case in this paper. This, of course, can not provide information about spatially heterogeneous biofilm structures, but it can demonstrate the interaction of the four components of the system, antibiotics, substrates, active and inert biomass, which is the primary purpose of this study.

For our simulations, at the substratum (x = 0) homogeneous Neumann conditions are specified for all components, and Dirichlet conditions at the other end of the domain (x = L). The boundary values for the dissolved substrates are positive and describe unrestricted availability in the bulk. For the biomass components, the density is kept at 0 for x = L. Thus, X +Y < M_{max} is guaranteed by Proposition 3.1. Initially (at t = 0) no inert biomaterial is in the system, and a homogeneous biofilm is considered with a constant distribution of active biomass close to the maximum value. This agrees in very good approximation with the results obtained in [8] for a single-substrate/single-species biofilm. This setup describes the disinfection process after a homogeneous biofilm with thickness L_{f} has developed. The initial conditions prescribed for the dissolved substrates B and C are not critical, as noted above since they reach a “quasiequilibrium” instantaneously, which is described by the analytical solution (3.4) for B. Thus, the initial and boundary conditions read

\begin{matrix}B(x, 0) = B_{0}, C(x, 0) = C_{0}, Y (x, 0) = 0 \\ X(x, 0) = \begin{cases}X_{0} & if x ≤ L_{f} \\ 0 & if x > L_{f}, \end{cases} \end{matrix}  (3.9)

\begin{matrix}B_{x}(0, t) = 0, C_{x}(0, t) = 0, X_{x}(0, t)) = 0, Y_{x}(0, t) = 0 \\ B(L, t) = B_{0}, C(L, t) = C_{0}, X(L, t) = 0, Y (L, t) = 0 , \end{matrix}  (3.10)

where X_{0} < M_{max} is a positive constant and L_{f} is the initial thickness of the biofilm. Thus one has Ω_{1}(0) = (L_{f}, L] and Ω_{2}(0) = [0, L_{f}].

Since the characteristic time scales of biofilm formation/disinfection (i.e. the evolution equations for X and Y) and for nutrient consumption (i.e. the evolution equations for B and C) differ by several orders of magnitude, two different approaches are made for the time integration of (3.1): An explicit 4th order Runge-Kutta method is deployed for the slower X and Y, while an implicit backward Euler scheme is used for the faster processes describing B and C. Finite volumes are used for space discretization. The nonlinear system in the calculation of B and C is solved with a restarted Newton method.

Illustration of Model Behavior

To illustrate the behavior of model (3.1), we carry out several simulations with changing model parameters. In the focus of our interest is the mixing of viable and dead biomass which is governed by the spatial biomass spreading mechanism at the core of the evolution equations for X and Y. Biofilm disinfection is an interaction of several processes, such as nutrient utilization and production of new biomass, antibiotic consumption, transfer of living biomass into dead biomass, and transport of soluble substrates in the biofilm by diffusion. Therefore, it can be expected that the qualitative behavior of the model depends strongly on how these processes balance. In the sequel, it will be investigated how the performance of antibiotic disinfection depends on the environmental conditions, using numerical simulations with different model parameters. The parameters to be varied are the bulk antibiotic concentration B_{0}, the initial biofilm thickness L_{f}, the inital biomass density X_{0}, and the supply of nutrients expressed by the bulk nutrient concentration C_{0}. The other model parameters are kept constant throughout all simulations in this section. The model parameters are taken from [16, 18], in order to allow a comparison of results. The variations of these default parameters are given in the text with the simulations. In all cases in this section λ = L_{F} /L = 0.91 at the beginning. Thus a thin concentration boundary layer with the same relative thickness for all cases is prescribed

Basic process description: As soon as antibiotic B is added to the system and transported from the boundary to the biofilm, disinfection starts and inert biomass is produced. The depth to which the antibiotics can penetrate into the biofilm and eradicate active biomass depends on factors such as utilization rates, diffusion coefficients, and bulk antibiotic concentration. In any case, disinfection starts at the biofilm/bulk interface and propagates to the substratum. Due to decay of antibiotics by viable biomass, the disinfection process is slowed down inside the biofilm. On the other hand, as long as oxygen C and living bacteria X are available, new biomass is produced. The biomass newly formed inside the biofilm pushes the inert biomass at the interface and the biofilm may even grow. Some typical simulations are shown in Figures 2, 3, and 4.

Variation of antibiotic bulk concentration: With increasing B_{0}, the production of inert biomass accelerates due to enhanced availability of antibiotics. In cases of very low B_{0}, the thickness of the biofilm may increase due to formation of new active biomass. This is shown in Figure 2. In particular, for the lowest test bulk antibiotic concentration (a), virtually no disinfection takes place due to antibiotic limitation. In the case (b) with increased B_{0}, inert biomass forms at the biofilm/bulk interface. However, due to the decay of antibiotics at the interface, they get limited in the interior of the biofilm, where no disinfection takes place for a long time (several days). In the case (c) with further increased B_{0}, the active biomass at the biofilm/bulk interface is immediately converted into dead biomass, the antibiotic does not penetrate into the biofilm fast enough. Thus, viable biomass exists at the substratum for some time. Only in the case (d) of a very high bulk antibiotic concentration, the entire biofilm is converted into inert biomass within a short time (less than half a day). In the presented simulations, L_{f} = 4.55 · 10^{−4}m, C_{0} = 0.035gm^{−3}, X_{0} = 0.95gm^{−3}. As the simulations demonstrate, there is no instantaneous spreading of inert biomass inside the biofilm but rather the propagation of an inert biomass wave. This shows that the spatio-temporal biomass spreading mechanism in our model is able to keep populations in a biofilm separated if the reaction terms allow this.

Figure 2: Production of inert biomass for different antibiotic bulk concentrations B_{0} in time. Shown is inert biomass density Y relative to maximum biomass density M_{Max}: (a) B_{0} = 0.1gm^{−3}, (b) B_{0} = 10gm^{−3}, (c) B_{0} = 100gm^{−3}, (d) B_{0} = 500gm^{−3}

This is one of the crucial properties of a reliable biofilm model. The profiles of X for a fixed t are the same as those shown in [16].

Variation of initial biofilm thickness: In a second experiment, we investigate the disinfection process in dependence of the initial biofilm thickness. Figure 3 shows simulations for B_{0} = 100gm^{−3}, C_{0} = 0.035gm^{−3}, X_{0} = 0.95gm^{−3}. Only for thin biofilms, i.e. case (a), the whole biofilm is disinfected quickly. With increasing initial biofilm thickness, it becomes more and more difficult for the antibiotics to penetrate the whole biofilm and to transform active biomass into inert biomass. In particular for the thickest example (c), over a long period of time virtually no disinfection takes place at the substratum. The time which is needed for full disinfection depends nonlinearly on the initial biofilm thickness.

Variation of initial biomass density and nutrient supply: In the last experiment we shall investigate how the the model (3.1) depends on the bulk oxygen concentration and on the initial biomass density. The first factor actively controls the growth of the biofilm. The second one is of interest since the underlying biofilm model (2.1) has the property that spatial spreading of biomass only takes place when the total biomass density approaches its maximum value Figure 3: Disinfection of biofilms with different inital thicknesses L_{f}. Plotted is the inert biomass density Y relative to the maximum biomass density M_{max}:

(a) L_{f} = 9.1 · 10^{−5}m, (b) L_{f} = 4.55 · 10^{−4}m, (c) L_{f} = 6.36 · 10^{−4}m

M_{max}. Thus, in this last illustration, it will be investigated how the model describes the growth of the binary biofilm under consideration. As an example, Figure 4, compares the active biomass X and the total biomass M = X + Y for different values of X_{0} and C_{0}. The bulk antibiotic concentration in these simulations is B_{0} = 10gm^{−3}, the initial biofilm thickness is L_{f} = 9.1·10^{−5}m. In neither of the cases with increased C_{0} in the bulk the oxygen limitation inside the biofilm could be overcome. Thus, a noteworthy amount of new biomass is produced only close to the interface of Ω_{1} and Ω_{2}. In case (b) with a higher initial biomass density, the new biomass leads immediately to an increase of the biofilm thickness. In case (a) with a lower inital biomass density, it propagates into the biofilm leading to a compression of the total biomass by pushing both the viable and the inert population. The growth of the biofilm starts only when the total biomass density approaches M_{max} at the interface. Again this is in agreement with the postulation that lead to the formulation of the underlying biofilm formation model (2.1) in [8]. In the cases (c) and (d) with lower C_{0}, much less new biomass is produced and the biofilm thickness remains virtually constant over the time interval shown. Under lower bulk nutrient concentrations, the active biomass density is smaller for higher C_{0}, i.e. disinfection is faster. Thus, increased active biomass production counteracts the conversion into inert biomass.

The qualitative model behavior in all the simulations shown above agrees Figure 4: Biomass growth and disinfection for various C_{0} and X_{0}. Plotted are the total biomass density M = X +Y and the viable biomass density X, relative to M_{max}: (a) C_{0} = 0.1gm^{−3}, X_{0} = 0.85M_{0}, (b) C_{0} = 0.1gm^{−3}, X_{0} = 0.95M_{max},

(c) C_{0} = 0.0035gm^{−3}, X_{0} = 0.85M_{max}, (d) C_{0} = 0.0035gm^{−3}, X_{0} = 0.95M_{max}

with the results published in [16, 17, 18]. Therefore, it may be concluded that the spatial biomass spreading mechanism introduced in (3.1) is able to reliably simulate binary biofilms like the system of viable and inert biomass in this study.

Growth vs. Decay of Active Biomass

Based on a frozen steady state assumption for slowly growing biomass, criterion (3.8) was derived, which allows an a priori statement whether at initial time net production of viable biomass takes place or whether the disinfection process is faster than biomass growth. This criterion was based on the assumption of a constant initial distribution of biomass and on the analytic steady state solution (3.4) for solutes. In this section the question will be investigated whether this criterion applies to transient processes as well. Thus, the question is, under which circumstances can the disinfection number D be used as a parameter to a priori decide about the qualitaive behaviour of the model over time, that is, extinction or growth of active biomass in the system.

For this purpose many computer simulations were run covering a broad range of model parameter sets, and hence values of D. The results are illustrated in Figures 5 and 6. Shown is the total active and inert biomass, relative to the active biomass at starting time. That is, ∫_{0}^{L} X(x, t)dt/(X_{0}L_{f}) and 136 A transient density-dependent diffusion-reaction model EJDE/Conf/10 Figure 5: Total active and inert biomass in the biofilm, relative to the total biomass at initial time for several disinfection numbers D.

Figure 6: Figure 5 continued ∫_{0}^{L} Y (x, t)dt/(X_{0}L_{f}).

For D ≫ 1, the rate of disinfection is much bigger than the production rate of new viable biomass. The total active biomass shrinks and the total inert biomass increases. Eventually this leads to total disinfection, that is, no active biomass remains. In the opposite case D ≪ 1, the production of new biomass outweighs conversion into dead material. That is, the total viable biomass grows. This leads to an increase of the biofilm thickness. For very low D, the disinfection rate appears negligibly small. For D ≈ 1 it is observed that the total amount of active biomass remains almost constant over a long period of time, while the inert biomass may increase. That is, the criterion (3.8) for disinfection applies to the transient case as well.

With one exception, all simulations presented here were carried out with parameters chosen such that C(λ) ≪ k. This was one of the assumptions made for the derivation of the disinfection criterion (3.8). Only in Figure 6b with D = 1 this condition was not met. While in Figure 6a (also with D = 1) the total amount of viable biomass remains unchanged in the course of the simulation, in Figure 6b disinfection starts. This effect is generally observed: When C(λ) ≪ k does not hold, the disinfection process is faster than for the same D number with lower bulk concentrations. Indeed, disinfection might start already for D < 1. The reason for this effect is that the reaction rate \frac{ξ_{1}c}{k+c} of the Monod kinetics is bounded while the first order reaction ξ_{1}c/k is not bounded. Therefore, Monod kinetics lead to a lower biomass production and D underestimates the disinfection process. This can be explained by introducing a disinfection number D_{Monod} based on the solution corresponding to the nonlinear reaction term in (3.1) in the same manner as D is based on the solution of the first order equation (3.6). Then the following holds.

Proposition 3.2 The Monod-based disinfection number D_{Monod} and the disinfection number D are related by

D_{Monod} ≥ D

Outline of proof: Let C be the solution of the steady state diffusion reaction model with first order kinetics and C̃ the solution of the equation with Monod kinetics. The following steps are taken:

i) The Monod disinfection number is the ratio of the integral I_{b} of the decay terms and the integral of the production term in (3.1):

D_{Monod} = I_{b}/I_{c,M}

with

I_{c,M} = ξ_{1} ∫_{0}^{L_{f}} \frac{C̃}{k + C̃} dx = \frac{D_{C,2}ξ_{1}}{γ} ∫_{0}^{L_{f}} C̃′′dx = \frac{D_{C,2}ξ_{1}}{γ} [C̃′(L_{f}) − C̃′(0)]

Substituting the boundary and interface conditions into this expression yields

I_{c,M} = \frac{D_{C,1}ξ_{1}}{γ(L − L_{f})}[C_{0} − C̃(L_{f})]

In the case of first order kinetics one obtains by the same argument

D = \frac{I_{b}}{I_{c}} with I_{c} = \frac{D_{C,1}ξ_{1}}{γ(L − L_{f})}[C_{0} − C(L_{f})]

I_{b} is a constant which is independent of C or C̃. ii) Using the interface conditions at x = L_{f} and the Dirichlet condition at x = L, one obtains smooth two-point boundary value problems for C and C̃ in the interval [0, L_{f}]

0 = D_{C,2}C′′ − \frac{γC}{k} , C′(0) = 0, \frac{C(L_{f})}{L − L_{f}} + τ_{c}C′(L_{f}) = \frac{C_{0}}{L − L_{f}}  (3.11)
0 = D_{C,2}C̃′′ − \frac{γC̃}{k + C̃} , C̃′(0) = 0, \frac{C̃(L_{f})}{L − L_{f}} + τ_{c}C̃′(L_{f}) = \frac{C_{0}}{L − L_{f}} (3.12)

¿From \frac{U}{k+U} ≤ \frac{U}{k} for U ≥ 0 it follows that

D_{C,2}C′′ − \frac{γC}{k + C} ≥ 0

With Theorem 3.4.1 of [3] one obtains C(x) ≤ C̃(x) for all x ∈ [0, L_{f}], in particular C(L_{f}) ≤ C̃(L_{f}). iii) Combining i) and ii) yields the assertion. □

Note that the evaluation of D_{Monod} requires the quadrature of a function of the solution of (3.12). Since no closed form for C̃ is available, it can not easily be used as an a priori defined parameter in order to determine if disinfection or biofilm growth takes place. However, if the first order reaction is not a valid approximation of the Monod reaction, we still have that D > 1 will lead to disinfection, even if the reverse criterion that D < 1 leads to growth of the biofilm does not hold anymore.

4 Conclusion

A fully continuous biofilm formation mechanism which had been introduced in [8] as a mono-species/mono-substrate prototype model was applied to a biofilm system with two species. The system describes antibiotic disinfection of biofilms, where in addition to viable biomass inert biomass must be taken into account. This simple system has the adavantage that it was extensively studied experimentally and numerically before. Thus, it offers a good opportunity to compare the results of the new model to previously published results.

The model is a system of diffusion-reaction equations which comprises degeneracy as well as fast diffusion. This is the reason why a full analytical treatment of the model is currently not possible. However, a statement about the existence, uniqueness, and boundedness of the model solution could be derived. The qualitative dynamic behaviour of the model was investigated by computer simulations. Although the simulations in the present study were restricted to the one-dimensional case in order to allow for a comparison with classical inherently one-dimensional model studies, the model is able to describe the general three-dimensional case of spatially heterogeneous biofilms, as well. This will be demonstrated in a forthcoming article. The new biofilm formation model was found capable of describing the disinfection process of bacterial biofilms and the interaction of viable and inert biomass, such as local coexistence or extinction. Moreover, it was possible to establish an a priori criterion for the extinction of active biomass, based on a dimensionless number combining model and system parameters.

The spatial biomass spreading mechanism adapted for the biofilm disinfection process was based only on the assumption that inert biomass locally coexists with viable biomass and that both fractions are moved together. This is the case for other biofilm processes as well, like the production of extracellular polymeric substances (EPS). Therefore, the model suggested here is not restricted to antibiotic disinfection but can be applied to further situations as well.

Acknowledgement: The authors are grateful to Professor Sebastian Walcher for having many fruitful discussions with us.

References

[1] Anderl J.N., Franklin M.J., Stewart P.S., Role of antibiotic penetration limitation in Klebsiella pneumoniae biofilm resistance to ampicillin and ciprofloxacin, Antimicrobial Agents & Chemotherapy, 44(7):1818-1824, 2000.

[2] Bryers J.D., Drummond F., Local macromolecule diffusion coefficients in structurally non-uniform bacterial biofilms using fluorescence recovery after photobleaching (FRAP). Biotechnology & Bioengineering 60(4):462-473, 1998.

[3] Carl S., Heikkilä S., Nonlinear Differential Equations in Ordered Spaces, Chapman & Hall CRC, Boca Raton - London - Washingtom DC, 2000.

[4] Costerton J.W., Stewart PS, Greenberg E.P., Bacterial biofilms: A common cause of persistent infections, Science, Vol 284, No 5418, Issue of 21 May, 1999.

[5] Dockery J.D., Keener J.P., A mathematical model for quorum sensing in Pseudomonas aeruginosa, Bulletin of Mathematical Biology, 63:95-116, 2001.

[6] Dockery J., Klapper I., Finger formation in biofilm layers, SIAM J. Applied Mathematics 62(3):853-869, 2001.

[7] Dodds M.G., Grobe K.J., Stewart P.S., Modeling Biofilm Antimicrobial Resistance, Biotechnology & Bioengineering, 68(4):456-465, 2000.

[8] Eberl H.J., Parker D.F., van Loosdrecht M.C.M., A new deterministic spatio-temporal continuum model for biofilm development, J. of Theoretical Medicine, 3(3):161-176, 2001.

[9] Efendiev M.A., Eberl H.J., Zelik S.V., Existence and Longtime behavior of solutions of a nonlinear reaction-diffusion system arising in the modeling of biofilms, Nonlinear Diffusive Systems and Related Topics, RIMS Kyoto, pp 49-71, 2002.

[10] Hermanowicz S.W., A simple 2D biofilm model yields a variety of morphological features, Mathematical Biosciences 169:1-14, 2001.

[11] Kissel J.C., McCarty P.L., Street R.L., Numerical simualtions of mixedculture biofilms, J. of Environmental Engineering, 110(2):393-411, 1984.

[12] Kreft J.-U., Booth G., Wimpenny J.W.T., BacSim, a simulator for individual-based modelling of bacterial colony growth, Microbiology, 144:3275-3287, 1998.

[13] Mehl M., Ein interdisziplinärer Ansatz zur dreidimensionalen numerischen Simulation von Strömng, Stofftransport und Wachstum in Biofilmsystemen auf der Mikroskala, PhD-Thesis Munich Univ. Technol., Faculty of Computer Science, 2001 (in German).

[14] Noguera D.R., Pizarro G., Stahl D.A., Rittmann B.E., Simulation of multispecies biofilm development in three dimensions, Water Sci. & Technology, 39(7):123-130, 1999.

[15] Picioreanu C., van Loosdrecht M.C.M., Heijnen J.J. A new combined differential-discrete cellular automaton approach for biofilm modeling: Application for growth in gel beads, Biotechnology & Bioengineering, 57(6):718–731, 1998.

[16] Stewart P.S., Biofilm accumulation model that predicts antibiotic resistance of pseudomonas aerginosa biofilms, Antimicrobial Agents & Chemotherapy, 38(5):1052–1058, 1994.

[17] Stewart P.S., Raquepas J.B., Implications of reaction-diffusion theory for the disinfection of microbial biofilms by reactive antimicrobial agents, Chemical Engineering Science, 50(19):3099–3104, 1995.

[18] Stewart P.S., Theoretical aspects of antibiotic diffusion into microbial miofilms, Antimicrobial Agents & Chemotherapy, 40(11):2517–2522, 1996.

[19] Wanner O., Gujer W., A multispecies biofilm model, Biotechnology and Bioengineering, 28:314–328, 1986.

[20] Wimpenny J.W.T., Colasanti R., A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models, FEMS Microbiology Ecology, 22:1-16, 1997.

Hermann J. Eberl Institute of Biomathematics and Biometry GSF – National Research Centre for Environment and Health Pf. 1129, D-85758 Neuherberg, Germany Dept. Mathematics and Statistics, University of Guelph Guelph, On, N1G2W1, Canada e-mail: eberl@gsf.de and heberl@uoguelph.ca fon: +1(519) 824 4120 ext. 52622

Messoud A. Efendiev Institute of Mathematics A, University of Stuttgart Pfaffenwaldring 57, D-70569 Stuttgart, Germany e-mail: efendiev@mathematik.uni-stuttgart.de