\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small Sixth Mississippi State Conference on Differential Equations and Computational Simulations, {\em Electronic Journal of Differential Equations}, Conference 15 (2007), pp. 77--95.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2007 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \setcounter{page}{77} \title[\hfilneg EJDE-2006/Conf/15\hfil A finite difference scheme] {A finite difference scheme for a degenerated diffusion equation arising in microbial ecology} \author[H. J. Eberl, L. Demaret \hfil EJDE/Conf/15 \hfilneg] {Hermann J. Eberl, Laurent Demaret} % in alphabetical order \address{Hermann J. Eberl \newline Dept. Mathematics and Statistics \\ University of Guelph, Guelph, On, Canada} \email{heberl@uoguelph.ca} \address{Laurent Demaret \newline Inst. Biomathematics and Biometry \\ GSF - National Research Centre for Environment and Health \\ Neuherberg, Germany} \email{laurent.demaret@gsf.de} \thanks{Published February 28, 2007.} \subjclass[2000]{35K65, 65M06, 92C17} \keywords{Finite differences; nonlinear diffusion; non-local representation; \hfill\break\indent non-standard discretisation; numerical simulation; biofilm} \begin{abstract} A finite difference scheme is presented for a density-dependent diffusion equation that arises in the mathematical modelling of bacterial biofilms. The peculiarity of the underlying model is that it shows degeneracy as the dependent variable vanishes, as well as a singularity as the dependent variable approaches its a priori known upper bound. The first property leads to a finite speed of interface propagation if the initial data have compact support, while the second one introduces counter-acting super diffusion. This squeezing property of this model leads to steep gradients at the interface. Moving interface problems of this kind are known to be problematic for classical numerical methods and introduce non-physical and non-mathematical solutions. The proposed method is developed to address this observation. The central idea is a non-local (in time) representation of the diffusion operator. It can be shown that the proposed method is free of oscillations at the interface, that the discrete interface satisfies a discrete version of the continuous interface condition and that the effect of interface smearing is quantitatively small. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Corollary}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{remark}[theorem]{Remark} \section{Introduction} \subsection{Biological background} The mathematical model in the focus of this study describes the formation and growth of bacterial biofilms. These are communities of migroorganisms that live embedded in a self-produced layer of extracellular polymeric substances ({\it EPS}). The EPS offers protection to the bacteria in the biofilm, which develop behavioral patterns that differ from those of planktonic bacteria \cite{Costerton}, \cite{Flemming}. As a consequence, harmful biofilm bacteria are more difficult to remove, mechanically or by antimicrobial agents, than suspended bacteria \cite{Stewart}, \cite{Watnick}. Biofilms occur naturally on surfaces and interfaces in aquatic systems wherever enough nutrients are available to sustain microbial growth \cite{Flemming}. Harmful occurrences of biofilms include bacterial infections, dental plaque and associated diseases, biocorrosion of drinking water pipes or industrial facilities, and food spoilage. On the other hand, the sorption properties and enhanced mechanical stability of biofilms are beneficially used by environmental engineers, {\it e.g.} in wastewater treatment, soil remediation, and groundwater protection \cite{Wanner:2006}. While the term biofilm indicates a homogeneous film-like layer, biofilms on the meso-scale ($10 \mu m \sim 1mm $, the actual biofilm scale) in reality can develop in structurally highly complicated architectures. The morphology of a biofilm is triggered by environmental conditions, such as the availability of required substrates like nutrients and, in the case of aerobic biofilms, oxygen \cite{vanL:1995}, \cite{vanL:1997}, \cite{Wanner:2006}, \cite{Xu:1998}. Well-known are the so-called mushroom or pillar shaped biofilm morphologies. These are clusters of biofilm colonies, trenched by water-filled pores and channels, cf. the confocal laser scanning micrograph in Figure \ref{CLSM} and the schematic in Figure \ref{schematic}. Biofilm morphologies of this type are observed in nutrient limited regimes. Initially, at low population density, the bacterial population develops unhindered. As the population grows, demand of food increases and inner-species competition for food sets in when nutrients become depleted. Large colonies have an advantage allowing them to grow faster towards the nutrient source and thus out-compete smaller colonies that stay behind in their development. As these biofilms grow under nutrient limitation they typically develop slowly. On the other hand, biofilms developing in a nutrient rich, un-limited environment tend to grow quickly and form homogeneous, compact, flat yet thick structures. In order to capture these features, spatially resolving techniques are required for biofilm research. In addition to modern non-invasive microscopy methods in experimental studies, mathematical models and computer simulations for theoretical studies become more and more accepted \cite{Wanner:2006}, \cite{Wilderer:2002}. \begin{figure}[ht] \begin{center} %\epsfxsize=10truecm\epsffile{ \includegraphics[width=0.8\textwidth]{figures/frame3gs} \end{center} \caption{Microscopy snapshot of a very young spatially heterogeneous biofilm of Listeria monocytogenes and Pseudomonas putida developing in a cluster-and-channel morphology. The CLSM figure was made available by Heidi Schraft, Lakehead University, Thunder Bay, On, Canada. The original colors were modified for better greyscale printing.} \label{CLSM} \end{figure} \subsection{Mathematical models and numerical simulation} Traditional biofilm models were formulated as one-dimensional free-boundary value problems, based on the assumption that newly produced biomass is converted into new biofilm volume. This leads to an increasing thickness of the biofilm and the speed of propagation of the biofilm/liquid interface normal to the substratum can be calculated from the production terms by integration over the biofilm thickness \cite{Wanner:2006}. As this concept is inherently one-dimensional, it cannot predict spatially heterogeneous biofilm morphologies. To fill this gap, a variety of models for spatially heterogeneous biofilms has been proposed in recent years, ranging from stochastic individual based models to cellular automata models to deterministic continuum models \cite{Pici:2003}. The underlying mathematical principles of these models are quite different and, therefore, the mathematical and computational challenges vary from model to model. Nevertheless, they all show the same qualitative behavior in predicting the development of biofilm morphologies in response to substrate limitation, cf. \cite{Pici:2003} and the references therein for a selection and overview. The model for biofilm formation that we consider was introduced originally in \cite{Eberl:2001} and studied in a small series of papers analytically, numerically and in applications \cite{Duvnjak:2004,Eberl:2003,Eberl:2004, Eberl:2006,Efendiev:2005,Khassehkhan:2006}. It is a quasi-linear diffusion equation for the dependent variable biomass density $u$, with two interacting causes of degeneracy. One is degeneracy as in the porous medium equation, i.e. the density dependent diffusion coefficient $D(u)$ vanishes as the dependent variable vanishes, $D(0)=0$; the other one is a singularity in the density-dependent diffusion coefficient as the dependent variable approaches its a priori known maximum density, $u\to 1$. In the current study we propose a numerical discretisation scheme that is able to handle these two simultaneously occurring effects and to render important qualitative features of the continuous model. The actual biofilm is the region in which the biomass density is positive. The transition from the surrounding aqueous phase to the biofilm is sharp, that is the biomass gradients are steep at the biofilm/liquid interface. Since the biofilm grows, this interface is not stationary but its location changes over time. The speed and direction of the interface propagation is not a priori known but an implicit consequence of the solution of the equation; furthermore, as two or more biofilm colonies merge into a bigger one the interfaces collide and become dissolved. Interface problems of this type are known to pose difficulties for many numerical methods, often leading to either un-physical and un-mathematical oscillations in the solution or smearing of the sharp interface. The underlying reason is that discretisation schemes often assume more regularity than the initial-boundary value problem offers. This effect can be observed already in very simple examples that possess classical smooth solutions, e.g. the heat equation with a discontinuous interface in the initial data \cite{Hundsdorfer}. A variety of numerical methods were used in previous studies of the biofilm model. In \cite{Eberl:2001,Eberl:2003,Eberl:2004}, based on a time-scale argument, the faster (in terms of characteristic time scales) substrate equations were solved implicitly while the slower biomass equation is solved with explicit schemes, e.g. the adaptive 4th/5th order Runge-Kutta-Fehlberg method, cf\cite{Epperson}. The maximum feasible time-step in this approach depends on the solution of the equation and becomes very small as the biofilm grows and maximum biomass density $u\to 1$ is approached somewhere. A very different idea was followed in \cite{Duvnjak:2004}, where the degenerated equation was transformed into new dependent variables such that the spatial operator becomes the Laplacian and all non-linear effects appear in the time-derivative. A fully implicit method for this equation was devised that has only a very mild time-step constraint. This method handles the diffusion singularity effect very well, but it was found that in some cases it emphasises interface oscillations, i.e. is not able to correctly deal with the porous medium degeneracy. In \cite{Khassehkhan:2006} the model equation was re-formulated in a weak form in the moving frame and subsequently discretised by a moving grid finite element method that explicitly tracks the biofilm/liquid interface, following an idea that was originally developed in \cite{Baines} for degenerate parabolic equations. This method describes the porous medium degeneracy very well. It is fully adaptive in time and space (with explicit time integration) and builds its biofilm grid adaptively. However, re-meshing and grid-refinement is required as new biomass is produced and the biofilm structure grows, in order to avoid finite elements becoming too large. This requires interpolation and becomes inefficient in higher dimensions. In this paper we will present a numerical method for the biofilm model that is based on a non-local (in time) representation of the nonlinear density dependent diffusive biomass flux. This leads to a semi-linear method that inherits the mild time-step constraint of the implicit method in \cite{Duvnjak:2004} and the advantages of the explicit method. In every time step the solution of only one linear algebraic system is required. In particular the method deals well with the moving interface, preserves the positivity of the solutions and is free of oscillations. \section{The continuous model}\label{sec_modeldesc} \begin{figure}[ht] \begin{center} \includegraphics[width=0.85\textwidth]{figures/schematic} \end{center} \caption{The domain $\Omega \subset \mathbb{R}^2$ with liquid region $\Omega_1(t)=\{x \in \Omega: u(t,x)=0\}$ and biofilm region $\Omega_2(t)=\{ x \in \Omega: u(t,x)>0\}$. As $u$ changes with time, the interface between $\Omega_1$ and $\Omega_2$ starts moving.} \label{schematic} \end{figure} We consider the biofilm as a continuum and study a density-dependent diffusion-reaction equation that describes the development of a spatially structured bacterial biofilm. The model is formulated as an evolution equation over the domain $\Omega \subset R^d$, $d=1,2,3$, in terms of the biomass density $u$, normalised with respect to the maximum biomass density. Under this {\it ansatz}, bacteria and EPS are considered together in one fraction. This is the classical approach in biofilm modelling \cite{Wanner:2006}. The region $\Omega_1(t):=\{x \in \Omega: u(t,x)=0\}$ describes the surrounding aquatic environment (bulk liquid, channels and pores of a biofilm) without biomass, while $\Omega_2(t)=\{ x \in \Omega: u(t,x)>0\}$ is the actual biofilm (cf. Figure \ref{schematic}). The model under consideration was introduced in \cite{Eberl:2001}. In the notation used here it reads \begin{equation}\label{model} u_t = \nabla_x (D(u) \cdot \nabla_x u) + ku \end{equation} where \begin{equation}\label{diffcoeff} D(u) = \delta {u^b \over (1-u)^a},\quad a,b \geq 1\gg \delta >0 \end{equation} together with Dirichlet, Neumann, or mixed boundary conditions. In (\ref{model}) the quantity $k$ describes the production rate. For our purpose it is sufficient to assume that it is bounded. If $k$ is a positive constant (\ref{model}) corresponds to a biofilm system in which nutrients are not limited. In this case we expect a homogeneous biofilm morphology for large $t$. In the case of nutrient limitations, $k(t,x)$ is a non-constant function. We shall present simulations for both scenarios below. Some analytical results about this model can be found in \cite{Duvnjak:2004}, \cite{Efendiev:2005}; biofilm applications of this type of equation and its generalisations were discussed in \cite{Eberl:2001}, \cite{Eberl:2003}, \cite{Eberl:2004}. The long term behavior of the initial-boundary value problem associated with (\ref{model}) was studied in \cite{Efendiev:2005}. It was shown that for initial data $u_0$ with \[ u_0\in L^\infty (\Omega), \quad F(u_0) \in H_0 ^1 (\Omega),\quad u_0\ge 0,\quad\Vert u_0 \Vert _{L^\infty} (\Omega)<1 \] where \[ F(u)=\int_0 ^u \frac{v^\alpha}{(1-v)^\beta} dv,\quad 0\le u<1 \] there exists a unique solution $u$ of (\ref{model}) in the class of functions \begin{gather*} u\in L^\infty (R_+ \times \Omega)\cap \, C([0,\infty),L^2 (\Omega))) \\ F(u)\in L^\infty(R_+,H^1(\Omega)\cap C([0,\infty),L^2(\Omega)) \\ 0\le u(t,x)\le 1 \quad \Vert u \Vert_{L^\infty (R_+ \times \Omega)}<1 \end{gather*} If the boundary conditions are homogeneous Neumann conditions everywhere, then the solution $u$ reaches its maximum density $1$ almost everywhere in finite time. Physically this corresponds to the situation where the biofilm growths unlimited in a closed vessel from where no biomass can escape. If Dirichlet conditions $u=0$ are specified somewhere on the boundary of $\Omega$, then the solution will remain below the maximum density for all $t>0$ almost everywhere. Further results include the existence of a global attractor \cite{Efendiev:2005} and the construction of a Lyapunov functional \cite{Duvnjak:2004}. An essential property of the solutions of (\ref{model}) is that initial data with compact support imply solutions with compact support and that the solution is indeed bounded by the maximum biomass density. The first effect is due to the porous medium degeneracy $u^b$ in (\ref{diffcoeff}), while the latter is due to the singularity $(1-u)^{-a}$ in (\ref{diffcoeff}). It is the production term $ku$ together with this second effect that drives the spatial spreading of $\Omega_2(t)$. This is counteracted by the degeneracy as $u=0$ at the interface. As a consequence, $u$ squeezes in $\Omega_2(t)$ and, in the absence of loss terms, approaches its maximum value 1. Hence, the interaction of both non-linear diffusion effects with the growth term in (\ref{model}) is needed to describe spatial biomass spreading. It leads to very steep gradients of $u$ at the moving interface between $\Omega_1(t)$ and $\Omega_2(t)$. If $\Omega_2(t)$ is initially not connected but consists of several ``sub-domains'' (each of which connected), these will eventually join if $k$ is positive and large enough, due to the expansion of each individual sub-domain. In biofilm applications this describes the merging of two initially isolated bacterial colonies. Thus, if two $\Omega_1$-$\Omega_2$ interfaces collide, they are dissolved and the sub-domains become connected. In \cite{Duvnjak:2004} a semi-discretisation in time was discussed that aimed at overcoming the obvious numerical difficulties arising when $u \to 1$. After transformation of the dependent variable $u \mapsto v:=F(u) $, the nonlinearities appear in the time-derivative, while the spatial spreading effect is described by the linear diffusion operator. Backward Euler discretisation leads to an elliptic boundary value problem in every time-step. Convergence results showed that the time-step restriction of this model is not critical, {\it i.e.} the maximum time-step size permitted for convergence is larger than the characteristic time-scale of the biofilm formation process, i.e. for given spatial discretisation it behaves like $\mathcal{O}(1/k)$. However, it was found that the numerical discretisation of the spatial terms with standard Finite Element or Finite Difference methods introduces oscillations with negative values at the biofilm/liquid interface. Local grid refinement can postpone and dampen this effect but not avoid it. In the current work we aim at a generally applicable numerical scheme for (\ref{model}) and model extensions that is free of this non-physical effect. \section{Numerical Method} \subsection{Finite Difference Scheme For The Continuous Problem} We introduce a finite difference scheme for (\ref{model}) on a regular grid based on a first order difference approximation for the time derivative and the usual second order spatial difference for self-adjoint diffusion operators, {\it cf} \cite{Morton:1994}. The key idea of the method is to represent the nonlinear diffusive flux $D(u)\cdot \nabla_x u$ in (\ref{model}) and the reaction term $k \cdot u$ non-local in time: The diffusion coefficient is evaluated in time level $t$ while the gradient is evaluated at the new time level $t+\Delta t$, {\it i.e.} one has \[ D(u)\cdot \nabla_x u \approx D\bigl(u(t,\cdot)\bigr)\cdot \nabla_x u(t+\Delta t, \cdot) \] and similarly for the reaction term \[ k\cdot u \approx k(t,\cdot) \cdot u(t +\Delta t,\cdot) \] Finite difference methods for (ordinary or partial) differential equations utilising such a non-local representation of non-linear terms are called {\it non-standard finite difference schemes} \cite{Lubuma:2001}. Despite the low order convergence of the chosen derivative approximations, such a discretisation is known to be optimal among linear finite difference schemes for certain problems under the aspect of positivity and monotonicity (cf. \cite{Hundsdorfer} for details). In the context of non-standard finite difference schemes this strategy of low order discretisation of derivatives is routinely used for the construction of discretisation schemes and referred to as {\it Rule 1} in \cite{Mickens:2000}. In the sequel, we denote by $u_j^n$ the numerical approximation of the exact solution $u(t^n, x_j)$ in the grid points $x_j \in R^d$. Here $j \in J$ can be a multi-index or a grid ordering in the multi-dimensional case $d>1$; then $J$ is an appropriate index set. The vector $U^n=(u_j^n)_{j\in J} \in R^{| J|}$ is the (ordered) vector, the coefficients of which are the numerical approximations of the solution at time level $n$. Employing a regular grid with spatial step size $\Delta x$ one obtains for the discretised diffusion operator \begin{equation}\label{diffdisc} \nabla_x \left(D(u) \nabla_x u\right) \dot= {1 \over 2 \Delta x^2}\sum_{k\in N_j} (D(u_k^n) +D(u_j^n))(u_k^{n+1} - u_j^{n+1}) \end{equation} where $N_j \subset J$ is the index set pointing at the direct neighbors of $x_j$ on the grid. For $d=1$ it has two elements, for $d=2$ there are four and for $d=3$ there are six. The non-standard finite difference scheme that is derived under the above assumptions reads \begin{equation}\label{imnsfd} U^{n+1} - U^n = \Delta t \mathcal{D}(U^n) \cdot U^{n+1} + \Delta t \mathcal{K} U^{n+1} \end{equation} where $\mathcal{K}=diag(k_j^n)$ is the diagonal matrix the coefficients of which are the net growth rates $k(t^n,x_j)$. The matrix $\mathcal{D}(U^n)$ is the banded matrix obtained by second order standard finite difference discretisation of the diffusion operator evaluated at the previous time-step $t^n$, according to (\ref{diffdisc}). System (\ref{imnsfd}) can be re-written in matrix-vector form as \begin{equation}\label{nsfd} U^{n+1} = \left[\mathcal{I} - \Delta t \mathcal{D}(U^n) - \Delta t \mathcal{K}\right]^{-1} U^n \end{equation} if the inverse matrix exists. This is the case if the time-step $\Delta t$ is chosen such that $1/\Delta t$ is not an eigenvalue of the matrix $\mathcal{D}(U^n) + \mathcal{K}$. Then, the numerical solution $U^{n+1}$ is unique. Hence, a sufficient condition for the existence of a unique solutions is the restriction of the time-step size $t^n \to t^n+\Delta t =: t^{n+1}$ by \begin{equation}\label{dtrestrict} \Delta t \cdot \max_i k_i^n <1 \end{equation} This ensures that $I-\Delta t \mathcal{K}$ is diagonal with strictly positive entries. Then in (\ref{nsfd}) the matrix $\mathcal{ I} - \Delta t \mathcal{D}(U^n) - \Delta t \mathcal{K}$ is diagonally dominant. Inequality (\ref{dtrestrict}) poses a time-step restriction for the numerical method. From the application point of view this is a weak restriction, as $1/k$ is the characteristic time scale of biomass production. Accordingly, time-steps larger than $1/k$ would be at the expense of inaccurate description of the growth process. The same time stepping condition (\ref{dtrestrict}) was obtained for the transformation method in \cite{Duvnjak:2004}. If $k$ is not a constant, (\ref{dtrestrict}) leads to a variable time-stepping strategy. For a given $U^n$ with $0\leq U^n <1$ (where all vector inequalities are understood coefficient-wise), $U^{n+1}$ according to (\ref{nsfd}) is a continuous function of $\Delta t$. Hence, for $\Delta t$ small enough $U^n <1$ implies $U^{n+1}<1$. A more quantitative statement will be presented in the next section. The non-standard finite difference scheme (\ref{nsfd}) is a system of explicit non-linear difference equations, {\it i.e.} it can be understood as a discrete dynamic system. In the numerical simulation the solution of one coupled linear system is required in every time-step. The matrix is sparse, banded and diagonally dominant. In two- and three-dimensional simulations, the system can be solved iteratively, e.g. by the stabilised bi-conjugate gradient method, which is described in \cite{Saad}; due to diagonal dominance, pre-conditioning is not required for convergence. In the one-dimensional case the system is tridiagonal and the Thomas algorithm (e.g. \cite{Epperson}) can be used. These are the methods that we employ in the simulations below. In the next section we will show that the numerical solution according to (\ref{nsfd}) inherits the following important qualitative properties of the continuous model: \begin{itemize} \item positivity \item boundedness $U^n < 1$ of the solution \item finite speed of propagation of the discrete liquid/biofilm interface \item the solution of (\ref{nsfd}) satisfies a discrete interface condition that corresponds to the interface condition for the continuous model (\ref{model}). \item the numerical interface is sharp, i.e. only weak interface smearing takes place \item the solution is monotonous at the interface and merging of two colonies is well-posed \end{itemize} \subsection{Properties Of The Scheme} \begin{lemma} [Positivity] \label{posLemma} Let (\ref{dtrestrict}) be always satisfied. If $0\leq U^0 <1$ then $U^{n}\geq 0$ \end{lemma} \begin{proof} We use an argument from the theory of $M$-matrices (cf. \cite{Hackbusch}) and mathematical induction: If $\Delta t$ is chosen according to (\ref{dtrestrict}) then the matrix in (\ref{nsfd}) is diagonally dominant with positive diagonal elements and non-positive off-diagonals. Hence, its inverse is positive and, thus, positivity of $U^n$ implies positivity of $U^{n+1}$ according to iteration (\ref{nsfd}). \end{proof} \begin{remark} \label{rmk3.2} \rm An alternate proof can be derived based on the following observation: In the time step $t^n \to t^{n+1}$ method (\ref{imnsfd}) can be understood as an implicit Euler method for the linear system of ordinary differential equations \begin{equation}\label{semidiscrete} \dot v = \left(\mathcal{D}(U^n)+\mathcal{K}\right) v, \quad v(t^n)=U^n, \quad t\in [t^n,t^{n+1}] \end{equation} where $v$ is a time-continuous vector valued function. Thus, method (\ref{imnsfd}) is a piecewise linear problem and in every time-step linear theory applies, which can be found {\it e.g.} in \cite{Hundsdorfer}. Thus in every time-step positivity is ensured if the matrix in (\ref{semidiscrete}) has only non-negative off-diagonal entries and if the diagonal is bounded from below, {\it i.e.} if there is an $\alpha_n>0$ such that $\left(\mathcal{D}(U^n)+\mathcal{K}\right)_{jj} > -\alpha_n$. That this is the case follows from a recursive argument: It is easily verified that these conditions are satisfied for $n=0$. If they hold for one $n$ then it follows from a simple calculation that it also holds for $n+1$, due to the assumption that $\Delta t$ is small enough to warrant $U^n<1$ and the definition of $\mathcal{ D}(U)$. Note that both methods of proof for Lemma \ref{posLemma} use the same property of the system matrix in (\ref{nsfd}), and thus are essentially equivalent. \end{remark} The M-matrix property is also used to derive the following sufficient (but not necessary) time-step condition for boundedness of the solution. \begin{lemma} [Boundedness] \label{bounded} The choice of a time-step $\Delta t$ such that (\ref{dtrestrict}) holds and \[ \Delta t < \min_i {1 -u_i^n \over \sum_j[\mathcal{D}(U^n) + \mathcal{ K}]_{ij}} \] guarantees that the numerical solution obeys the upper bound, i.e. $U^{n+1} \leq 1$ if $U^n <1$. \end{lemma} \begin{proof} We introduce $W^n:=e-U^n$ where $e:=(1,\dots ,1)^T$. Hence, $U^n\leq 1$ is equivalent to $W^n\geq 0$. Then (\ref{nsfd}) gives \begin{equation}\label{wnsfd} \left[\mathcal{I} - \Delta t \mathcal{D}(U^n) - \Delta t \mathcal{ K}\right](e-W^{n+1})=(e-W^n) \end{equation} and, therefore, \begin{equation} \left[\mathcal{I} - \Delta t \mathcal{D}(U^n) - \Delta t \mathcal{ K}\right] W^{n+1} =W^n- \Delta t \left[\mathcal{D}(U^n) + \mathcal{ K}\right]e \end{equation} If $W^n>0$, positivity of $W^{n+1}$ is ensured if $\Delta t$ is chosen such that the right hand side of (\ref{wnsfd}) is non-negative. This is certainly the case if \begin{equation}\label{dtrestrict2} \Delta t < \min_i {1 -u_i^n \over \sum_j[\mathcal{D}(U^n) + \mathcal{ K}]_{ij}} \end{equation} \end{proof} It is easy to verify that criterion (\ref{dtrestrict2}) can be more strict than (\ref{dtrestrict}). It can be relaxed by the observation that it is sufficient but not necessary. Indeed, in the numerical simulations conducted below, this time-step criterion was implemented as an emergency strategy to guarantee the upper bound if other time-step control mechanisms fail, but never activated. Note that if in one coefficient $U^n$ reaches 1, (\ref{dtrestrict}) implies the termination of the simulation. This is in accordance with the continuous problem (\ref{model}), since it was shown in \cite{Efendiev:2005} that the solution of (\ref{model}) can reach 1 in finite time in dependence of boundary conditions and reaction rates. To formulate the next result we introduce the following concept. \noindent{\bf Definition.} We denote by $\Omega_1^n$ a discrete grid analogy of $\Omega_1(t^n)$ as \[ \Omega_1^n = \{ x_j : u_j^n=0 \} \] An interior point of $\Omega_1^n$ is a point in $\Omega_1^n$ whose direct neighbors are in $\Omega_1^n$ as well, i.e. a point $x_j \in \Omega_1^n$ such that $x_i \in \Omega_1^n$ for all $i \in N_j$. Similar definitions can be made for the biofilm region $\Omega_2$ Note that this definition is based on the numerical solution $u_j^n$. Thus it does not imply that $\Omega_1^n$ is contained in $\Omega_1(t^n)$ or {\it vice versa}. For now, the discrete interface at time level $t^n$ is defined by grid points of $\Omega_1^n$ with a direct neighbor in $\Omega_2^n$. The following result states that the discrete interface travels at most with finite speed. In particular it states that in one time-step it moves across at most one grid cell. \begin{lemma} [Finite speed of propagation] \label{finspeedLemma} Let \eqref{dtrestrict} hold. If $x_j$ is an interior point of $\Omega_1^n$ then $x_j \in \Omega_1^{n+1}$. \end{lemma} \begin{proof} If $x_j\in \Omega_1^n$ is interior then the corresponding $j$th row of the matrix $\mathcal{D}(U^n)$ contains only zeros. Thus one has from (\ref{nsfd}) \[ (1-\Delta t k_j^{n}) u_j^{n+1} = u_j^n= 0 \] with (\ref{dtrestrict}) it follows that $u_j^{n+1}=0$ and thus $x_j \in \Omega_1^{n+1}$. \end{proof} In the one-dimensional case $d=1$ an explicit expression can be derived for the speed of interface propagation in the continuous model (\ref{model}). An interface between $\Omega_1(t)$ and $\Omega_2(t)$ can be parameterised locally by a curve $x^\ast(t)$. The speed of the interface is $\dot x^\ast(t)$. Due to continuity and $u(x^{\ast}(t),t)=0$ we obtain from (\ref{model}) \begin{equation}\label{intcond} \dot x^\ast(t) = -[u_t/u_x]_{x^\ast(t)} \end{equation} where the right hand side is to be understood in the sense of a limit of the quotient of the one-sided derivatives as one approaches the interface. The theory developed in \cite{Efendiev:2005} implies that this is finite. We show that the numerical solution satisfies a corresponding discrete interface condition if the derivatives in (\ref{intcond}) are replaced by the usual difference quotients. Without loss of generality, we restrict ourselves to the case of an expanding biofilm (i.e. $\Omega_2$ increases) with an interface with positive speed. We define the location $x_{j^{\ast,n}}$ of a discrete interface at time $t^n$ as follows: $x_{j^{\ast,n}-1} \in \Omega_2^n$, $x_{j^{\ast,n}} \in \Omega_2^n$, $x_{j^{\ast,n}+1} \in \Omega_1^n$ and $x_{j^{\ast,n}+2}$ an interior point of $\Omega_1^n$. That is, $j^{\ast,n}$ is the index of the last biofilm point before the liquid/biofilm interface at time $t^n$. \begin{lemma} [Moving interface condition] The following discrete interface condition holds \begin{equation}\label{dintcond} { x_{j^{\ast,n+1}} - x_{j^{\ast,n}} \over \Delta x} = - { u_{j^{\ast,n+1}}^{n+1} - u_{j^{\ast,n+1}}^{n} \over u_{j^{\ast,n+1}+1}^{n+1} - u_{j^{\ast,n+1}}^{n+1}} \end{equation} \end{lemma} \begin{proof} For brevity of the notation let $i:=j^{\ast,n}$. We obtain from the $(i+1)$th coefficient of (\ref{nsfd}) and the hypothesis on the position of the interface \begin{equation}\label{interface} -{\Delta t \over 2} D(u_i^n) u_i^{n+1} +\Big(1 +{\Delta t \over 2} D(u_i^n) - \Delta t k_{i+1}^{n}\Big) u_{i+1}^{n+1}=0 \end{equation} and thus $u_{i+1}^{n+1} >0$. With the previous Lemma we have that $j^{\ast,n+1}=i+1=j^{\ast,n}+1$. Thus $u_{j^{\ast,n+1}}^{n}=0$. Furthermore we have by hypothesis $u_{j^{\ast,n+1}+1}^{n+1}=0$. The assertion follows by substituting these observations into (\ref{dintcond}). \end{proof} Note that (\ref{dintcond}) is a discrete version of (\ref{intcond}) after multiplying with $\Delta x/\Delta t$. The proof showed that it could have been formulated in a stronger version that essentially states that the speed of propagation of the discrete interface is $\Delta x/ \Delta t$. This is a consequence of the smearing around the interface that is introduced by the spatial discretisation of the diffusion operator. However, from (\ref{interface}) it follows \[ {\Delta t \over 2} D(u_i^n)\left[u_{i+1}^{n+1} - u_i^{n+1}\right] + \left[1 - \Delta t k_{i+1}^{n}\right] u_{i+1}^{n+1} =0 \] and thus with (\ref{dtrestrict}) and the positivity of (\ref{nsfd}) we obtain $0< u_{i+1}^{n+1} < u_i^{n+1}$. Moreover, after re-arranging terms we have \[ u_{i+1}^{n+1} = \mathcal{O}\left((u_i^n)^b u_i^{n+1}\right) \] with the definition (\ref{diffcoeff}), since $u_i^n$ small at the interface implies $(1-u_i^n) \gg 0$. In biofilm applications one has $b$ in the range 2 through 6. Thus we may expect the interface smearing effect to be quantitatively small. This will be demonstrated in the numerical simulations in the next section. The last theoretical result in this section shows that the numerical solution at time $t^{n+1}$ in a grid point in $\Omega_1^n$ is bounded from above by the values of the solution in neighboring grid points. This is of particular interest in the case where the grid point changes its membership from $\Omega_1^n$ to $\Omega_2^{n+1}$ as $t^n \to t^{n+1}$, i.e. is passed by a biofilm/water interface. \begin{lemma} [Merging of two colonies] \label{mergeLemma} Let $x_i \in \Omega_1^n$ and assume that (\ref{dtrestrict}) holds. Then $u_i^{n+1}$ is bounded by the values of $u$ in the neighboring nodes, i.e. \[ 0\leq u_i^{n+1} \leq \max_{j\in N_i}\{u_j^{n+1} \} \] \end{lemma} \begin{proof} We have $x_i^n \in \Omega_1^n \Longrightarrow u_i^n =0 \Longrightarrow D(u_i^n)=0$. Thus (\ref{imnsfd}) and (\ref{diffdisc}) imply \[ u_i^{n+1} = {\Delta t \over \Delta x^2}\sum_{j\in N_i} D(u_j^n)(u_j^{n+1} - u_i^{n+1}) + \Delta t k_i^{n} u_i^{n+1} \] and, hence, \begin{align*} u_i^{n+1} &={{\Delta t \over \Delta x^2}\sum_{j\in N_i} D(u_j^n)u_j^{n+1} \over 1+ {\Delta t \over \Delta x^2} \sum_{j\in N_i} D(u_j^n) - \Delta t k_i^{n} } \\ &\leq {{\Delta t \over \Delta x^2}\sum_{j\in N_i} D(u_j^n) \cdot \max_{j\in N_i} u_j^{n+1} \over 1+ {\Delta t \over \Delta x^2} \sum_{j\in N_i} D(u_j^n) - \Delta t k_i^{n} } \\ &\leq \max_{j\in N_i} u_j^{n+1} \end{align*} The last inequality follows from (\ref{dtrestrict}) with $1-\Delta t k_j^n>0$. From Lemma \ref{posLemma}, we have $0\leq u_i^{n+1}$. \end{proof} Thus, if $x_i \in \Omega_1^n$ separates two approaching interfaces then Lemma \ref{mergeLemma} implies that the colonies merge smoothly without forming bulges or inducing oscillations. In the case of a single moving interface, Lemma \ref{mergeLemma} implies monotonicity, i.e. the absence of non physical oscillations. Note that Lemma \ref{mergeLemma} includes the previous Lemma \ref{finspeedLemma} for the trivial case where $x_i$ is interior to $\Omega_1^n$. \section{Numerical Simulations}\label{simulation_sec} We present three sets of simulations of (\ref{model}) with the nonstandard finite difference scheme (\ref{nsfd}). The first one is one-dimensional in space. While this does not fully capture the features of biofilm formation and growth it allows to investigate the behavior of the numerical scheme. In the second part of this section we carry out some two-dimensional numerical illustrations of biofilm growth for constant net growth rate $k$. In both cases we had for the parameters of the dimensionless equation $a=b=4$, $k=0.1$, $\delta=10^{-8}$. In the third simulation study we apply the discretisation scheme to a more complicated system, where biofilm formation is controlled by two dissolved substrates. All simulations were carried out using double accuracy arithmetics. Codes were implemented in Fortran 95 on Linux based personal computers. In all simulations the time-step was variable, as \begin{equation}\label{dtchoice} \Delta t = \min\Big\{1/k, \min_j {M \Delta x^2 \over D(u_j^n)} , 0.1 \Big\} \end{equation} The first condition is according to (\ref{dtrestrict}), the last one to keep the time-step bounded for the sake of accuracy. This was also the reason to introduce the second condition, where $M$ is a user-defined constant. It renders the typical $\Delta t / \Delta x^2$ dependency for parabolic problems and allows comparability of results for different choices of $\Delta x$. In our simulations, we use $M=40$; this choice of time-step is well above the limitations of the explicit method. That said, this time-step constraint was introduced as a trade-off between between fast and accurate computation. The condition (\ref{dtrestrict2}) was implemented as an emergency brake, such that it becomes active only if (\ref{dtchoice}) would lead to $\quad \max U^{n+1}>1$. It never became activated. \subsection{One-dimensional simulations} \begin{figure} \begin{center} \includegraphics[width=0.488\textwidth]{figures/nsfd1} \includegraphics[width=0.488\textwidth]{figures/nsfd2} \end{center} \caption{Simulation of biofilm formation in time: the inoculum at $t=0$ consist of (right) three sinusoidal colonies of different height and (left) two colonies with constant densities of different height.} \label{nsfdfig} \end{figure} The results of two numerical simulations of (\ref{model}) with the finite difference method presented here are shown in Figure \ref{nsfdfig}. In case the (I), the initial data were specified by the function \[ u(x,0)=\max( 0, -0.8\sin(7\pi x)(1-x^4)). \] In the case (II), the initial data were specified by \[ u(x,0)=\begin{cases} 0.8, & 0.3\leq x \leq 0.4 \\ 0.9, & 0.5 \leq x\leq 0.6, \\ 0, & \mbox{elsewhere}. \end{cases} \] In the case (III), finally one sinusoidal colony was placed around the centre of the interval $(0,1)$. This scenario is included to study the dependence of the interface on the spatial resolution. That is, in the first and third case the initial data are continuous but not differentiable at the interface, while in the second case they are only piecewise continuous. Both scenarios are known to be problematic if treated with standard techniques. It is easily verified that in the first case the biomass is organised in three colonies originally. Figure \ref{nsfdfig} shows the simulations of cases I and II. They were carried out with a spatial step-size $\Delta x = 1/200$. In the third scenario the dependency of the simulation results on the spatial resolution was tested. To this end, the location of the actual interface between $\Omega_1(t)$ and $\Omega_2(t)$ for different step-sizes $\Delta x=1/N,\quad N=50, 100, 200, 400, 800$ is plotted in the $x$-$t$ plane. The second panel in Fig. \ref{interfacefig} shows the solution surface for various $\Delta x$. In Figures \ref{nsfdfig} and \ref{interfacefig} (right panel), the solution surfaces $u(t,x)$ are plotted where $| u | > 10^{-16}$ for the sake of readability, assuming that $10^{-16}$ is a reasonably good approximation of $0$ for our purpose. \begin{figure}[t] \includegraphics[width=0.48\textwidth]{figures/interfc} \includegraphics[width=0.48\textwidth]{figures/nsfdsc3} \caption{left: Location of the left interface between liquid and biofilm in the third scenario for different choices of $\Delta x=1/N, \quad N=50, 100, 200, 400, 800$; right: solution of scenario 3 for various choices of $N$.} \label{interfacefig} \end{figure} As predicted by the theory outlined above, the new numerical method gives non-negative and oscillation-free solutions. The simulations also confirm the above statement that the diffusive smearing around the interface is quantitatively small, {\it i.e.} negligible. In both first scenarios the colonies of the inoculum compress initially, that is the biomass density $u$ increases while only very little spatial expansion is observed. As $u$ approaches 1, $\Omega_2$ expands exponentially due to the first order growth term in (\ref{model}) and the squeezing property. In both cases the colonies eventually merge. In the first scenario, the first initial colony is bigger than the second one, which is bigger than the third one. Accordingly the first and second colony merge first. The actual interface was defined somewhat arbitrary as the location at which $u$ becomes bigger than $10^{-5}$ (as an approximation of $0$). Both figures show a very good agreement. It can be seen that the computed interface converges as $\Delta x$ becomes smaller. Even for the coarsest discretisation the deviations are accurate within $\Delta x$. This again confirms the above comment that the smearing effect at the interface is small. \subsection{A two-dimensional illustration: Formation of a homogeneous biofilm} \begin{figure}[ht] \begin{center} \includegraphics[width=0.48\textwidth]{figures/2d1} \includegraphics[width=0.48\textwidth]{figures/2d12} \\ \includegraphics[width=0.48\textwidth]{figures/2d20} \includegraphics[width=0.48\textwidth]{figures/2d28} \\ \includegraphics[width=0.48\textwidth]{figures/2d36} \includegraphics[width=0.48\textwidth]{figures/2d42} \\ \includegraphics[width=0.48\textwidth]{figures/2d48} \includegraphics[width=0.48\textwidth]{figures/2d59} \end{center} \caption{Formation of a homogeneous, thick biofilm from originally heterogeneously distributed biomass. The substratum $y=0$ is at the top of each picture. Shown are snap shots at $t=0.01, 1.22, 2.1, 2.98, 3.85, 4.45, 5.06, 6.16$ (top left to bottom right). Biomass density $u$ is coded in a linear greyscale: black corresponds to $u=0$, white to $u=1$.} \label{Fig2D} \end{figure} The simplified biofilm model (\ref{model}) with positive $k={\rm const}$ describes biofilm formation under non-limiting conditions. This can be the case for nutrient rich environments, e.g. bacterial growth in foods or relatively thin biofilms (in the initial stages) in wastewater treatment plants. For many other biofilms this is not a realistic assumption because often either nutrients or oxygen or ph-value or other controlling substances are not constant inside the biofilm. Typically, the thicker a biofilm is, the more limited become nutrients in the deeper layers. Often, decay due to lysis will even become dominant, leading to negative local net growth rates. This spatially heterogeneous situation will be addressed in the next section. Since we want to focus on the spatial discretisation of models of this kind first, we nevertheless accept the case $k={\rm const}$ for an illustrative example. Under such non-limiting conditions it is expected (cf. \cite{Eberl:2001}, \cite{Wanner:2006}) that even starting from a heterogeneous initial distribution of biomass, a homogeneous biofilm will form eventually and that the biomass distribution within $\Omega_2$ will be more or less homogeneous. This will be verified in the following simulations. Initially biomass is randomly distributed in 15 different locations across the substratum, {\it i.e.} the surface on which the biofilm forms. That is, the biomass is placed initially in the first layer of grid cells. The biomass density in each of these locations is chosen randomly between 0 and 1. The computational domain is the interval $[0,1] \times [0,0.3]$. We specify Neumann boundary conditions at the substratum $y=0.3$ and on the boundaries $x=0$, $x=1$. The first assumptions ensures that no biomass leaves the system across the solid surface, while the second assumption allows us to consider the computational domain as part of a larger system (periodic conditions would do this as well). At $y=0$ we specify homogeneous Dirichlet conditions. The results of a simulation are shown in Figure \ref{Fig2D}. The rate at which the biofilm colonies grow depends initially on the biomass density at $t=0$. Colonies with larger $u(x,0)$ grow faster and merge earlier with their neighbors. Eventually all colonies merge and form a compact film as expected. As it was already observed in the one-dimensional simulations, the biomass density inside the biofilm approaches the maximum density rapidly everywhere, due the homogeneous growth rate. \subsection{A two-dimensional illustration: Formation of a cluster-and-channel biofilm}\label{mushroom_sec} We investigate now the more general case where the local growth of biomass is controlled by the availability of required substances. As an example we choose a heterotrophic biofilm that is limited by two dissolved substrates, nutrients and oxygen. This is a standard example of biofilm systems in wastewater engineering and was defined by the International Water Association's taskgroup on biofilm modeling as a first benchmark system (BM1) for biofilm models, cf. \cite{Morgenroth:2004}, \cite{Wanner:2006}. The local net growth rate $k(t,x)$ in this example is given by \begin{equation}\label{reacs} k(t,x)= k_1 {c_1(t,x) \over k_{s,1} + c_1(t,x)} {c_2(t,x) \over k_{s,2} + c_2(t,x)} - k_3 \end{equation} where the positive constant $k_1$ is the maximum specific growth rate and the positive constant $k_3$ is the lysis rate describing biomass deactivation. The positive parameters $k_{s,1}$ and $k_{s,2}$ are the Monod half saturation constants. By $c_1$ and $c_2$ we denote the concentrations of the dissolved substrates. They are governed by the semi-linear diffusion-reaction equations \begin{equation}\label{c1eq} c_{1,t} = D_1 \Delta_x c_1 - k_4 u {c_1 \over k_{s,1} + c_1} {c_2 \over k_{s,2} + c_2} \end{equation} and \begin{equation}\label{c2eq} c_{2,t} = D_2 \Delta_x c_2 - k_5 u {c_1 \over k_{s,1} + c_1} {c_2 \over k_{s,2} + c_2} \end{equation} where the consumption rates $k_4$ and $k_5$ are essentially the maximum specific growth rate divided by the respective yield factors, which quantify how much substrate is needed to produce one unit mass of microbial biomass. The positive parameters $D_1$ and $D_2$ are the diffusion coefficients of both substrates. For small molecules these are constants, attaining the same values in the biofilm and in the surrounding bulk phase. For our simulation we choose the same reaction and diffusion parameters as in the benchmark problem specification in \cite{Morgenroth:2004}. The model is completed by initial and boundary conditions for $c_1$ and $c_2$. We choose the computational domain to be rectangular of size $L_1 \times L_2$ and specify the mixed boundary conditions for $j=1,2$ \[ c_{j}=c_{\infty,j} \quad\mbox{for}\quad x_2=L_2 \] and \[ {\partial c_j \over \partial n} =0 \quad\mbox{elsewhere on}\quad \partial \Omega \] The initial conditions for $c_{1,2}$ are \[ c_j(x,0)= c_{\infty,j} \] For the biomass density $u$ we specify the same initial and boundary conditions as in the previous example. The purpose of this example is to show that the numerical scheme is able to simulate the formation of spatially heterogeneous of mushroom-shaped cluster-and-channel biofilm architectures. It is well understood in the biofilm literature, cf \cite{Wanner:2006}, that a good predictor for the surface heterogeneity of a biofilm is how biomass growth terms relate to substrate availability. If substrate does not become limited, a flat, spatially homogeneous biofilm is eventually obtained, as in the previous example. If substrate in the biofilm becomes limited, growth of biomass slows down. Microbial inner-species competition for nutrients leads to irregular biofilm morphologies. Eventually larger biofilm colonies closer to the food source become dominating and smaller colonies stay back in their development. With our choice of parameters from \cite{Morgenroth:2004}, the growth terms are fixed. Hence, we use substrate availability to control the biofilm structure. It is quantified in terms of the diffusion coefficients (also fixed from \cite{Morgenroth:2004}) and the environmental conditions, in particular the Dirichlet values $c_{\infty,j}$ and the diffusion length $L_2$. We choose the bulk concentrations in the same range as the Monod half saturation concentrations, $c_{\infty,1}=k_{s,1}$ and $c_{\infty,2}=2k_{s,2}$. The maximum principle for parabolic equations implies that these values are upper estimates for the concentrations $c_1$ and $c_2$ and that inside the domain and, hence, in the biofilm lower concentration values will be obtained. Thus, our choice of environmental concentrations ensures that substrates are not available in abundance. The system height $L_2$ is the distance over which substrates need to be transported by diffusion to reach the bacteria at the substratum. The smaller this value is, the more compact a biofilm develops \cite{Wanner:2006}. We choose $L_2=L_1=500\mu m$, i.e. carry out our simulations in a square domain. \begin{figure}[ht] \begin{center} \includegraphics[width=0.47\textwidth]{figures/BM1_001} \includegraphics[width=0.47\textwidth]{figures/BM1_005} \\ \includegraphics[width=0.47\textwidth]{figures/BM1_007} \includegraphics[width=0.47\textwidth]{figures/BM1_010} \\ \includegraphics[width=0.47\textwidth]{figures/BM1_020} \includegraphics[width=0.47\textwidth]{figures/BM1_040} \end{center} \caption{Formation of a cluster-and-channel biofilm morphology (top left to bottom right): Shown are the biofilm/liquid interface and the limiting oxygen concentration $c_2(t,x)$ in time $t= 0.02, 4.03, 6.05, 9.08, 19.12, 39.13 d$. The food source is located at the right boundary $x_2=L_2$.} \label{BM1_sim_fig} \end{figure} \begin{figure*} \begin{center} \includegraphics[width=0.48\textwidth]{figures/BM1_080} \includegraphics[width=0.48\textwidth]{figures/BM1_120} \\ \includegraphics[width=0.48\textwidth]{figures/BM1_160} \includegraphics[width=0.48\textwidth]{figures/BM1_200} \\ \includegraphics[width=0.48\textwidth]{figures/BM1_230} \includegraphics[width=0.48\textwidth]{figures/BM1_250} \end{center} \caption{Figure \ref{BM1_sim_fig} continued for $t= 79.13, 119.13, 159.13, 199.13, 229.20, 249.35 d$} \end{figure*} For the numerical solution of (\ref{c1eq}) and (\ref{c2eq}) we use a standard finite volume scheme on the same grid on which (\ref{model}) is discretised, 2nd order in space and 1st order in time. The simulation results of this example are shown in Figure \ref{BM1_sim_fig}; plotted are the biofilm/liquid interface and the concentration field of oxygen, $c_2$. The nutrient concentration field $c_1$ is qualitative similar, albeit on a higher quantitative level. Initially the substratum is inoculated by few small bacterial colonies. As time increases, first a smooth, almost homogeneous, smooth biofilm layer develops which than changes in a cluster-and-channel morphology as substrates become limited. It should be noted that due to the lysis term in (\ref{reacs}) and due to substrate limitations the active biomass density remains distinctively separated from the maximum biomass density, i.e. $u<1$ clearly. Thus, the simulation results agree with our a priori expectations on model behavior. \section{Conclusion} A finite difference scheme was developed for a parabolic evolution equation with two distinct nonlinear effects in the density-dependent diffusion coefficient. One is degeneracy as in the porous medium equation, while the other one is a singularity as the dependent variable approaches its {\it a priori} known maximum value. In the numerical method, the nonlinearity of the diffusion operator was handled in a non-local representation in time. It was shown that the numerical method renders the essential qualitative features of the solution of the continuous problem. In particular it is free of non physical oscillations that often are observed in problems of this kind, it shows the finite speed of propagation of interfaces between the regions $\Omega_1(t)$ where $u=0$ and $\Omega_2(t)$ where $u>0$, and it guarantees that the upper bound of the solution. The nonstandard finite difference scheme for the proto-type biofilm model (\ref{model}) is based on a non-local in time representation of the nonlinear diffusive flux. The method was constructed in a manner that allows a straightforward application to more complicated biofilm systems with several particulate substances as in section \ref{simulation_sec}\ref{mushroom_sec} and dissolved substrates, {\it i.e.} for models as studied in {\it e.g.} \cite{Eberl:2003}, \cite{Eberl:2004}. Although only one- and two-dimensional simulations have been carried out here for the illustration of the new method, an application to more realistic three-dimensional biofilm descriptions is possible. In fact, the method (\ref{nsfd}) itself is formulated in the general setup and the analytical results were derived for the three-dimensional case. This makes the new method attractive for further studies in mathematical modeling of biofilms. \subsection*{Acknowledgments} This study was supported by the Volkswagen Stiftung, Hannover, Germany with a grant in {\it Interdisciplinary Environmental Research}, by the NSERC Canada with a {\it Discovery Grant} and by the {\it Advanced Foods and Materials Network of Centres of Excellence (AFMnet)} with a project grant. The original draft version of this manuscript was written while the first author visited the GSF Research Centre in Neuherberg, Germany. The authors also wish to thank Messoud Efendiev for many stimulating discussions. The CLSM micrograph was made available by Heidi Schraft, Dept. Biology, Lakehead University, Thunder Bay, ON. \thebibliography{00} \bibitem{Lubuma:2001} Anguelov R, Lubuma J. M. S. Contributions to the mathematics of the nonstandard finite difference method and its applications, {\it Num. Meth. PDE} 17:518-543, 2001 \bibitem{Baines} Baines M. J., Hubbard M. E., Jimack P. K., A moving mesh finite element algorithm for the adaptive solution of time dependent PDEs with moving boundaries, {\it Appl. Num. Math.}, 54:450-469, 2005 \bibitem{Costerton} Costerton JW, Lewandowsky Z, Caldwell DE, Korber DR, Lappin-Scott HM, Microbial Biofilms, {\it Annu. Rev. Microbiol.} 49:711-745, 1995 \bibitem{Duvnjak:2004} Duvnjak A, Eberl H. J. Time-discretisation of a degenerate reaction-diffusion equation arising in biofilm modeling, {\it El. Trans Num. Analysis}, 23:15-38, 2006 \bibitem{Eberl:2001} Eberl H. J., Parker D. F., van Loosdrecht M. C. M. A new Deterministic Spatio-Temporal Continuum Model For Biofilm Development, {\it J. Theor. Medicine}, 3(3), 2001 \bibitem{Eberl:2003} Eberl H. J., Efendiev M. A. A Transient Density Dependent Diffusion-Reaction Model for the Limitation of Antibiotic Penetration in Biofilms. {\it El. J. Diff Eq}. CS10:123-142, 2003 \bibitem{Eberl:2004} Eberl H. J., A deterministic continuum model for the formation of EPS in heterogeneous biofilm architectures, {\it Proc. Biofilms 2004}, Las Vegas, 2004 \bibitem{Eberl:2006} Eberl H. J., Schraft H., A diffusion-reaction model of a mixed culture biofilm arising in food safety studies, {\it accepted} \bibitem{Efendiev:2005} 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. {\it RIMS Kyoto}, 1258:49-71, 2002 \bibitem{Epperson} Epperson J. F., {\it An Introduction to Numerical Methods and Analysis}, Wiley \& Sons, 2002 \bibitem{Flemming} Flemming H. C., Biofilme -- das Leben am Rande der Wasserphase, {\it Nachr. Chemie}, 48:442-447, 2000 \bibitem{Hackbusch} Hackbusch W., {\it Theorie und Numerik elliptischer Differentialgleichungen}, Teubner, Stuttgart, 1986 \bibitem{Hundsdorfer} Hundsdorfer W., Verweer J. G., {\it Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations}, Springer, 2003 \bibitem{Khassehkhan:2006} Khassekhan H., Eberl H. J., Interface tracking for a non-linear degenerated diffusion-reaction equation describing biofilm formation, {\it accepted} \bibitem{vanL:1995} Van Loosdrecht M. C. M., Eikelboom D., Gjaltema A., Mulder A., Tijhuis L., Heijnen J. J., Biofilm Structures, {\it Wat. Sci. Tech.} 32(8):35-43, 1995 \bibitem{vanL:1997} Van Loosdrecht M. C. M. , Picioreanu C., Heijnen J. J., A more unifying hypothesis for the structure of microbial biofilms, {\it FEMS Microb. Ecol.}, 24:181–183, 1997 \bibitem{Mickens:2000} Mickens R. E. Nonstandard finite difference schemes, {\it in} Mickens R. E. (ed), {\it Applications of nonstandard finite difference schemes}, World Scientific, Singapore, 2000 \bibitem{Morgenroth:2004} Morgenroth E., Eberl H. J., Van Loosdrecht M.C. M., Noguera D. R., Picioreanu C, Rittmann BE, Schwarz AO, Wanner O. Results from the single species benchmark problem (BM1), {\it Wat. Sci. Tech.}, 49(11/12):145-154, 2004 \bibitem{Morton:1994} Morton K. W., Mayers D. F., {\it Numerical solution of partial differential equations}, Cambridge Univ. Press, 1994 \bibitem{Pici:2003} Picioreanu C., Van Loosdrecht M. C. M., Use of mathematical modelling to study biofilm development and morphology. In: Lens P et al (eds). {\it Biofilms in Medicine, Industry and Environmental Biotechnology – Characteristics, Analysis and Control}, IWA Publishing. pp 413-437, 2003 \bibitem{Saad} Saad Y., {\it Iterative Methods for Sparse Linear Systems}, 2nd edition, 2000 \bibitem{Stewart} Stewart P. S., Multicellular nature of biofilm protection from antimicrobial agents, in: McBain A et al (eds), {\it Biofilm Communities: Order from Chaos}, BioLine, 2003 \bibitem{Wanner:2006} Wanner O., Eberl H. J., Van Loosdrecht M. C. M., Morgenroth E., Noguera D. R., Picioreanu C, Rittmann BE. {\it Mathematical Modeling of Biofilms}, IWA Publishing, 2006 \bibitem{Watnick} Watnick P., Kolter R., Biofilm -- City of Microbes (Minireview), {\it J. Bacteriol.}, 182(10):2675-2679, 2000 \bibitem{Wilderer:2002} Wilderer P. A., Bungartz H. J., Lemmer H., Wagner M., Keller J., Wuertz S., Modern Scientific Methods and their Potential in Wastewater Science and Technology, {\it Water Res.} 36:370-393, 2002 \bibitem{Xu:1998} Xu K. D., Stewart P. S., Xia F., Huang C.-T., McFeters G. A., Spatial physiological heterogeneity in Pseudomonas aeruginosa biofilm is determined by oxygen availability, {\it Appl. and Env. Microbiol.}, 64(10):4035-4039, 1998 \endthebibliography \end{document}