Free boundary problem for the role of planktonic cells in biofilm formation and development

The dynamics of biofilm lifecycle are deeply influenced by the surrounding environment and the interactions between sessile and planktonic phenotypes. Bacterial biofilms typically develop in three distinct stages: attachment of cells to a surface, growth of cells into colonies, and detachment of cells from the colony into the surrounding medium. The attachment of planktonic cells plays a prominent role in the initial phase of biofilm lifecycle as it initiates the colony formation. During the maturation stage, biofilms harbor numerous microenvironments which lead to metabolic heterogeneity. Such microniches provide conditions suitable for the growth of new species, which are present in the bulk liquid as planktonic cells and can penetrate the porous biofilm matrix. We present a 1D continuum model on the interaction of sessile and planktonic phenotypes in biofilm lifestyle which considers both the initial attachment and colonization phenomena. The model is formulated as a hyperbolic-elliptic free boundary value problem with vanishing initial value. Hyperbolic equations reproduce the transport and growth of sessile species, while elliptic equations model the diffusion and conversion of planktonic cells and dissolved substrates. The attachment is modelled as a continuous, deterministic process which depends on the concentrations of the attaching species. The growth of new species is modelled through a reaction term in the hyperbolic equations which depends on the concentration of planktonic species within the biofilm. Existence and uniqueness of solutions are discussed and proved for the attachment regime. Finally, some numerical examples show that the proposed model correctly reproduces the growth of new species within the biofilm and overcomes the ecological restrictions characterizing the Wanner-Gujer type models.


Introduction
In recent years, the study of how the sessile and planktonic phenotypes interact in biofilm lifestyle has become a theme of intense interest and scrutiny [1]. Biofilms are microbial assemblies which commonly develop attached to abiotic or biotic surfaces. They are characterized by a solid matrix of extracellular polymeric substance (EPS) in which microorganisms are embedded [2]. The biofilm dynamics are deeply influenced by microbial mass exchanges between biofilm and the surrounding environment, which involve both the sessile and planktonic biomasses. The biofilm formation is initiated by pioneer microbial cells in planktonic form, which attach to a solid support through an initial attachment process. Such cells switch their mode of growth from planktonic to sessile and constitute the first sessile microbial colony [3], which develops and expands over time as a result of the microbial metabolic growth. Meanwhile, large EPS production by sessile cells confers high density and compactness to the aggregate and protects it from external agents. During the maturation stage, the high density induces large spatial gradients in biofilm properties, leading to numerous microenvironments and extremely heterogeneous microbial distributions. Specifically, new biological conditions arising within the biofilm can promote the phenomenon of microbial invasion: motile planktonic cells colonize the aggregate by penetrating the biofilm matrix, and proliferate as new sessile biomass where ideal conditions for their metabolic activity occur [4]. This means that the number of microbial species constituting the biofilm can increase over time, since microbial species initially not present can join the biofilm when new metabolic microniches arise. Furthermore, external shear forces, nutrients depletion and biomass decay lead to the detachment of cells from the biofilm colony into the surrounding medium [5]. Lastly, in the final stage of the biofilm lifecycle, microbial dispersal phenomena can occur: as a result of habitat decay (resource depletion and cell competition for space), planktonic cells, known as dispersed cells, are released in the surrounding environment, migrate to new surfaces and subsequently constitute new biofilm aggregates [6].
Despite the high amount of mathematical works on multispecies biofilms growth developed in the framework of the Wanner and Gujer model [7] or as multidimensional partial differential equation models [8; 9; 10; 11; 12], most of them completely neglect the attachment process in the initial phase of biofilm formation, since the initial data that prescribe location, size, and composition of colonies at the onset of the simulations are arbitrarily assigned. This strongly affects the biofilm development and maturation as highlighted by a recent work [13] where the attachment has been incorporated as a discrete stochastic process in a density-dependent diffusion-reaction model for cellulolytic biofilms. Furthermore, the Wanner-Gujer type models [7; 14; 15; 16] can lead, in some cases, to ecological restrictions on the number of species constituting the biofilm [17]. Indeed, they are characterized by a restriction on the number of species that can inhabit the biofilm under the detachment regime: that is if a species is not initially present within the biofilm on the support, it will be washed out from the system. The free boundary problem introduced in this work is intended to overcome these limitations by considering the initial biofilm formation mediated by planktonic cells as well as the colonization process. In particular, we present a one-dimensional continuous model considering two state variables representing the planktonic and sessile phenotypes and reproducing the transition from the former to the other in the biofilm lifecycle. The underlying model is a coupled hyperbolic-elliptic free boundary value problem with nonlocal effects. The attachment is modelled as a continuous, deterministic process which depends on the concentrations of the attaching species in the bulk liquid [16]. The colonization process which results in the establishment of new species in sessile form is modelled by considering an additional reaction term in the hyperbolic equations, which depends on the concentration of planktonic species within the biofilm [18]. The concentration of the planktonic species within the biofilm is governed by elliptic partial differential equations which describe their diffusion from the bulk liquid within the biofilm. A reaction term is considered to account for the conversion of the planktonic phenotype into the sessile mode of growth.
The work is organized as follows. Section 2 introduces the mathematical background for the attachment process in the initial phase of multispecies biofilm formation, in the framework of the Wanner-Gujer approach to biofilm modelling [16]. The free boundary is constituted by the biofilm thickness and it is assumed to be initially zero. The growth of the attaching species is governed by nonlinear hyperbolic partial differential equations. The free boundary is governed by a first order differential equation that depends on attachment, detachment and biomass growth velocity. It is recalled that the free boundary velocity is greater than the characteristic velocity of the mentioned hyperbolic system during the first instants of biofilm formation. As a consequence, the free boundary is a space-like line. The initial-boundary conditions for the microbial concentrations are assigned on this line and they are equal to the relative abundance of the species in the biomass attached to the biofilm-bulk liquid interface. The free boundary value problem is completed by a system of semi-linear elliptic partial differential equations that governs the quasi-static diffusion of substrates. In Section 3, a numerical experiment shows that the free boundary problem introduced in [16] needs to be generalized to eliminate any restriction on the number of species inhabiting the mature biofilm as described in [17]. Section 4 introduces the new free boundary problem which accounts for both the initial phase of biofilm formation and the diffusion and colonization of planktonic species within the biofilm. Section 5 introduces the integral version of the differential free boundary problem provided in Section 4, which is derived by adopting characteristics coordinates. An existence and uniqueness theorem of solutions is shown in Section 5 in the class of continuous functions. The proposed model is also solved numerically to simulate the biofilm evolution during biologically relevant conditions and provides interesting insights towards quantitative understanding of biofilm dynamics and ecology. Numerical results are reported in Section 6. Finally, the conclusions of the work are outlined in Section 7.

