1. Introduction

Biofilms are dense accumulations of microbial cells on biotic or abiotic surfaces (called substrata) in aqueous environments. Once the microbial cells become sessile, they produce extracellular polymeric substances (EPS) that protect them against antibiotic attacks and mechanical washout [29]. Due to the sorption properties and enhanced mechanical stability of biofilms, they are beneficially used in wastewater treatment, soil remediation and groundwater protection [43]. On the other hand, biofilm formation and detachment can have very disadvantageous effects and lead to serious infections in the human body, biocorrosion of drinking water pipes or industrial facilities [11], contamination in food processing plants [27,30], etc. Biofilm formation is characterized by the balance of attachment, growth and detachment or dispersal processes [23,41]. Among these phenomena, there is a growing interest in the study of the latter, i.e. the release of microbial cells from the biofilm into the aqueous environment. Commonly, by detachment one refers to cell losses into the aqueous environment that are caused by external forces. Usually these external forces are shear forces due to bulk flow hydrodynamics [32]. These detachment losses can manifest themselves as sloughing or erosion of cells from the outer layers of the biofilm. By dispersal we refer to cell losses that can be internally triggered, e.g. by enzyme-mediated breakdown of the biofilm matrix [6], production of surfactants which loosen cells from the biofilm [8]; or externally triggered, e.g. by changes in nutrient availability [23], production of free radicals [4], or controlled by quorum sensing systems [34,37,45]. Dispersed cells can originate from inner layers of the biofilm and can contribute to downstream colonization, and thus eventually result in pipe obstruction, bacterial infection (biomedical implants), or increased microbial contamination in food processing plants [40].

Numerous mathematical models of biofilms have been proposed in the literature, focusing on different time and length scales and processes, and utilizing different mathematical concepts, from agent based to continuum mechanistic and from stochastic to deterministic. All dynamic models of biofilm formation include growth processes, including their dependence on nutrients. On the other hand, detachment and dispersal processes are often neglected, or treated in a rather qualitative manner, e.g. by coupling the detachment rate to biofilm thickness [12,43,44] or geometrical properties of the biofilm structure [46], without accounting for the processes that induce cell loss. No biofilm model is known that includes several or all known detachment and dispersal mechanisms. Rather each modeling study focuses on one particular trigger. A number of biofilm growth models have been proposed that include a physics based description of sloughing or erosion due to external forces, cf [1] for a simple one-dimensional model set in the Wanner-Gujer framework, and [31] for a 2D cellular automaton model. Only few papers have been published that use models of internally triggered biofilm dispersal: A cellular automaton model for nutrient limited dispersal was presented in [22,23], a cellular automaton model for detachment caused by enzymatic breakdown of the EPS matrix was presented in [47]. A partial differential equation based model for cell dispersal triggered by quorum sensing was presented in [18] and investigated there in first computer simulations. The question of well-posedness of the solutions of this model remained open. To give an answer to this question will be the focus of our current paper.

The model in [18] is an extension of a prototype biofilm growth model, which was originally introduced in [13] and did not include any detachment or dispersal processes. The biofilm is characterized in terms of the volume fraction that cells and EPS locally occupy. This is described by a highly non-linear diffusion-reaction equation for biomass, with two non-Fickian effects: (a) the diffusion operator degenerates like the porous medium equation for vanishing biomass densities and (b) it blows up if the local cell density approaches its maximum value. These effects ensure that the biofilm/water interface spreads at a finite speed and that the maximum biomass density is never exceeded. The biofilm expands spatially if the local cell density fills up the available volume while it does not spread notably if there is space locally available to accommodate new cells. In the prototype biofilm growth model, this biomass model is coupled with an additional non-degenerate diffusion-reaction equation for the nutrient that limits biomass growth. In [18] this model was extended to account also for the concentration of the quorum sensing molecules that trigger dispersal, and for an equation that describes the dispersed cells.

Previous extensions of this prototype biofilm growth model that included quorum sensing effects focused on the role of bulk hydrodynamics to facilitate non-local up-regulation due to advective transport [21], and on quorum sensing controlled EPS production [20] as a mechanism to switch from a colonization mode, in which resources are primarily invested in proliferation, to a protected mode of growth, in which EPS is produced, e.g. to mechanically stabilize the biofilm. In these models, down-and up-regulated cells are treated as two different cell fractions. For the model in [21] existence and uniqueness were derived in [38], whereby special features of the model assumptions could be used that do not hold in other, seemingly similar applications of the same biofilm modeling framework. In the current model, up-regulation is implicit. Rather than distinguishing between down-and up-regulated cells, we distinguish between cells that are sessile in the biofilm and motile ones, which after up-regulation disperse from the colony. The structure of the quorum sensing induced dispersal model is different from the previous multi-component biofilm models [10,24,38,40], in which one biomass type is described by a degenerate diffusion-reaction equation, and the other one by a semi-linear one.

To prove existence and uniqueness of solutions and continuous dependence of solutions on initial data we will use here ideas applied in [17] for the mono-species model and in [3] for a scalar degenerate reaction-diffusion equation of porous-medium type.

---

2. Mathematical model

We analyze the mathematical model of quorum sensing induced detachment in biofilms which was proposed in [18]. In this model we consider particulate biomass, i.e. biofilm bacteria, and suspended biomass, i.e. dispersed cells. As is common in many biofilm modeling studies, cf [44], we do not explicitly track EPS but subsume them in the biomass. The growth of both biomass fractions depends on a growth limiting nutrient. The bacterial cells have the ability to produce quorum sensing signal molecules which can trigger a switch from the sessile to the suspended mode of growth. The model is formulated for the dependent variables local density of particulate biomass (biofilm), $\tilde M$, the concentration of dispersed cells, $\tilde N$, the concentration of a nutrient, $\tilde C$, and the quorum sensing molecule or autoinducer concentration $\tilde A$. The independent variables are time $\tilde t$ and the spatial coordinate $\tilde x \in \Omega \subset \mathbb{R}^n, n=2,3$; $\Omega$ is here the model domain that will be made more precise later. The model in dimensional form reads

$\begin{split} \partial_{\tilde t}\tilde M &= \nabla_{\tilde x}\cdot(\tilde D_M(M)\nabla_{\tilde x} \tilde M) + \tilde \mu \frac{\tilde C}{\tilde k_1+\tilde C}\tilde M - \tilde k_2 \tilde M - \tilde \eta_1 \left(\frac{\tilde A^m}{\tilde \tau^m + \tilde A^m}\right)\tilde M, \\ \partial_{\tilde t} \tilde N &= \tilde d_1\Delta_{\tilde x} \tilde N+ \tilde \mu \frac{\tilde C}{\tilde k_1+ \tilde C}\tilde N - \tilde k_2 \tilde N + \tilde \eta_1 \left(\frac{\tilde A^m}{\tilde \tau^m + \tilde A^m}\right)\tilde M, \\ \partial_{\tilde t}\tilde C &= \tilde d_2\Delta_{\tilde x} \tilde C - \tilde \sigma \frac{\tilde C}{\tilde k_1+\tilde C}\left(\tilde M+ \tilde N\right),\\ \partial_{\tilde t}\tilde A &= \tilde d_3\Delta_{\tilde x} \tilde A - \tilde \lambda \tilde A + \left[\tilde \alpha + \tilde \beta\frac{\tilde A^m}{\tilde \tau^m+\tilde A^m}\right]\left(\tilde M+ \tilde N\right). \\ \end{split}$

(1)

In model (1) we have intentionally omitted the re-attachment terms which were originally included in the model [18], because the simulations in [18] show that only a negligible amount of dispersed cells gets re-attached. Hence the inclusion or exclusion of the re-attachment terms will not make a significant difference from an application point of view, but it simplifies the mathematical analysis.

In the first equation of model (1), the diffusion term describes spatial expansion of the biofilm. This is a volume filling problem: as long as locally space is available to accommodate new cells the biofilm expands very slowly or not at all. When the maximum biomass density $\tilde M_\infty$ that can be locally accommodated is nearly reached the biofilm starts to expand quickly. Following [13], this is described here by nonlinear, density dependent diffusion, with coefficient

$\tilde D_M(M) = \tilde d\frac{M^a}{(1-M)^b} , \mbox{where} \;\;\;a>1, b>1, \tilde d > 0.$

(2)

Here $M:=\tilde M/\tilde M_\infty$ is the volume fraction occupied by the biomass. This implies that we must have $0\leq M<1$. The diffusion coefficient $\tilde D(M)$ vanishes when the biomass density approaches zero and blows up when it tends to its maximum value. Note that the prototype biofilm growth model of [13] is a special case for this model in the absence of quorum sensing and/or dispersal effects, i.e. if $\tilde \alpha=\tilde \beta=0$ (and initially $\tilde A \equiv 0 \equiv \tilde N$) or if $\tilde \eta_1=0$. The spatial expansion of the biofilm is thus driven by biomass accumulation. The polynomial degeneracy $M^a$, well known from the porous medium equation, guarantees that spatial spreading is negligible for low values of $M$ and yields the separation of biofilm and liquid phase. For $0\ll M \approx 1$, the equation shows a super-diffusion effect, since the diffusion coefficient possesses a singularity at $M=1$ and blows up. In the prototype biofilm growth model of [13], this ensures the maximal bound for the biomass density which is a physical limitation as the number of cells that fits into a unit volume is bounded [16]. In particular, the blow up of the diffusion coefficient guarantees that $M<1$ if homogeneous Dirichlet conditions are specified for $M$ on some parts of the boundary of the domain [16]. Since the production of biomass depends on the availability of nutrients, the upper bound on $M$ cannot be guaranteed by the growth terms alone. The degeneracy $M^a$ alone does not yield this maximum bound for the cell density, while the singularity $(1-M)^{-b}$ does not guarantee the separation of biofilm and liquid region by a sharp interface.

