A Class of Mathematical Models Describing Processes in Spatially Heterogeneous Biofilm Communities

We present a class of deterministic continuum models for spatially heterogeneous biofilm communities. The prototype is a single-species biofilm growth model, which is formulated as a highly non-linear system of reaction-diffusion equations for the biomass density and the concentration of the growth controlling substrate. While the substrate concentration satisfies a standard semi-linear reaction-diffusion equation the equation for the biomass density comprises two non-linear diffusion effects: a porous medium-type degeneracy and super diffusion. When further biofilm processes are taken into account equations for several substrates and multiple biomass components have to be included in the model. The structure of these multi-component extensions is essentially different from the mono-species case, since the diffusion operator for the biomass components depends on the total biomass in the system and the equations are strongly coupled. We present the prototype biofilm growth model and give an overview of its multi-component extensions. Moreover, we summarize analytical results that were obtained for these models. Keywords-biofilm; antibiotic disinfection; probiotic control; quorum-sensing; degenerate reaction-diffusion systems;


I. INTRODUCTION
The dominant mode of microbial life in aquatic ecosystems are biofilm communities rather than planktonic cultures ( [1]).Biofilms are dense aggregations of microbial cells encased in a slimy extracellular matrix forming on biotic or abiotic surfaces (called substrata) in aqueous surroundings.Such multicellular communities are a very successful life form and able to tolerate harmful environmental impacts that would eradicate free The author is supported by the ERC Advanced Grant FPT-246775 NUMERIWAVES.
floating individual cells ( [3], [17]).Whenever environmental conditions allow for bacterial growth, microbial cells can attach to a substratum and switch to a sessile life form.They start to grow and divide and produce a gel-like layer of extracellular polymeric substances (EPS) often forming complex spatial structures.The selfproduced EPS yields protection and allows survival in hostile environments.For instance, the mechanisms of antibiotic resistance in biofilm cultures are essentially different from those of free swimming cells, which makes it difficult to eradicate bacterial biofilm infections.The EPS retards diffusion of antibiotics and the antibiotic agents fail to penetrate into the inner cores of the biofilm ( [3], [17], [4]).
Biofilms play a significant role in various fields.They are beneficially used in environmental engineering technologies for groundwater protection and wastewater treatment.However, in most occurrences biofilm formations have negative effects.If they form on implants and natural surfaces in the human body they can provoke bacterial infections such as dental caries and otitis media ( [3]).Biofilm contamination can lead to health risks in food processing environments, and biofouling of industrial equipment or ships can cause severe economic defects for the industry ( [4], [21]).
Mathematical models of biofilms have been studied for several decades.They range from traditional onedimensional models that describe biofilms as homogeneous flat layers, to more recent two-and threedimensional biofilm models that account for the spatial heterogeneity of biofilm communities.A variety of mathematical modeling concepts has been suggested, including discrete stochastic particle based models and deterministic continuum models, that are based on the description of the mechanical properties of biofilms ( [7], [21]).We are concerned with the latter, where biofilm and liquid surroundings are assumed to be continua, and its time evolution is governed by systems of deterministic PDEs.The first continuum model [22] was a one-dimensional biofilm growth model and essentially based on the assumption that biofilms are homogeneous flat layers.Such models serve well for engineering applications on the macro-scale (larger than 1cm) are, however, not capable to predict the often highly irregular spatial structure of microbial populations and the behavior of biofilms on the meso-scale (between 50µm and 1mm), the actual length scale of mature biofilms ( [7]).Biofilms can form mushroom-like caps and contain clusters and channels, where substrates can circulate.Cells in different regions of the biofilm live in diverse micro-environments and exhibit differing behavior ( [3]).
To capture the spatial heterogeneity of biofilms a higher dimensional biofilm growth model was proposed in [6], which is based on the interpretation of a biofilm as a continuous, spatially structured microbial population.The essential difficulty is to model the spatial spreading mechanism of biomass and to reproduce the growth characteristics of biofilms that have been observed in experiments ( [6]): • Biofilm and aqueous surroundings are separated by a sharp interface.• The biomass density is bounded by a known maximum value.• Spatial spreading only takes place where the local biomass density approaches values close to its maximum possible value, while it does not occur in regions where the biomass density is low.
The mathematical model [6] is formulated as system of highly non-linear reaction-diffusion equations for the biomass density and the concentration of a growth limiting nutrient, and is the prototype of the biofilm models we consider.While the substrate concentration satisfies a standard semi-linear reaction-diffusion equation the governing equation for the biomass density exhibits two non-linear diffusion effects.The biomass diffusion coefficient degenerates like the porous medium equation and shows super diffusion, which causes difficulties in the mathematical analysis of the model.It was shown by numerical experiments that the model is capable of predicting the heterogeneous spatial structure of biofilms and is in good agreement with experimental findings ( [6]).In [10] and [9] the model was studied analytically.In particular, the well-posedness of the model and the existence of a compact global attractor was shown.
The prototype single-species single-substrate model was extended to model biofilms which consist of several types of biomass and account for multiple dissolved substrates.The model introduced in [4] describes the diffusive resistance of biofilms against the penetration by antibiotics.In [13] an amensalistic biofilm control system was modelled, where a beneficial biofilm controls the growth of a pathogenic biofilm.The structure of these multi-species models differs essentially from the mono-species model, and the analytical results for the prototype model could not all be carried over to the more involved multi-species case.In both articles, existence proofs for the solutions were given, and numerical studies presented, but the question of uniqueness of solutions remained unanswered in [4] and [13].Recently, another multi-component biofilm model was proposed in [11], which describes quorum-sensing in growing biofilm communities.Quorum-sensing is a cell-cell communication mechanism used by bacteria to coordinate behavior in groups.The model behavior was studied by numerical experiments in [11] and [21], analytical questions were addressed in [21].Compared to the multicomponent biofilm models [4] and [13], the particularity of the quorum-sensing model is, that adding the governing equations for the involved biomass components we recover exactly the monospecies biofilm growth model.Taking advantage of the known results for the prototype model the existence and uniqueness of solutions and the continuous dependence of solutions on initial data could be established in [21].It is the first uniqueness result for multi-species reaction-diffusion models of biofilms that extend the single-species model [6].
In Section II we introduce the prototype biofilm growth model.Multicomponent extensions are addressed in Section III.In Section IV we give an overview of the analytical results obtained for these models.For numerical simulations and further details we refer to [6], [10], [13], [4], [21] and [8].In [8] also extensions in other directions of the prototype biofilm growth model are discussed, e.g., models that take the effect of hydrodynamics into account.

