Qualitative analysis of the moving boundary problem for a biofilm reactor model

Accepted Manuscript Qualitative analysis of the moving boundary problem for a biofilm reactor model

B. D’Acunto, L. Frunzo, M.R. Mattei

PII: DOI: Reference:

S0022-247X(16)00132-3 http://dx.doi.org/10.1016/j.jmaa.2016.02.008 YJMAA 20181

To appear in:

Journal of Mathematical Analysis and Applications

Received date:

24 July 2015

Please cite this article in press as: B. D’Acunto et al., Qualitative analysis of the moving boundary problem for a biofilm reactor model, J. Math. Anal. Appl. (2016), http://dx.doi.org/10.1016/j.jmaa.2016.02.008

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

J Math Anal App Journal of Mathematical Analysis and Applications 00 (2016) 1–16

Qualitative analysis of the moving boundary problem for a biofilm reactor model B. D’Acuntoa , L. Frunzoa,∗, M.R. Matteia a Department

of Mathematics and Applications ”Renato Caccioppoli”, University of Naples Federico II, Complesso Monte Sant’Angelo, 80124 Naples, Italy

Abstract The work presents the qualitative analysis of the free boundary value problem related to a biofilm reactor model. In the framework of continuum approach to mathematical modelling of biofilm growth, the problem consists of a system of nonlinear hyperbolic partial differential equations governing the microbial species growth, a system of semilinear elliptic partial differential equations describing the substrate trends, and a system of nonlinear differential equations for the mass balance in the reactor. The free boundary evolution is governed by a differential equation that also accounts for detachment and attachment. The main result is a uniqueness and existence theorem. By using the method of characteristics, the original differential system is converted to Volterra integral equations and then the fixed point theorem is used. The work is completed with numerical simulations describing the free boundary evolution and mainly focused on attachment/detachment process. Keywords: Multispecies biofilms, Biofilm reactor, Free boundary value problems, Nonlinear hyperbolic partial differential equations.

1. Introduction The term biofilm refers to a form of microbial ecosystem constituted by accumulations of microorganisms, associated to a solid surface or phase inter-phase and embedded in a self-produced primarily polysaccharide matrix [1]. The importance and attractiveness of biofilm systems is undisputed. Biofilms are found in extremely varied environments and are characterized by an inherent robustness and flexibility which have been exploited for many applications, such as bioremediation, water purification and wastewater treatment, biofuels and electricity production. Nearly all biofilm communities in nature comprise a variety of microbial species, which function as a cooperative consortium in a relatively complex and coordinate fashion [2]. Generally, the metabolic activities of microorganisms constituting the biofilm determine the formation of various ecological niches due to the local variations in the biofilm environment [4]. Each species fills in a specific ecological niche depending on its metabolism and morphology [5]. It is precisely the presence of multiple species within the biofilm that makes it particularly suitable for the treatment of wastewaters, usually diverse in composition and fluctuate in component concentration. Biofilms have long been considered as a valid alternative technology for wastewater treatment. Over the last decades, their application has considerably increased due to the need of meeting the more stringent requirements for clean water sources and mainly to the significant expansion in the understanding of biofilm properties and formation ∗ Corresponding

author Email addresses: [email protected] (B. D’Acunto), [email protected] (L. Frunzo), [email protected] (M.R. Mattei)

1

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

2

mechanism [6, 7]. Biofilm reactors are characterized by some unique properties when compared to suspended growth systems: bacterial cells in biofilms are more resistant to antimicrobial substances and heavy metals[8]; their residence time is independent of the fluid phase flow rate so they are not subject to washout; biofilm reactors operate at higher biomass density, volumetric productivity and performance stability. Biofilm structure, composition and activity are crucial parameters for the success in operating biofilm reactors [9]. For example, biofilm density and thickness directly affect the efficiency of conversion processes as they define the achievable biomass in the reactor. In addition, microbial diversity is directly related to substrate uptake and to the feasibility of treating contemporary in the same reactor various wastewater constituents. Mathematical modelling represents a valuable tool in understanding the biofilm microenvironment and ultimately elucidate the macroscale process performance. In the context of engineering practice, biofilm reactor models should account for the multiple physical, chemical and biological processes occurring at different time scales within the biofilm and should be able to describe the fate of wastewater constituents. Specific attention should be paid to the evaluation of biofilm composition in terms of relative mass proportion due to microbial competition; concentration of the particulate components in the bulk liquid with respect to impact on system performance and sludge production; spatial profiles of any number of dissolved components in the biofilm; removal rates and effluent concentrations of the dissolved components; biofilm development in terms of thickness. Specifically, the mathematical modelling of a biofilm reactor usually takes into account two different compartments: the bulk liquid where the dissolved component concentrations vary according to the inlet and outlet flow and the flux into/from the biofilm, and the biofilm itself which grows and consumes the substrates provided by the bulk liquid. The two compartments are interconnected through a flux [10, 11]. The biofilm growth can be modeled by using a continuum or discrete approach [12, 13, 14, 15]. 1D continuum models have been widely used to predict the whole process dynamics of biofilm reactors as they have been recognized as suitable tools for most engineering applications [16, 17]. Indeed, nowadays many wastewater treatment plant simulators include a mathematical biofilm module typically based on a 1D approach for biofilm representation (see Boltz et al. [16] for a review). In this work, we refer to the mathematical model derived by coupling the reactor mass balance for substrates with a full one-dimensional Wanner & Gujer type biofilm model as presented in [10, 11]. The resulting mathematical problem is constituted by a system of nonlinear hyperbolic partial differential equations which describes biomass dynamics, a system of semi-linear parabolic partial differential equations for substrate distribution in biofilm, a nonlinear ordinary differential equation for biofilm thickness and a system of nonlinear ordinary differential equations for what concerns the bulk liquid phase. All equations are mutually connected and give rise to a free boundary value problem which is essentially hyperbolic. As for the majority of 1D biofilm models, the previous model has been extensively studied through numerical simulations as a rigorous analysis of the dynamic behavior is usually not easily accessible due the complexity of the system, [11]. Indeed, even when the equations governing the reactor dynamics are not considered, there are only few contributions on the analysis of the Wanner & Gujer type biofilm model, e.g. [18, 19]. In this context, a qualitative analysis of the dynamic behavior of the system biofilm-reactor is presented. We assume that the number n of species into the biofilm is arbitrary as well as the number m of substrates. In addition, also the number of dissolved components in the bulk compartment is arbitrary. We use the characteristic-lines to convert all equations, both partial and ordinary differential equations, to integral equations, mainly of Volterra type. Then, the fixed point theorem is used to obtain uniqueness and existence theorem of solutions to the free boundary value problem. Apart from being interesting for itself, this analysis constitutes a discriminator element for the validity of numerical findings, mainly in relation to the uniqueness result. The work is organized as follows. In Section 2 the engineering problem related to the biofilm reactor is introduced and the main variables defined. Then, the governing equations are presented and discussed. In Section 3 the free boundary problem is stated and the integral equations derived. Section 4 is devoted to the uniqueness and existence theorem. In Section 5 we discuss the detachment phenomenon, neglected in the previous sections. A few numerical simulations are developed in Section 6. The numerical simulations show the free boundary profiles, distinguishing between attachment and detachment problem. Moreover, biofilm evolution is described in terms of microbial species volume fractions and substrate trends.