In model (1), diffusion of $\tilde N, \tilde C, \tilde A$ is Fickian with constant coefficients $\tilde d_{1,2,3}$. This is a simplification of [18], where it was assumed that the diffusion coefficients are smaller in the biofilm than in the aqueous phase, but well within one order of magnitude.

The reaction terms in (1) describe the following processes that are the same as in the model of [18]:

• Growth of sessile and dispersed cells $\tilde M$ and $\tilde N$ is controlled by the local availability of nutrients. This is described by standard Monod kinetics, cf [44], in the first and second equation of (1), where $\tilde k_1$ is the half saturation concentration and $\tilde \mu_1$ is the maximum specific growth rate. The growth kinetics is assumed to be the same for both sessile and dispersed biomass.

• Natural cell death occurs at the same rate $\tilde k_2$ for $\tilde M$ and $\tilde N$.

• Dispersal, i.e. the transition of bacteria from the sessile state $\tilde M$ into the suspended motile state $\tilde N$, is controlled by the local concentration of the quorum sensing molecule $\tilde A$. For small values of $\tilde A \ll \tilde \tau$ this transition is nearly zero (i.e. the biofilm is primarily down-regulated), but increases for large concentration values $\tilde A \gg \tilde \tau$, and levels off at a maximum rate $\tilde \eta_1$ if $\tilde A$ is large (i.e. the biofilm is primarily up-regulated). This transition between both states is described by a Hill function, where $\tilde \tau$ is the induction threshold and the exponent $m>1$.

• Consumption of nutrients is proportional to biomass growth, i.e. described by the same Monod kinetics. The maximum consumption rate is defined as $\tilde \sigma_1=\tilde \mu_1/Y$ where $Y$ is the yield coefficient.

• The signal molecules $\tilde A$ in the fourth equation of (1) are produced at a base rate $\tilde \alpha$ if the local signal concentration is small, $\tilde A \ll \tilde \tau$, and at the increased rate $\tilde \alpha + \tilde \beta$ if it is large, $\tilde A \gg \tilde \tau$. The transition between these two regimes is modeled by the same Hill function that was used to describe the triggering of dispersal, cf also [18]. According to [19] we assume $\tilde \beta =10\tilde \alpha$. Furthermore we have included an abiotic decay term, at rate $\tilde \lambda$, for the signaling molecule, which was neglected in [18], but is included in many other quorum sensing models, e.g. [20,26,38,42].

To non-dimensionalize the model we express the biofilm biomass density $\tilde M$ by the volume fraction $M$ as described above and further introduce the new dependent variables

$N:=\frac{\tilde N}{\tilde M_\infty}, C:=\frac{\tilde C}{\tilde C_\infty}, , A:=\frac{\tilde A}{\tilde \tau},$

where $\tilde C_\infty$ is the bulk substrate concentration that enters the model via boundary conditions (see below). The dimensionless independent variables are introduced as

$t := \tilde \mu \tilde t, x :=\frac{\tilde x}{\tilde L},$

where $\tilde L$ is a characteristic length scale of the domain, e.g. its length in case of rectangular $\Omega$. In dimensionless form, the model becomes then

$\begin{split} \partial_{t}M &= \nabla\cdot(D_M(M)\nabla M) + \frac{C}{k_1+C}M - k_2 M - \eta_1 \ \left(\frac{A^m}{1 + A^m}\right)M, \\ \partial_{t}N &= d_1\Delta N+ \frac{C}{k_1+C}N - k_2 N + \eta_1 \ \left(\frac{A^m}{1 + A^m}\right)M, \\ \partial_{t}C &= d_2\Delta C - \frac{\sigma \ C}{k_1+C}(M+N),\\ \partial_{t}A &= d_3\Delta A - \lambda A + \left[\alpha + \beta\frac{A^m}{1+A^m}\right](M+N). \\ \end{split}$

(3)

Here the spatial derivative operators $\nabla$ and $\Delta$ are now with respect to the dimensionless independent variable $x$. The parameters in (3) are $d_i :=\frac{\tilde d_i}{\tilde \mu \tilde L^2}, i=1,2,3$, $d :=\frac{\tilde d}{\tilde \mu \tilde L^2}$, $D_M(M):=\frac{d}{\tilde d} \tilde D(\tilde M/\tilde M_\infty)$, $k_1:=\frac{\tilde k_1}{\tilde C_{\infty}}$, $\alpha=\frac{\tilde \alpha \tilde M_{\infty}}{\tilde \tau \tilde \mu}$, $\beta=\frac{\tilde \beta \tilde M_{\infty}}{\tilde \tau \tilde \mu}$, $\eta_1=\frac{\tilde \eta_1}{\tilde \mu}$, ${\sigma} = \frac{\tilde \sigma \tilde M_{\infty}}{\tilde \mu \tilde \tau}$ and $\lambda := \frac{\tilde \lambda}{\tilde \mu}$. The full list of the parameters and their values in dimensionless form that we use in our simulations later on is found in table 1.

In this model, the actual biofilm is the region where sessile biomass is present, $\Omega_2(t):=\{x \in \Omega : M(t,x)>0\}$, whereas the surrounding aqueous phase corresponds to the region where sessile biomass is absent, $\Omega_1(t):=\{x \in \Omega: M(t,x)= 0\}$, cf. Figure 1. Since $M$ changes with time due to growth and dispersal, also these two regions change. They are separated by the biofilm/water interface, $\Gamma(t) := \partial \Omega_2(t) \setminus \partial \Omega$. Neither $\Omega_1(t)$ nor $\Omega_2(t)$ need to be connected domains. In fact, $\Omega_2(t)$ will in general consist of several colonies that are separated from each other by water-filled channels and voids. If a biofilm is contained in the inner region of $\Omega$, away from its boundary, it is often called microbial floc (biofilm without substratum). This situation plays a major role in biological wastewater treatment.

It remains to specify initial and boundary values for the biomass fraction $M$ and the concentrations of detached cells $N$, nutrient $C$ and quorum sensing molecule $A$ to complement the model. This will be made precise in the following section.

---

3. Analysis of the quorum sensing dispersal model

---

3.1. Preliminaries

For technical reasons we study the model in the auxiliary form

$\begin{split} \partial_{t}M &= \nabla\cdot(D_M(M)\nabla M) + \underbrace{\frac{C}{k_1+C}M - k_2 M - \eta_1 \ \left(\frac{|A|^m}{1 + |A|^m}\right)M} _{=g(M,N,C,A)},\\ \partial_{t}N &= d_1\Delta N + \underbrace{\frac{C}{k_1+C}N - k_2 N + \eta_1 \ \left(\frac{|A|^m}{1 + |A|^m}\right)M }_{=f_1(M,N,C,A)},\\ \partial_{t}C &= d_2\Delta C - \underbrace{\frac{\sigma \ C}{k_1+C}(M+N)}_{=f_2(M,N,C,A)},\\ \partial_{t}A &= d_3\Delta A - \underbrace{\lambda |A| + \left[\alpha + \beta\frac{|A|^m}{1+|A|^m}\right](M+N)}_{=f_3(M,N,C,A)}, \end{split}$

(4)

with the initial and boundary conditions

$M|_{\partial \Omega} = 0,\;\;\;\; N|_{\partial \Omega}=0, \;\;\;\; C|_{\partial \Omega}=C_{\infty}, \;\;\;\;A|_{\partial \Omega} = 0, \\ M|_{t=0} = M_0,\;\;\;\; N|_{t=0}=N_{0},\;\;\;\; C|_{t=0}=C_{0}, \;\;\;\; A|_{t=0}=A_{0},$

(5)

where $C_{\infty}=1$, and the functions $M_0,N_0,C_0,A_0:\Omega\rightarrow\mathbb{R}$ are non-negative and bounded. Moreover, we assume that

$\Vert M_0\Vert_{L^{\infty}(\Omega)} < 1-\rho,$

(6)

for some $\rho \in (0,1)$, and $M_0,N_0,C_0$ and $A_0$ satisfy the compatibility conditions.

We point out that non-negative solutions of (3) solve (4) and vice versa, i.e. after the non-negativity is shown, the absolute value $|.|$ can be removed from (4) to obtain the original model (3). Constant Dirichlet boundary conditions are imposed on the nutrient concentration $C$ reflecting a constant unlimited nutrient supply at the boundary of the considered domain, while homogeneous Dirichlet boundary conditions are assumed for $N$ and $A$ enforcing the removal of dispersed cells and quorum sensing signal molecules. Similarly, homogeneous Dirichlet boundary conditions are assumed for the biomass fraction $M$. This describes the situation of a growing biofilm in the interior of the considered domain, away from the boundary. These specific boundary conditions are primarily chosen for convenience regarding the analysis, albeit boundary conditions of mixed type are often considered more appropriate in applications. Typically, Dirichlet boundary conditions are prescribed on some parts of the boundary while Neumann or Robin boundary conditions are specified on the other parts. In particular, the substratum on which the biofilm grows is impermeable for all dependent variables, which can be modelled by homogenous Neumann boundary conditions. An extension of our results for the Dirichlet problem to different boundary conditions can be achieved with the same strategy that was used in [17] for the prototype single-species/single-substrate biofilm model.