Background
A free boundary approach was introduced in [16] for modelling the initial phase of the multispecies biofilm formation and growth in the framework of Wanner and Guyer model [7]. In this context, denoting by X i pz, tq the concentration of the generic bacterial species i, the one-dimensional multispecies biofilm growth is governed by the following system of nonlinear hyperbolic partial differential equations where upz, tq denotes the velocity of the microbial mass, r M,i the specific growth rate, and ρ i the constant density. In addition, the substratum is assumed to be placed at z " 0. The function r M,i depends on X " pX 1 , ..., X n q, and substrates S j , j " 1, ..., m, as well r M,i " r M,i pXpz, tq, Spz, tqq, S " pS 1 , ..., S m q.
upr, tq is governed by the following equation: where Lptq represents the biofilm thickness. The substrate diffusion is governed by semi-linear parabolic partial differential equations that are usually considered in quasi-static conditions [14] D j B 2 S j Bz 2 " r S,j pXpz, tq, Spz, tqq, j " 1, ..., m, (2.4) where the functions r S,j denote the conversion rate of substrate j and D j the diffusion coefficients assumed constant. The biofilm thickness Lptq represents the free boundary of the mathematical problem. Its evolution is governed by the following ordinary differential equation [7; 16; 19; 20],

9
Lptq " upL, tq`σ a pψ˚q´σ d pLq, (2.5) where σ a denotes the attachment velocity of biomass from bulk liquid to biofilm and σ d the detachment velocity of biomass from biofilm to bulk liquid. The function σ a depends linearly on the concentrations ψi , i " 1, ..., n, ψ˚" pψ1 , ..., ψnq, of the microbial species in planktonic form present in the bulk liquid [7; 14; 16]. According to the experimental evidence, the ability of colonizing a clean surface is a feature of few microbial species, which are able to switch from their planktonic state, attach to the surface and start to secrete a polymeric matrix anchoring the cells to each other and to the surface. Even the formation of a single layer of cells can lead to a change on the electrostatic nature and mechanical properties of the surface, that can facilitate the attachment of new species that were initially unable to colonize the clean surface. According to [16], this is taken into account by considering in the formulation of the attachment flux σ a " ř n i"1 v a,i ψi {ρ i different attachment velocities v a,i for the single microbial species living in the liquid environment. Such velocities can be assigned constant or can be considered as functions of the environmental conditions affecting biofilm growth, that is substrate concentrations, biofilm composition itself, electrostatic and mechanical properties of the surface.
The function σ d is usually assumed to be proportional to L 2 : σ d " δL 2 , [21], where δ depends on the mechanical properties of the biofilm. In the initial phase of biofilm formation, where Lp0q " 0, the c(t 0 ,t) attachment is the prevailing process and σ d is very small, since so is L 2 . Therefore, it is σ a´σd ą 0 and the free boundary velocity is greater than the characteristic velocity, 9 Lptq ą upL, tq. The free boundary is a space-like line, as illustrated in Fig. 1. In the same figure, the characteristic-like lines of system (2.1) are also depicted. These lines, z " cpt 0 , tq, are defined by the differential initial value problem Bc Bt pt 0 , tq " upcpt 0 , tq, tq, cpt 0 , t 0 q " Lpt 0 q.
For mature biofilms the free boundary L becomes large, the detachment is the prevailing process, it is σ a´σd ă 0 and the free boundary is a time-like line, Fig. 2. The free boundary value problem (2.1)-(2.4) was discussed in [16] under the following initialboundary conditions: Lp0q " 0. (2.9) In equations (2.7), X i,0 ptq is the relative abundance of the species i in the biomass attached to the biofilm-bulk liquid interface [22]. More precisely, X i,0 ptq can be evaluated as where σ a,i " v a,i ψi ptq denotes the attachment flux of the single species i and σ a " ř n i"1 v a,i ψi ptq the total attachment flux. According to (2.10), the concentration of the microbial species at the biofilmliquid interface X i pLptq, tq for a multispecies biofilm growing under attachment regime, depends on both the concentrations of the same species in planktonic form in the bulk liquid and their attachment propensity. Note that, when all the microbial species in the bulk liquid are characterized by the same attachment velocity, equations (2.10) reduces to that is the volume fraction of the microbial species i at the biofilm-bulk liquid interface assumes the same value of the volume fraction within the bulk liquid. This reproduces the case of a biofilm that will be initially constituted by all microbial species inhabiting the surrounding liquid environment. However, going on with time the biofilm composition is affected by other factors such as substrate availability, specific microbial growth rate, detachment flux. For what concerns substrate diffusion, the first boundary condition (2.8) is the no flux condition at substratum. The functions Sj ptq in the second boundary condition (2.8) are prescribed functions in general.