2

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

3

2. Biofilm reactor model A mathematical model intended to describe the biofilm development and substrate consumption in a biofilm reactor is studied. Following the common practice in biofilm modelling, the latter is assumed to be constituted by two compartments: a completely mixed bulk phase, which supplies the substrates necessary for bacterial growth, and the biofilm itself, which is assumed to cover homogeneously the support surfaces. The two compartments exchange solutes by diffusion. The main processes involved in the biofilm reactor are governed by the following nonlinear differential equations where the notations in Table 1 have been used: ∂ ∂Xi + (uXi ) = ρi rM,i (z, t, X, S), i = 1, ..., n, 0 ≤ z ≤ L(t), t > 0, ∂t ∂z

(1)

∂u  = rM,i , 0 < z ≤ L(t), t ≥ 0, ∂z i=1

(2)

˙ L(t) = u(L(t), t) + σa (t) − σd (L(t)), t > 0,

(3)

∂ 2 Sj ∂Sj − Dj = rS,j (z, t, X, S), 0 < z < L(t), j = 1, ..., m, t > 0, ∂t ∂z 2

(4)

∂Sj (L(t), t) + Q(Sjin − Sj∗ (t)), j = 1, ..., m, t > 0. V S˙ j∗ = −ADj ∂z

(5)

n

n Xi (z, t) = ρi fi ρi fi (z, t) u(z, t) m Sj (z, t) rM,i (z, t,X,S) L(t) σa (t) σd (L(t)) Dj rS,j (z, t, X, S) V Sj∗ (t) Sjin A Q

number of species concentration of microorganism i, X = (X1 , ..., Xn ) constant density n volume fraction of microbial species i, i=1 fi = 1 velocity of biomass number of substrates concentration of substrate j, S = (S1 , ..., Sm ) specific growth rate of species i biofilm thickness, free boundary attachment biomass flux from bulk liquid to biofilm detachment biomass flux from biofilm to bulk liquid diffusion coefficient of substrate j production rate of substrate j volume of the reactor concentration of substrate j in the bulk liquid concentration of substrate j in the influent total biofilm area flow rate through the reactor Table 1. Notations

Equation (1) derives from local mass balance. It was first presented in [12] and introduced in the general form above in [20]. Such equation regulates the dynamics of the various types of microbial biomass constituting the biofilm whose spreading has been modeled as an advective transport mechanism. The concentration of the various microbial species varies on time and space according to the bioconversion processes and the advective transport caused by biofilm spreading. Note that equation (1) has been generally used to describe the dynamics of not only the active microbial species but of all the particulate components, including inert biomass and EPS, which constitute the biofilm. In particular, the inert has been usually treated as a fictitious microbial species whose growth derives from the decay 3

4

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

n

and inactivation of the active biomass. Equation (2) is obtained from (1) by summing on i and using i=1 fi = 1. It regulates the expansion velocity of the microbial mass, assumed equal for all the microbial species. Equation (3) describes the free boundary evolution. It derives from global mass balance by taking into account the attachment and detachment flux, [20]. Equations (4) and (5) describe the dynamics of dissolved compounds within the biofilm and the bulk liquid respectively. In particular, Equation (5) follows from [10] and it basically represents the mass balance for substrate concentration in the bulk liquid, which has been modeled as a Continuously Stirred Tank Reactor of volume V continuously fed and withdrawn at the same flow rate Q. This equation accounts also for the influent and outflow fluxes to or from the biofilm and provides as result the time evolution of the substrate concentrations in the bulk liquid. Equation (4) is the commonly used reaction-diffusion equation which governs substrate diffusion from bulk liquid into the biofilm and its further consumption or production due to microbial activity. This equation is usually considered at pseudo-steady state conditions as the time-scale of substrate mass transport and reaction is orders of magnitude shorter than that of biomass spreading [21] −Dj

∂ 2 Sj = rS,j (z, X, S), 0 < z < L(t), j = 1, ..., m, t ≥ 0. ∂z 2

(6)

3. Statement of the free boundary value problem This section is devoted to the following special free boundary value problem, obtained by associating the initial and boundary conditions to the equations introduced in the previous section:  ∂Xi ∂ ∂t + ∂z (uXi ) = ρi rM,i (z, t, X, S), 0 ≤ z ≤ L(t), t > 0, (7) Xi (z, 0) = ϕi (z), i = 1, ..., n, 0 ≤ z ≤ L(0), 





n = i=1 rM,i , 0 < z ≤ L(t), t ≥ 0, u(0, t) = 0, t ≥ 0,  ˙ L(t) = u(L(t), t), t > 0, L(0) = L0 , ∂u ∂z

(8) (9)

∂2S

−Dj ∂z2j = rS,j (z, X, S), 0 < z < L(t), j = 1, ..., m, t ≥ 0, ∂Sj ∗ ∂z (0, t) = 0, Sj (L(t), t) = Sj (t), j = 1, ..., m, t ≥ 0,

(10)

V S˙ j∗ = −ADj ∂zj (L(t), t) + Q(Sjin − Sj∗ (t)), j = 1, ..., m, t > 0, Sj∗ (0) = Sjin , j = 1, ..., m.

(11)

∂S