Here and in the sequel, we use the following notations, $Q_T := (0,T]\times\Omega$ for some $T>0$, $Q= \mathbb{R}^+\times\Omega$ and

$\Phi(M) := \int_0^M D_M(s)ds= \int_0^Md\frac{s^a}{(1-s)^b}ds \;\;\;\text{for} \ 0\leq M<1 .$

(7)

Definition 3.1. We call $(M,N,C,A)$ a solution of system (4) on $[0,T]$ with the initial and boundary data (5), if the functions

$M,N,C,A \in C([0,T];L^1(\Omega)) \cap L^{\infty}(Q_T)$

and satisfy (4) in distributional sense.

More precisely, if $M$ is a solution of system (4), then

$\int_\Omega (M\varphi)\big|_{s=t} -\int_\Omega M_0\varphi\big|_{s=0}- \int_{Q_t}(M\partial_s\varphi+\Phi(M)\Delta\varphi) =\int_{Q_t} g\left(M,N,C,A\right)\varphi,$

(8)

for all $t\in[0,T]$ and $\varphi\in C^{2}(\overline{Q_T})$ such that $\varphi\geq 0$ in $Q_T$ and $\varphi|_{\partial\Omega}=0.$ Similarly, for $C$ we have

$\begin{split} &\int_\Omega(C\varphi)\big|_{s=t} -\int_\Omega C_0\varphi\big|_{s=0}- \int_{Q_t}(C\partial_s\varphi+C\Delta\varphi) +\int_0^t\int_{\partial\Omega}C_\infty \partial_\nu \varphi\\ =&\int_{Q_t} f_2\left(M,N,C,A\right)\varphi, \end{split}$

(9)

for all $t\in[0,T]$ and $\varphi\in C^2(\overline{Q_T})$ such that $\varphi\geq 0$ in $Q_T$ and $\varphi|_{\partial\Omega}=0$. As usual, $\partial_\nu$ denotes the outward unit normal derivative at the boundary. The identities for $N$ and $A$ are accordingly.

---

3.2. Existence for smooth initial data

We consider smooth non-degenerate approximations for system (4) and show that their solutions converge to the solution of the degenerate problem (4). The ideas are based on the proof developed for scalar degenerate reaction-diffusion equations of porous medium type in [2], the solution theory in [17] for the single species biofilm model and the ideas applied in [10,24,38] and [39] for multi-species biofilm models.

For small $\varepsilon >0$, we define

$\begin{split} D_{\varepsilon}(M) &:= \left\{ \begin{array}{ll} d\varepsilon^a&\mbox{if }\;\; M < 0 , \\ d\frac{(M+\varepsilon)^a}{(1-M)^b}&\mbox{if }\;\; 0\leq M\leq 1- \varepsilon ,\\ d\frac{1}{\varepsilon^b}&\mbox{if}\;\; M\geq 1-\varepsilon, \end{array} \right.\\ \Phi_{\varepsilon}(M) &:= \int_0^M D_{\varepsilon}(s)ds, \end{split}$

(10)

and denote the solutions of the regular auxiliary systems

$\begin{split} \partial_{t}M^{\varepsilon} &=\nabla\cdot(D_{\varepsilon}(M^{\varepsilon})\nabla M^{\varepsilon}) + g(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon}),\\ \partial_{t}N^{\varepsilon} &= d_1 \Delta N^{\varepsilon} + f_1(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon}),\\ \partial_{t}C^{\varepsilon} &= d_2 \Delta C^{\varepsilon} + f_2(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon}),\\ \partial_{t}A^{\varepsilon} &= d_3 \Delta A^{\varepsilon} + f_3(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon}). \end{split}$

(11)

where $D_M$ is replaced by $D_\varepsilon$, by $(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon})$.

Lemma 3.2. Let the boundary and initial data be non-negative and smooth, $M_0$, $N_0$, $A_0 \in C_0^{\infty}(\Omega)$, $C_0 \in C^{\infty}(\overline{\Omega})$ such that $C_0|_{\partial\Omega}=1$ and $\Vert M_0\Vert_{L^{\infty}(\Omega)} < 1-\rho$ for some $\rho \in (0,1)$. Then, there exist unique solutions $(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon}) \in C^{1,2}(\overline{Q_T})$ of (4), that are non-negative and uniformly bounded w.r.t. $\varepsilon>0$.

Proof. By the classical theory for quasilinear parabolic equations there exist unique solutions $(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon}) \in C^{1,2}(\overline{Q_T})$ of (4) (cf. [25]). We use comparison theorems for quasilinear parabolic equations (e.g., see [2] or [5]) for each equation separately to show the non negativity and uniform boundedness of the solutions.

All components of the solution take non-negative values at the boundary, and the initial data are non-negative. Moreover, we observe that

$g(0,N,C,A)=0=f_2(M,N,0,A),$

and, hence, zero is a subsolution for $M^\varepsilon$ and $C^\varepsilon$. Since $M^\varepsilon$ is non-negative and

$f_1(M,0,C,A)\geq 0\;\; \text{for} \;\;\ M\geq 0,$

we conclude the non-negativity of $N^\varepsilon$. Finally, we observe that

$f_3(M,N,C,0)\geq 0 \;\;\text{for}\;\; \ M\geq 0, N\geq 0,$

which implies that zero is a subsolution for $A^\varepsilon$. Consequently, the parabolic comparison principle implies that $M^{\varepsilon},N^{\varepsilon},C^{\varepsilon}$ and $A^{\varepsilon}$ are non-negative.

To show the uniform boundedness of $M^\varepsilon$ we introduce the time-independent barrier function $M_\theta=1+\theta$, where $\theta$ is a solution of the elliptic problem

$\begin{split} \Delta\theta &= -1 \;\;\;\; \text{in} \ \Omega, \\ \theta|_{\partial\Omega} &= 0. \end{split}$

(12)

The maximum principle for elliptic equations [33] implies that $\theta\geq 0$ in $\Omega$, and

$1\leq M_{\theta}(x)\leq 1+ R_1, \ \ x\in\Omega$

for some constant $R_1\geq 0$. Moreover, we observe that

$M_0 = M^{\varepsilon}|_{t=0} \leq M_{\theta}|_{t=0} , 0=M^{\varepsilon}|_{\partial\Omega}\leq M_{\theta}|_{\partial\Omega}$

and since $M_\theta$ is time-independent, it follows that

$\begin{split} & \partial_t M_{\theta} - \nabla\cdot (D_{\varepsilon}(M_{\theta})\nabla M_{\theta}) - \frac{C^{\varepsilon}}{k_1 + C^{\varepsilon}}M_{\theta} + k_2M_{\theta} +\eta_1\left(\frac{|A^{\varepsilon}|^{n}}{1 + |A^{\varepsilon}|^{ n}}\right)M_{\theta} \\ & = 0+ \frac{d}{\varepsilon^b}- \frac{C^{\varepsilon}}{k_1 + C^{\varepsilon}}M_{\theta} + k_2M_{\theta} + \eta_1\left(\frac{|A^{\varepsilon}|^{n}}{1 + |A^{\varepsilon}|^{n}}\right)M_{\theta} \\ & \geq \frac{d}{\varepsilon^b} - M_{\theta} + k_2M_{\theta} \geq \frac{d}{\varepsilon^b} - (1+R_1) \geq 0\\ & = \partial_tM^{\varepsilon} - \nabla\cdot(D_{\varepsilon}(M^{\varepsilon})\nabla M^{\varepsilon}) -g(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon})\ , \end{split}$

(13)

for all sufficiently small $\varepsilon > 0$. Consequently, by the parabolic comparison principle there exists $\varepsilon'_0>0$ such that the solutions $M^\varepsilon$ are uniformly bounded for all $0<\varepsilon\leq \varepsilon'_0$.

Next, we show the uniform boundedness of the nutrient concentration $C^\varepsilon$ by defining $C_{max} : = \max\{||C_0||_{L^{\infty}(\Omega)}, 1\}$. Then

$\partial_tC_{max} - d_2\Delta C_{max} + \sigma\frac{C_{max}}{k_1 + C_{max}}(M^\varepsilon + N^\varepsilon)= \sigma\frac{C_{max}}{k_1 + C_{max}}(M^\varepsilon + N^\varepsilon)\geq 0,$

(14)

where we used the non-negativity of the sessile biomass $M^{\varepsilon}$ and the dispersed cells $N^{\varepsilon}$. This shows that $C_{max}$ is an upper solution, and hence, by the parabolic comparison principle $C^\varepsilon$ is bounded by $C_{max}.$ To prove the uniform boundedness of $N^\varepsilon$ we denote by $\hat N$ the solution of the initial value problem

$\partial_t \hat N =d_1\Delta \hat N + \hat N -k_2 \hat N+\eta_1 (1+R_1),\\ \nonumber \hat N|_{\partial\Omega}=0,\\ \hat N|_{t=0}=N_0,$

(15)

where $1+R_1$ is the upper bound for the sessile biomass fraction. We observe that $\hat N$ is non-negative, it satisfies $\hat N\in L^\infty(Q_T)$ and