II. PROTOTYPE BIOFILM GROWTH MODEL
The multi-dimensional biofilm growth model proposed in [6] is formulated as a non-linear reaction-diffusion system for the biomass density and the concentration of the growth controlling nutrient in a bounded domain Ω ⊂ R n , n = 1, 2, 3, where the boundary of the domain ∂Ω is piecewise smooth.In dimensionless form the substrate concentration S is scaled with respect to the bulk concentration, and the biomass density is normalized with respect to the maximal bound for the cell density.Consequently, the dependent model variable M represents the volume fraction occupied by biomass.The EPS is implicitly taken into account, in the sense that the biomass volume fraction M describes the sum of biomass and EPS assuming that their volume ratio is constant.Both unknown functions depend on the spatial variable x ∈ Ω and time t ≥ 0, and satisfy the parabolic system where the constants d, d S and k 2 are positive, the constants k 1 , k 3 and k 4 are non-negative.Furthermore, ∂ t denotes the partial derivative with respect to time t > 0, ∆ the Laplace operator and the gradient with respect to the spatial variable x ∈ Ω and • the inner product in R n .The solid region occupied by the biofilm as well as the liquid surroundings are assumed to be continua.The actual biofilm is described by the region and the liquid area by The substratum, on which the biofilm grows, is part of the boundary ∂Ω.Biomass is produced due to the consumption of nutrients, which is described the Monod reaction functions Natural cell death is also included in the model and given by the lysis rate k 4 in the equation for the biomass fraction.
While the nutrient is dissolved in the domain and the substrate concentration satisfies a standard semilinear reaction-diffusion equation, the spatial spreading of biomass is determined by the density-dependent diffusion coefficient The biomass motility constant d is small compared to the diffusion coefficient d S of the dissolved substrate, which reflects that the cells are to some extent immobilized in the EPS matrix.Accumulation of biomass leads to spatial expansion of the biofilm.We observe that the biomass diffusion coefficient vanishes when the total biomass approaches zero and blows up when the biomass density tends to its maximum value.The polynomial degeneracy M a is well-known from the porous medium equation and guarantees that spatial spreading is negligible for low values of M .Moreover, it yields the separation of biofilm and liquid phase, that is, a finite speed of interface propagation.Spreading of biomass takes place when and where the biomass fraction takes values close to its maximal value.When M = 1 instantaneous spreading occurs, which is known as the effect of super diffusion.This singularity at M = 1 of the biomass diffusion coefficient ensures the maximal bound for the biomass density.
It was shown in numerical experiments that the model ( 1) is in good agreement with experimental findings and is capable to reproduce the irregular, heterogeneous spatial structure of biofilms observed on the mesoscale ([6], [10]).More precisely, the simulations show that the biofilm develops a rather regular, homogeneous structure if nutrients are nowhere limited in the system.On the other hand, when the nutrient supply is not symmetric and nutrients become limited the colonies grow in the direction of higher nutrient concentrations, which can lead to cluster and channel morphologies and mushroomshaped architectures.