The functions ϕi (z) in (7)2 represent the initial concentrations of the microbial species. According to the biological problem they are assigned positive. Condition (8)2 means that no biomass flux is assumed at the biofilm support z = 0. Condition (9)2 assigns the initial biofilm thickness. The first boundary condition in (10)2 prescribes zero substrate flux at the support z = 0. The second condition assigns substrate values on the moving boundary z = L(t). Note that the solution Sj∗ (t) of equation (11)1 is used as boundary condition for equation (10)1 and incorporates in the model the influence that the bulk liquid phase exerts on the biofilm and viceversa. Condition (11)2 assumes the concentration of substrates in the influent as initial condition. In the case of microbial products the inlet concentration is set to zero. Detachment flux, neglected in equation (9)1 , will be studied in Section 5. In the following, the differential problem (7)-(11) will be converted to an integral system (21), (25), (27)-(30) to be used for the uniqueness and existence theorem. Consider the characteristic-like lines z = c(z0 , t) defined by ∂c (z , t) = u(c(z0 , t), t), c(z0 , 0) = z0 , z0 ∈ [0, L0 ], t > 0. ∂t 0

(12)

Two results follow immediately from (12). First, L(t) = c(L0 , t), 4

(13)

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

5

because of (9). In addition, consider the characteristic for z0 = 0 ∂c (0, t) = u(c(0, t), t), c(0, 0) = 0. ∂t Since u(0, t) = 0, as stated in (8)2 , it is apparent that c(0, t) = 0.

(14)

d Xi (c(z0 , t), t) = ρi rM,i (c(z0 , t), t, X(c(z0 , t), t), S(c(z0 , t), t)) dt −Xi (c(z0 , t), t)g(c(z0 , t), t, X(c(z0 , t), t), S(c(z0 , t), t)), z0 ∈ [0, L0 ],

(15)

Use (12) in (7)1 and consider (8)2

where g=

n 

(16)

rM,i .

i=1

The following initial conditions are associated to (15) Xi (c(z0 , 0), 0) = Xi (z0 , 0) = ϕi (z0 ), z0 ∈ [0, L0 ],

(17)

because of (7)2 and (12). Now, let us introduce the following positions x(z0 , t) = X(c(z0 , t), t), s(z0 , t) = S(c(z0 , t), t),

(18)

x˙ i (z0 , t) = Fi (c(z0 , t), t, x(z0 , t), s(z0 , t)), z0 ∈ [0, L0 ],

(19)

Fi = ρi rM,i − xi g.

(20)

and rewrite (15) as follows

where Integrating (19) yields the following system of integral equations  t xi (z0 , t) = ϕi (z0 ) + Fi (c(z0 , τ ), τ, x(z0 , τ ), s(z0 , τ ))d τ, z0 ∈ [0, L0 ],

(21)

0

which incorporate the initial conditions (17). Let us discuss (10). Integrating (10)1 and considering (10)2 leads to the following integral equations Sj (z, t) = Dj−1





L

η

dη 0

z

rS,j (ζ, X(ζ, t), S(ζ, t))dζ + Sj∗ (t), 0 < z < L(t),

(22)

which are equivalent to problem (10). When z = c(z0 , t), equation (22) are rewritten as Sj (c(z0 , t), t) = Dj−1





c(L0 ,t)

η

dη c(z0 ,t)

0

rS,j (ζ, X(ζ, t), S(ζ, t))dζ + Sj∗ (t), 0 < z0 < L0 ,

(23)

where (13) was used. Consider the change of variables η = c(η0 , t), ζ = c(ζ0 , t), 0 < η0 , ζ0 < L0 ,

(24)

in (23) and obtain sj (z0 , t) = Dj−1



L0 z0

∂c (η , t)dη0 ∂η0 0



η0 0

rS,j (c(ζ0 , t), x(ζ0 , t), s(ζ0 , t)) 5

∂c (ζ , t)dζ0 ∂ζ0 0

6

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

+Sj∗ (t), 0 < z0 < L0 ,

(25)

where notation (18)2 was used. From (22) it follows that ∂Sj (z, t) = −Dj−1 ∂z and therefore



z 0

rS,j (ζ, X(ζ, t), S(ζ, t))d ζ,

 c(z0 ,t) ∂Sj (c(z0 , t), t) = −Dj−1 rS,j (ζ, X(ζ, t), S(ζ, t))d ζ ∂z 0  z0 dc = −Dj−1 rS,j (c(ζ0 , t), x(ζ0 , t), s(ζ0 , t)) dζ0 dζ 0 0

because of (24). Hence, ∂Sj ∂Sj (L(t), t) = (c(L0 , t), t) = −Dj−1 ∂z ∂z



L0 0

rS,j (c(ζ0 , t), x(ζ0 , t), s(ζ0 , t))

dc dζ . dζ0 0

(26)

Considering (26) in (24)1 yields S˙ j∗ (t) = (A/V )



L0 0

rS,j (c(ζ0 , t), x(ζ0 , t), s(ζ0 , t))

dc dζ + (Q/V )(Sjin − Sj∗ (t)), dζ0 0

which explicitly shows the dependence of the substrate concentrations in the bulk liquid from the production and consumption rates characterizing the biofilm activity. Note also that the substrate dynamics in the bulk liquid are affected by the operational parameters of the biofilm reactor, such as the total biofilm area and the dilution rate which corresponds to the ratio Q/V . Integrating the last equation with respect to time yields Sj∗ (t) =



t

0

 exp(−Q(t − τ )/V )dτ

L0 0

rS,j (c(ζ0 , τ ), x(ζ0 , τ ), s(ζ0 , τ ))

dc dζ dζ0 0

+Sjin , j = 1, ..., m, t > 0.

(27)

Consider equation (8)1 rewritten by using notations (16), (18) ∂c ∂u(c(z0 , t), t) = g(c(z0 , t), t, x(z0 , t)), s(z0 , t)) (z , t). ∂z0 ∂z0 0 Hence,

 u(c(z0 , t), t) =

z0 0

g(c(ζ0 , t), t, x(ζ0 , t)), s(ζ0 , t))

(28)

∂c (ζ , t)dζ0 , ∂ζ0 0