$\partial_t \hat N -d_1\Delta \hat N - \frac{C^\varepsilon}{k_1+C^\varepsilon}\hat N +k_2 \hat N-\eta_1 \frac{|A^\varepsilon|^m}{1+|A^\varepsilon|^m}M^\varepsilon\\ \nonumber \geq \ \partial_t \hat N -d_1\Delta \hat N - \hat N +k_2 \hat N-\eta_1 (1+R_1)=0\\ \nonumber =\ \partial_t N^\varepsilon -d_1\Delta N^\varepsilon - \frac{C^\varepsilon}{k_1+C^\varepsilon}N^\varepsilon +k_2 N^\varepsilon-\eta_1 \frac{|A^\varepsilon|^m}{1+|A^\varepsilon|^m}M^\varepsilon.$

(16)

Consequently, $\hat N$ is an upper solution for the dispersed cells $N^\varepsilon.$

Finally, we prove the uniform boundedness of the quorum sensing signal molecule concentration by showing that there exists a constant $A_{max}\geq \|A_0\|_{L^\infty (\Omega)}$ which is an upper solution for $A^{\varepsilon}$. It satisfies $A_{max}|_{\partial\Omega} \geq 0 = A^{\varepsilon}|_{\partial\Omega}, \ A_{max}|_{t=0} \geq A_0 =A^{\varepsilon}|_{t=0}$ and

$\begin{split} \partial_tA_{max} - d_3\Delta A_{max}&+ \lambda A_{max} - \left[\alpha + \beta\frac{|A_{max}|^m }{1 + |A_{max}|^m }\right](M^{\varepsilon}+N^{\varepsilon}) \\ & = \lambda A_{max} - \left[\alpha + \beta\frac{|A_{max}|^m }{1 + |A_{max}|^m }\right](M^{\varepsilon}+N^{\varepsilon}) \\ &\geq \lambda A_{max} - (\alpha + \beta)(\vert M^\varepsilon\vert + \vert N^\varepsilon\vert)\\ & \geq \lambda A_{max} - (\alpha + \beta)(1+R_1+\hat N) \geq 0, \end{split}$

(17)

where we use the boundedness of $M^\varepsilon$ and $N^\varepsilon$, and $\frac{|A_{max}|^m }{1 + |A_{max}|^m } \leq 1$. This shows that if the constant $A_{max}$ is sufficiently large, it is an upper solution, and by the parabolic comparison principle, it is a uniform upper bound for $A^{\varepsilon}$.

In the following Lemma, we improve the upper bound on the sessile biomass density, in particular, we show that the singularity for $M=1$ is not attained.

Lemma 3.3. Under the hypothesis of Lemma 3.2, there exist $\delta>0$ and $\varepsilon_0>0$ such that $M^{\varepsilon}\leq 1-\delta$ in $Q_T$ for all $0<\varepsilon<\varepsilon_0$.

Proof. In order to improve the upper estimate on $M$ we construct a suitable barrier function and consider the elliptic problem

$\begin{split} \Delta\theta &= -c_1 \text{in} \ \Omega , \\ \theta|_{\partial\Omega} &= c_2 . \end{split}$

(18)

The constants $c_1$ and $c_2$ are defined by

$\begin{split} c_1:=& \Vert g(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon})\Vert_{L^\infty(Q_T)},\\ c_2:=& \Vert\Phi_{\varepsilon}(M_0)\Vert_{L^\infty(\Omega)}, \end{split}$

(19)

for $0<\varepsilon<\varepsilon'_0$, where $\varepsilon_0'$ was defined in the proof of Lemma 3.2, and

$\Phi_\varepsilon(M_0) = \int_0^{M_0} \frac{(s+ \varepsilon)^a}{(1-s)^b}ds \;\;\; \text{for} \;\;\;\ 0\leq M_0<1-\varepsilon.$

(20)

We remark that for sufficiently small $\varepsilon_1<\varepsilon_0'$, the constants $c_1$ and $c_2$ can be chosen uniform for all $0<\varepsilon<\varepsilon_1$. Moreover, the solution $\theta$ of (18) is bounded on $\Omega$, and by the maximum principle for elliptic problems [33] it follows that $\theta \geq c_2$ in $\Omega$.

For $0<\varepsilon < \varepsilon_1$ we define $Z_\varepsilon := \Phi ^{-1}_{\varepsilon}(\theta)$ and observe that

$\partial_tZ_\varepsilon - \Delta(\Phi_\varepsilon(Z_\varepsilon))=c_1=\Vert g(M^{\varepsilon},N^{\varepsilon},C^{\varepsilon},A^{\varepsilon})\Vert_{L^\infty(Q_T)} \geq \partial_tM^{\varepsilon} - \Delta(\Phi_\varepsilon(M^\varepsilon))$

in $Q_T$. Moreover, the boundary conditions imply that

$Z_\varepsilon|_{\partial\Omega} = \Phi^{-1}_{\varepsilon}(\theta)|_{\partial\Omega} = \Phi^{-1}_{\varepsilon}(c_2)\geq M^\varepsilon|_{\partial\Omega}=0,$

and the initial data satisfies

$Z_\varepsilon|_{t=0} = \Phi^{-1}_{\varepsilon}(\theta)|_{t=0} \geq \Phi^{-1}_{\varepsilon}(c_2) \geq \Phi^{-1}_{\varepsilon}(\Phi_\varepsilon(M_0))=M_0 ,$

(21)

where we used the monotonicity of the function $\Phi^{-1}_{\varepsilon}$. Consequently, the function $Z_\varepsilon$ is an upper solution for the sessile biomass $M^\varepsilon$. Using the fact that $\theta$ is bounded in $\Omega$ and that $\Phi_{\varepsilon}(M^\varepsilon)=d\frac{M^\varepsilon}{\varepsilon^b}$ for $M^\varepsilon\geq 1-\varepsilon$, we conclude by the parabolic comparison principle that there exist $0<\varepsilon_0\leq \varepsilon_1$ and $\delta \in (0,1)$ such that $M^\varepsilon\leq Z_\varepsilon= \Phi^{-1}_\varepsilon(\theta) < 1-\delta$ for all $0<\varepsilon < \varepsilon_0$.

Lemma 3.4. Under the hypotheses of Lemma 3.2, the solutions $M^\varepsilon$ of (11) satisfy

$\int_{Q_T} D_\varepsilon(M^\varepsilon)\Big(\partial_t M^\varepsilon\Big)^2 + \sup\limits_{t \in[0,T]} \int_{\Omega} \big | \nabla \Phi_\varepsilon(M^\varepsilon) \big |^2 \leq c, \\ \nonumber \ \int_{Q_T}\Big(|\partial_t N^\varepsilon|^2+|\partial_t C^\varepsilon|^2+|\partial_t A^\varepsilon|^2 \Big) + \sup\limits_{t\in[0,T]}\int_{\Omega}\Big(|\nabla N^\varepsilon|^2 + | \nabla C^\varepsilon|^2+|\nabla A^\varepsilon|^2 \Big)\leq c,$

(22)

for some constant $c> 0.$

Proof. We first multiply the second equation in (11) by $\partial_s N^\varepsilon$ and integrate,

$\int_\tau^t\int_\Omega |\partial_s N^\varepsilon|^2 =-\frac{d_1}{2} \int_\tau^t\int_\Omega \partial_s |\nabla N^\varepsilon|^2 +\int_\tau^t \int_\Omega\partial_s N^\varepsilon f_1(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon),$

(23)

where $0\leq \tau\leq t\leq T.$ Using Young's inequality leads to the estimate

$\int_\tau^t\int_\Omega |\partial_s N^\varepsilon|^2 +\frac{d_1}{2} \int_\Omega |\nabla N^\varepsilon|^2\Big|_{s=t}\\ \nonumber =\frac{d_1}{2} \int_\Omega |\nabla N^\varepsilon|^2\Big|_{s=\tau} +\int_\tau^t\int_\Omega\partial_s N^\varepsilon f_1(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)\\ \leq \frac{d_1}{2} \int_\Omega |\nabla N^\varepsilon|^2\Big|_{s=\tau} +\xi \int_\tau^t\int_\Omega |\partial_s N^\varepsilon|^2 +C_\xi \int_\tau^t\int_\Omega | f_1(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)|^2,$

(24)

for small $\xi>0$, and some constant $C_{\xi}$. Setting $\tau=0$ we obtain

$(1-\xi) \int_{Q_t} |\partial_s N^\varepsilon|^2 +\frac{d_1}{2} \int_\Omega |\nabla N^\varepsilon|^2\Big|_{s=t}\\ \nonumber \leq \frac{d_1}{2} \int_\Omega |\nabla N_0|^2 +C_\xi \int_{Q_t} | f_1(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)|^2,$

(25)

for small $\xi>0$, and some constant $C_{\xi}$. Lemma 3.2 and the continuity of $f_1$ now imply that

$\int_{Q_T} |\partial_t N^\varepsilon|^2\leq C', \sup\limits_{t\in[0,T]}\int_{\Omega}|\nabla N^\varepsilon|^2 \leq C',$

(26)

for some constant $C'\geq 0.$ The corresponding bounds for the solutions $C^\varepsilon$ and $A^\varepsilon$ can be obtained in the same way.

To derive the estimates for the biomass fraction $M^\varepsilon$ let

$G(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon):= \int_0^{M^\varepsilon}D_\varepsilon(\zeta)g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)d\zeta.$

(27)

Multiplying the first equation of (11) by $\partial_s(\Phi_\varepsilon(M^\varepsilon))$ and integrating over $\Omega$ and from $\tau$ to $t$ we obtain