III. MULTICOMPONENT BIOFILM MODELS
The prototype biofilm growth model (1) was extended to incorporate further biofilm processes.It requires to distinguish different types of biomass and dissolved substrates and to include governing equations for these multiple biomass fractions and dissolved substrates in the model.The pattern of the multi-component biofilm models is essentially different from the prototype model, the equations for the biomass components are strongly coupled through the diffusion operators.In this section we discuss multi-component biofilm models, that were proposed and studied in [4], [5], [11], [13], [21].

A. Antibiotic Disinfection of Biofilms
The first multi-species multi-substrate generalization of the prototype model ( 1) was suggested in [5].In [4] existence results for the solutions were established and numerical simulations presented.The model describes a growing biofilm community and its disinfection by antimicrobial agents.Bacteria in biofilm populations are better protected than free floating cells and behave essentially different under antibiotic treatment.The EPS retards diffusion of antimicrobial agents into the biofilm, cells in the outer layers are attacked first while bacteria in the inner cores are well protected and continue to grow.
The dependent model variables are:

S nutrient concentration B concentration of the antimicrobial agent X volume fraction occupied by active biomass Y volume fraction occupied by inert biomass
The dissolved nutrient S controls the growth of the biomass, and the antimicrobial agent B regulates the disinfection process.As previously, the EPS is implicitly taken into account, and the total biomass fraction M := X +Y is normalized with respect to the maximum bound for the cell density.In dimensionless form the model is represented by the parabolic system with non-negative and bounded initial and boundary data where we use the same notations as in (1).Apart from the diffusion of the dissolved substrates and the death, growth and spatial spreading of biomass the disinfection mechanism is included in the model.During this process antibiotic agents are consumed and active biomass is directly converted into inert biomass, which is determined by the disinfection parameters ζ 1 and ζ 2 .Like in the mono-species model, the production of active biomass due to the consumption of nutrients is described by Monod reaction functions.In the absence of antimicrobial agents and inert biomass, the model reduces to the single species biofilm growth model (1).
In [4] numerical simulations were presented to illustrate the model behavior and to analyze the efficiency of different disinfection strategies.The numerical experiments show that cells in the outer layers of the biofilm are attacked first while cells in the inner scores remain protected and survive longer.