where (8)2 and (14) have been considered. Substitute the result above in equation (12)1 and integrate with respect to time with initial condition (12)2  t  z0 ∂c c(z0 , t) = z0 + dτ g(c(ζ0 , τ ), τ, x(ζ0 , τ )), s(ζ0 , τ )) (ζ0 , τ )dζ0 . (29) ∂ζ 0 0 0 In addition, ∂c(z0 , t) =1+ ∂z0



t 0

g(c(z0 , τ ), τ, x(z0 , τ )), s(z0 , τ ))

∂c (z , τ )dτ, ∂z0 0

(30)

which follows easily from (29). The last equation has been introduced as ∂c/∂z0 constitutes an additional unknown function which is involved in all the integral equations. The integral system (21),(25), (27)-(30) will be discussed in the next section. Note that the equation for the biofilm thickness L is included as specified in (13). 6

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

7

4. Uniqueness and existence theorem In order to obtain more compact equations, let us introduce the vector x∗ as follows ∗ , s), x∗ = (x, c, ∂c/∂z0 , S1∗ , ..., Sm

so that

(31)

x∗i = xi , i = 1, ..., n, x∗i+1 = c, x∗i+2 = ∂c/∂z0 , x∗i+2+j = Sj∗ , x∗i+2+m+j = sj , j = 1, ..., m.

Therefore, equation (21) writes x∗i (z0 , t) = ϕi (z0 ) +



t 0

Fi (τ, x∗ (z0 , τ ))d τ, i = 1, ..., n, z0 ∈ [0, L0 ].

(32)

Equation (29) writes x∗n+1 (z0 , t)

 = z0 +



t

z0

dτ 0

0

Fn+1 (τ, x∗ (ζ0 , τ ))dζ0 , z0 ∈ [0, L0 ],

where Fn+1 (τ, x∗ (ζ0 , τ )) = g(c(ζ0 , τ ), τ, x(ζ0 , τ )), s(ζ0 , τ )) Equation (30) writes x∗n+2 (z0 , t) = 1 +



t 0

∂c (ζ , τ ). ∂ζ0 0

Fn+2 (τ, x∗ (z0 , τ ))dτ, z0 ∈ [0, L0 ],

where Fn+2 (τ, x∗ (z0 , τ )) = g(c(z0 , τ ), τ, x(z0 , τ )), s(z0 , τ ))

(33)

(34)

∂c (z , τ ). ∂z0 0

Equation (27) writes x∗n+2+j (t) = Sjin +





t

L0

dτ 0

0

Fn+2+j (t, τ, x∗ (ζ0 , t))dζ0 , j = 1, ..., m,

where Fn+2+j (t, τ, x∗ (ζ0 , t)) = e−Q(t−τ )/V rS,j (c(ζ0 , t), x(ζ0 , t), s(ζ0 , t)) Equation (25) writes  +

z0

dη0

∂c (ζ , τ ). ∂ζ0 0

x∗n+2+m+j (z0 , t) = x∗n+2+j (t)



L0

(35)

η0 0

Fn+2+m+j (x∗ (ζ0 , t))dζ0 , j = 1, ..., m, 0 < z0 < L0 .

where Fn+2+m+j (x∗ (ζ0 , t)) = Dj−1

∂c ∂c rS,j (x∗n+1 (ζ0 , t), x(ζ0 , t), s(ζ0 , t)) . ∂η0 ∂ζ0

Now, we can show the following theorem. Theorem 1. Suppose that: (i) x∗k (z0 , t) ∈ C 0 ([0, L0 ] × [0, T1 ]), L0 > 0, T1 > 0, k = 1, ..., n + 2 + 2m; (ii) ϕi (z0 ) ∈ C 0 ([0, L0 ]), i = 1, ..., n; (iii) Sjin = positive constant, j = 1, ..., m; (iv) |x∗i − ϕi | ≤ hi , i = 1, ..., n; |x∗n+1 − z0 | ≤ hn+1 ; 1 ≤ x∗n+2 ≤ 1 + hn+2 ; |x∗n+2+j − Sjin | ≤ hn+2+j , j = 1, ..., 2m, where hk = constant > 0; 7

(36)

8

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

(v) the function g, defined in (16), is essentially positive; (vi) the initial biofilm thickness is relatively small; (vii) Fk continuous and bounded Mk = max |Fk |, k = 1, ..., n + 2 + 2m, when (z0 , t) ∈ [0, L0 ] × [0, T1 ] and the functions x∗k satisfy the assumptions (i)-(iv). Then, integral system (32)-(36) has a unique solution x∗k ∈ C 0 ([0, L0 ] × [0, T ]), k = 1, ..., n + 2 + 2m, where   h1 hn+1 hn+2 hn+3 hn+2+m , ..., , , , ..., . T = min T1 , M1 L0 Mn+1 Mn+2 L0 Mn+3 L0 Mn+2+m P ROOF OF T HEOREM 1. Denote by Ω the space of continuous functions x∗i (z0 , t), z0 ∈ [0, L0 ], t ∈ [0, T ], and let hn+2+m+j = hn+2+j + L20 Mn+2+m+j , j = 1, ..., m. Introduce the norm ||x∗ || =

n+2+2m 

max exp(−γt)|x∗k |, Ω

k=1