$\int_\tau^t\int_\Omega D_\varepsilon(M^\varepsilon)\Big(\partial_s M^\varepsilon\Big)^2= \int_\tau^t\int_\Omega\partial_s (\Phi_\varepsilon(M^\varepsilon))\partial_s M^\varepsilon \\ \nonumber = \int_\tau^t\int_\Omega\nabla\cdot (D_\varepsilon(M^\varepsilon)\nabla M^\varepsilon)\partial_s(\Phi_\varepsilon(M^\varepsilon)) +\int_\tau^t\int_\Omega\partial_s(\Phi_\varepsilon(M^\varepsilon))g(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon) \\ \nonumber =- \int_\tau^t\int_\Omega D_\varepsilon(M^\varepsilon)\nabla M^\varepsilon\cdot \partial_s\nabla (\Phi_\varepsilon(M^\varepsilon)) +\int_\tau^t\int_\Omega\partial_s(\Phi_\varepsilon(M^\varepsilon))g(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)\\ \nonumber =-\frac{1}{2} \int_\tau^t\int_\Omega \partial_s | D_\varepsilon(M^\varepsilon)\nabla M^\varepsilon|^2 +\int_\tau^t\int_\Omega\partial_s(\Phi_\varepsilon(M^\varepsilon))g(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon),$

(28)

where we used integration by parts in the second step. Consequently, using the identity

$\int_\Omega G(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)\Big|_{s=t} -\int_\Omega G(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)\Big|_{s=\tau}\\ \nonumber =\int_\tau^t\int_\Omega \partial_s \big(G(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)\big)\\ \nonumber =\int_\tau^t\int_\Omega \int_0^{M^\varepsilon} D_\varepsilon(\zeta)\partial_s \big(g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)\big) +\int_\tau^t\int_\Omega \partial_s M^\varepsilon D_\varepsilon(M^\varepsilon)g(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)\\ \nonumber =\int_\tau^t\int_\Omega \int_0^{M^\varepsilon}D_\varepsilon(\zeta)\Big( \partial_s N \partial_N g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)+ \partial_s C \partial_C g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon) \\ \nonumber +\partial_s A \partial_A g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)\Big) +\int_\tau^t\int_\Omega\partial_s(\Phi_\varepsilon(M^\varepsilon))g(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon),$

(29)

it follows that

$\int_\tau^t\int_\Omega D_\varepsilon(M^\varepsilon)\Big(\partial_t M^\varepsilon\Big)^2 +\frac{1}{2} \int_\Omega |D_\varepsilon(M^\varepsilon)\nabla M^\varepsilon|^2\Big|_{s=t} \\ \nonumber =\frac{1}{2} \int_\Omega |D_\varepsilon(M^\varepsilon)\nabla M^\varepsilon|^2\Big|_{s=\tau} +\int_\Omega G(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)\Big|_{s=t} -\int_\Omega G(M^\varepsilon,N^\varepsilon,C^\varepsilon,A^\varepsilon)\Big|_{s=\tau}\\ \nonumber +\int_\tau^t\int_\Omega \int_0^{M^\varepsilon}D_\varepsilon(\zeta)\Big( \partial_s N \partial_N g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)+ \partial_s C \partial_C g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon) \\ \nonumber +\partial_s A \partial_A g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)\Big),$

(30)

where $\partial_N, \partial_C,\partial_A$ denote the partial derivatives w.r.t. the variables $N,C$ and $A,$ respectively. Lemma 3.2 and Lemma 3.3 imply that there exists $\varepsilon_0>0$ and $\delta\in(0,1)$ such that the approximate solutions $M^\varepsilon$ satisfy $M^\varepsilon<1-\delta$ in $Q_T$ for all $0<\varepsilon<\varepsilon_0$, where $\delta$ is independent of $\varepsilon$. Consequently, $D_\varepsilon(M^\varepsilon)$ is positive and uniformly bounded from above by a constant, which is independent of $\varepsilon$,

$d\varepsilon^a\leq D_\varepsilon(M^\varepsilon(t,x))= d\frac{(M^\varepsilon(t,x)+\varepsilon)^a}{(1-M^\varepsilon(t,x))^b} \leq d\frac{(1-\delta+\varepsilon)^a}{(1-(1-\delta))^b}\leq\frac{d}{\delta^b}, (t,x)\in\ Q_T,$

for all $0<\varepsilon<\varepsilon_1:=\min\{\delta,\varepsilon_0\}$. The last integral can therefore be estimated by

$\int_\tau^t\int_\Omega \int_0^{M^\varepsilon}D_\varepsilon(\zeta)\big( \partial_s N \partial_2 g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)+ \partial_s C \partial_3 g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon) \\ \nonumber +\partial_s A \partial_4 g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)\big)\\ \nonumber \leq \frac{d}{\delta^b}\int_\tau^t\int_\Omega \Big( (|\partial_s N|^2 +|\partial_s C|^2 +|\partial_s A |^2) +\int_0^{M^\varepsilon} \sum\limits_{i=2}^4|\partial_i g(\zeta,N^\varepsilon,C^\varepsilon,A^\varepsilon)|^2\Big).$

(31)

By the above estimates (26), Lemma 3.2, Lemma 3.3 and setting $\tau=0$, we conclude that

$\begin{align*} & \int_{Q_T} D_\varepsilon(M^\varepsilon)\Big(\partial_t M^\varepsilon\Big)^2\leq C', \sup\limits_{t\in[0,T]} \int_\Omega |D_\varepsilon(M^\varepsilon)\nabla M^\varepsilon|^2 \leq C', \end{align*}$

for some constant $C'\geq 0$.

Lemma 3.5. Under the hypotheses of Lemma 3.2 the approximate solutions of (11) converge to a solution of the degenerate system (4) as $\varepsilon>0$ tends to zero.

Proof. By Lemma 3.3 the solutions $M^\varepsilon$ are uniformly bounded by $1-\delta$, and the following estimate for the diffusion coefficient was shown in the proof of Lemma 3.4,

$d\varepsilon^a\leq D_\varepsilon(M^\varepsilon(t,x))\leq\frac{d}{\delta^b}, (t,x)\in\ Q_T,$

for all $0<\varepsilon<\varepsilon_1$. Moreover, if $\varepsilon<\varepsilon_1$ we observe that

$D_M(M^\varepsilon)\leq D_\varepsilon(M^\varepsilon)\leq\frac{d}{\delta^b},$

and consequently,

$\Phi(M^\varepsilon(t,x))\leq \Phi_\varepsilon(M^\varepsilon(t,x))\leq (1-\delta)\frac{d}{\delta^b}, (t,x)\in\ Q_T.$

Lemma 3.3 and Lemma 3.4 further imply that

$\int_{Q_T}\Big(\partial_t( \Phi(M^\varepsilon))\Big)^2 =\int_{Q_T} \Big( D_M(M^\varepsilon)\partial_t M^\varepsilon\Big)^2\leq \int_{Q_T} \Big( D_\varepsilon(M^\varepsilon)\partial_t M^\varepsilon\Big)^2\\ \nonumber \leq \|D_\varepsilon(M^\varepsilon)\|_{L^\infty(Q_T)} \int_{Q_T} D_\varepsilon(M^\varepsilon)\Big(\partial_t M^\varepsilon\Big)^2\leq c\frac{d}{\delta^b}, \\ \nonumber \sup\limits_{t\in[0,T]}\int_{\Omega}\big|\nabla \Phi(M^\varepsilon)\big|^2 = \sup\limits_{t\in[0,T]}\int_{\Omega}\big|D_M(M^\varepsilon)\nabla M^\varepsilon\big|^2 \leq \sup\limits_{t\in[0,T]}\int_{\Omega}\big|D_\varepsilon(M^\varepsilon)\nabla M^\varepsilon\big|^2\leq c,$

(32)

for some constant $c\geq 0.$ This shows that the family $\Psi^\varepsilon:=\Phi(M^\varepsilon),\ 0<\varepsilon<\varepsilon_1,$ is uniformly bounded in $W=\{u\in L^\infty(0,T;H^1(\Omega))\,|\,\partial_t u\in L^2(0,T;L^2(\Omega))\},$ which is compactly embedded into $C([0,T];L^2(\Omega))$ by Aubin-Lions' Lemma (e.g., see [7], Theorem Ⅱ.1.5). Consequently, there exists $\Psi\in C([0,T];L^2(\Omega))$ and a sequence $\varepsilon_n$ tending to zero as $n\rightarrow\infty$ such that $\Psi^{\varepsilon_n}\rightarrow \Psi$ in $C([0,T];L^2(\Omega)).$ This implies that $M^{\varepsilon_n}=\Phi^{-1}(\Psi^{\varepsilon_n})\rightarrow M$ and $\Phi(M^{\varepsilon_n})\rightarrow \Phi(M)$ in $C([0,T];L^2(\Omega)).$

Moreover, by Lemma 3.4 the approximate solutions $N^\varepsilon,C^\varepsilon$ and $A^\varepsilon$ are uniformly bounded in $W$, which implies that there exist $N,C,A\in C([0,T];L^2(\Omega))$ and a sequence $\varepsilon_n$ tending to zero as $n\rightarrow\infty$ such that

$N^{\varepsilon_n}\rightarrow N,\ C^{\varepsilon_n}\rightarrow C, \ A^{\varepsilon_n}\rightarrow A \;\;\;\text{in}\;\;\; C([0,T];L^2(\Omega)).$