B. Amensalistic Biofilm Control System
The model of an amensalistic biofilm control system [13] extends the single-species probiotic model [14] and possesses a similar structure as the model of antibiotic disinfection.In [13] existence results for the solutions were established and numerical simulations presented.
The model describes how a beneficial biofilm controls the growth of a pathogenic biofilm community by alternating the environmental conditions.The probiotic biofilm modifies the local concentration of protonated lactic acids, which decreases the pH concentration and deteriorates the growth conditions for the pathogens, while the controlling bacteria are more tolerant to these changes.The dependent model variables are: C concentration of protonated lactic acids P concentration of hydrogen ions X volume fraction occupied by pathogens Y volume fraction occupied by probiotics Z volume fraction occupied by inert biomass In dimensionless form the model is represented by the parabolic system with non-negative and bounded initial and boundary data where we use the same notations as in (1).The constants d C , d P , α 1 , α 2 , α Inert probiotics and pathogens are not distinguished in the model.As previously, the EPS is implicitly taken into account, and the total biomass fraction M := X +Y +Z is normalized with respect to the maximum bound for the cell density.Protonated lactic acids C are produced by both bacterial species until a saturation level is reached.
The hydrogen ion concentration P is related to the local pH value by pH = − log P. It increases, facilitated by the protonated lactic acids, until a threshold value is archived.
The growth and inhibition functions ψ 1 and ψ 2 are piecewise linear such that they are positive if C and P are small, and negative if C or P becomes large.Between the growth and inhibition range there is an extended neutral range.More precisely, the functions ψ i are given by where h i 1 and h i 2 are defined as and the constants ζ 1 j and ζ 2 j , j = 1, . . ., 4, are positive with For the probiotic strategy to be effective we require that ζ 2 j ≥ ζ 1 j , j = 1, . . ., 4 ( [13]).
The mechanism of probiotic control is different from traditional antibiotic control strategies of biofilms, where the inner layers of the film are protected by the outer layers and the antibiotics fail to fully penetrate the biofilm.The numerical experiments for the probiotic biofilm model in [13] show that pathogens in the core of the biofilm, close to the substratum, are eradicated first.

C. Quorum-Sensing in Patchy Biofilm Communities
A model for quorum-sensing in growing biofilm communities was proposed and studied by numerical simulations in [11].It extends the prototype biofilm growth model ( 1) and combines it with the model for quorumsensing in planktonic cultures in [16].Analytical aspects of the quorum-sensing model [11] were addressed in [21].
Quorum-sensing is a cell-cell communication mechanism used by bacteria to coordinate gene expression and behavior in groups.Bacteria constantly produce low amounts of signaling molecules that are released into the environment.Accumulation of autoinducers triggers a response by the cells and since the producing cells respond to their own signals the molecules are also called autoinducers ( [16], [12]).When the concentration of autoinducers locally passes a certain threshold, the cells are rapidly induced, and switch from a so-called downregulated to an up-regulated state.In an up-regulated state they typically produce the signaling molecule at a highly increased rate ( [11]).
The dependent model variables are: S concentration of growth controlling substrate A concentration of autoinducers X volume fraction occupied by down-regulated cells Y volume fraction occupied by up-regulated cells In dimensionless form the model is represented by the parabolic system with non-negative and bounded initial and boundary data where we use the same notations as in (1) and |•| denotes the absolute value.The constants d A and γ are positive, m ≥ 1, and α, β and k 5 are non-negative.Moreover, we require that α + β > γ.
Apart from the constants in the prototype model (1) the constants in (4) have the following meaning: The total biomass density M = X + Y is normalized with respect to the maximum bound for the cell density and the EPS is implicitly taken into account.Assuming that induction switches the cells between down-and upregulated states without changing their growth behavior we can assume that the spatial spreading of both biomass fractions is described by the same diffusion operator.The biomass motility constant d is small compared to the diffusion coefficients d S and d A of the dissolved substrates.Like in the mono-species biofilm growth model, biomass is produced due to the consumption of nutrients, which is described by the Monod reaction functions in the equations for the biomass fractions X and Y .Natural cell death is included in the model and determined by the lysis rate k 4 .
If we do not distinguish between down-regulated and up-regulated cells in the model ( 4), we recover the prototype biofilm growth model (1) for the total biomass M = X + Y and the growth controlling nutrient S.
The autoinducer concentration A is normalized with respect to the threshold concentration for induction.Down-regulated cells produce the signaling molecule at rate α, while up-regulated cells produce it at the increased rate α + β, where β is one order of magnitude larger than α.Due to an increase of the autoinducer concentration A, down-regulated cells are converted into up-regulated cells at rate k 5 A m .In applications, typical values for the degree of polymerization are 2 < m < 3 ( [11], [21]).Up-regulated cells are converted back into down-regulated cells at constant rate k 5 .If the molecule concentration A < 1 the latter effect dominates, if A > 1 up-regulation is super-linear.Moreover, abiotic decay of signaling molecules is taken into account in the model and determined by the constant rate γ.
The numerical simulations for the model in [21] indicate that quorum-sensing in spatially structured biofilm populations does not only depend on the local cell density of the population but also on mass transfer effects.Namely, the location of the cell colonies relative to each other and on the prescribed boundary conditions for the substrates.