Criticism
As outlined in [16], the model for the initial biofilm formation, summarized in the previous section, should be generalized to include the possibility that new attaching bacterial species can move downward within the biofilm matrix and colonize the regions where the conditions for their growth are optimal. An example, referred to as Case 1, could help to better understand the question. To discuss this special problem, an equivalent expression will be used for equations (2.1), where X i is replaced by the volume fraction f i defined by subjected to the constraint For equations above, conditions (2.7) are replaced by Let us consider a three species and substrate biofilm n " 3, m " 3 growing under time-dependent conditions. In particular, the model simulates the case of a biofilm growing in a liquid environment initially inhabited by species ψ1 and ψ2 and continuously fed with substrates S 1 and S 2 . At time t " t 1 ą 0, a third species ψ3 is supposed to be fed into the system Species ψ1 and ψ2 start to attach at t " 0 while the third at t " t 1 ą 0 Functions f i,0 ptq can be derived from equations (2.10), (3.6) and (3.7). Species f 1 and f 2 grow on substrate S 1 and S 2 , respectively. Species f 1 by consuming substrate S 1 produces S 3 , which is uptaken by f 3 . All species are supposed to grow only in sessile form, and the reactor is considered as an infinite reserve of substrates and planktonic species (Sj ptq " Sj ,0 ). The reaction terms r M,i and r S,j in equations (3.3) and (2.4) are modelled by using Monod type kinetics and are expressed as The values of the kinetic parameters and boundary conditions used in the numerical simulations are reported in Table 1. All the sessile species are supposed to have the same density ρ i " ρ, i " 1, ..., n.
The simulation time adopted for the numerical experiment is t " 10 d. We are aware that such simulation time will cover both the initial biofilm formation and the maturation phase where the detachment will be prevalent on the attachment flux. This choice is justified by the fact that we were interested in showing also the mature biofilm configuration, which is achieved under detachment regime. Fig. 3 shows the free boundary evolution and the characteristic line cpt 1 , tq starting from pLpt 1 q, t 1 q up to 0.8d simulation time. Figs. 4 and 5 show the biofilm composition and substrate trends within the biofilm over time, under attachment and detachment regimes respectively. The third species begins to adhere to the biofilm-bulk liquid interface at t " t 1 . Substrates S 1 and S 2 are consumed within biofilm by species f 1 and f 2 . As a consequence, favorable conditions for f 3 growth occurs within the inner biofilm region due to S 3 production. According to the uniqueness and  Parameter Definition Unit Value µmax,1 Maximum specific growth rate for f1 d´1 0.4 µmax,2 Maximum specific growth rate for f2 d´1 1.5 µmax, 3 Maximum specific growth rate for f3 d´1 0.5 K1 Half saturation constant for f1 on S1 g m´3 1 K2 Half saturation constant for f2 on S2 g m´3 20 K3 Half saturation constant for f3 on S3 g m´3 1 Y1 Yield of f1 on S1´´0.4 Y2 Yield of f2 on S2´´0.9 Y3 Yield of f1 on S1´´0.9 D1 Diffusion coefficient of S1 in biofilm Biofilm density g m´3 5000 δ Biomass shear constant m´1 d´1 2000 S1 ,0 S1 concentration in the bulk liquid g m´3 100 S2 ,0 S2 concentration in the bulk liquid g m´3 100 S3 ,0 S3 concentration in the bulk liquid g m´3 0 ψ1 ,0 ψ1 concentration in the bulk liquid g m´3 100 ψ2 ,0 ψ2 concentration in the bulk liquid g m´3 100 ψ3 ,0 ψ3 concentration in the bulk liquid g m´3  existence theorem provided in [16], the third species is confined within the region z ą cpt 1 , tq and does not colonize the region z ă cpt 1 , tq where there are favorable conditions for its growth f 3 pz, tq " " 0, 0 ď z ă cpt 1 , tq, 0 ď t ă t 1 , ą 0, z ě cpt 1 , tq, t ě t 1 . (3.12) As shown in Figs. 4 and 5, the third species is unable to penetrate the inner biofilm region. Under detachment regime, its concentration tends to zero and it is completely washed out from the biofilm system. This is a well-known behavior of the model [7], as stated in [17]. The model introduced in this work is intended to eliminate this great limitation.