We can now pass to the limit $\varepsilon\rightarrow 0$ in the distributional formulation of the degenerate system (3) using the uniform boundedness of the approximate solutions and the continuous embedding $C([0,T];L^2(\Omega))\hookrightarrow C([0,T];L^1(\Omega)),$ and conclude that the limits $M,N,C,A$ are solutions of the degenerate problem.

---

3.3. Uniqueness and well-posedness for general initial data

Lemma 3.6. Let the hypotheses of Lemma 3.2 be satisfied. If $(M,N,C,A)$ and $(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A})$ are two solutions corresponding to initial data $(M_0,N_0,C_0,A_0)$ and $(\widetilde{M_0},\widetilde{N_0},\widetilde{C_0},\widetilde{A_0})$, respectively, then

$\begin{split} \Vert M(T)-\widetilde{M}(T)\Vert_{L^1(\Omega)} - \Vert M_0-\widetilde{M_0}\Vert_{L^1(\Omega)} &\leq \int^T_0\int_\Omega |g_0(t,x)|dxdt,\\ \Vert N(T)-\widetilde{N}(T)\Vert_{L^1(\Omega)} - \Vert N_0-\widetilde{N_0}\Vert_{L^1(\Omega)} &\leq \int^T_0\int_\Omega |h_1(t,x)|dxdt,\\ \Vert C(T)-\widetilde{C}(T)\Vert_{L^1(\Omega)} - \Vert C_0-\widetilde{C_0}\Vert_{L^1(\Omega)} &\leq \int^T_0\int_\Omega |h_2(t,x)|dxdt,\\ \Vert A(T)-\widetilde{A}(T)\Vert_{L^1(\Omega)} - \Vert A_0-\widetilde{A_0}\Vert_{L^1(\Omega)} &\leq \int^T_0\int_\Omega |h_3(t,x)|dxdt, \end{split}$

(33)

where the functions $g_0$ and $h_i$ are defined as

$\begin{split} g_0 &:= g(M, N, C, A) - g(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A}),\\ h_i &:= f_i(M, N, C, A) - f_i(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A}), \end{split}$

(34)

for $i=1,2,3$.

Proof. The estimates immediately follow from Lemma 3.3 in [17].

Lemma 3.7. Let $\mathbb{B}$ be a bounded subset of $\ \mathbb{R}_+^4$. Then for all $(M,N,C,A)$ and $(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A}) \in \mathbb{B}$ we have

$\begin{split} &| g(M,N,C,A) - g(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A}) | + \sum\limits_{i=1}^3 | f_i(M,N,C,A) -f_i(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A})\\ \leq &\ c \left(| M- \widetilde{M} | + | N- \widetilde{N} | + | C- \widetilde{C} | + | A- \widetilde{A} | \right), \end{split}$

(35)

for some constant $c\geq 0$.

Proof. Let $(M,N,C,A),(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A}) \in \mathbb{B}$. For the function $f_2$ we obtain,

$\begin{split} &| f_2(M,N,C,A) - f_2(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A})| \\ =& -\sigma\Big| \left[\left(\frac{C}{k_1+C}\right)(M+N) - \left(\frac{\widetilde{C}}{k_1+\widetilde{C}}\right)(\widetilde{M}+\widetilde{N}) \right]\Big| \\ \leq&\ r_1\Big (| M-\widetilde{M}| + | N-\widetilde{N} | \Big ), \end{split}$

(36)

for some constant $r_1\geq 0$. To show that the functions $f_1,f_3$ and $g$ satisfy (35) we observe, that

$\begin{split} &A^m_1X_1 - A^m_2X_2 = A^m_1(X_1-X_2) + X_2(A^m_1-A^m_2)\\ =&\ A^m_1(X_1-X_2) + \nu X_2(A_1-A_2)\int_0^1(sA_1 + (1-s)A_2)^{m-1}ds, \end{split}$

(37)

which implies that $|A^m_1X_1 - A^m_2X_2| \leq r|(X_1-X_2) + X_2(A_1-A_2)|$ for some constants $\nu$ and $r\geq 0$ and $(X_1,A_1), (X_2,A_2)$ in bounded subsets of $\mathbb{R}_+^2$.

Applying this to $f_1$, we obtain

$\begin{split} &| f_1(M,N,C,A) - f_1(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A})| \\ =&\ \bigg| \left[\left(\frac{C}{k_1+C} - k_2\right)N + \eta_1\left(\frac{A^m}{1+A^m}\right)M\right] \\ &\ - \left[\left(\frac{\widetilde{C}}{k_1+\widetilde{C}} - k_2\right)\widetilde{N} + \eta_1\left(\frac{\widetilde{A}^m}{1+\widetilde{A}^m}\right)\widetilde{M} \right]\bigg| \\ \leq &\ | 1-k_2\vert | N-\widetilde{N}| + |\eta_1| \Big [| M-\widetilde{M} | + | M(A-\widetilde{A})|\Big ] \\ \leq& \ r_2(| N-\widetilde{N}| + | M-\widetilde{M}| + | A-\widetilde{A}|) \end{split}$

(38)

for some constant $r_2\geq 0$. Similarly, for $f_3$ we obtain

$\begin{split} & | f_3(M,N,C,A) - f_3(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A})| = \ \bigg|\left[-\lambda A + \left[\alpha + \beta \frac{A^m}{1+A^m}\right](M+N)\right] \\ &\ - \bigg[-\lambda \widetilde{A} + \bigg[\alpha + \beta \frac{\widetilde{A}^m}{1+\widetilde{A}^m}\bigg](\widetilde{M}+\widetilde{N})\bigg]\bigg|\\ \leq&\ \lambda | A - \widetilde{A} | + \alpha \Big[| M-\widetilde{M}| + | N-\widetilde{N}|\Big] \\ &\ + \beta \Big [| M-\widetilde{M}| + | \widetilde{M}| | A-\widetilde{A}|\Big] + \beta [| N-\widetilde{N}| + |\widetilde{N}| | A-\widetilde{A}|\Big ] \\ \leq&\ r_3(| A-\widetilde{A}| + | M-\widetilde{M}| + | N-\widetilde{N}|) \end{split}$

(39)

for some constant $r_3\geq 0$. Finally, for $g$ we obtain

$\begin{split} &| g(M,N,C,A) - g(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A})| \\ =&\ \Big|\left[\left(\frac{C}{k_1+C} - k_2\right)M - \eta_1\left(\frac{A^m}{1+A^m}\right)M \right] \\ &\ - \left[\left(\frac{\widetilde{C}}{k_1+\widetilde{C}} - k_2\right)\widetilde{M} - \eta_1\left(\frac{\widetilde{A}^m}{1+\widetilde{A}^m}\right)\widetilde{M} \right]\Big| \\ \leq&\ | 1 -k_2 | | M-\widetilde{M}| + \eta_1 \Big[| M - \widetilde{M}| + | \widetilde{M}|| A-\widetilde{A}|\Big] \\ \leq &\ r_4| M - \widetilde{M}| + r_2| A-\widetilde{A}| \end{split}$

(40)

for some constant $r_4\geq 0$. Combining equations (36), (38), (39) and (40) gives

$\begin{split} &| g(M,N,C,A) - g(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A}) | + \sum\limits_{i=1}^3 | f_i(M,N,C,A) -f_i(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A})|\\ \leq& \ r^\prime (| M- \widetilde{M} | + | N- \widetilde{N} | + | C- \widetilde{C} | + | A- \widetilde{A} |) \end{split}$

(41)

for some constant $r^\prime\geq 0$.

Theorem 3.8. For every $T>0$ and initial data $(M_0,N_0,C_0,A_0)$ in $C(\overline{\Omega})$ such that

$\begin{align*} M_0 &\geq 0, \Vert M_0 \Vert_{L^\infty(\Omega)} < 1- \rho, N_0\geq 0, C_0 \geq 0, A_0 \geq 0,\\ C_0|_{\partial\Omega}&=1, \ \, M_0|_{\partial\Omega}= N_0|_{\partial\Omega}= A_0|_{\partial\Omega}=0, \end{align*}$

for some $\rho \in (0,1)$, there exists a unique solution $(M, N, C, A)$ on $[0,T]$ of model (4). The functions $M, N, C$ and $A$ are non-negative, belong to $L^\infty(Q_T)$ and there exists $\delta \in (0,1)$ such that $M$ satisfies $\Vert M \Vert _{L^\infty(Q_T)} < 1- \delta$.

Proof. We first assume the initial data are smooth and satisfy the hypotheses of Lemma 3.2.

If $(M,N,C,A)$ and $(\widetilde{M},\widetilde{N},\widetilde{C},\widetilde{A})$ are two solutions corresponding to initial data $(M_0,N_0,C_0,A_0)$, and $(\widetilde{M}_0,\widetilde{N}_0,\widetilde{C}_0,\widetilde{A}_0)$ respectively, then Lemma 3.6 and the estimate in Lemma 3.7 imply that

$F(T)-F(0) \leq c\int^T_0 F(s) ds,$

where

$F(t) := \Vert M(t) - \widetilde{M}(t)\Vert_{L^1(\Omega)} + \Vert N(t) - \widetilde{N}(t)\Vert_{L^1(\Omega)} \\ \nonumber + \Vert C(t) - \widetilde{C}(t)\Vert_{L^1(\Omega)} + \Vert A(t) - \widetilde{A}(t)\Vert_{L^1(\Omega)}.$

(42)

By Gronwall's Lemma we conclude that

$F(T) \leq F(0)e^{cT},$

(43)

