Time adaptive numerical solution of a highly degenerate diffusion-reaction biofilm model based on regularisation

TL;DR

This paper develops a regularisation-based, adaptive numerical method for solving a complex degenerate diffusion-reaction biofilm model, enabling accurate simulation of biofilm dynamics with moving interfaces and gradient blow-up.

Contribution

It introduces a regularisation approach combined with adaptive integration to effectively solve a highly degenerate PDE system in biofilm modeling.

Findings

The method accurately captures biofilm formation dynamics.

Adaptive time-stepping improves computational efficiency.

Signal diffusion impacts quorum sensing induction.

Abstract

We consider a quasilinear degenerate diffusion-reaction system that describes biofilm formation. The model exhibits two non-linear effects: a power law degeneracy as one of the dependent variables vanishes and a super diffusion singularity as it approaches unity. Biologically relevant solutions are characterized by a moving interface and gradient blow-up there. Discretisation of the PDE in space by a standard Finite Volume scheme leads to a singular system of ordinary differential equations. We show that regularisation of this system allows the application of error controlled adaptive integration techniques to solve the underlying PDE. This overcomes the major limitation of existing methods for this type of problem which work with fixed time-steps. We apply the resulting numerical method to study the effect of signal diffusion in the aqueous phase on quorum sensing induction in a biofilm colony.

Figures

Tables

Table 1. Table 1: Least square error between 𝐔 ϵ superscript 𝐔 italic-ϵ \mathbf{U}^{\epsilon} and the the non-regularised numerical solution with ϵ = 0 italic-ϵ 0 \epsilon=0 , E 0 ϵ := 1 256 × 256 ​ ‖ U ϵ − U 0 ‖ 2 assign superscript subscript 𝐸 0 italic-ϵ 1 256 256 subscript norm superscript U italic-ϵ superscript U 0 2 E_{0}^{\epsilon}:=\frac{1}{256\times 256}\parallel\textbf{U}^{\epsilon}-\textbf{U}^{0}\parallel_{2} , and the error between regularised and Barenblatt solutions, E ϵ := 1 256 × 256 ​ ‖ 𝐔 − 𝐔 ϵ ‖ 2 assign superscript 𝐸 italic-ϵ 1 256 256 subscript norm 𝐔 superscript 𝐔 italic-ϵ 2 E^{\epsilon}:=\frac{1}{256\times 256}\parallel\mathbf{U}-\mathbf{U}^{\epsilon}\parallel_{2} , at t = 1 𝑡 1 t=1 . The cell resolution is N = M = 256 𝑁 𝑀 256 N=M=256 and the tolerance of the ROW method is set at T ​ O ​ L = 1 ​ e − 7 𝑇 𝑂 𝐿 1 𝑒 7 TOL=1e-7 .

$ϵ$	$E_0^{ϵ}$	$E^{ϵ}$
${10}^{-3}$	$0.5456777671\times {10}^{-5}$	$0.2670409320\times {10}^{-6}$
${10}^{-4}$	$0.5458369969\times {10}^{-6}$	$0.2332513375\times {10}^{-6}$
${10}^{-5}$	$0.5483556720\times {10}^{-7}$	$0.2302420952\times {10}^{-6}$
${10}^{-6}$	$0.5489546560\times {10}^{-8}$	$0.2299452153\times {10}^{-6}$
${10}^{-7}$	$0.5511348599\times {10}^{-9}$	$0.2299155176\times {10}^{-6}$

Table 2. Table 2: Error between the Barenblatt solution ( 36 ) and the numerical solution with ϵ = 0 italic-ϵ 0 \epsilon=0 of PME ( 37 ), E 0 N = 1 N 2 ​ ‖ 𝐔 − 𝐔 0 ‖ 2 superscript subscript 𝐸 0 𝑁 1 superscript 𝑁 2 subscript norm 𝐔 superscript 𝐔 0 2 E_{0}^{N}=\frac{1}{N^{2}}\parallel\mathbf{U}-\mathbf{U}^{0}\parallel_{2} , for different number of grid cells at t = 1 𝑡 1 t=1 .

$N=M$	$E_0^N$
$32$	$1.0314434\times {10}^{-3}$
$64$	$1.9763426\times {10}^{-4}$
$128$	$8.7552260\times {10}^{-5}$
$256$	$2.2991221\times {10}^{-5}$
$512$	$6.8993723\times {10}^{-6}$
$1024$	$2.1955107\times {10}^{-6}$

Table 3. Table 3: Dimensionless model parameters of system ( 1 ), taken from [ 37 ] .

Parameter	Symbol	Value
Decay rate	$k$	$0.67$
Monod half saturation	$K_U$	$0.13$
substrate uptake rate	${\nu }_U$	$530$
Nutrient diffusion coefficient	$d_c$	$33$
biomass motility coefficient	$\delta$	${10}^{-8}$

Table 4. Table 4: The error between two consecutive solutions for grids with 2 κ × 2 κ superscript 2 𝜅 superscript 2 𝜅 2^{\kappa}\times 2^{\kappa} and 2 κ − 1 × 2 κ − 1 superscript 2 𝜅 1 superscript 2 𝜅 1 2^{\kappa-1}\times 2^{\kappa-1} cell resolution at t = 2 𝑡 2 t=2 (before merging) and Elapsed CPU Time (ECPUT).

$\kappa$	${E_U^{\kappa }|}_{t=2}$	${E_C^{\kappa }|}_{t=2}$	${ECPUT|}_{t=2}$
5	$0.2818501\times {10}^{-2}$	$0.1102250\times {10}^{-2}$	$0.2403\times {10}^2$
6	$0.5379241\times {10}^{-3}$	$0.1528744\times {10}^{-3}$	$0.1524\times {10}^3$
7	$0.1908778\times {10}^{-3}$	$0.1745382\times {10}^{-4}$	$0.2111\times {10}^4$
8	$0.1352093\times {10}^{-3}$	$0.1209376\times {10}^{-4}$	$0.3424\times {10}^5$

Table 5. Table 5: The error between two consecutive solutions for grids with 2 κ × 2 κ superscript 2 𝜅 superscript 2 𝜅 2^{\kappa}\times 2^{\kappa} and 2 κ − 1 × 2 κ − 1 superscript 2 𝜅 1 superscript 2 𝜅 1 2^{\kappa-1}\times 2^{\kappa-1} cell resolution at t = 6 𝑡 6 t=6 (after merging) and Elapsed CPU Time (ECPUT).

$\kappa$	${E_U^{\kappa }|}_{t=6}$	${E_C^{\kappa }|}_{t=6}$	${ECPUT|}_{t=6}$
5	$0.2238254\times {10}^{-2}$	$0.2687217\times {10}^{-3}$	$0.5431\times {10}^2$
6	$0.5709526\times {10}^{-3}$	$0.2166271\times {10}^{-4}$	$0.4148\times {10}^4$
7	$0.1972021\times {10}^{-3}$	$0.2606688\times {10}^{-5}$	$0.6121\times {10}^4$
8	$0.1095095\times {10}^{-3}$	$0.1136270\times {10}^{-5}$	$0.1161\times {10}^6$

Table 6. Table 6: Dimensionless model parameters of system ( 38 ).

Parameter	Symbol	Value	Source
Biomass decay rate	$k$	$0.67$	$RE$
Monod half saturation	$K_U$	$0.13$	$RE$
substrate uptake rate	${\nu }_U$	$530$	$RE$
Nutrient diffusion coefficient	$d_c$	$33$	$RE$
biomass motility coefficient	$\delta$	${10}^{-8}$	$RE$
AHL diffusion coefficient	$d_s$	$16.5$	$-$
dimerization exponent	$m$	$2.5$	$-$
AHL production rate	$\alpha$	$4500$	$-$
AHL production rate	$\beta$	$45000$	$-$
AHL decay rate	$\psi$	$0.02$	$-$
Reference: RE=[37]

Table 7. Table 7: The values of the coefficients of fitted quadratic function for T i subscript 𝑇 𝑖 T_{i} and M t ​ o ​ t ​ a ​ l , i subscript 𝑀 𝑡 𝑜 𝑡 𝑎 𝑙 𝑖 M_{total,i} and the error between the exact and approximated data.

$T_i\left(L\right)$	${\alpha }_1$	${\beta }_1$	${\gamma }_1$	$RMSErro{r}_T$
$T_1\left(L\right)$	$-0.331332$	$2.04645$	$7.98685$	$0.014949$
$T_2\left(L\right)$	$-0.343759$	$2.10267$	$7.95014$	$0.011887$
$T_3\left(L\right)$	$-0.339661$	$2.09279$	$7.98219$	$0.0138615$
$T_4\left(L\right)$	$-0.253235$	$1.86745$	$8.26481$	$0.0175622$

Equations