A. Prototype Biofilm Growth Model
A solution theory for the prototype model (1) was developed in [10].In particular, the well-posedness was established and the existence of a compact global attractor was shown.
We recall the main results in [10] regarding the wellposedness of the model (1).The following theorem yields the existence and regularity results for the solutions (Theorem 3.1, [10]).
Theorem 1.We assume the initial data satisfies where the function Then, there exists a solution (S, M ) satisfying System (1) in the sense of distributions, and the solution belongs to the class Moreover, it was shown that the solutions are L 1 (Ω)-Lipschitz continuous with respect to initial data, which implies its uniqueness.The following result recalls Theorem 3.2 in [10].Proposition 1.Let (S, M ) and ( S, M ) be two solutions of System (1) corresponding to initial data (S 0 , M 0 ), ( S 0 , M 0 ) respectively, and the initial data satisfy the assumptions of the previous theorem.Then, the following estimate holds for t ≥ 0 and some constant c ≥ 0.
We shortly indicate the main ideas of the proofs, for all details and further results we refer to [10].The solutions for the degenerate system (1) are obtained as limits of the solutions of smooth regular approximations.More precisely, the non-degenerate auxiliary systems for the single-species model are given by where the diffusion coefficient in the equation for the biomass fraction in (1) is replaced by the regularized function for small > 0.
First, the non-negativity of the approximate solutions (S , M ) is shown and uniform a-priori L ∞ (Ω)-bounds are established for all sufficiently small > 0. This is archived by comparison principles and the construction of appropriate barrier functions (Proposition 1, [10]).Once a-priori bounds for the solutions are known the existence of solutions of the regular auxiliary systems (5) follows by standard arguments ( [15]).
We further remark that the fast diffusion effect prevents the biomass density from attaining the singular value.More precisely, if the initial data satisfies M 0 L ∞ (Ω) = 1 − δ for some 0 < δ < 1, then, there exists 0 < η < 1 such that the solutions of the nondegenerate approximations (5) satisfy for all sufficiently small > 0 (Proposition 6, [10]).
Having established uniform a-priori bounds we can pass to the limit tends to zero in the distributional formulation of the equation for the nutrient concentration in (5) and obtain a solution S of the original degenerate problem (1).In the governing equation for the biomass density the passage to the limit in the reaction terms is also immediate, the difficulty is to justify the limit for the biomass diffusion terms (see the proof of Theorem 3.1, [10]).
The Lipschitz continuity with respect to initial data in L 1 (Ω)-norm can be deduced from Kato's inequality and is shown in Proposition 2 and Theorem 3.2, [10].