where γ is a positive constant that will be fixed later on. Assume Fk Lipschitz continuous with respect to x∗k , k = 1, ...n + 2 + 2m, ˜ ∗ )| ≤ λi |Fi (τ, t, x∗ ) − Fi (τ, t, x

n+2+2m 

|x∗k − x ˜∗k |, i = 1, ...n + 2 + 2m.

k=1

Consider the map y = Ax∗ , yk = right hand side of equations (32)-(36) and assume g ≥ 0. Hypothesis (v), jointly with x∗n+2 ≥ 1, implies Fn+2 ≥ 0. Let us prove that A maps Ω into itself. Indeed, |yi − ϕi | ≤ Mi T ≤ hi , i = 1, ..., n, |yn+1 − z0 | ≤ Mn+1 T L0 ≤ hn+1 , 1 ≤ yn+2 ≤ 1 + Mn+2 T ≤ 1 + hn+2 , |yn+2+j − Sjin | ≤ Mn+2+j T L0 ≤ hn+2+j , j = 1, ..., m, |yn+2+m+j − Sjin | ≤ |yn+2+j − Sjin | + Mn+2+m+j L20 ≤ hn+2+m+j , j = 1, ..., m. ˜ ∗ ∈ Ω and let y ˜ = A˜ Consider x x∗ . We easily obtain ˜ ∗ ||, i = 1, ..., n, |yi − y˜i | exp(−γt) ≤ (λi /γ)||x∗ − x ˜ ∗ ||, |yn+1 − y˜n+1 | exp(−γt) ≤ (λn+1 L0 /γ)||x∗ − x ˜ ∗ ||, |yn+2 − y˜n+2 | exp(−γt) ≤ (λn+2 /γ)||x∗ − x ˜ ∗ ||, i = n + 3, ..., n + 2 + m, |yi − y˜i | exp(−γt) ≤ (λi /γ)||x∗ − x ˜ ∗ ||, i = n + 2 + m + 1, ..., n + 2 + 2m. |yi − y˜i | exp(−γt) ≤ λi L20 ||x∗ − x Hence,

˜ ∗ || ≤ Λ||x∗ − x ˜ ∗ ||, ||y∗ − y

where 1 Λ1 = γ



Λ = Λ1 + Λ2 , n 

λi + λn+1 L0 + λn+2 +

n+2+m 

i=1

i=n+3

8

λi

, Λ2 = L20

n+2+2m  i=n+2+m+1

λi .

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

9

Note that the positive constant γ was not specified earlier. Therefore, it can be chosen very large so that Λ1 becomes very small, say < ε, ∀ε > 0. In addition, if the parameter L0 is small enough such that  L20

≤ (1 − ε)

n+2+2m 

−1 λi

,

i=n+2+m+1

then Λ2 ≤ 1 − ε. Overall this leads to Λ < 1. Remark 1. Let us discuss the hypotheses of Theorem 1. The constants Sjin , j = 1, ..., m, represent the initial concentrations in the bioreactor and so they are essentially positive. The function g is essentially positive according to [12, 20]. Finally, the initial biofilm thickness L0 used in these models is very very small [10, 11]. It should be noted that the theorem provides the functions xi (z0 , t), i = 1, ..., n, and sj (z0 , t), j = 1, ..., m. On the other hand the original problem required the the functions Xi (z, t) and Sj (z, t). This question is discussed in the following corollary. Corollary 1. Under the same hypotheses as Theorem 1, Xi , i = 1, ..., n and Sj , j = 1, ..., m can be determined as functions of z, t. P ROOF OF C OROLLARY 1. The function z = c(z0 , t) can be inverted and gives z0 = c−1 (z, t). Therefore, xi (z0 , t) = xi (c−1 (z, t), t) = Xi (z, t), sj (z0 , t) = sj (c−1 (z, t), t) = Sj (z, t). 5. Detachment and attachment The exchange of particulate species between biofilm and bulk liquid, namely the attachment and detachment flux, plays a crucial role in biofilm technology as it contributes to defining the microbial ecology of biofilm ecosystems. In this work, the attachment process has been considered as the mechanical aggregation of cell clusters on biofilm surface. The attachment as diffusion of species from the bulk liquid to the biofilm requires a consistent modification of the original model and it has been recently discussed in [22]. The modelling approach introduced in [22] allows us to hold the model hyperbolic and so, no boundary condition is needed for the microbial biomass on the free boundary. In this contribution, the attachment process is analytically described by the positive function σa (t). It plays a relevant role mostly in the initial phase of biofilm formation and it is usually neglected for mature biofilms, where the detachment flux is prevalent. Moreover, the mathematical modelling of the attachment process leads to a free boundary value problem which is completely different from the one discussed in Sec. 3 and will be mentioned at the end of this section. Therefore, we focus our attention on the detachment process which was neglected in the free boundary value problem stated in Sec. 3. Considering the detachment flux yields the following equation for the free boundary  ˙ L(t) = u(L(t), t) − σd (L(t)), t > 0, (37) L(0) = L0 . As shown in the previous sections, the method of characteristics has been effectively applied for σd = 0. Therefore, in this section, we propose to transform suitably equation (37) by writing it in terms of characteristics. In this way all reasonings of Sections 3 and 4 can be used. Denote by y0 = y0 (t) the initial point of the characteristic that assumes the value L(t) at the time t, fig. 1, y0 = y0 (t), 0 ≤ y0 (t) ≤ L0 ,

(38)

L(t) = c(y0 (t), t).

(39)

9

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

10

Figure 1. Characteristic from y0 (t)

By using the notations above, equation (37)1 can be rewritten as c(y ˙ 0 (t), t) = u(c(y0 (t), t), t) − σd (c(y0 (t), t)), t > 0.

(40)

Of course, the function y0 (t) is unknown, as L(t). However, by considering the reasonings developed in Sections 3 and 4, it is apparent that finding L(t) is equivalent to finding y0 (t). Consider the equation (12) for characteristics and note that it is equivalent to  t u(c(z0 , τ ), τ )dτ, 0 ≤ z0 ≤ L0 , t > 0. c(z0 , t) = z0 + 0

For z0 = y0 (t) we obtain  c(y0 (t), t) = y0 (t) +

0

t

u(c(y0 (t), τ ), τ )dτ, 0 ≤ y0 (t) ≤ L0 , t > 0.

Differentiating with respect to t yields  c(y ˙ 0 (t), t) = y˙ 0 (t) + u(c(y0 (t), t), t) + y˙ 0 (t)   = y˙ 0 (t) 1 +

t 0

t 0

∂c ∂u (c(y0 (t), τ ), τ ) (y (t), τ )dτ ∂z ∂z0 0

∂u (c(y0 (t), τ ), τ )dτ ∂z0

 + u(c(y0 (t), t), t).

Hence, by considering (30) and (28) c(y ˙ 0 (t), t) = y˙ 0 (t)

∂c (y (t), t) + u(c(y0 (t), t), t), 0 ≤ y0 (t) ≤ L0 , t > 0. ∂z0 0

(41)

Finally, from (40), (41) we derive the desired differential equation for the unknown function y0 (t) y˙ 0 (t)

∂c (y (t), t) = −σd (c(y0 (t), t)), 0 ≤ y0 (t) ≤ L0 , t > 0. ∂z0 0

(42)

Note that the last equation can be solved with respect to y˙ 0 (t) since ∂c/∂z0 ≥ 1. Also, equation (42) suggests efficient numerical methods. 10

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

Consider the attachment process. Equation for free boundary writes  ˙ L(t) = u(L(t), t) + σa (t), t > 0, L(0) = L0 ,

11

(43)

where σa (t) is a known positive function of time. In this situation the free boundary z = L(t) is a space-like line, fig. 2. It is well known that the equation for the initial conditions is replaced by the following, [12],

Figure 2. Free boundary space-like line and a characteristic

Xi (L(t), t) = ψi (t), i = 1, ..., n, t > 0.

(44)

The last equation expresses that the biomass concentrations at biofilm boundary L(t) are the same as the bulk liquid. Condition (44) is mainly used to describe the initial phase of biofilm formation. In addition, the free boundary that arises is completely different from the ones discussed previously. It can be analyzed as in [23, 24]. 6. Numerical simulations This section is devoted to numerical simulations for a specific biofilm reactor system described in details in the following. The numerical analysis has been developed by using the method of characteristics and an original software has been properly set-up. The main objective is to investigate the effect that attachment and detachment processes exert separately or contextually on biofilm system evolution, which has been evaluated in terms of biofilm thickness, microbial species distribution and substrate trends. The primary parameters that have been varied in the numerical experiments are the attachment σa (t) and the detachment flux σd (L(t)), as summarized in Table 2. In particular, the first simulation experiment examines the dynamics of a biofilm reactor system subjected to a continuous and constant on time σa (t). The second numerical experiment investigates the case of a biofilm system continuously eroded at a variable rate σd (L(t)) = λL2 (t), where λ denotes the erosion parameter. In the last numerical simulation the concomitant effects exerted by the attachment and detachment processes on biofilm dynamics have been evaluated. Simulation 1 Simulation 2 Simulation 3

Attachment + +

Detachment + +

Table 2. Overview of the different numerical simulations performed in this study.

All numerical simulations presented in this section are related to a microbial symbiosis for methane production. Two microbial species are considered: fermentative bacteria X1 using organic acid as substrate, and hydrogenotrophic methanogens X2 growing on dioxide carbon and hydrogen. The decay of these microbial 11

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

12

species generates inert material, which is treated as an additional particulate component X3 . Four reacting dissolved components are simultaneously considered: organic acid S1 , which is continuously provided from the bulk liquid, acetate S2 , dioxide carbon S3 and hydrogen S4 , which on turns are produced within the biofilm by fermentative bacteria and subsequently consumed in part by methanogenic bacteria. For the specific biological case analyzed in this work, the net specific biomass growth rates in equations (1) are given by S1 f − kd,1 f1 , kS,1 + S1 1

(45)

S3 S4 f2 − kd,2 f2 , kS,3 + S3 kS,4 + S4

(46)

rM,1 = μmax,1 rM,2 = μmax,2

rM,3 = kd,1 f1 + kd,2 f2 ,

(47)

where μmax,i is the maximum specific growth rate for species i, kS,i is the half-saturation constant for substrate Si and kd,i is the decay coefficient for species i. The conversion rates in equations (6) characterizing microbial activity are expressed as rS,1 = −

S1 μmax,1 X , Y1 kS,1 + S1 1

(48)

μmax,1 S1 X , Y1 kS,1 + S1 1

(49)

rS,3 = g3 (1 − Y1 )

μmax,1 μmax,2 S1 S3 S4 X − (β1 − Y2 ) X , Y1 kS,1 + S1 1 Y2 kS,3 + S3 kS,4 + S4 2

(50)

rS,4 = g4 (1 − Y1 )

μmax,1 μmax,2 S1 S3 S4 X − (β2 − Y2 ) X , Y1 kS,1 + S1 1 Y2 kS,3 + S3 kS,4 + S4 2

(51)

rS,2 = g2 (1 − Y1 )

being Yi the yield coefficient of species i, gi and βi stoichiometric coefficients. For all the simulations, the initial biofilm thickness has been set to 2μm. The values of the kinetics, stoichiometric and operational parameters of the biofilm reactor system used for numerical simulations are reported in Table 3. The initial microbial distribution and the substrate concentration in the inlet flow are stated in Table 3 as well. Model outcomes for each simulation experiment have been summarized in Figures 3-5. More precisely, simulation results have been reported for each investigated biofilm reactor system in terms of biofilm thickness evolution over time (Figures 3-5(A1)), microbial distribution (Figures 3-5(A2)) and substrate trends (Figures 3-5(A3)) within biofilm after 15 days of simulation time. Simulation 1 examines the case of a biofilm reactor system subjected to a constant attachment flux over time. The attachment process exerts a positive effect on biofilm growth: biofilm thickness is substantially higher than the one obtained in absence of attachment and due only to biomass expansion in z direction, as shown in Figure 3(A1). However, the two curves show the same behavior. Indeed, in absence of ”negative processes” and according to equation (3), biofilm thickness increases indefinitely over time; for this reason only the first 6 days of simulation time have been reported in Figure 3(A1). For what concerns microbial species distribution and substrate trends, it is possible to note that after 5 days simulation time, the fermentative bacteria predominate all over the biofilm and the methanogens are confined at the bottom of the biofilm (Figure 3(A2)). The organic acid shows a fully penetrated profile; the carbon dioxide and the acetate have a higher concentration in the inner part of the biofilm, while the hydrogen is almost zero all over the biofilm (Figure 3(A3)).

12

13

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16 Parameter Initial Biofilm thickness Initial Volume Fraction of f1 Initial Volume Fraction of f2 Initial Volume Fraction of f3 Initial concentration of S1 Initial concentration of S2 Initial concentration of S3 Initial concentration of S4 Reactor Volume Biofilm surface Inlet and outlet flow rate Inlet organic acid concentration Inlet acetate concentration Inlet carbon dioxide concentration Inlet hydrogen concentration Maximum growth rate of X1 Maximum growth rate of X2 Half saturation constant of X1 for S1 Half saturation constant of X2 for S3 Half saturation constant of X2 for S4 Yield of X1 on organic acid Yield of X2 on organic acid Microbial decay constant of X1 Microbial decay constant of X2 Stoichiometric coefficient for S2 production Stoichiometric coefficient for S3 production Stoichiometric coefficient for S4 production Stoichiometric coefficient Stoichiometric coefficient Erosion parameter Attachment flux

Symbol L0 f1 (z, 0) f2 (z, 0) f3 (z, 0) S1 (z, 0) S2 (z, 0) S3 (z, 0) S4 (z, 0) V S Q S1in S2in S3in S4in μmax,1 μmax,2 KS,1 KS,3 KS,4 Y1 Y2 kd,1 kd,2 g2 g3 g4 β1 β2 λ σa

Unit μm – – – – – – – m3 m2 m3 /d mgl−1 mgl−1 mgl−1 mgl−1 d−1 d−1 mgl−1 mgl−1 mgl−1 gbiomass /gsubstrate gbiomass /gsubstrate d−1 d−1 – – – – – mmd−1

Value 2 0.8 0.2 0.0 0.0 0.0 0.0 0.0 102 10 103 1 0 0 0 0.75 0.35 10−1 10−2 10−3 0.3 0.3 5 · 10−3 5 · 10−3 0.663 0.685 0.045 0.016 0.5 6 · 103 5 · 10−7

Table 3. Parameter values and initial conditions used for numerical simulations

(A1) 6

5

Time [d]

4

3

2

1

0

0

10

20

30

40 50 60 Biofilm thickness [μm]

(A2)

1

Methanogens

0

20

40 60 Space [μm]

80

90

100

(A3)

1.5 1

13

Fermentative

0.5 0

Substrate concentration trends [mgl−1]

Biofilm Volume Fraction

1.5

70

80

100

0.5 0

0

20

40 60 Space [μm]

80

100

Figure 3. Free boundary evolution over time (A1), microbial species distribution (A2) and substrate trends (A3) after 5 days simulation time of a biofilm reactor system experiencing a constant attachment flux. (A1): square line - free boundary; continuous line - characteristic line z = (L0 , t). (A3): continuous line - S1 ; dash dot line - S2 ; dot line - S3 ; circle line - S4 .

14

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

In Simulation 2, the biofilm is continuously detached at a variable rate σd (L(t)), until reaching a steady state. The inclusion of the detachment term results in a shorter biofilm thickness than the one corresponding to the characteristic line z(L0 , t) and in the achievement of a steady-state value of 80 μm after 15 days (Figure 4(A1)). The gap between the two curves is visibly small at the beginning of the simulation and starts to increase going on with time until the biofilm thickness approaches the steady-state condition. This is due to the fact that the detachment rate is a function of L2 (t). After 10 days of simulation time, the curves start to show a divergent behavior: the biofilm thickness reaches the steady-state while the characteristic line z(L0 , t) continues to increase over time. In Figure 4(A2) the microbial distribution within the biofilm after 15 days of simulation time is shown. The fermentative bacteria are found to predominate all over the biofilm, while the methanogens develop mostly in the inner layer of the biofilm where they found more appropriate growth conditions. Substrate profiles (Figure 4(A3)) reflect the microbial composition: the organic acid concentration decreases from the free boundary to the substratum reaching zero in the inner part of the biofilm, while the carbon dioxide and the acetate concentrations are higher in the inner layer of the biofilm. The hydrogen is produced by the fermentative bacteria and successively consumed by the methanogens at the bottom of the biofilm. In the last simulation experiment, both the attachment and detachment fluxes are considered in equation (3). At the beginning of the simulation, the attachment flux is higher than the detachment rate; therefore, the function σ = σa − σd is positive and the overall effect consists in a longer biofilm thickness than the characteristic line z(L0 , t). Going on with simulation time, the detachment rate increases until assuming the same value of the attachment flux. (A1) 15

Time [d]

10

5

0

0

10

20

30

40 50 60 Biofilm thickness [μm]

(A2)

1

Methanogens Fermentative

0.5 0

Substrate concentration trends [mgl−1]

Biofilm Volume Fraction

1.5

0

20

40 60 Space [μm]

80

100

70

80

90

100

(A3)

1.5 1 0.5 0

0

20

40 60 Space [μm]

80

100

Figure 4. Free boundary evolution over time (A1), microbial species distribution (A2) and substrate trends (A3) at steady-state of a biofilm reactor system experiencing a continuous detachment. (A1): square line - free boundary; continuous line - characteristic line z = (L0 , t). (A3): continuous line - S1 ; dash-dot line - S2 ; dot line - S3 ; circle line - S4 .

14

15

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

At this time, the function σ = σa − σd becomes zero and the two curves, free boundary and characteristic line z(L0 , t), intersect. In the following instants, the detachment rate starts to be higher than the attachment flux and the function σ = σa − σd assumes negative values. This leads to a biofilm thickness shorter than the line z(L0 , t) (Figure 5(A1)), and similarly to the previous case, the biofilm approaches a steady-state thickness of 80 μm after 15 days. In addition, as shown in Figure 5(A2,A3), microbial distribution and substrate trends resemble the ones obtained in the previous numerical experiment. Note that the numerical outputs of the last simulation experiment, confirm that the attachment process plays a crucial role mostly in the initial phase of biofilm formation. Conversely, the detachment affects mostly the dynamics of mature biofilms. This result highlights that when both detachment and attachment are simultaneously considered, the prevailing process determines the type of free boundary problem. (A1) 15

Time [d]

10

5

0

0

10

20

30

40 50 60 Biofilm thickness [μm]

(A2)

1

Substrate concentration trends [mgl−1]

Biofilm Volume Fraction

1.5 Methanogens Fermentative

0.5 0

0

20

40 60 Space [μm]

80

100

70

80

90

100

(A3)

1.5 1 0.5 0

0

20

40 60 Space [μm]

80

100

Figure 5. Free boundary evolution over time (A1), microbial species distribution (A2) and substrate trends (A3) at steady-state of a biofilm reactor system subjected to both attachment and detachment processes. (A1): square line - free boundary; continuous line - characteristic line z = (L0 , t). (A3): continuous line - S1 ; dash dotted line - S2 ; dotted line - S3 ; circle line - S4 .

7. Conclusion In this work a qualitative analysis of the free boundary problem related to a biofilm reactor model was presented. The model takes into account the dynamics of both biofilm and bulk liquid compartments which are connected by the substrate flux. The main result was an existence and uniqueness theorem to the solutions of the free boundary value problem. This result is interesting for itself due to the relatively few contributions on the topic in literature. In addition, the 15

D’Acunto et al. / Journal of Mathematical Analysis and Applications 00 (2016) 1–16

16

work provides consistency to the model mainly for the existence of solutions and covers the gap of a solid base for the numerical integration especially for uniqueness result. A valuable extension of the model may be related to the phenomenon of microbial colonization of an already constituted biofilm by new species, [22]. The presented model takes into account the influence that the bulk liquid exerts on the biofilm in terms of substrate flux. By a similar mechanism, the presence and concentration of microbial species in the bulk liquid can affect the biofilm ecology and lead to the penetration of new species. Also in this case, the link between the two compartments is constituted by the flux which should be incorporated in mass balance equations. References [1] H. J. E BERL , D. F. PARKER , M. VAN L OOSDRECHT, A new deterministic spatiotemporal continuum model for biofilm development, Comput. Math. Methods Med., 3 (2001), pp. 161–175. [2] F. V ILLA , F. C APPITELLI, Plant-derived bioactive compounds at sub-lethal concentrations: towards smart biocide-free antibiofilm strategies, Phytochem. Rev., 12 (2013), pp. 245–254. [3] M. E. DAVEY, G. A. O’T OOLE , M. VAN L OOSDRECHT, Microbial biofilms: from ecology to molecular genetics, Microbiol. Mol. Bio. Rev., 64 (2000), pp. 847–867. [4] B. H ALAN , K. B UEHLER , A. S CHMID, Biofilms as living catalysts in continuous chemical syntheses, Trends Biotechnol., 30 (2012), pp. 453– 465. [5] N. Q URESHI , B. A. A NNOUS , T. C. E ZEJI , P. K ARCHER , I. S. M ADDOX, Biofilm reactors for industrial bioconversion processes: employing potential of enhanced reaction rates, Microb. Cell Fact., 4 (2005), pp. 24-45. [6] M. R. M ATTEI , B. D’ACUNTO , G. E SPOSITO , L. F RUNZO , F. P IROZZI, Mathematical modeling of competition and coexistence of sulfatereducing bacteria, acetogens, and methanogens in multispecies biofilms, Desalin. Water Treat., 55 (2015), pp. 740–748. [7] M. R. M ATTEI , L. F RUNZO , B. D’ACUNTO , G. E SPOSITO , F. P IROZZI, Modelling microbial population dynamics in multispecies biofilms including Anammox bacteria, Ecol. Model., 304 (2015), pp. 44–58. [8] F. V ILLA , S. V ILLA , A. G ELAIN , F. C APPITELLI, Sub-lethal activity of small molecules from natural sources and their synthetic derivatives against biofilm forming nosocomial pathogens, Curr. Top. Med. Chem., 13 (2013), pp. 3184–3204. [9] K. C. C HENG , A. D EMIRCI , J. M. C ATCHMARK, Advances in biofilm reactors for production of value-added products, Appl. Microbiol. Biotechnol., 87 (2010), pp. 445–456. [10] A. M A Sˇ I C´ , J. B ENGTSSON , M. C HRISTENSSON, Measuring and modeling the oxygen profile in a nitrifying Moving Bed Biofilm Reactor, Math. Biosci., 227 (2010), pp. 1–11. [11] A. M A Sˇ I C´ , H.J. E BERL, A Modeling and Simulation Study of the Role of Suspended Microbial Populations in Nitrification in a Biofilm Reactor, Bull. Math. Biol., 76 (2014), pp. 27–58. [12] O. WANNER , W. G UJER, A multispecies biofilm model, Biotechnol Bioeng., 28 (1986), pp. 314–328. [13] E. A LPKVIST, I. K LAPPER, A multidimensional multispecies continuum model for heterogeneous biofilm development, Bull. Math. Biol., 69 (2007), pp. 765–789. [14] J. B. X AVIER , C. P ICIOREANU , M. C. M. VAN L OOSDRECHT, A modelling study of the activity and structure of biofilms in biological reactors, Biofilms, 1 (2004), pp. 377–391. [15] C. P ICIOREANU , J. U. K REFT, M. C. M. VAN L OOSDRECHT, Particle-based multidimensional multispecies biofilm model, J. Appl. Environ. Microbiol., 70 (2004), pp. 3024–3040. [16] J. P. B OLTZ , E. M ORGENROTH , D. S EN, Mathematical modelling of biofilms and biofilm reactors for engineering design, Water Sci. Technol., 62 (2010), pp. 1821–1836. [17] J. P. B OLTZ , E. M ORGENROTH , D. B ROCKMANN , C. B OTT, W. J. G ELLNER , P. A. VANROLLEGHEM, Systematic evaluation of biofilm models for engineering practice: components and critical assumptions, Water Sci. Technol., 64 (2011), pp. 930–944. [18] B. S ZOMOLAY, Analysis of a moving boundary value problem arising in biofilm modelling, Math. Method Appl. Sci., 31 (2008), pp. 1835– 1859. [19] B. S ZOMOLAY, I. K LAPPER , M. D INDOS, Analysis of Adaptive Response to Dosing Protocols for Biofilm Control, SIAM J. Appl. Math., 70(8) (2010), pp. 3175–3202. [20] B. D’ACUNTO , L. F RUNZO, Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models, Math. Comput. Model., 43 (2011), pp. 1596–1606. [21] J. B. X AVIER , C. P ICIOREANU , M. VAN L OOSDRECHT, A framework for multidimensional modelling of activity and structure of multispecies biofilms, Environ Microbiol., bf 7 (2005), pp. 1085–1103. [22] B. D’ACUNTO , L. F RUNZO , I. K LAPPER , M. R. M ATTEI, Modeling multispecies biofilms including new bacterial species invasion, Math. Biosci., 259 (2015), pp. 20–26. [23] B. D’ACUNTO , L. F RUNZO, Free boundary problem for an initial cell layer in multispecies biofilm formation, Appl. Math. Lett., 25 (2012), pp. 20–26. [24] B. D’ACUNTO , G. E SPOSITO , L. F RUNZO , M. R. M ATTEI , F. P IROZZI Analysis and simulations of the initial phase in multispecies biofilm formation. Comm. Appl. Ind. Math., 2013, pp. 1–23, DOI: 10.1685/journal.caim.448.

16