Statement of the free boundary problem
This section presents the free boundary value problem for biofilm growth which considers its initial formation (σ a´σd ą 0, Lp0q " 0) and the diffusion and colonization of the planktonic species within the biofilm. It is a generalization of the problem discussed in [14; 16] and it is obtained by developing some ideas introduced in [18]. Specifically, an additional state variable is considered, Ψ i , Ψ " pΨ 1 , ..., Ψ n q, which represents the concentration of the planktonic species within the biofilm. The differential mass balance equations (2.1) are modified by adding a growth rate term r i that takes into account for the colonizing bacterial species, and further equations are introduced for the diffusion of the planktonic species. The resulting model is able to overcome the criticism outlined in Sec. 3, as shown through the simple examples reported in Sec. 6.
The biofilm growth is governed by the following equations BX i Bt`B Bz puX i q " ρ i r M,i pX, Sq`ρ i r i pΨ, Sq, 0 ď z ď Lptq, t ą 0, i " 1, ..., n, (4.1) Lptq " upLptq, tq`σ a pψ˚q, t ą 0, Lp0q " 0, BS j Bz p0, tq " 0, S j pL, tqq " Sj ptq, t ą 0, j " 1, ..., m. (4.7) The diffusion of the colonizing species within the biofilm is governed by semi-linear parabolic partial differential equations that are considered in quasi-static conditionś where r Ψ,i indicates the conversion rate due to the switch from planktonic to sessile mode of growth and D Ψ,i is the diffusivity coefficient of the planktonic species within the biofilm. Diffusion equations for Ψ are considered in quasi-static conditions for the same reason as S. Equations (4.8) are integrated with the following Neumann-Dirichlet boundary conditions BΨ i Bz p0, tq " 0, Ψ i pL, tq " ψi ptq, t ą 0, i " 1, ..., n, (4.9) where the no flux boundary conditions on the support are evident and the Dirichlet boundary conditions state that the values of the planktonic species on the free boundary are the same as in the bulk liquid.
Note that equation (4.3) refers to the initial phase of biofilm formation, when the detachment flux σ d is negligible compared to σ a . The free boundary Lptq is a space-like line and equation (4.2) provides the initial conditions for the microbial species in sessile form on the free boundary. Conversely, during the maturation stage of biofilm growth the detachment flux is predominant and the free boundary is represented by a time-like line as stated in Sec. 2. The free boundary value problem referring to the mature phase of biofilm growth and considering the interaction between the planktonic and sessile phenotype through the colonization process has been investigated both qualitatively and numerically in [23; 24].
The functions introduced in equations (5.4)-(5.7) are defined below Σpt 0 q " An existence and uniqueness theorem for the integral problem (5.4)-(5.9) can be proved in the space of the continuous functions as generalization of the results in [16].

Numerical applications
Numerical simulations have been performed to test the behavior of the model formulated in Sec. 4. Specifically, we have considered the same biofilm system of Sec. 3 composed of 3 microbial species and 3 dissolved substrates. The planktonic species present in the bulk liquid are able to initiate the biofilm formation through the attachment process, penetrate the biofilm matrix once constituted and establish where the most appropriate growth conditions are found. We have explored two ideal biological situations. In the first case the species ψ3 is not initially present in the bulk liquid but it arrives at time t " t 1 and starts to attach to the external surface of the biofilm as well as penetrate the biofilm matrix. In the second case, the species ψ3 is not able to attach to the biofilm surface (v a,3 " 0) but it can establish in sessile form through the colonization process. These biological situations will be referred to as Case 2 and Case 3. The reaction terms r M,i and r S,j in equations (4.1) and (4.6) have been adopted according to (3.10) and (3.11). The values for Ψ i on the free boundary have been set according to (3.6) and (3.7). The values for S j on the free boundary are reported in Table 1. The reaction terms concerning the colonization process in equations (4.1) and (4.8) are modelled using Monod type kinetics and are expressed as where k col,1 , k col,2 , k col,3 are the maximum colonization rates of motile species, and Y Ψ,1 , Y Ψ,2 , Y Ψ,3 are the yields of the sessile species on planktonic ones. The values of such kinetic parameters and the diffusion coefficients for Ψ i are reported in

Case 2: Attachment and colonization of microbial species ψ3
The results reported in Figs. 6 and 7 highlight model capability to reproduce both the attachment and colonization phenomena that strongly affect biofilm lifecycle. In particular, it is possible to notice that during the initial phase of biofilm formation, the biofilm undergoes the same development illustrated in Sec. 3 and reported in Fig. 7. However, at time t " 0.50 d it is visible that the volume fraction of the third species f 3 is slightly positive even in the region z ă cpt 1 , tq due to the colonization phenomenon ( Fig. 6(A2)). Going on with the simulation time, f 3 increases all over the biofilm leading to a higher biofilm thickness at time t " 1 d (Fig. 7(A1)). At the final simulation time t " 10 d and under detachment regime, the biofilm is constituted by all the species inhabiting the bulk liquid ( Fig.  7(A2)) conversely to the numerical results reported in Sec. 3 where the complete washout of species f 3   has been observed. The different biofilm stratification affects substrate trends as it is possible to notice that at the final simulation time, S 3 concentration is much lower when compared to the numerical example of pure attachment regime.
Case 3: Pure colonization of microbial species ψ3 Figs. 8 and 9 illustrate the biofilm development and substrate trends when the species ψ3 is not able to attach to the biofilm surface, but it can penetrate the biofilm matrix and establish in sessile form. According to Figs. 8-9, numerical results reveal that for all simulation times the biofilm thickness is smaller compared to the pure attachment case. This contributes to have different substrate trends within the biofilm (Figs. 8-9(B1-B2)). In terms of biomass distribution, it is possible to notice that the third species grow in sessile form in the inner layers of the biofilm where there is the highest S 3 concentration. The biomass stratification and subtrate trends at the final simulation time resemble the one achieved for Case 2. Such results highlight an important feature of the model: the attachment and colonization phenomena are both dependent on the planktonic cells present in the bulk liquid. They can occur simultaneously reproducing the case of planktonic cells able to attach to the surface and penetrate the biofilm matrix. Conversely, the planktonic cells can be characterized by a certain motility which drives them to the biofilm region where there are the most appropriate conditions for their growth.

Conclusion
The proposed model comprehensively describes the transition from planktonic to sessile phenotype which governs the biofilm dynamics. This allows to properly reproduce the evolution of biofilms starting from the initial formation and including the establishment and growth of new species. The criticism of Wanner and Gujer type models, discussed in [17], is here emphasized through a numerical example. Such models are not able to properly describe the growth of microbial species which do not participate in the initial biofilm formation attaching later to a pre-existing aggregate. The presented model is able to overcome this issue as it considers both the initial attachment phase and the growth of new sessile species within the biofilm mediated by the invasion process. The modelling of the initial phase of biofilm formation allows to describe the biofilm growth without arbitrarily fixing the initial composition of the biofilm. The existence and uniqueness of solutions is proved in the case of attachment regime. Numerical examples are provided to show model capability to reproduce the different stages of biofilm growth as affected by the planktonic phenotype. Future work may be related to the role of biofilm porosity on planktonic species diffusion and the qualitative analysis under detachment regime.