B. Multicomponent Biofilm Models
The pattern of the multi-component biofilm models is essentially different from the prototype model, the equations are strongly coupled through the diffusion operators and the analytical results for the single-species model could not all be carried over.In [4] and [13] the behavior of solutions was studied in numerical simulations and the existence of solutions was established, but the question of uniqueness of solutions remained unanswered in both cases.The first uniqueness result for multispecies models was obtained in [21].
1) Antibiotics and Probiotics Model: The following existence result for solutions of the antibiotics model (2) was shown in [4] (Theorem 2.3).
Theorem 2. We assume the functions B r and S r are nonnegative and belong to the class L ∞ (∂Ω).Moreover, if the initial data X 0 , Y 0 , S 0 , B 0 are non-negative, belong to L ∞ (Ω) and satisfy then, there exists a global solution of the antibiotics model, the functions S, B, X and Y belong to the space L ∞ (Ω × (0, ∞)), are non-negative and satisfy System (2) in distributional sense.
It was shown in [4] that the limit of solutions of non-degenerate approximations yields a solution of the degenerate problem (2).The smooth regular approximations for the antibiotics model are obtained from the system of equations ( 2) by replacing the diffusion coefficient D(M ) in the equations for the biomass components by the regularized function D (M ), The non-negativity and uniform boundedness of the solutions (S , B , X , Y ) of the auxiliary systems can be deduced from comparison principles and by constructing suitable barrier functions.Once a-priori L ∞ (Ω)-bounds are known, the existence of solutions of the regular approximations follows by standard arguments ( [15]).
Using the uniform a-priori bounds for the solutions of the non-degenerate approximations we can pass to the limit tends to zero in the distributional formulation of the equations for the dissolved substrates.The difficulty is to justify the limit for the diffusion terms in the governing equations for the biomass fractions.Since the equations are strongly coupled, the proof requires different arguments than the ones applied for the single-species model in [10].For all details and the proof of Theorem 2 we refer to [4].
The structure of the probiotics model ( 3) is similar to the structure of the antibiotics model (2) and the existence of solutions was shown by similar arguments (Theorem 3.3, [13]).Theorem 3. We assume the functions C r and P r belong to the class L ∞ (∂Ω) and satisfy Moreover, if the initial data C 0 , P 0 , X 0 , Y 0 , Z 0 are nonnegative, belong to L ∞ (Ω) and satisfy then, there exists a global solution of the probiotics model, the functions C, P, X, Y and Z belong to L ∞ (Ω × (0, ∞)), are non-negative and satisfy System (3) in distributional sense.
For the proof of Theorem 3 we refer to [13].Further analytical results were not obtained for the models (2) and (3).In particular, the question of uniqueness of solutions remained unanswered in [4] and [13].
2) Quorum-Sensing Model: A different approach was developed for the quorum-sensing model ( 4) in [21], which led to a uniqueness result for the solutions.The following theorem states the well-posedness of the model (Theorem 3.5 and Theorem 3.11, [21]).
Then, there exists a unique global solution of the quorum-sensing model (4), the functions A, S, X and Y are non-negative and satisfy System (4) in distributional sense.
The particularity of the quorum-sensing model ( 4) is that adding the equations for the biomass fractions X and Y we recover exactly the prototype biofilm growth model (1) for the total biomass M = X + Y and the growth limiting nutrient S. Hence, we may regard M and S as known functions, and the problem reduces to show the well-posedness of the semi-linear degenerate problem for the biomass fraction X and the signaling molecule concentration A, Like for the antibiotics and probiotics model the existence of solutions was established by considering non-degenerate approximations, where the diffusion coefficient D(M ) in the equations for the biomass fractions in (4) is replaced by the regularized function D (M ), Moreover, using the results obtained for the solutions of the single-species model in [10], further regularity results for the solutions could be established, which led to a uniqueness result.For further details and the proofs we refer to [21].

V. CONCLUDING REMARKS
The existence results in the previous section are formulated assuming homogeneous Dirichlet boundary conditions for the biomass components.This situation resembles growing biofilms without substratum, which are commonly called microbial flocs.Such bacterial aggregates enclosed in an EPS matrix are used in the industry for wastewater treatment and also occur in natural settings ( [18]).Boundary conditions of mixed type are, however, often more appropriate in applications.Typically, Dirichlet conditions are specified on some part of the boundary, while Neumann or Robin conditions are imposed on the other parts.In particular, the substratum, on which the biofilm grows is impermeable for all dependent variables, which is described by homogeneous Neumann boundary values.All existence proofs of the previous section carry over to these more general situations as long as homogeneous Dirichlet conditions are imposed for the biomass fractions on one part of the boundary (for details see Theorem 4.1, [10]).
On the other hand, if homogeneous Neumann conditions are assumed for all biomass components and constant Dirichlet conditions for the nutrient concentration, which reflects the situation that no biomass can leave the system and nutrients are constantly added, it was shown that the biomass density reaches the singular value in finite time (Proposition 7 and Proposition 8, [10]).
We finally remark that the solution theory for the single-species model was also extended to less regular initial data.It suffices to assume that the total biomass initially fulfils M 0 ∈ L ∞ (Ω) and 0 ≤ M 0 ≤ 1.If the biomass components satisfy homogeneous Dirichlet conditions on one part of the boundary, the biomass spreading mechanism is strong enough to prevent the solution from attaining the singular value (see Theorem 3.5, [10]).

d A diffusion coefficient of autoinducers k 5
up-regulation rate α autoinducer production rate of down-regulated cells β increased autoinducer production rate of up-regulated cells γ abiotic decay rate of autoinducers m polymerization exponent The additional constants ζ 1 , ζ 2 and d B in (2) are positive and have the following meaning: d B diffusion coefficient of antibiotics ζ 1 antibiotics consumption rate ζ 2 inert biomass production rate