$$\left\{\begin{array}[]{lllll}\frac{\partial u}{\partial t}=\nabla(D(u)\nabla u)+\frac{c}{K_{U}+c}u-ku,\\ \frac{\partial c}{\partial t}=d_{c}\Delta c-\frac{\nu_{U}c}{K_{U}+c}u,\end{array}\right.$$

$$D(u)=\delta \frac{u^{\alpha }}{(1-u{)}^{\beta }},\alpha ,\beta \ge 1,d_c\gg \delta >0.$$

$$\left\{\begin{array}[]{lllll}&&\partial_{n}u=\partial_{n}c=0~{}\ \mbox{at}\ x=0,x=L~{}~{}~{}\mbox{and}~{}~{}~{}y=0,\\ &&u=0,~{}~{}~{}~{}c+\lambda\partial_{n}c=1~{}\ \mbox{at}\ y=H,\end{array}\right.$$

$$\left\{\begin{array}[]{lllll}1.~{}c_{0}\in L^{\infty}(\Omega)\cap H^{1}(\Omega),~{}~{}0\leq c_{0}\leq 1,~{}c_{0}\mid_{\partial\Omega}=1,\\ 2.~{}u_{0}\in L^{\infty}(\Omega),~{}F(u_{0})\in H_{0}^{1}(\Omega),\\ 3.~{}u_{0}\geq 0,~{}\|u_{0}\|_{\small{L^{\infty}}}<1,\end{array}\right.$$

$$F(u)=F_{\infty }(u):={\int }_0^u\frac{v^{\alpha }}{(1-v{)}^{\beta }}dv,0\le u<1.$$

$$\left\{\begin{array}[]{lllll}1.~{}u,c\in L^{\infty}(\mathbb{R_{+}}\times\Omega)\cap C([0,\infty),L^{2}(\Omega)),\\ 2.~{}c,F(u)\in L^{\infty}(\mathbb{R_{+}},H^{1}(\Omega))\cap C([0,\infty),L^{2}(\Omega)),\\ 3.~{}0\leq u(t,x,y),c(t,x,y)\leq 1,~{}\|u\|_{\small{L^{\infty}(R_{+}\times\Omega})}<1.\end{array}\right.$$

$$\Vert c_1(t)-c_2(t){\Vert }_{L^1(\Omega )}+\Vert u_1(t)-u_2(t){\Vert }_{L^1(\Omega )}$$

$$\le e^{K}t\left(\Vert c_1(0)-c_2(0){\Vert }_{L^1(\Omega )}+\Vert u_1(0)-u_2(0){\Vert }_{L^1(\Omega )}\right),$$

$$\frac{d}{dt}{\int }_{v_{i,j}}udxdy={\int }_{\partial v_{i,j}}{J}_{n}ds+{\int }_{v_{i,j}}R(c)udxdy,i=1,...,N,j=1,...,M$$

$$U_{i,j}(t):=u(t,x_i,y_j)\approx u\left(t,\left(i-\frac{1}{2}\right)\Delta x,\left(j-\frac{1}{2}\right)\Delta x\right)$$

$$C_{i,j}(t):=c(t,x_i,y_j)\approx c\left((t,\left(i-\frac{1}{2}\right)\Delta x,\left(j-\frac{1}{2}\right)\Delta x\right).$$

$$\frac{d}{dt}{U}_{i,j}=\frac{1}{\Delta x}\left(J_{i+\frac{1}{2},j}+J_{i-\frac{1}{2},j}+J_{i,j+\frac{1}{2}}+J_{i,j-\frac{1}{2}}\right)+R_{i,j}{U}_{i,j},$$

$$J_{i+\frac{1}{2},j}=\left\{\begin{array}[]{ll}\frac{1}{2\Delta x}\big{(}D(U_{i+1,j})+D(U_{i,j})\big{)}(U_{i+1,j}-U_{i,j})&\mbox{for}\quad i<N,\\ -\frac{2}{\Delta x}D(0)U_{N,j}&\mbox{for}\quad i=N,\end{array}\right.$$

$$J_{i-\frac{1}{2},j}=\left\{\begin{array}[]{ll}0&\mbox{for}\quad i=1,\\ \frac{1}{2\Delta x}\big{(}D(U_{i,j})+D(U_{i-1,j})\big{)}(U_{i-1,j}-U_{i,j})&\mbox{for}\quad i>1,\end{array}\right.$$

$$J_{i,j+\frac{1}{2}}=\left\{\begin{array}[]{ll}\frac{1}{2\Delta x}\big{(}D(U_{i,j+1})+D(U_{i,j})\big{)}(U_{i,j+1}-U_{i,j})&\mbox{for}\quad j<M,\\ 0&\mbox{for}\quad j=M,\end{array}\right.$$

$$J_{i,j-\frac{1}{2}}=\left\{\begin{array}[]{ll}0&\mbox{for}\quad j=1,\\ \frac{1}{2\Delta x}\big{(}D(U_{i,j-1})+D(U_{i,j})\big{)}(U_{i,j-1}-U_{i,j}),&\mbox{for}\quad j>1.\end{array}\right.$$

$$\frac{d}{dt}{C}_{i,j}=\frac{1}{\Delta x}\left({\hat{J}}_{i+\frac{1}{2},j}+{\hat{J}}_{i-\frac{1}{2},j}+{\hat{J}}_{i,j+\frac{1}{2}}+{\hat{J}}_{i,j-\frac{1}{2}}\right)-{\hat{R}}_{i,j}{U}_{i,j}$$

$${\hat{R}}_{i,j}=\frac{{\nu }_U{C}_{i,j}}{K_U+C_{i,j}}$$

$$\hat{J}_{i+\frac{1}{2},j}=\left\{\begin{array}[]{ll}\frac{d_{c}}{\Delta x}(C_{i+1,j}-C_{i,j})&\mbox{for}\quad i<N,\\ \frac{d_{c}}{\Delta x}\left(\frac{2\Delta x}{2\lambda+\Delta x}-C_{i,j}\left(1+\frac{\triangle x-2\lambda}{\Delta x+2\lambda}\right)\right)&\mbox{for}\quad i=N,\end{array}\right.$$

$$\hat{J}_{i-\frac{1}{2},j}=\left\{\begin{array}[]{ll}0&\mbox{for}\quad i=1,\\ \frac{d_{c}}{\Delta x}(C_{i-1,j}-C_{i,j})&\mbox{for}\quad i>1,\end{array}\right.$$

$$\hat{J}_{i,j+\frac{1}{2}}=\left\{\begin{array}[]{ll}\frac{d_{c}}{\Delta x}(C_{i,j+1}-C_{i,j})&\mbox{for}\quad j<M,\\ 0&\mbox{for}\quad j=M,\end{array}\right.$$

$$\hat{J}_{i,j-\frac{1}{2}}=\left\{\begin{array}[]{ll}0&\mbox{for}\quad j=1,\\ \frac{d_{c}}{\Delta x}(C_{i,j-1}-C_{i,j}),&\mbox{for}\quad j>1.\end{array}\right.$$

$$\pi :\left\{1,...,N\right\}\times \left\{1,...,M\right\}\to \left\{1,...,NM\right\},(i,j)\mapsto p=(i-1)N+j$$

$$\left\{\begin{array}[]{lllll}\frac{d\textbf{U}}{dt}=\mathfrak{D}(\textbf{U})\textbf{U}+\mathfrak{R}_{U}(\textbf{C})\textbf{U}\\ \frac{d\textbf{C}}{dt}=\mathfrak{L}\textbf{C}-\mathfrak{R}_{C}(\textbf{C})\textbf{U}+\textbf{b}.\end{array}\right.$$

$$D_{\epsilon}(u)=\left\{\begin{array}[]{ll}\delta\epsilon^{\alpha}&u<0\\ \delta\frac{(u+\epsilon)^{\alpha}}{(1-u)^{\beta}}&0\leq u\leq 1-\epsilon\\ \delta\epsilon^{-\beta}&u\geq 1-\epsilon.\end{array}\right.,$$

$$\left\{\begin{array}[]{lll}\frac{d\textbf{U}^{\epsilon}}{dt}&=&\mathfrak{D_{\epsilon}}(\textbf{U}^{\epsilon})\textbf{U}^{\epsilon}+\mathfrak{R}_{U}(\textbf{C}^{\epsilon})\textbf{U}^{\epsilon}\\ \frac{d\textbf{C}^{\epsilon}}{dt}&=&\mathfrak{L}\textbf{C}^{\epsilon}-\mathfrak{R}_{C}(\textbf{C}^{\epsilon})\textbf{U}^{\epsilon}+\textbf{b},\end{array}\right.$$

$$\left\{\begin{array}[]{lllll}\frac{\partial u_{\epsilon}}{\partial t}=\nabla(D_{\epsilon}(u_{\epsilon})\nabla u_{\epsilon})+\frac{c_{\epsilon}}{K_{U}+c_{\epsilon}}u_{\epsilon}-ku_{\epsilon},\\ \frac{\partial c_{\epsilon}}{\partial t}=d_{c}\Delta c_{\epsilon}-\frac{\nu_{U}c_{\epsilon}}{K_{U}+c_{\epsilon}}u_{\epsilon},\end{array}\right.$$

$$1\le {\text{U}}^{\theta }\le 1+\text{q}$$

$${\text{U}}^{ϵ}(0)\le {\text{U}}^{\theta },$$

$$\frac{d{\text{U}}^{\theta }}{dt}$$

$$\frac{d{\mathbf{V}}^{ϵ}}{dt}={\mathfrak{D}}_{ϵ}({\mathbf{V}}^{ϵ}){\mathbf{V}}^{ϵ}+\mathbf{Q}+\mathbf{b},{\mathbf{V}}^{ϵ}(0)=\rho 1$$

Full text

Discretisation of the PDE in space by a standard Finite Volume scheme leads to a singular system of ordinary differential equations. We show that regularisation of this system allows the application of error controlled adaptive integration techniques to solve the underlying PDE. This overcomes the major limitation of existing methods for this type of problem which work with fixed time-steps.

We apply the resulting numerical method to study the effect of signal diffusion in the aqueous phase on quorum sensing induction in a biofilm colony.

Keywords: Biofilm, Degenerate diffusion-reaction equation, Quorum sensing, Regularization, Semi-discretization, Time adaptivity

MSC: primary: 35K65, 65M08; secondary: 68U20,92D25

1 Introduction

Bacterial biofilms are microbial communities on immersed surfaces, embedded in layers of self-produced extracellular polymeric substances (EPS), which protect the sessile cells against mechanical or chemical washout [6, 18, 20, 27, 30, 53]. Biofilms are prevalent in natural, industrial and hospital settings. Depending on the context, they can be harmful or beneficial. For example biofilms can lead to corrosion problems in fresh water pipes, and oil pipelines [44]; foremost cause of failure of medical implants are bacterial infections caused by biofilm formation [2, 42]; biofilms can lead to crop disease in plants [7]; dental plaque is an oral biofilm on teeth that causes dental disease [1]. On the other hand, the adsorption and absorption properties and enhanced mechanical stability of biofilms make them advantageous to environmental engineering technologies, e.g. for waste water treatment, remediation of contaminated soil, and elimination of petroleum oil from contaminated oceans or marine systems [51].

Although the term biofilm suggests a homogeneous film-like layer, biofilms on the meso-scale ($10\mu m\sim 1mm$, the actual biofilm scale) in reality are spatially heterogeneous assemblages of colonies which may merge as they grow and expand [28]. The architecture of a biofilm depends on some biological factors such as maximum cell density, specific growth rate, and local nutrient availability. In environments with unlimited amount of substrate, biofilms tend to grow quickly and form homogeneous, compact, thick layer. On the other hand, biofilms in nutrient limited regimes tend to form in patchy heterogeneous structures [9].

Several mathematical models have been proposed in the literature to describe the growth of spatially heterogeneous biofilms, ranging from stochastic individual based models to cellular automata models to deterministic continuum models, cf [26, 28, 50, 51] and the references that they cite. 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 behaviour in predicting the development of biofilm morphologies in response to substrate limitation [51].

An early prototype biofilm growth model has been introduced in [9], and was extended later to account for further biofilm aspects and processes, such as quorum sensing, response to antibiotics, internally triggered dispersal, EPS production, and multispecies systems, e.g. in [11, 12, 14, 15, 16, 24, 29, 38]. The model has been derived both from the view point of biofilms as spatially structured populations [25] and as a fluid [15, 34], thus reflecting the ecological-mechanical duality of biofilms. In its basic form as a biofilm growth model, it consists of a density-dependent diffusion-reaction equation for biomass that is coupled with a semi-linear diffusion-reaction equation for growth limiting nutrients. The biomass equation shows two interacting non-linear diffusion effects: (i) the diffusion coefficient vanishes where the biomass vanishes as in the porous medium equation, and (ii) a super-diffusion singularity as the dependent variable in the diffusion coefficient approaches the known maximum density. The interplay of both effects assures that the solution of the biomass equation is bounded by the maximum cell density [13], that spatial expansion of the biofilm only takes place if locally no space is available to produce new biomass (i.e. a volume filling effect), and that initial data with compact support have solutions with compact support [9]. The latter property is a characteristic property of the porous medium equation and related problems. In our application this expresses itself as a sharp biofilm/water interface, at which the biomass gradient blows up. Since biofilm colonies grow, their interfaces with the aqueous phase are not stationary but change over time. Eventually neighbouring colonies can merge into a bigger colony, in which case their interfaces merge and dissolve.

Interface problems of this type pose difficulties for numerical and analytical methods. The steep biomass gradients at the interface can lead to interface smearing or spurious oscillations in numerical solutions. The singularity in the density-dependent diffusion coefficient leads to a blow up of model coefficients and forces the numerical method to work with very small time-steps. Several methods have been proposed in the literature for the numerical treatment of this highly non-linear degenerate model.

In [9] a hybrid time-integration strategy was applied. The slower biomass processes were solved explicitly with a Runge-Kutta-Fehlberg method, whereas the faster nutrient processes were solved implicitly. This required inefficiently small time steps for the biomass equation when the biomass density approaches maximum density somewhere. A different idea was explored in [8] in which the degenerate equation is transformed into new dependent variables such that the spatial operator converted into the Laplacian operator and all non-linear effects appeared in the time-derivative. A fully implicit, fixed time-step method for the new equation could handle the diffusion singularity effect very well, but emphasized interface oscillations. In [22] a fully adaptive in space and time method was proposed in the 1D case, based on a weak moving frame formulation. This method explicitly tracked the biofilm/water interface and was able to deal well with the porous medium degeneracy but had problems with the super-diffusion singularity. In [23] the transformation idea of [8] was combined with the moving frame approach of [22]. This method worked very well in 1D but in the 2D case becomes cumbersome; moreover, it was not obvious how it could be adapted for multi-species biofilm systems.

A semi-implicit method that is based on a non-local (in time) representation of the non-linear density dependent biomass flux is described in [10, 43]. This method can be efficiently parallelised for shared memory architectures [33] and has been used for several extensions and applications of the underlying biofilm model [11, 12, 14, 15, 16, 24, 37]. This low order method has shown to be robust and to handle the steep biomass gradients at the biofilm/water interface reasonably well. This method has also been studied and used for a variant of the biofilm model with a slightly different biomass diffusion coefficient in [34]. Minor adaptations and extensions of this approach were suggested in [29, 32], however without major documented improvement. A related fully-implicit first order time integration method was proposed in [21] for a cellulosic biofilm system in which nutrients are stationary and do not diffuse. It uses a fixed-point iteration at each time-step to solve the resulting non-linear system. This method reduces to the semi-implicit approach of [10] as a special case if only one iteration step is carried out. In [3] the implicit trapezoidal rule was used for time integration in the spatially one-dimensional case and explicit variants were proposed in [31, 45], however, without demonstrated gain. All these methods use fixed time steps and because of the low regularity of the solution they do not possess error control properties for time integration. Choosing the time step such that an acceptable trade-off between accuracy and compute time is achieved, and that a breakdown of the method near the super-diffusion singularity is avoided, requires some user experience.

In this paper we want to overcome the constraints around time integration. Our problem at hand has two sources of stiffness. One stems from the disparity of characteristic time scales of substrate diffusion and uptake on one hand, and biomass growth on the other hand [9]. The other one comes into play when the biomass density approaches the super-diffusion singularity. For stiff ordinary differential equations with sufficiently smooth solutions, a host of such methods exists. Due to the singularity in the biomass diffusion coefficient and the resulting missing regularity of solutions, it is not a priori clear whether these can be applied for the time integration of the biofilm equation. Furthermore, it will be important that the resulting numerical scheme is able to deal with sharp biomass gradients at the biofilm/water interface without introducing spurious oscillation or extensive interface smearing. To address these questions, and to propose a fully adaptive method for time integration of the biofilm model is our main objective.

We will demonstrate the utility of this method by studying a question that arises in the context of quorum sensing in biofilms. Bacterial cells produce and secrete chemical signals to communicate with each other, and measure the concentration of the signal molecules in their environment. This mechanism is known as quorum sensing. Once a critical signal concentration is reached, changes in gene expressions are induced and bacteria collectively synchronize their behaviour [16, 19]. Since the producer cells respond to their own signal these are sometimes referred to as autoinducers (e.g., acyl-homoserine lactones (AHL) in Gram-negative bacteria). In most bacterial autoinducer systems, the autoinducer synthase gene is upregulated, i.e a positive feedback stimulus is induced, leading to production of AHL molecules at an increased rate [16].

The name quorum sensing stems from an interpretation of autoinduction as a mechanism to estimate cell density. However, in environments that are not completely mixed, such as spatially structured biofilm systems, the autoinducer concentration is affected by transport of chemical signals in the aqueous environment. In fact they might use the autoinducer concentration as an estimate of any of the factors that can affect the concentration of this molecule, such as diffusion limitation [19]. This has lead to the interpretation of quorum sensing as diffusion sensing and a mechanism to survey the environment [40]. An important question to ask is then how the environment affects the time to autoinduction in a biofilm colony, and the size of the colony that is required for autoinduction, i.e. for the signal concentration to exceed induction threshold.

In a recent study [46] it was investigated experimentally using a synthetic biofilm matrix and mathematically using a simple diffusion model, how much autoinducer needs to be added instantaneously in a colony center to achieve concentrations in the biofilm above induction threshold for different colony sizes. Our goal is to build on this using a model that accounts also for biofilm growth and dynamic autoinducer production, and that also considers diffusion of signal molecules into the environment. Related, older studies are [4, 5], where the relationship between biofilm size and induction was studied in a one-dimensional setting, using a Wanner-Gujer type biofilm growth model. The 1D setting in that study, however, is limiting in that it does not allow to investigate effects of signal diffusion parallel to the substratum, and thus effects of neighboring colonies or confined domains. How convective transport of signal molecules due to bulk flow hydrodynamics affects washout of autoinducers and the onset of induction in narrow flow channels, on the other hand, was simulated in [15, 48]; the presence of background flow added another complexity that does not allow to hone in on aspects of diffusion sensing.

Our paper is organized as follows: In Section 2 we introduce the prototype biofilm model and obtain its semi-discrete version by applying a Finite Volume Method to discretize the partial differential equations in space. In Section 3 we analyze the semi-discrete model; in particular we show, using regularization and the method of super and sub solutions, that the solution of the resulting ODE never attains the singularity and is thus sufficiently smooth to apply standard numerical methods for time integration, such as embedded Rosenbrock-Wanner methods, which we will use. In order to demonstrate that the resulting numerical method is able to handle gradient blow-up at the biofilm/water interface, we apply in section 4 the numerical scheme to a simpler problem (without super-diffusion singularity), for which an exact solution is known, namely the porous medium equation with linear growth. In section 5 we investigate the behaviour of our numerical method for the full prototype biofilm growth model. In particular we give a grid refinement study and explore quantitatively the dependence of the numerical method on the regularization parameter that was introduced for analytical treatment. In section 6, finally we present an illustrative application of the method in the context of quorum sensing. Concluding remarks are provided in Section 7.

2 Mathematical Model

2.1 Governing equations

The mathematical model studied in this paper was originally introduced in [9]. It is formulated as a density-dependent degenerate diffusion-reaction equation over domain $\Omega\subset\mathbb{R}^{2}$. The dependent variables are the volume fraction occupied by biomass, $u$, and the concentration of a growth limiting nutrient, $c$. Domain $\Omega$ is divided into region $\Omega_{1}(t)=\left\{(x,y)\in\Omega\subset\mathbb{R}^{2}:u(t,x,y)=0\right\}$ that describes the aquatic phase (bulk liquid, channels and pores of a biofilm) without biomass, and region $\Omega_{2}(t)=\left\{(x,y)\in\Omega\subset\mathbb{R}^{2}:u(t,x,y)>0\right\}$, which is the actual biofilm with positive density, cf. Figure 1.

The independent variables $t\geq 0$ and $(x,y)\in\Omega$ denote here time and spatial location, respectively. Both regions are separated by the biofilm/water interface $\Gamma(t)=\partial\Omega_{1}\cap\partial\Omega_{2}$ which changes as the biofilm grows, i.e. as $u$ changes.

As is usual in mathematical models of biofilms, the EPS that is produced by the bacteria is subsumed in the biofilm volume fraction in this prototype biofilm model. An extension of this model that accounts for EPS explicitly has been proposed, for example, in [16].

The model in non-dimensionalised form reads [9]:

[TABLE]

$D(u)$ is a density-dependent diffusion coefficient defined as:

[TABLE]

Diffusion coefficient $D(u)$ shows two non-linear effects: (i) a power law degeneracy as in the porous medium equation, i.e. $D(u)$ vanishes as $u$ vanishes and (ii) a super diffusion singularity as $u$ approaches to unity. The porous medium degeneracy, $u^{\alpha}$, guarantees that the biofilm does not spread notably if the biomass density is small, $u\ll 1$, and it is also responsible for the formation of a sharp interface between biofilm and surrounding liquid, i.e. initial data with compact support lead to solutions with compact support. The second effect (ii) enforces the solution to be bounded by unity [10, 51]. This is counteracted by the degeneracy as $u=0$ at the interface. Consequently, $u$ squeezes in the biofilm region and approaches its maximum value 1. Hence, the interaction of both non-linear diffusion effects with the growth term is needed to describe spatial biomass spreading [10]. Diffusion of the dissolved nutrient is assumed to be Fickian, i.e. $d_{c}$ is a constant. Biomass spreading is much slower than diffusion of dissolved substrate, thus the biomass motility coefficient, $\delta$, is several orders of magnitude smaller than the substrate diffusion coefficient.

In (1), parameter $k$ is the cell lysis rate, $K_{U}$ the Monod half saturation concentration, and $\nu_{U}$ the maximum substrate uptake rate.

To study the mathematical model (1) we restrict ourselves to the rectangular domain $\Omega=[0,L]\times[0,H]$. The substratum, on which biofilm colonies form is the bottom boundary, $y=0$. We assume the substratum is impermeable to biomass and dissolved substrate. At the lateral boundaries, $x=0$ and $x=L$, we assume a symmetry boundary condition for both dependent variables, which allows us to view the domain as a part of a continuously repeating larger domain. At the top boundary, $y=H$, we pose a homogeneous Dirichlet condition for the biomass and a Robin condition for the nutrient which reflects that the substrate is added to the system through this segment of the domain boundary. Thus the imposed boundary conditions on domain $\Omega=[0,L]\times[0,H]$ are defined as:

[TABLE]

Here $\lambda$ can be understood as the external concentration boundary layer thickness, and $\partial_{n}$ denotes the outward normal derivative. The concentration boundary layer is introduced to mimic the convective contribution of external bulk flow to substrate supply, i.e. it is related to the bulk flow velocity, in the sense that a small bulk flow velocity implies a thick concentration boundary layer, while a thin concentration boundary layer represents fast bulk flow [11].

For the reader’s convenience we recall the main result of [13] that studies the longtime behaviour of PDE model (1), without proof:

Theorem 2.1** (Existence).**

Let the initial data $(u_{0},c_{0})$ satisfy the conditions

[TABLE]

where

[TABLE]

Then, there exists a solution $(u(t),c(t))$ of degenerate problem (1) (in the sense of distributions) belonging to the following class:

[TABLE]

Theorem 2.2** (Uniqueness).**

Let $(u_{1}(t),c_{1}(t))$ and $(u_{2}(t),c_{2}(t))$ be two solutions of (1) belonging to the class (6). Then, the following estimate is valid:

[TABLE]

for some positive constant $K$ that depends on reaction parameters. In particular, the solution of (1) is unique in the class (6).

An extension of these results to more general boundary conditions is described in [13].

2.2 Spatial discretization

We introduce a uniform grid of size $N\times M$ for the rectangular domain $[0,L]\times[0,H]$. Integrating the first equations of (1) over each grid cell and using the Divergence Theorem yields

[TABLE]

where $v_{i,j}$ denotes the domain of the cell with grid index $(i,j)$, $J_{n}=D(u)\partial_{n}u$ denotes the outward normal flux across the grid cell boundary, and $R(c)=c/(K_{U}+c)-k$ stands for the reaction terms.

To evaluate the area integrals in (2.2), we evaluate the dependent variables at the center of the grid cells,

[TABLE]

and similarly for the nutrient concentration

[TABLE]

for $i=1,...,N$ and $j=1,...,M$ with $\Delta x=L/N=H/M$, and we approximate the integrals by the midpoint rule. Similarly, the line integral in (2.2) is evaluated by considering every edge of the grid cell separately, using the midpoint rule. To this end, the diffusion coefficient $D(u)$ in the midpoint of the cell edge is approximated by arithmetic averaging from the neighbouring grid cell center points, and the derivative of $u$ across the cell edge by a central finite difference. We get then for the biomass density in the grid cell center the ordinary differential equation

[TABLE]

where $R_{ij}=C_{i,j}/(K_{U}+C_{i,j})-k$ and for the fluxes we have, accounting for the boundary conditions,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The spatial discretization of the equation for the nutrient follows the same principle, with the major difference being that at the top of the domain we have a Robin boundary condition instead of homogeneous Dirichlet condition. That $d_{c}$ is constant simplifies the flux expressions. We have then

[TABLE]

with

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

By introducing the lexicographical grid ordering

[TABLE]

and the vector notation $\textbf{U}=(U_{1},...,U_{NM})$, $\textbf{C}=(C_{1},...,C_{NM})$ with $U_{p}:=U_{\pi(i,j)}=U_{i,j}$, $C_{p}:=C_{\pi(i,j)}=C_{i,j}$ for $i=1,...,N$, $j=1,...,M$, we arrive at the coupled system of $2\cdot N\cdot M$ ordinary differential equations

[TABLE]

Remark 2.3**.**

By construction the $NM\times NM$ matrices $\mathfrak{D}(U),\mathfrak{L}$ are symmetric, and weakly diagonally dominant with non-positive main diagonals and non-negative off-diagonals. They contain the spatial derivative terms. They can efficiently be stored in sparse diagonal format with 4 off-diagonals, with offsets $\pm 1,\pm M$. The matrices $\mathfrak{R}_{U}(\textbf{C}),\mathfrak{R}_{C}(\textbf{C})$ are diagonal matrices which contain the reaction terms; in particular the $p$th entry in $\mathfrak{R}_{C}(\textbf{C})$ vanishes if $C_{p}=0$. The vector b contains contributions from the Robin boundary conditions in (18); its entries are zero for all grid cells $(i,j)$ with $i<N$, and $b_{\pi(i,j)}=\frac{d_{c}}{\Delta x}\frac{2}{2\lambda+\Delta x}>0$ for grid points with $i=N$.

Remark 2.4**.**

The discretization can easily be extended to more general boundary conditions, such as Robin conditions everywhere, in analogy to (18).

3 The regularized semi-discrete system

Due to the singularity in the biomass diffusion coefficient (2) standard arguments from the theory for ordinary differential equations, such as the Picard-Lindelöf theorem and the tangent condition cannot be readily applied to study the well-posedness of the ODE system (23) and to prove that the biomass volume fraction U is indeed bounded by unity, as is required for physical reasons and for agreement with the underlying PDE (1), as per Theorem 2.1. To overcome this limitation we regularize the semi-discretised system. Following [13], we introduce the regularised density-dependent biomass diffusion coefficient

[TABLE]

where the continuous extension of $D_{\epsilon}$ for negative $u$ is for technical reasons. We show lateron, that indeed $u\geq 0$ for the solutions of our model. We define the following regularized version of ODE system (23),

[TABLE]

where matrix $\mathfrak{D_{\epsilon}}(\cdot)$ is defined like $\mathfrak{D}(\cdot)$, only with $D(u)$ in the flux terms (12)-(15) replaced by $D_{\epsilon}(u)$ as defined in (24).

Remark 3.1**.**

The regularised semi-discretised system (25) is a standard, 2nd order (in space) convergent, finite difference approximation of the regularised partial differential equation system

[TABLE]

the solutions of which are smooth. It was shown in [13] that the solutions of the initial value problem of (26) converge to the solution of the corresponding initial value problem of (1) as $\epsilon\rightarrow 0$, and that $u_{\epsilon}\leq 1-\xi$ for some $\xi>0$ and sufficiently small $\epsilon$.

Going forward, we define vector inequalities component wise, i.e for two vectors $\textbf{U}=(U_{1},...,U_{n})^{T}$ and $\textbf{V}=(V_{1},...,V_{n})^{T}$ the inequality $\textbf{U}\leq\textbf{V}$ means $U_{p}\leq V_{p}$ for all $p=1,...,n$, and accordingly for strict inequalities. Furthermore, by $\mathbf{1}\in\mathbb{R}^{n}$ we denote the unity vector $(1,...,1)^{T}$, and by $\mathbf{0}\in\mathbb{R}^{n}$ we denote the vector $(0,...,0)^{T}$.

We show that the initial value problem of (25) with initial data for $\textbf{U}^{\epsilon}$ and $\textbf{C}^{\epsilon}$ such that $0\leq\textbf{U}^{\epsilon}(0)<\textbf{1}$ and $0\leq\textbf{C}^{\epsilon}(0)\leq\textbf{1}$ posseses a unique solution that is non-negative and bounded.

Proposition 3.2**.**

Let $0<\epsilon\ll 1$. Suppose the initial values $\textbf{U}^{\epsilon}(0)$ and $\textbf{C}^{\epsilon}(0)$ of the regularized ODE system (25) are non-negative and satisfy $\textbf{U}^{\epsilon}(0)\leq(1-\rho)\mathbf{1}$ with $\rho\in(0,1)$ and $\textbf{C}^{\epsilon}(0)\leq\mathbf{1}$. Then there exists a unique solution $(\textbf{U}^{\epsilon},\textbf{C}^{\epsilon})$ of the regularized problem (25), which is non-negative and bounded by a constant for sufficiently small values of $\epsilon$. Also, an upper bound on the solution exists that does not depend on the regularization parameter $\epsilon$.

Proof.

The semi-discrete system (25) satisfies the Lipschitz condition in the non-negative cone. Therefore, the initial value problem has a unique solution. Non-negativity of this solution follows with standard arguments, such as the tangent condition (cf [49]), from the properties of matrices $\mathfrak{R}_{C},\mathfrak{R}_{U},\mathfrak{L},\mathfrak{D}_{\epsilon}$ as per Remark 2.3. That $\textbf{C}^{\epsilon}\leq\mathbf{1}$ follows with the non-negativity of $\mathbf{U}^{\epsilon}$ from the definitions of $\mathfrak{R}_{C},\mathfrak{L},\textbf{b}$ again from the tangent condition.

To show the boundedness of $\textbf{U}^{\epsilon}$ we introduce the barrier function $\textbf{U}^{\theta}=\mathbf{1}+\vec{\theta}$, where the vector $\vec{\theta}$ is the solution of the linear system $\mathcal{A}\vec{\theta}=\mathbf{1}$ with $\mathcal{A}=-\frac{1}{\delta\epsilon^{\alpha}}\mathfrak{D}_{\epsilon}(0)$. Note that this matrix $\mathcal{A}$ does not depend on $\epsilon$, whence also $\vec{\theta}$ is independent of $\epsilon$.

By construction, $\mathcal{A}$ is invertible and weakly diagonally dominant with positive main diagonal entries and negative off-diagonal entries. Therefore, it is an M-matrix [17]. Hence its inverse exists and is non-negative, which implies the non-negativity of $\vec{\theta}.$ Thus,

[TABLE]

for some positive vector q. Since $\vec{\theta}$ does not depend on $\epsilon$, neither does q. The initial data obviously satisfy

[TABLE]

and for sufficiently small values of $\epsilon$ we have

[TABLE]

Note that $\mathcal{A}$ is a weakly diagonally dominant matrix such that $\mathcal{A}\mathbf{1}\geq 0$. Consequently, $\textbf{U}^{\theta}$ is an upper solution of (25) for sufficiently small $\epsilon$, wherefore then $\textbf{U}^{\theta}\geq\textbf{U}^{\epsilon}$, by the comparison theorem for quasimonotonic systems [49]. Thus, inequality (29) implies that there exists $\hat{\epsilon}$ such that the solutions $\textbf{U}^{\epsilon}$ are uniformly bounded for all $0<\epsilon\leq\hat{\epsilon}.$ ∎

We will use this result to improve the upper bound on $\mathbf{U}^{\epsilon}$. In particular we will show that $\mathbf{U}^{\epsilon}<\mathbf{1}$ for small enough $\epsilon$.

According to Proposition 3.2, there is a constant $Q$ independent of $\epsilon$, such that $Q\geq\|\mathfrak{R}_{U}(\mathbf{C}^{\epsilon})\mathbf{U}^{\epsilon}\|_{\infty}$. As before, we assume that there is a $0<\rho<1$ such that $\mathbf{U}^{\epsilon}(0)\leq\rho\mathbf{1}$. Let us denote by $\textbf{V}^{\epsilon}$ the solution of

[TABLE]

where vector $\mathbf{Q}=Q\mathbf{1}$ and vector $\mathbf{b}$ is a boundary correction that one obtains if imposing the inhomogeneous Dirichlet boundary condition $u=\rho$ at $y=H$ instead of the homogeneous one, defined in analogy with equation (23) and Remark 2.3. This ensures that at the boundary the $\textbf{V}^{\epsilon}$ bounds $\textbf{U}^{\epsilon}$ from above.

We have the following relationship between $\mathbf{U}^{\epsilon}$ and $\mathbf{V}^{\epsilon}$:

Lemma 3.3**.**

Let $\mathbf{V}^{\epsilon}$ be the solution of (30) and $\mathbf{U}^{\epsilon}$ be the solution of (25) with initial data $\mathbf{U}^{\epsilon}(0)\leq\rho\mathbf{1}$. Then $\mathbf{V}^{\epsilon}(t)>\mathbf{U}^{\epsilon}(t)\geq\mathbf{0}$ for all $t>0$.

Proof.

This follows with standard comparison theorems for essentially positive ordinary differential equation systems, cf [49], applied to the first equation in (25). By introducing the defect $\mathfrak{P}:=\frac{d}{dt}-\mathfrak{D}_{\epsilon}(.)-\mathfrak{R}_{U}(\cdot)$, we have $\mathfrak{P}\mathbf{V}_{\epsilon}(t)>\mathfrak{P}\mathbf{U}_{\epsilon}(t)\equiv 0=\mathfrak{P}\textbf{0}$, which gives $\mathbf{V}^{\epsilon}(t)>\mathbf{U}^{\epsilon}(t)\geq\mathbf{0}$. ∎

For each vector index $p\in\left\{1,...,NM\right\}$ we denote by $(i_{p},j_{p})$ the corresponding location on the $N\times M$ grid, i.e.

[TABLE]

where $\pi^{-1}$ is the inverse of the grid ordering map (22).

The following observation allows us to reduce the analysis of the 2D system of $NM$ ordinary differential equations to a 1D problem consisting of $M$ ordinary differential equations, which will simplify notation.

Lemma 3.4**.**

The coefficients $V_{\epsilon,p}$ of the solution $\mathbf{V}^{\epsilon}$ of (30) satisfy $V_{\epsilon,p}(t)=V_{\epsilon,j_{p}}(t)$ for $t>0$ and $p=1,...,NM$.

Proof.

This follows from the symmetry of boundary conditions of (30), the spatially homogeneous initial data and that the source term is a constant in (30), and from the uniqueness of solutions of (30). ∎

With this observation and the definition of $\pi$ in (22), we have $V_{\epsilon,\pi(i,j)}=V_{\epsilon,\pi(\hat{i},j)}=V_{\epsilon,j}$ for all $\hat{i}=1,...,i-1,i+1,...,N$. Thus the ordinary differential equations determining the functions $V_{\epsilon,j}(t)$ reduce to

[TABLE]

For given $\epsilon,0<\epsilon<1$, we furthermore define $\mathbf{Z}^{\epsilon}:=(Z_{\epsilon,1},...,Z_{\epsilon,M})^{T}\in\mathbb{R}^{M}$ as the solution of system of linear algebraic equations

[TABLE]

Lemma 3.5**.**

Let $\mathbf{V}^{\epsilon}(t)$ be the solution of (30) and $\mathbf{Z}^{\epsilon}$ be the solution of the linear system (33). If $\rho\leq 1-\epsilon$, then $\mathbf{V}^{\epsilon}(t)\leq\mathbf{Z}^{\epsilon}$ for all $t>0$.

Proof.

We use the usual invariance theorem, cf. [49], to show that the set $\mathcal{Z_{\epsilon}}=\prod\limits_{i=1}^{M}[\rho,Z_{\epsilon,i}]$ is positively invariant under the differential equation (32): In the equation for $V_{\epsilon,j}$ substitute $V_{\epsilon,j}=\rho$ and assume $\rho\leq V_{\epsilon,j\pm 1}\leq Z_{\epsilon,j\pm 1}$. Then $\frac{dV_{\epsilon,j}}{dt}\geq 0$, wherefore $V_{\epsilon,j}(t)\geq\rho$. Similarly, in the equation for $V_{\epsilon,j}$ substitute $V_{\epsilon,j}=Z_{\epsilon,j}$ and assume $\rho\leq V_{\epsilon,j\pm 1}\leq Z_{\epsilon,j\pm 1}$. Using $D_{\epsilon}(V_{\epsilon,j\pm 1})\leq\delta\epsilon^{-\beta}$ and the definition of $Z_{\epsilon,j}$ in (33) gives $\frac{dV_{\epsilon,j}}{dt}\leq 0$, wherefore $V_{\epsilon,j}(t)\leq Z_{\epsilon,j}$. Thus, the solution of the initial value problem (32) remains in $\mathcal{Z}_{\epsilon}$ which proves the assertion. ∎

Lemma 3.6**.**

Let $\mathbf{Z}^{\epsilon}$ be the solution of the linear system (33). For small enough $\epsilon$ and $\Delta x$, the inequality $\mathbf{Z}^{\epsilon}<\mathbf{1}$ holds.

Proof.

By straightforward calculation we have

[TABLE]

for $j=1,...,M-1$, where $k_{1}=0$ and $k_{j+1}:=k_{j}+j$. Therefore, $Z_{\epsilon,1}>...>Z_{\epsilon,M}$. Finally we have from the last equation of (33)

[TABLE]

For small enough $\epsilon$ and $\Delta x$, the second term on the right hand side is dominated by the first term, i.e. $1>Z_{\epsilon,1}>...>Z_{\epsilon,M}$. ∎

Theorem 3.7**.**

Let $\mathbf{U}^{\epsilon}(t)$ be a solution of (25) with $\mathbf{U}^{\epsilon}(0)<\rho\mathbf{1}$ for some $\rho$ with $0<\rho<1$. Then for sufficiently small $\epsilon$ and $\Delta x$ it holds $\mathbf{U}^{\epsilon}(t)<\mathbf{1}$.

Proof.

Proposition 3.2 and Lemmas 3.3-3.6 imply directly that for small enough $\epsilon$ and $\Delta x$ the following inequalities hold

[TABLE]

where $\hat{\mathbf{Z}}^{\epsilon}=(\hat{Z}_{\epsilon,1},...,\hat{Z}_{\epsilon,NM})^{T}$ with $\hat{Z}_{\epsilon,\pi(i,j)}:=Z_{\epsilon,j}$ for $1\leq i\leq N$, $1\leq j\leq M$. ∎

Theorem 3.8**.**

Let $\mathbf{U}(0)\leq\rho\mathbf{1}$ for some $0<\rho<1$, and let $T>0$. For $\epsilon>0$, let $(\mathbf{U}^{\epsilon}(t),\mathbf{C}^{\epsilon}(t))$ be the solution of (25) with $\mathbf{U}^{\epsilon}(0)=\mathbf{U}(0),\mathbf{C}^{\epsilon}(0)=\mathbf{C}(0)$. Then $(\mathbf{U}^{\epsilon}(t),\mathbf{C}^{\epsilon}(t))$, $0\leq t\leq T$ converges to the solution of (23) as $\epsilon\rightarrow 0$.

Proof.

We have shown that $\textbf{U}^{\epsilon}(t)$ and $\textbf{C}^{\epsilon}(t)$ are bounded for small values of $\epsilon$ and $\Delta x$. Boundedness of $\textbf{U}^{\epsilon}(t)$ and $\textbf{C}^{\epsilon}(t)$ implies that there is a sequence $\epsilon_{n}$ that goes to zero when $n\rightarrow\infty$ such that $\textbf{U}^{\epsilon_{n}}(t)$ and $\textbf{C}^{\epsilon_{n}}(t)$ converge to some $\textbf{U}(t)$ and $\textbf{C}(t)$ component wise. We can now pass to the limit $\epsilon\rightarrow 0$ in the regular semi-discrete system (25) and claim that $(\textbf{U},\textbf{C})$ are the desired solutions of degenerate semi-discrete problem (23). In order to prove this hypothesis for U it is sufficient to verify that $\mathfrak{D}_{\epsilon}(\textbf{U}^{\epsilon})\textbf{U}^{\epsilon}\rightarrow\mathfrak{D}(\textbf{U})\textbf{U}$ if $\epsilon\rightarrow 0$. We have

[TABLE]

Boundedness of $\textbf{U}^{\epsilon}(t)$ and continuity of $D_{\epsilon}(u)$ imply that each term in the right hand side of inequality (34) tends to zero if $\epsilon\rightarrow 0$. This proves that $\textbf{U}(t)$ is indeed the solution of the biomass equation of the degenerate semi-discrete system (23). In the same manner we can also show that $\textbf{C}(t)$ is the solution of the substrate equation of (23).

∎

4 Numerical method

4.1 Time integration

In order to obtain a numerical approximation of (1) we will solve the regularized semi-discrete problem, (25) to approximate (26), and pass $\epsilon$ to 0. The regularized system (25) satisfies a Lipschitz condition and has bounded solutions. This suggests the application of standard ODE solvers.

In typical biofilm applications, different time scales for the substrate equation and the biomass equation induce stiffness, which can be exacerbated by non-linear diffusion effects if the biomass approaches unity somewhere. We use a time adaptive, error controlled, embedded Rosenbrock-Wanner method (ROW), more specifically ROS3PL, a third order method with 4 stages [39]. Rosenbrock-Wanner methods require the solutions of linear systems. In our case the coefficient matrix of these systems are sparse and non-symmetric. We use the stabilised bi-conjugate gradient method [47] to solve them. More specifically we use a routine from the SPARSKIT library [41], that was prepared for parallel execution using OpenMP in [33].

4.2 Spatial discretisation and gradient blow-up at the interface

The ordinary differential equation (23) was derived by applying a space discretisation to the underlying partial differential equation (1). Since at the biofilm/water interface the solutions of (1) are known to have blow-up of the gradient, the question arises to which extent our spatial discretisation introduces smearing around the interface. We investigate this, applying our method to a related test problem with known exact solution, that has the same phenomenon. Our test problem is the Porous Medium Equation with linear source term [36],

[TABLE]

which admits the self-similar Barenblatt solution [36]

[TABLE]

with $\tau=\frac{1}{km}e^{kmt}$ and $k_{0}^{2}=r_{0}^{2}(\frac{1}{km}e^{kmt_{0}})^{-(m+1)^{-1}}$, where $r_{0}$ is the radius of the initial spherical colony and $t_{0}$ is the initial time. Notation $u_{+}$ is defined as $u_{+}:=\max(u,0)$. This solution is induced by a Dirac delta function as initial data. For our numerical test we chose as initial data function (36) evaluated at time $t_{0}=0.1$ and $r_{0}=0.1$. The values of $m$ and $k$ used in our simulations are chosen as $m=4,~{}k=3$.

For our test simulations, we pose the homogeneous Neumann boundary condition on the boundaries and initially one semi-spherical colony is placed in the center of the domain. Note that (36) and the solution of the biofilm model satisfy this boundary condition as well as the homogeneous Dirichlet condition up until the moving interface reaches the boundary somewhere. We terminated the simulations before this occurs.

After regularisation, (35) becomes

[TABLE]

in analogy to the regularization in (24).

We solve (37) on a uniform grid of size $256\times 256$ for several choices of $\epsilon$ and compare the results against the exact solution (36) and the numerical solution obtained with the choice $\epsilon=0$, i.e. the solution of (35) in Table 1. Reported are the 2-norm differences between the regularised numerical solutions and the Barenblatt solution, $E^{\epsilon}:=\frac{1}{256\times 256}\parallel\mathbf{U}-\mathbf{U}^{\epsilon}\parallel_{2}$, and between $\mathbf{U}^{\epsilon}$ and $\mathbf{U}^{0}=\mathbf{U}^{\epsilon}\mid_{\epsilon=0}$, i.e $E_{0}^{\epsilon}:=\frac{1}{256\times 256}\parallel\mathbf{U}^{\epsilon}-\mathbf{U}^{0}\parallel_{2}$. We observe that with decreasing $\epsilon$ the error $E_{0}^{\epsilon}$ decreases to zero, i.e. convergence of the regularisation is confirmed. The error $E^{\epsilon}$ converges as well, but to a positive value. This suggests that error due to spatial discretisation is larger than the regularisation error.

Results of a computational convergence test that computes the least square error between the numerical and exact solutions of equation (35), $E_{0}^{N}=\frac{1}{N^{2}}\parallel\mathbf{U}-\mathbf{U}^{0}\parallel_{2}$, for different numbers of grid cells are given in Table 2. The error is steadily decreasing as the grid is refined, suggesting convergence of the method.

To illustrate spatio-temporal behaviour, the results of the numerical simulations with $\epsilon=0$ at different times are shown in Figure 2. The support of the solution expands radially, in agreement with the property of the porous medium equation, that initial data with compact support lead to solution with compact support (finite speed of interface propagation). The solution attains its maximum in the center of the domain and decreases toward the interface, where it indicates blow-up of the gradient, as expected from (35) for $m>1$. Interface smearing effects that are introduced by the spatial discretisation are rather mild, and the numerical solution is free of spurious oscillations. The numerical errors occur primarily at the interface, whereas agreement between the exact and the numerical solution is very good in the interior of the region where $U>0$ and the solutions are smooth, cf. Figure 3. But also here we see that the errors at the interface decrease as the grid resolution is increased.

Figure 4 shows the time step chosen by the ROW method. Initially it remains constant at 0.003, at approximately $t=0.56$, it drops to 0.0015 where it remains for the remainder of the simulation, up to $t=1$.

These tests show that the spatial discretisation employed here produces results of acceptable accuracy for highly degenerate diffusion reaction equations.

5 Numerical examples and validation

In section 5.1 we show an illustrative simulation of biofilm model (1) with boundary conditions (3) in growth limited and transport limited regimes with the numerical method discussed in section 4. Subsequently, in section 5.2 we investigate the behaviour of the numerical solution with respect to grid step size $\Delta x$, regularization parameter $\epsilon$ and the tolerance required of the Rosenbrock-Wanner method. The model parameters that we use in these simulations are given in Table 3.

5.1 Illustrative Simulation

In a first illustration of the method we consider a case where initially the substratum is inoculated by six semi spherical colonies, which are symmetrically arranged around the center of the substratum on the bottom boundary. The spacing between the three colonies to the left and to the right is equidistant, but the spacing between the two inner most colonies is slightly larger. Substrate is added from the top. The domain $[0,1]\times[0,1]$ is square, uniformly discretised by a mesh with $N\times M=256\times 256$ grid cells. The simulations shown here are carried out with $\epsilon=0$, i.e, the non-regularised, degenerate problem is solved. For the time integration of the spatially discretised system (23), we set the tolerance of the ROW method at $TOL=1e-7$.

Snapshots of the simulation are shown in Figure 5. Biofilm colonies are color-coded by biomass density, the grey-scale isolines are drawn to depict the substrate concentration $C$.

As the graphical results show, the numerical solution is positive and bounded by a value less than one as we have shown mathematically in section 3. Furthermore, the numerical simulation preserves the symmetry of the solution that was imposed by symmetric initial conditions.

Initially the individual colonies grow and expand but are still separated. At $t=3$ the colonies have merged into one larger colony left of center, and one right of center, but with a gap still in the center where the colonies initially were spaced further apart. The colonies continue to grow and stratify. At $t=6$ the gap is closed and the biofilm stratifies. At $t=12$ finally, one almost homogeneous layer of biofilm is formed. Throughout the simulation, the substrate concentration appears stratified, parallel to the substratum. While the substrate concentration decreases toward the substratum due to consumption in the biofilm, it does not become severely limiting.

In Figure 5 we also show how the time step $h(t)$ changes over the progress of the simulation. Initially it increases, indicating that the initial guess for the time step could have been chosen larger. After a short plateau phase it decreases in an oscillatory fashion and attains a plateau phase at time $t=3$ which is occasionally disrupted by intermittent drops in time-step from which the method recovers quickly. This shows not only the adaptability of the numerical method but also indicates that the adaptivity of the method allows a simulation with time steps that usually are much larger than the smallest occurring time step.

It was shown in [35] that the biofilm structure depends on the ratio of biomass growth and nutrient supply. For low nutrient supply, biofilms form patchy structures while for high nutrient availability more compact structures are obtained. We have shown in Figure 5 that for the values of parameters which we used for numerical simulations, biofilm colonies expand and form a homogeneous thick layer. Nutrient supply can be controled by the bulk concentration and by the diffusion parameters, namely the diffusion coefficient $d_{c}$ and the diffusion length $H$. In our non-dimensional formulation, decreasing the bulk concentration is equivalent to increasing the half saturation concentration $K_{U}$. In order to test our method for biofilms under nutrient limitations, we change $K_{U}$ from $0.13$ to $0.65$ and decrease $d_{c}$ from $33$ to $3.3$.

Snapshots of the simulation are shown in Figure 6. Biofilm colonies are color-coded by biomass density. It is observed in Figure 6 that for lower nutrient availability, the biofilm has rough structure and bacterial colonies do not merge together even at $t=350$. They grow upwards, where nutrients are coming from, rather than spreading horizontally. The simulation in Figure 6 also shows the hollow inside the biofilm which emerges at $t=150$ as consequence of nutrient limitation.

5.2 Validation of the Numerical Method

5.2.1 Convergence of the regularization

We have seen in Section 4 for the porous medium equation that solving the regularised problem and passing $\epsilon$ to 0 gives a solution of the degenerate problem. We investigate this here for the full prototype biofilm growth model (1) by solving (25), on a uniform grid of size $256\times 256$. In order to determine the rate of convergence approximately we consider degenerate and regular semi-discrete problems (23) and (25) with a similar initial setting as in Section 5.1.

We define the error between solutions of degenerate and regular semi-discrete problems (23) and (25) as $E^{\epsilon}_{U}=\frac{1}{256\times 256}\parallel\mathbf{U}^{0}-\mathbf{U}^{\epsilon}\parallel_{2}$ and $E^{\epsilon}_{C}=\frac{1}{256\times 256}\parallel\mathbf{C}^{0}-\mathbf{C}^{\epsilon}\parallel_{2}$ and compute them at $t=6$. This is the time at which colonies are merged into a larger colony and there is not any gap between them.

The values of $E^{\epsilon}_{U}$ and $E^{\epsilon}_{C}$ are shown in Figure 7 for different values of $\epsilon$. The logarithmic scale is chosen for both axes to show the results more clear.

As the graphical results in Figure 7 show the error between solutions of regular and degenerate semi-discrete problems is decreasing by decreasing the value of $\epsilon$. For $\mathbf{U}$ the error is initially decreasing, at approximately $\epsilon=10^{-5}$ it becomes constant and does not show any considerable variation for $\epsilon<10^{-5}$ indicating the convergence of $\mathbf{U}^{\epsilon}$ to $\mathbf{U}^{0}$. For $\mathbf{C}$ the error is initially decreasing, at approximately $\epsilon=10^{-7}$ it reaches a plateau and remains constant for $\epsilon<10^{-7}$ confirming the convergence of $\mathbf{C}^{\epsilon}$ to $\mathbf{C}^{0}$. From the results shown in Figure 7 we conclude that $\mathbf{U}^{\epsilon}$ converges faster than $\mathbf{C}^{\epsilon}$.

5.2.2 Grid refinement

For a grid refinement study of the full degenerate biofilm model, we consider the degenerate semi-discrete problem (23) with the same initial setting as Section 5.1 to account for spatial accuracy explicitly. We refine the grids with $2^{\kappa}$, where $\kappa$ is an integer and compute the least square norm of the difference between two subsequent solutions. The least square errors are defined as $E^{\kappa}_{U}=\frac{\parallel\mathbf{U}_{\kappa}-\mathbf{U}_{\kappa-1}\parallel_{2}}{2^{2(\kappa-1)}}$ and $E^{\kappa}_{C}=\frac{\parallel\mathbf{C}_{\kappa}-\mathbf{C}_{\kappa-1}\parallel_{2}}{2^{2(\kappa-1)}}$.

We compute the errors at two different times $t=2$ and $t=6$. At $t=2$ individual colonies have not merged into one larger colony yet and at $t=6$ the gap between colonies is closed and the biofilm stratifies. We also compute the Elapsed CPU Time (ECPUT) for each cell resolution.

Tables 4-5 show the least square error between two subsequent solutions at times $t=2$ and $t=6$ and also the elapsed CPU time. We observe a steady decrease in errors $E^{\kappa}_{U}$ and $E^{\kappa}_{C}$ at both times as the grid is refined, indicating the convergence of the method and its reliability to give tolerable smearing around the interface.

Since we have the density value at the center grids, to compute these data interpolation must be used which is a source of error so these data cannot be used to determine the convergence rate accurately.

5.2.3 Influence of grid refinement on interface location

In order to study in more detail the movement of the interface and its location for different choices of $\Delta x$ we change the initial setting. One colony is placed at the middle of the substratum. The domain $[0,1]\times[0,1]$ is square and the tolerance of the ROW method is set at $TOL=1e-7$.

To show the location of biofilm/liquid interface we pick the interface point located on the symmetry line of the initial colony, i.e. the line parallel to the $y$-axis with $x=0.5$. The location of the interface between biofilm and liquid in the $x$-$t$ plane is plotted in Figure 8 for various $\Delta x$. Its movement is described by a piecewise constant increasing function, which represents the discrete nature of the grid and that the interface location can be determined at most with a $\Delta x$ accuracy. As the grid is refined, the discrete jumps diminish and convergence of the interface location is observed as $\Delta x\rightarrow 0$.

5.2.4 Effect of method parameters on the time step

We have shown in Sections 5.2.2 that the accuracy of spatial discretisation increases if grids are refined. We have also shown in section 5.2.1 how fast the solution of the regularised semi-discrete problem is convergent to the solution of the degenerate semi-discrete problem. A third parameter, besides $\Delta x$ and $\epsilon$, that determines the accuracy of ROW method and controls the time-step size $h(t)$ is the tolerance of the ROW method. Here we study the effect of these three parameters on the time-step size. We consider the same initial setting as in Section 5.2.3.

The effect of cell resolution on the value of time-step $h(t)$ for the degenerate semi-discrete problem (23) is shown in Figure 9.A. For each cell resolution, initially the time-step size increases, indicating that initial guess could have been chosen larger. After a short plateau phase, $h(t)$ decreases in an oscillatory fashion. The time-step oscillates during the simulation for all spatial-step size except for $\Delta x=\frac{1}{256}$ in which the oscillation stops at approximately $t=2.96$ and $h(t)$ remains constant for the remainder of the simulation up to $t=6$. Figure 9.A also shows that finer grids require finer time steps.

The size of the time-step for the regular semi-discrete problem (25) with different values of $\epsilon$ including $\epsilon=0$ is plotted in Figure 9.B, for $M\times N=256\times 256$ and $TOL=1e-7$. The time step $h(t)$ is not strongly sensitive with respect to the value of $\epsilon$.

The time-step size $h(t)$ for the degenerate problem (23) for different values of tolerance of the ROW method is shown in figure 9.C. The time step is sensitive with respect to this parameter, and decreasing the tolerance, i.e. increasing the desired accuracy, in the ROW method decreases the value of time-step size. For all three values of tolerance, initially $h(t)$ increases but the deviation from the initial guess is considerable for $TOL=1e-6$. After the first short plateau phase, $h(t)$ decreases in oscillatory fashion for all case studies. For $TOL=1e-6$ this oscillation can be seen during the simulation up to $t=6$. For $TOL=1e-7$ oscillatory behaviour stops at approximately $t=2.96$ while this time for $TOL=1e-8$ is approximately $t=1.96$. We conclude that for smaller values of tolerance, $h(t)$ shows less oscillation and becomes stable faster.

5.3 Discussion of the Numerical Method vis-a-vis the Semi-Implicit Method with Fixed Time Steps

The spatial discretisation described in Section 2.2 is the same that underlies the semi-implict method of [10, 43] and related schemes, which have been predominantly used in applications. The resulting semi-discrete equation is a large, highly nonlinear and singular ordinary differential equation system. The difference between the methods lies in the approach to solving this initial value problem. We showed in Section 3, via a regularisation analysis, that the solutions of this system indeed possess sufficient smoothness to employ higher order methods, such as embedded Rosenbrock-Wanner methods (which we used; other choices would have been possible), despite the low regularity of the biologically relevant solutions of the underlying partial differential equation, which display gradient blow-up along an interface moving with finite speed. Compared to the semi-implicit method, which requires in each time step the solution of 2 sparse linear systems of size $(NM)^{2}$, the method that we used here is more expensive per time step, as it it requires the solution of 4 sparse linear systems of size $4(NM)^{2}$, and the evaluation of the Jacobian. On the other hand, it has better convergence properties, allowing larger time steps to achieve the same accuracy. In applications, the semi-implicit method and is relatives that use fix time-steps require sometimes extensive a priori explorations by the user to find for a particular simulation experiment a suitable time step that gives an acceptable effort/accuracy trade-off. In the method presented here, on the other hand, accuracy is guaranteed by an automatic, adaptive time-step choice.

6 Quorum sensing induction in isolated colonies

6.1 Mathematical model

The mathematical model that we propose to study quorum sensing induction in the growing biofilm is an extension of (1), that accounts for the autoinducer, or signal molecule concentration $s$. Following [14] for autoinducer production, the modified model reads

[TABLE]

where $D(u)$ is defined in (2) and $d_{c}$ and $d_{s}$ are the constant diffusion coefficients for nutrient and signal molecules respectively. In equation (38) autoinducers are produced at a base rate $\alpha$ if the local autoinducer concentrations is small relative to induction threshold, which is scaled here to be 1. If the autoinducer concentration exceeds the induction threshold, the production rate will increase to $\alpha+\beta$. The transition is described by a Hill function with exponent $m$. Quorum sensing parameters $\alpha,~{}\beta,~{}\psi$ and $m$ are chosen in the range of data in [16].

For our simulations we consider a rectangular computational domain $[0,L]\times[0,H]$. The boundary condition for biomass $u$ is the homogeneous Neumann condition everywhere. Note that, as long as the biofilm interface does not reach the domain, the biomass density satisfies both the homogeneous Dirichlet and Neumann conditions. For the concentrations $c$ and $s$ we pose the homogeneous Neumann condition at the lateral sides of the domain and at the substratum, and Robin condition at the top, where nutrients are added and autoinducer signals diffuse out of the domain.

Thus, the specific boundary condition for problem (38) is defined as:

[TABLE]

where $\lambda$ is the concentration boundary layer thickness.

Initially the system is inoculated by one semi-spherical colony placed at the center of the substratum on the bottom boundary. The simulation is terminated when the signal concentration everywhere in the domain exceed switching threshold, i.e. the first time at which $s(t,x,y)\geq 1$ for all $(x,y)\in\Omega$.

We investigate how the domain size affects the time to induction. We keep the height of the domain constant at $H=1$ and change the length of the computational domain in the range $1\leq L\leq 2.5$. For each $L$ we measure the following lumped output parameters:

• •

$T_{1}$: The first time at which $s(t,x,y)\geq 1$ for some $x\in\Omega$

• •

$T_{2}$: The first time at which the average signal concentration in the colony exceeds switching threshold, i.e. when $\int_{\Omega_{2}}s(t,x,y)dxdy\geq\int_{\Omega_{2}}dxdy$

• •

$T_{3}$: The first time at which $s(t,x,y)\geq 1$ for some $(x,y)\in\Omega_{2}(t)$

• •

$T_{4}$: The first time at which $s(t,x,y)\geq 1$ for all $(x,y)\in\Omega$

We also compute the total biomass in the system and its value at each time $T_{i},~{}i=1,...,4$ to describe the relation between biofilm density and quorum sensing induction. The value of total biomass and its value at each time $T_{i},~{}i=1,...,4$ are defined as

[TABLE]

The parameter values used in the simulations and their definition are shown in Table 6. The internal tolerance is set in all simulations to $TOL=1e-7$, spatial discretization is such that $\Delta x=1/256$.

6.2 Results

Figure 10 shows the autoinducer concentration $s$ and the interface between biofilm and surrounding liquid at $T_{1}$ for four different system sizes with $L=1,1.5,2,2.5$. In all cases the biofilm colony remains approximately semi-spherical, but we observe that the larger $L$, the larger is the biofilm colony at the time that the signal concentration reaches 1 first (see also Figure 11.C). This is a consequence of increased substrate supply as the domain length is increased, and hence the length of the boundary through which nutrient is added. The maximum signal concentration is attained in the center of the colony; it decreases from there toward the lateral boundaries of the domain and toward the top boundary. The concentration at the top boundary is lower, as this is where autoinducers are removed by diffusion due to the boundary conditions (39). The larger the system length $L$, the steeper are the autoinducer gradients and the lower is the autoinducer concentration at the boundaries. This is a consequence of signals only being produced in the colony. For the smallest system length $L=1$ the autoinducer concentration appears stratified across the domain, whereas in the the highest systems length $L=2.5$ it appears to be much larger in the colony than in the aqueous phase.

The value of parameters $T_{i},~{}i=1,...,4$ and also the value $M_{total,i}$ of total biomass at each $T_{i}$ for computational domains with length $L=1,1.25,1.5,1.75,2,2.5$ are shown in Figure 11. We observe in Figure 11.A that increasing the length of the computational domain increases the value of $T_{i}$’s at which QS induction occurs. In fact, by increasing $L$ more autoinducers transport away from the producer. Hence, more biomass is needed to have induction, therefore the time at which QS induction occurs increases. We also observe that times $T_{i},~{}i=1,2,3$ are very similar but there is a considerable difference between them and $T_{4}$ especially in the larger domains.

The values of total biomass at each time $T_{i},i=1,...,4$ are shown in Figure 11.B. For each time $T_{i},~{}i=1,...,4$, increasing the length of the domain increases the total amount of biomass. Increasing the domain length increases mass transfer thus more biomass is needed for differentiation. Mirroring the differences between $T_{1,2,3}$ and $T_{4}$, Figure 11.B shows considerable difference between total biomass at $T_{i},~{}i=1,2,3$ and its value at $T_{4}$. The time course of total biomass $M(t)$ for $L=1,1.5,2,2.5$ is plotted in Figure 11.C to show the effect of the length of the domain on the required biomass for QS induction vs time. We observe that initially the values of total biomass in all four systems are the same and deviation emerges at $T\simeq 2.2$. We also observe that $M(t)$ in the system with $L=2.5$ has the maximum value indicating the maximum mass transfer occurs in the system with largest domain. Furthermore, it can be seen that the value of $T_{4}$ which is in fact the time at which simulation stops is increasing by increasing the length of computational domain as we have shown in Figure 11.A.

The simulation results in Figure 11.A,B, suggest that $T_{i}$’s and $M_{total,i}$ change with respect to the system length quadratically. This is confirmed by fitting a quadratic function to our data points, $T_{i}(L)=\alpha_{1}L^{2}+\beta_{1}L+\gamma_{1}$ and $M_{\small{total},i}(L)=\alpha_{2}L^{2}+\beta_{2}L+\gamma_{2}$, which are also included in the plots. The coefficients $\alpha_{1},~{}\beta_{1},~{}\gamma_{1}$ and $\alpha_{2},~{}\beta_{2},~{}\gamma_{2}$ are given in Table 7. Also given is the root mean squared error defined as $RMSError_{T}=\sqrt{\frac{\sum_{i=1}^{n}(\hat{T}_{i}-T_{i})^{2}}{n}}$ for $T_{i}$ and $RMSError_{M}=\sqrt{\frac{\sum_{i=1}^{n}(\hat{M}_{i}-M_{i})^{2}}{n}}$ for $M_{total,i}$ in which $n$ is the number of collected data and $\hat{T}_{i}$ and $\hat{M}_{i}$ are the predicted values of $T_{i}$ and $M_{total,i}$. We observe that $RMSError_{M}<RMSError_{T}$ indicating the error between the predicted value of $M_{\small{total},i}$ and obtained data is less than the error between the predicted and real values of $T_{i}$.

6.3 Discussion of the results of QS induction in isolated colonies

In agreement with earlier studies [4, 16, 52] our results show that upregulation is rapid. In our simulations, the time between the signal concentration first reaching the switching threshold somewhere in the colony and until the switching threshold is exceeded in the colony everywhere is very short, about 1% of the induction time.

Our simulations suggest that the environment and the physical conditions there affect the time to induction for a single isolated colony. This is to a large extent determined by mass transfer of signals from the colony into the surrounding aqueous phase, and transport of signal molecules there. In our simulations we have neglected abiotic decay of signal molecules in the aqueous phase, but it is to conjecture that including such an effect would exacerbate this phenomenon. This suggests that upregulation of an individual colony is not a mere function of colony size, in accordance with [19]. A related, but not directly comparable study in this regard is [4], where a one-dimensional model was used to determine the minimum thickness required for a homogeneous biofilm layer to induce. In that study the aqueous phase was not explicitly considered, but the diffusive flux of signal molecules from the colony in direction of the bulk phase is specified as a boundary condition. This, in our setting would correspond to changing the depth of the cavity, which we kept constant. In contrast to this, we varied the lateral extension of the domain, showing that also horizontal mass transfer, orthogonal to the main diffusion direction matters. The quadratic dependence of induction time and biomass on the system size also suggest that diffusion is a dominant factor in the process.

On the other hand, if we understand the lateral homogeneous Neumann conditions as symmetry conditions instead of as no-flux conditions, then we can interpret our system also as an infinite array of equidistantly spaced identical colonies. Our system then models a biofilm community consisting of many individual colonies. System length $L$ is then the distance between the centers of 2 consecutive colonies, i.e. the smaller $L$, the more biomass in the biofilm. From the viewpoint of the entire biofilm community our results confirm the traditional interpretation of quorum sensing, i.e. induction takes place when the population reaches a certain size. We can reconcile this with the findings discussed above by the observation that the horizontal signal fluxes between two neighbouring colonies cancel each other, i.e. the net flux between neighbouring colonies is zero. Again this confirms the important role of signal diffusion.

7 Conclusion

Our primary objective was to study a time-adaptive, error-controlled numerical solution strategy for a non-linear degenerate diffusion system that arises in biofilm modelling. This system has two non-linear effects (i) porous medium degeneracy for $u=0$ and (ii) super-diffusion singularity for $u=1$, which make this problem difficult to treat numerically and analytically.

Biologically relevant solutions of this system are characterized by an interface that propagates with finite speed, along which the biomass gradient blows up. Due to the low reglarity of these solutions, low order time intergation methods have been used previously for the time integration of the semi-discrete system that is obtained after spatial discretisation. These methods mostly have been explicit or semi-implicit, with small enough fixed time step to avoid numerical difficulties (e.g. overshooting of the singularity), and without error control capabilites. The important question that we answered is whether we can also use a higher order time integration method to solve the semi-discrete approximation of this highly non-linear PDE. By regularisation we showed that the spatially discretised problem indeed satisfies a Lipschitz condition, i.e. has classical smooth solutions, which never reach the singular point. This allows the application of time-adaptive, error controlled time integration techniques, such as embedded Rosenbrock-Wanner methods, to simulate these highly nonlinear problems with singular and degenrate diffusion.

To demonstrate the utility of the numerical method we carried out a simulation experiment of quorum sensing induction of a biofilm colony in a protected niche. We found that induction time does not depend only on colony size, but also on the properties of the environment that affect transport of signal molecules, e.g niche size in our case.

Acknowledgements This study has been financially supported by the Natural Science and Engineering Research Council of Canada (NSERC): MG holds a PGS-D graduate scholarship, HJE a Discovery Grant, the equipment was purchased with a Research Tools and Infrastructure Grant (HJE).