for some constant $c\geq 0$, which implies the uniqueness of solutions corresponding to smooth initial data. For general initial data, let $M^n_0,N^n_0,C^n_0,A^n_0$ be an approximating sequence satisfying the hypothesis of Lemma 3.2 such that

$\Vert M^n_0 - M_0\Vert_{L^1(\Omega)} + \Vert N^n_0 - N_0\Vert_{L^1(\Omega)} + \Vert C^n_0 - C_0\Vert_{L^1(\Omega)} + \Vert A^n_0 - A_0\Vert_{L^1(\Omega)} \rightarrow 0$

as $n \rightarrow \infty$. The Lipschitz continuity of solutions in $L^1(\Omega)$ norm (43) implies that the corresponding solutions $M^n,N^n,C^n,A^n$ satisfy

$\begin{split} &\sup\limits_{t\in [0,T]} \Big\{\Vert M^n(t) - M^k(t)\Vert_{L^1(\Omega)} + \Vert N^n(t) - N^k(t)\Vert_{L^1(\Omega)}\\ & + \Vert C^n(t) - C^k(t)\Vert_{L^1(\Omega)} + \Vert A^n(t) - A^k(t)\Vert_{L^1(\Omega)}\Big\}\\ \leq&\ C_T\Big(\Vert M^n_0-M^k_0\Vert_{L^1(\Omega)} + \Vert N^n_0-N^k_0\Vert_{L^1(\Omega)} \\ & + \Vert C^n_0-C^k_0\Vert_{L^1(\Omega)} + \Vert A^n_0-A^k_0\Vert_{L^1(\Omega)}\Big). \end{split}$

(44)

for some constant $C_T \geq 0$ and all $n,k \ \in \ \mathbb{N}$. Consequently, $M^n,N^n,C^n,A^n$ form a Cauchy sequence in $C([0,T]; L^1(\Omega))$ and converge to the unique solution $M,N,C,A$ of (4).

---

4. Numerical simulations

In the previous section we established the well-posedness of the quorum sensing induced biofilm detachment model, however, we are currently unable to describe the solutions of the model qualitatively based on rigorous analytical arguments. Hence, we illustrate the model behaviour in computer simulations.

For the numerical solution of the model we use a straightforward extension to the problem at hand of the numerical method for the prototype biofilm model that is described in detail in [14,36]. For space discretisation this uses a Finite Volume method on a uniform grid, which uses Finite Difference approximations for the diffusive fluxes across grid cell edges. We extended this method to account for the new dependent variables $A$ and $N$ which are treated in the same manner as $C$. A simple first order semi-implicit method is used for time-stepping. In every time step, this requires the solution of a sparse linear system for each dependent variable. By construction of the method, the system matrices are at least weakly diagonally dominant. These linear systems are solved with the stabilized biconjugate gradient method [35]. The linear solver is prepared for parallel execution on multi-core and shared memory multiprocessor architectures using OpenMP as described in [28]. Simulations are terminated when the biofilm (or the microbial floc) reaches a set target size, or when a set maximum simulation time is reached. The model parameters used in our simulations are summarised in Table 1.

In the numerical experiments, we investigate the model behaviour under different boundary conditions reflecting biofilms or microbial flocs with particular emphasis on the internal structure of the colonies. We define the following output parameters:

• Relative variation: This is the standard deviation of the sessile biomass density from its mean in the biofilm,

$R(t) := \left[\int_{\Omega_2}\left( M(t,x) -\int_{\Omega_2}M(t,y)dy \right)^2dx\right]^{\frac{1}{2}}$

(45)

• Relative biofilm (floc) size: This is the size of the biofilm relative to the domain size,

$\omega (t) := \frac{1}{|\Omega |} \int_{\Omega_2(t)} dx$

(46)

• Sessile biomass in the biofilm:

$M_{tot} (t) := \int_{\Omega} M(t,x) dx$

(47)

• Dispersed cells:

$N_{tot} (t) := \int_{\Omega} N(t,x) dx$

(48)

• Average signal molecule concentration in the biofilm:

$A_{ave}(t) := \frac{\int_{\Omega_2(t)} A(t,x) dx}{ \int_{\Omega_2(t)} dx }.$

(49)

In these definitions, $\Omega$ is the fixed computational domain, and $\Omega_2(t)\subset \Omega$ is the biofilm subdomain where $M>0$, which evolves in time, as defined in Section $2$.

---

4.1. Microbial floc

In the first simulation experiments, we consider a biofilm without substratum i.e. a microbial floc. We restrict ourselves to a two-dimensional setting with rectangular computational domain $\Omega = [0,1]\times [0,1]$. The boundary conditions used are the Dirichlet boundary conditions (5). One small circular floc is placed in the centre of the domain which at time $t=0$ contains only sessile cells of biomass density $M_0=0.1$, while everywhere else in the domain $M_0\equiv 0$; this floc is circular and occupies 0.03% of the domain $\Omega$.

Initially, no signal molecule $A$ and no dispersed cells $N$ are assumed to be in the system, thus $A_0\equiv 0$ and $N_0\equiv 0$. The nutrient concentration in the interior at time $t=0$ is everywhere at the same level as the concentration at the boundary, $C_0\equiv 1$.

For the visualization presented in Figure 2, we used an autoinducer production rate $\alpha=30.7$ and a maximum dispersal rate $\eta_1=0.6$, every other parameter is as listed in Table 1. Immediately after the simulation starts, the floc is very small in size and the concentration of the signal molecule $A$ is very low. At the snapshots $t=5$, $t=10$ and $t=15$, the original circular shape of the floc is still visible, the nutrient is not strongly growth limiting and the floc expands. As the floc grows within the rectangular domain, it looses its circular shape and expands faster toward the sides of the domain rather than towards the corners. This is a boundary condition effect that is prompted by the smaller diffusion length for substrates from the middle of the lateral sides than from the corners. Overall, however, the colony remains in a rather compact shape, in analogy with biofilms which have been found to grow in compact, homogenous layers when nutrients are nowhere severely limited. With the increase in biomass density in the system, the signal concentration also increases. The autoinducer concentration attains its highest values inside the floc and, due to diffusion and the parabolic maximum principle, decreases from there towards the floc/water interface and the boundary of the domain where the signal concentration is kept at $A=0$. At the next snapshot $t=20$, we observe a visible decrease in biomass density in the centre of the floc, thus creating a hollow in the floc due to conversion of sessile cells into dispersed cells and their subsequent diffusion out of the domain. The last snapshot taken at $t=25$ shows that the microbial floc has started growing again both internally and on the surface; and consequently the signal concentration has started rising.

The 2-D structural representation of the biofilm shown in Figure 2 reveals that the bacterial cells leave from the inner core of the microfloc, but it does not illustrate the extent of the dispersal effect on the biofilm structure. Hence, we present a 1-D spatial representation of a typical biofilm dispersal event in Figure 3 whereby the biofilm is cut vertically to reveal the dispersal taking place in the inner core of the biofilm. The sessile biomass density and the concentration of the quorum sensing molecule are shown at selected times $t=0.0002,5,10,15,20,25$. From the top view of the microfloc, we observe that at times $t=0.0002$, $t=5$, and $t=10$ the microfloc is still in its initial growth phase, the sessile biomass density has increased leading to spatial spreading. The signal molecule concentration also rises until its level is high enough to induce cell dispersal, as seen at the next selected time $t=15$. The following two snapshots at times $t=20$ and $t=25$ show that the dispersal of cells from the microfloc creates voids in the floc and its depth increases over time (as seen from the top view).

Figures 2 and 3 illustrate that cell dispersal occurs from the inner core of the microfloc thereby creating hollowing structures, as reported in experimental studies, e.g. [9,22]. Furthermore, we will investigate the general extent of the hollow effect over a longer period of time through the lumped quantities $M_{tot},N_{tot}, A_{ave}$ and $R$ defined in equations (45)-(49).

The temporal plots are shown in Figure 4 where the autoinducer production rate is varied as $\alpha =92.0, 46.0, 30.7, 23.0, 18.4, 15.3, 13.1$. For $\alpha \geq 30.7$, we observe a rapid removal of biomass and signal molecule (see Figure 4c, e) due to the Dirichlet conditions prescribed on the boundary. After the first dispersal event, the bacterial cells that are left behind in the floc are too few to maintain the signal concentration at a level high enough to sustain cell dispersal, so the biomass in the floc starts increasing. This also leads again to an increase in autoinducer concentration until the next dispersal event is triggered. This pattern continues, thus producing an almost periodic pattern of discrete dispersal events, cf. Figure 4. For $\alpha < 30.7$, the floc grows larger with higher levels of biomass density and autoinducer concentration before cell dispersal is induced. Here, the dispersal appears continuously and the biomass density in the floc reaches a plateau. The simulations in Figure 4 illustrate that cell dispersal is taking place, depending on the value of $\alpha$ in a continuous or periodic pattern. The relative size of the biofilm colony, however, does not reflect this as seen in Figure 4b. It indicates that most of the dispersed cells are from the inner core of the floc.

The relative variation of the biomass density defined in equation (45) is the standard deviation of the sessile biomass density from its mean in the microfloc. This is evaluated and plotted in Figure 4d. So far, we established that cell dispersal occurs in the inner core of the biofilm, it creates voids (hollows) whose depth increases over time when viewed from the top and the floc size does not shrink as a result of dispersal. The importance of the variable $R$ is to investigate the spatiotemporal effect of cell dispersal, in order to analyze the extent of the voids over time. When the standard deviation is close to zero, the hollow in the floc is small (or narrow); on the other hand, if the standard deviation is large, the hollow is big (or wide). The plots in Figure 4d show that the size of the hollow decreases over time as the standard deviation from the mean decreases. This is more visible in smaller flocs (i.e. when $\alpha \geq 30.7$). Hence, after the voids are created by cell dispersal, the floc grows both inside and outside. Loss of biomass occurs in the inner core even though the total floc size does not shrink. For bigger flocs (i.e. when $\alpha < 30.7$), the dispersal occurs rather continuously. The floc hardly grows, especially in the inner core, and hence the created hollow becomes constant, seen in Figure 4d as $R$ approaches a plateau.

So far, we presented simulation experiments for a microfloc. We will compare this situation with a biofilm. Besides, we will investigate the effect of quorum sensing induced dispersal on merged colonies which has not been done in the previous study [18]. More specifically, we will analyze the behaviour of merging colonies in terms of biomass growth, cell dispersal and hollow structures.

---

4.2. Microbial biofilm

The simulated biofilm community consists of bacterial cells accumulating on a surface (substratum) surrounded by an aquatic region. The substratum is inoculated by two colonies of small pockets of sessile cells with biomass density $M_0=0.1$ in each colony.

The substratum forms the bottom boundary of the domain $\Omega$ and is impermeable to biomass, substrate and signal molecule. This is described by homogeneous Neumann boundary conditions imposed on the dependent variables $M,N,C,A$. These can be understood as symmetry boundary conditions which enable us to view the small simulation section as a part of a much larger system. At the top boundary, we impose homogenous Dirichlet conditions for the dependent variables $M$, $N$ and $A$ by setting $M=0$, $N=0$ and $A$; this enforces a diffusion gradient from the biofilm in the interior of the domain to the boundary and mimics removal of quorum sensing molecule and dispersed cells into the surrounding bulk phase, where they are negligible due to instantaneous dilution. An inhomogeneous Dirichlet condition is imposed for $C$ at the top boundary by setting $C$ to the bulk concentration value (i.e. $C= C_\infty=1$); this reflects that the growth limiting substrate is added to the system through the top boundary. For the left and right boundaries, we impose homogeneous Neumann boundary conditions for all the dependent variables $M,N,C,A$ which makes the left and right boundaries impermeable to biomass, substrate and signal molecules.

Initially, no dispersed cells and no autoinducers are in the system, and the concentration of nutrients is at bulk level, i.e. $C_0\equiv C_{\infty}=1, N_0=A_0\equiv 0$.

To investigate the effect of quorum sensing triggered dispersal on the spatial structure of the biofilm colony, we visualize the development, growth and cell dispersal of the biofilm in Figure 5. We used the autoinducer production rate $\alpha=30.7$ and a maximum dispersal rate $\eta_1=0.6$, which mimics the periodic dispersal events described in [18]. All other parameter values used are as listed in Table 1. We show the spatial distribution of the sessile biomass $M$ and the iso-lines of the autoinducer concentration $A$ within the first dispersal event (cycle); the selected time instances illustrate the biofilm growth and dispersal in response to autoinducer concentration.

After the simulation starts, the biomass $M$ starts growing and colonies gradually expand as shown in the snapshots at $t=0.002$ and $t=5$, but due to the removal of signal molecules from three out of four boundary segments of the domain, accumulation of signal molecules is not strong enough to induce dispersal. Expansion starts locally when and where the biomass density $M$ approaches its maximum value. As long as the nutrient substrate is not severely limited, the total biomass density $M$ is close to $1$ in the biofilm.

At the next shown time instance $t=10$, the biomass density $M$ inside the colonies has increased. The colonies which were initially placed apart are now drawing close and then merge, albeit the biomass density inside the colonies still remains below the maximum biomass density. The autoinducer concentration has risen, though not strong enough to induce cell dispersal. Moreover, the merged colonies now act as one colony, whereby the iso-lines of the signal concentration are no longer separated. The autoinducer concentration in the middle of the merged colony has increased and few cells started dispersing as shown at the time instance $t=11$.

At the next snapshot $t=12$, we notice a huge biomass loss resulting in a significant decrease in the biomass density of the merged colony. An important observation is the void region of the colonies, seen at the centre of the merged colonies and not at the centre of the individual colonies as in the study [18]. At the snapshot $t=13$, the biomass density in the merged colony has started to increase again.

Furthermore, we investigate the behaviour of the biofilm if cell dispersal occurs before merging of the colonies. The simulation setup here is the same as in Figure 5 except that the autoinducer production rate is set to $\alpha=92.0$ which relates to a small biofilm. Here, the biofilm grows less and merging of the colonies is delayed. Similar to the results in Figure 5, we notice that at the first two time instances $t=0.002$ and $t=5$, the colonies are still in the growth phase (see Figure 6). At the next two snapshots $t=10$ and $t=11$, we observe a significant biomass loss leading to hollow structures. After the dispersal event, the biofilm starts growing again albeit the colonies do not merge (see Figure 6f); the inability of the colonies to merge could be due to mass transfer. The signal concentration attains its highest values at the sides of the colony close to the neighbouring colonies, in accordance with the parabolic maximum principle, and hence more cells disperse from the sides which inhibits merging.

The extent of the hollow effect in the quorum sensing induced biofilm dispersal is investigated through the lumped quantities $M_{tot}, N_{tot}, A_{ave}, \omega$ and $R$, shown in Figure 7, by varying the value of $\alpha$ as $\alpha = 92.0, 46.0, 30.7, 23.0, 18.4, 15.3, 13.1$. The time evolution profile presented in Figure 7 for the biofilm is similar to the one for the microbial floc in Figure 4. Small values of $R$ indicate small hollows while large values of $R$ indicate large deviation of the bacterial density from the mean, i.e. large hollows. The deviations in $R$ show that the size of the hollows is not steady but rather decreases as new cells are produced during the biofilm intermittent growth phase. It is interesting that the size of the void in the biofilm changes as dispersal and growth takes place, which is not evident when dispersal happens continuously, whereby $R$ attains a near constant value.

In the simulations discussed above, we prescribed homogeneous Dirichlet boundary conditions for the autoinducer $A$. We compare these results with the situation where homogeneous Neumann boundary conditions are imposed for the autoinducer. This assumption is made solely to isolate and emphasise the effect of accumulating QS molecules vs their removal from the system, and not thought of as a simulation of a realistic physical environment. The comparison is presented in Figure 8 where we used an autoinducer production rate $\alpha = 18.4$ and maximum dispersal rate $\eta_1=0.6$; for the case of a microbial floc (see Figure 8a) with homogeneous Neumann boundary conditions imposed for $M, N, C, A$ on all the boundaries, and for the microbial biofilm (see Figure 8b) with homogeneous Neumann conditions imposed for $M, N, C, A$ on all boundaries except for the top boundary. The Neumann boundary condition models the case where the signal molecule cannot leave the domain and therefore accumulates faster. As a consequence, the onset of quorum sensing occurs much earlier than with Dirichlet boundary conditions. As a result of unhindered accumulation of signal molecules in the system, very high autoinducer concentrations are attained, hence quick up-regulation and cell dispersal occurs which reduces the density of sessile biomass in the floc. In the case of Dirichlet boundary conditions, the accumulation of the signal molecule is slower due to its early removal through the boundaries, hence causing a delay in up-regulation and cell dispersal.

---

5. Conclusion

In order to describe dispersal of cells from a biofilm colony into the aqueous environment, simple prototype biofilm growth models must be extended to include both the dispersing cells as well as the trigger that causes such detachment. In the case of quorum sensing induced dispersal these are autoinducer molecules that are produced in the biofilm colony. In total, this introduces two new dependent variables. The existence and uniqueness proofs for the underlying prototype model must be adapted and extended to account for this new complexity. Our particular model is based on a biofilm growth model that consists of a density-dependent diffusion-reaction equation for sessile biomass that is coupled with a semi-linear diffusion-reaction equation for nutrient, which are consumed in the biofilm. The extension of the model introduces additional semi-linear diffusion-reaction equations that describe quantities that are produced in the biofilm. Since the associated reaction terms have opposite signs than the nutrient terms, the approach to obtain estimates for the nutrient concentration does not carry over to the two new dependent variables and alternative arguments were developed. We formulated the analysis for the case of Dirichlet boundary conditions, but the results can be generalized to other situations with the same ideas that were used for the underlying prototype biofilm growth model. In numerical simulations we focus on the effect of quorum sensing controlled dispersal on the colony. The simulations suggest: (ⅰ) Depending on parameters, microflocs and biofilms do not shrink as a result of dispersal but hollow out, with lower biomass densities in the inner layers of the colonies. (ⅱ) After a rapid dispersal event the number of cells remaining in the colony drops which also leads to a drop in the amount of autoinducer molecules produced; if sufficient nutrients are available cells grow inside the colonies after the dispersal event, leading to increasing biomass there, i.e. a shrinkage of the inner hollow regions, until the next dispersal event. Thus, (ⅲ) quorum sensing induced hollowing of biofilm colonies is a dynamic feature, changing in size and depth over time.

---

Acknowledgments

We thank the referees for their careful reading of the manuscript and their helpful comments and remarks that greatly improved the writing of the paper. BOE acknowledges financial support received from the Ontario Ministry for Agriculture, Food and Rural Affairs through the HQP Scholarship program; HJE acknowledges financial support received from the Natural Science and Engineering Research Council of Canada through the Discovery Grant program.

---