Cette thèse s'intèresse à l'application d'un modèle mathématique unidimensionnel de formation et de croissance de biofilms multi-espèces. Le modèle se compose d'un système d'équations non linéaires aux dérivées partielles hyperboliques, décrivant la croissance d'espèces microbiennes dans le biofilm, et un système d'équations semi-linéaires aux dérivées partielles paraboliques, qui régit la diffusion de substrat de la phase aqueuse vers la matrice du biofilm. L'ensemble I would like to express my deep and sincere gratitude to all the wonderful people I met during these three years as a PhD student. Despite all the challenges and sleepless nights, this worthwhile experience really changed my life and made me growth both personally and scientifically. First of all, I would like to thank the Erasmus Mundus Joint Program ET eCoS 3 (Environmental Technologies for Contaminated Solids, Soils and Sediments) committee for giving me the opportunity of being part of this stimulating research program. I specially thank my promotor Prof. Giovanni Esposito and co-promotor Prof. Berardino D'Acunto for all their help and guidance along the way. With their constant encouragement they helped me overcoming the difficulties and taught me how to be open minded and always critically solve the problems. I would also like to thank my supervisor Prof. Francesco Pirozzi for his precious collaboration; his comments and suggestions have always driven me to do more

Abstract English

This dissertation relates to the applications of a one-dimensional mathematical model for multispecies biofilm formation and growth. The model consists of a system of nonlinear hyperbolic partial differential equations, describing the growth of microbial species in biofilms, and a system of semilinear parabolic partial differential equations, which governs substrate diffusion from the surrounding aqueous phase into the biofilm. Overall, this leads to a free boundary value problem, essentially hyperbolic.

In a first study, the analysis and simulations of the attachment phenomena in the initial phase of biofilm growth have been addressed.

The resulting mathematical problem has been discussed by using the method of characteristics and the fixed-point theorem has been used to obtain existence, uniqueness and properties of solutions.

A second aspect of the thesis deals with the analysis and prediction of population dynamics in multispecies biofilms for wastewater treatment. The model has been applied to simulate the bacterial competition and to evaluate the influence of substrate diffusion on microbial stratification for a nitrifying multispecies biofilm including Anammox bacteria and a sulfate-reducing biofilm. In both cases, specific kinetics equations have been introduced to describe biomass growth and substrate consumption. The method of characteristics has been used for numerical purposes and the mass conservation equation plays a crucial role in checking the accuracy of simulations. The simulation results reveal that the model is able to evaluate properly the effects that boundary conditions exert on bacterial competition.

Finally, the biofilm model has been extended to include the colonization phenomenon. The new model is able to take into account the invasion of new species diffusing from bulk liquid to biofilm, still based on a set of nonlinear hyperbolic partial differential equations for what concerns the growth process. Indeed, the biological invasion process of new species into the biofilm has been modeled by a system of nonlinear parabolic partial differential equations. The invasion model has been successfully applied to simulate the invasion of heterotrophic bacteria in a constituted autotrophic biofilm and viceversa. i

Italian

Il presente lavoro di ricerca riguarda le applicazioni di un modello monodimensionale relativo alla formazione e crescita di biofilm multispecie. Il modello in esame consiste in un sistema di equazioni differenziali a derivate parziali di tipo iperbolico, che descrive la crescita delle diverse specie microbiche all'interno del biofilm, ed un sistema di equazioni differenziali a derivate parziali paraboliche semilineari, che governa la diffusione dei substrati dalla fase acquosa all'interno del biofilm. Nel complesso, ciò comporta la formulazione di un problema a frontiera libera, essenzialmente iperbolico.

In un primo studio, sono state condotte analisi e simulazioni riguardo la fase iniziale di crescita del biofilm. Il problema matematico che ne deriva è stato discusso utilizzando il metodo delle caratteristiche ed il teorema del punto fisso è stato utilizzato per ottenere esistenza, unicità e proprietà delle soluzioni.

Un secondo aspetto della tesi riguarda l'analisi e predizione delle dinamiche di popolazione di biofilm multispecie applicati al trattamento delle acque. Il modello è stato applicato al caso di un biofilm nitrificante contente batteri Anammox ed un biofilm solfato-riduttore, al fine di simulare la competizione batterica e valutare l'influenza che la diffusione dei substrati esercita sulla stratificazione microbica. In entrambi i casi, specifiche equazioni cinetiche sono state introdotte per descrivere la crescita della biomassa e il consumo dei substrati. Il metodo delle caratteristiche è stato utilizzato per scopi numerici e l'equazione di conservazione della massa ha giocato un ruolo fondamentale nella valutazione dell'accuratezza delle soluzioni. I risultati delle simulazioni rivelano che il modello è in grado di valutare correttamente gli effetti che le condizioni al contorno esercitano sulla competizione microbica.

Infine, il modello sul biofilm è stato esteso al fine di includere il fenomeno di colonizzazione. Il nuovo modello è in grado di trattare l'invasione di nuove specie che diffondono dal bulk liquido verso il biofilm pur basandosi su un set di equazioni differenziali a derivate parziali di tipo iperbolico per quanto riguarda il processo di crescita microbica. In effetti, il processo di invasione di nuove specie all'interno del biofilm è stato modellato tramite un sistema di equazioni differenziali a derivate parziali paraboliche non lineari. Il modello sull'invasione è stato applicato con successo alla simulazione dell'invasione eterotrofa in un biofilm autotrofo già costituito e viceversa. conduit à un problème de valeur limite libre, essentiellement hyperbolique.

Dans une première étude, l'analyse et la simulation de la phase initiale de croissance du biofilm ont été examinées. Le problème mathématique résultant a été discuté en utilisant la méthode des caractéristiques et le théorème du point fixe a été utilisé pour déterminer l'existence et l'unicitè des solutions mathématiques.

Un deuxième aspect de la thèse porte sur l'analyse et la prévision de la dynamique des populations microbienne dans plusieurs types biofilms pour le traitement des eaux usées. Le modèle a été appliqué pour simuler la compétition bactérienne et évaluer l'influence de la diffusion du substrat sur la stratification microbienne des biofilms multi-espèces, en incluant les bactéries nitrifiantes, Anammox et bactéries sulfato-réductrices. Dans les deux cas, des spécifiques équations de cinétique ont été introduites pour décrire la croissance de la biomasse et la consommation du substrat. La méthode des caractéristiques a été utilisée à des fins numériques et l'équation de conservation de masse joue un rôle crucial pour vérifier l'exactitude des simulations. Les résultats des simulations montrent que le modèle est en mesure d'évaluer correctement les effets des conditions limites qui s'exercent sur la concurrence bactérienne.

Enfin, ce modèle a été étendu pour inclure le phénomène de colonisation microbienne. Le nouveau modèle est capable de prendre en compte l'invasion de nouvelles espèces en se basant sur un ensemble d'équations non linéaires aux dérivées partielles hyperboliques pour ce qui concerne le processus de croissance. De plus, le processus d'invasion biologique d'espèces nouvelles dans le biofilm a été modélisé par un système d'équations non linéaires aux dérivées partielles paraboliques. Ce modèle d'invasion a été appliqué avec succès pour simuler l'invasion des bactéries hétérotrophes dans les biofilms autotrophes.

Dutch

Dit proefschrift beschrijft de applicaties van een één-dimensionaal wiskundig model voor de vorming en groei van multi-species biofilms. Het model bestaat uit een systeem van niet-lineaire hyperbole partiële differentiaal vergelijkingen, welke de groei van microorganismen in biofilmen beschrijft, en een systeem van semi-lineaire parabool partiële differentiaal vergelijkingen, welke de substraat diffusie van de omliggende waterige fase naar de biofilm beschrijft. Samen leidt dit tot een vrij randvoorwaarden probleem, in essentie hyperbolisch.

Een eerste studie behandelde de analyse en simulatie van hechtingsfenomenen gedurende de initiële biofilmvorming. Het resulterende mathematisch probleem werd besproken aan de hand van de methode van karakteristieken en het vastepunt theorema is gebruikt om het bestaan, de uniekheid en de eigenschappen van de oplossing aan te tonen.

Het tweede deel van deze thesis behandelt de analyse en voorspelling van populatie dynamica van multi-species biofilmen voor afvalwaterzuivering. Het model is toegepast om de bacteriële competitie te simuleren en de invloed van substraat diffusie op microbiële laagvorming voor een nitrificerende multi-species biofilm, inclusief Anammox bacteriën en een sulfaat reducerende biofilm. In beide gevallen werden specifieke kinetische vergelijkingen gebruikt om de groei van biomassa en het substraatverbruik te beschrijven.

De methode van karakteristieken is gebruikt voor numerieke doeleinden en de wet van behoud van massa behoud speelt een belangrijke rol in het checken van de nauwkeurigheid van de simulaties. De simulaties laten zien dat het model in staat is om de effecten die de randvoorwaarden uitoefenen op de bacteriële competitie te evalueren.

Ten slotte is het biofilm model uitgebreid om kolonisatie fenomenen te berekenen. Het nieuwe model is in staat om rekening te houden met de invasie van een nieuw species diffunderend van de bulk oplossing naar de biofilm, nog steeds gebaseerd op een set van niet-lineaire hyperbool partiële differentiaal vergelijkingen voor het groeiproces. Inderdaad, het biologische invasie proces van een nieuw species naar de biofilm is gemodeleerd door een systeem van niet-lineaire parabolische partiëel differentiaal vergelijkingen. Het invasie model is succesvol toegepast om de invasie van heterotrofe bacteriën in een autotrofe biofilm en vice versa te simuleren.

List of Figures

Schematic representation of biofilm model classification adopted

in this review and based on biomass representation and dimensionality: a) 1D continuum models; b) multidimensional continuum models; c) Cellular Automata; d) Individual based Models (for multidimensional models only the 2D representation has been reported). . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Spreading mechanism adopted in [START_REF] Wanner | A multispecies biofilm model[END_REF] and [START_REF] Kissel | Numerical simulation of mixedculture biofilm[END_REF] (figure adapted from [START_REF] Be Rittmann | Development and experimental evaluation of a steady-state, multispecies biofilm model[END_REF]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 2D continuum models: a) spreading mechanism adopted by Eberl et al. [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF]; b) spreading mechanism adopted by Dockery and Klapper [START_REF] Klapper | Finger formation in biofilm layers[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Schematic representation of the spreading rules adopted by a)

Picioreanu et al. [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF][START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF]; b) Laspidou and Rittmann [START_REF] Laspidou | Modeling the development of biofilm density including active bacteria, inert biomass, and extracellular polymeric substances[END_REF][START_REF] Laspidou | Evaluating trends in biofilm density using the umcca model[END_REF]. Figure adapted from [START_REF] Cs Laspidou | Cellular-automata and individual-based approaches for the modeling of biofilm structures: Pros and cons[END_REF][START_REF] Tang | An improved cellular automaton method to model multispecies biofilms[END_REF]. . . . . . . . . . . . . . . . . . . . 2.5 Spreading mechanisms adopted by: a), c) Xavier et al. [START_REF] Xavier | A framework for multidimensional modelling of activity and structure of multispecies biofilms[END_REF]; b) Lardon et al. [START_REF] La Lardon | idynomics: next-generation individual-based modelling of biofilms[END_REF].

Research context

Biofilm heterogeneity and complex functioning offer research perspectives in different fields. Biofilms have been successfully appplied to wastewater treatment thanks to their ability to remove unwanted compounds from wastewater. Their structure and activities have been investigated using a broad variety of microscopic, physico-chemical and molecular biological techniques. Parallel to these investigations more and more complex mathematical models and simulations have been developed to describe the growth, structures and interactions of biofilms [START_REF] Wimpenny | Heterogeneity in biofilms[END_REF].

The analysis and modeling of a biofilm system are a scientific challenge and require in many cases interdisciplinary cooperation. Mathematical modeling of biofilm systems and in particular of their application to wastewater treatment is a widely accepted tool to improve the understanding of the fundamental mechanisms regulating biofilm formation and performance, to formulate and validate hypothesis and operating strategies, to predict system's behavior under different conditions without the cost, time, and risk of building an experimental prototype [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF].

The main elements of biofilm literature, mostly related to this study are briefly presented in the following sections.

Biofilm definition

Biofilms are a form of microbial ecosystems constituted by accumulations of microorganisms, irreversibly associated to a solid surface or phase interphase and embedded in a self-produced primarily polysaccharide matrix [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF]. This matrix protects the cells against mechanical washout, facilitates communication among them through biochemical signals and offers resistance to antimicrobials and biocides.

Bacteria living in a biofilm are not randomly distributed but they live in distinct niches. They benefit from interspecies cooperation [START_REF] Elias | Multi-species biofilms: living with friendly neighbors[END_REF] and show more resistance to toxic substances such as antibiotics, chlorine and detergents thanks 1 to diffusion barriers [START_REF] Jw Costerton | Biofilms, the customized microniche[END_REF]. Based on a mechanistic hypothesis, the spatial organization of the microorganisms in biofilms can be seen as the result of differences in local substrate availability. Neverthless, recent studies have shown that bacterial signaling and other interactions of that type, such as the gene expression, often affect the regulatory activities in the cells resulting in coordinated performances and in the regulation of the organization of microbial communities [START_REF] Klapper | Mathematical description of microbial biofilms[END_REF][START_REF] Tolker-Nielsen | Spatial organization of microbial biofilm communities[END_REF]. Therefore, biofilms are not just simple aggregations of selfish individuals in competition for the nutrients they need, neither barely multicellular organisms governed by hormone-based communication control: they rather act as highly differentiated and multicultural communities similar to our own cities [START_REF] Davey | Microbial biofilms: from ecology to molecular genetics[END_REF][START_REF] Watnick | Biofilm, city of microbes[END_REF]. Nearly all biofilms communities in nature comprise a variety of microorganisms embedded in a exopolymeric matrix and functioning as a cooperative consortium, in a relatively complex and coordinate manner [START_REF] Davey | Microbial biofilms: from ecology to molecular genetics[END_REF].

Harmful and beneficial biofilms

The importance and attractiveness of biofilm systems is undisputed. They are found in extremely varied environments, ranging from water distribution systems and wastewater treatment plants to stream beds, ship hulls, and teeth surfaces [START_REF] Ee Gonzo | Continuum heterogeneous biofilm model-a simple and accurate method for effectiveness factor determination[END_REF]. Biofilms play both beneficial and detrimental roles in their environment depending on whether their formation is controlled or unintentional [START_REF] Mašić | Investigation of a biofilm reactor model with suspended biomass[END_REF]. By corroding pipes, reducing heat transfer, contaminating drinking water or infecting medical implants, biofilms can affect production quality or constitute a risk for human health [START_REF] Coufort | Cohesion and detachment in biofilm systems for different electron acceptor and donors[END_REF]. Conversely, an extensive use of biofilms is made within the field of environmental biotechnology, i.e. self-purification of water, groundwater protection, soil remediation and pollution treatment. Biofilms have been widely used to treat wastewaters in attached growth systems since the end of the 19th century.

Attached and suspended growth systems are based on the same biological metabolic processes to remove carbon and nutrients from wastewater, but there are some inherent differences that provide several advantages and some challenges to the application of biofilm-based processes. The suspended growth systems are characterized by the presence of biological flocs which are removed by gravity settling and whose shape makes dissolved substrates available to all of the cells. However, suspended flocs forming biomasses are often washed out from the system and hence they experience low microbial diversity function [START_REF] Ras | Extracellular polymeric substances diversity of biofilms grown under contrasted environmental conditions[END_REF]. In attached growth systems, the treatment process can be operated at higher biomass concentration in the reactor, with no need of settlers which are used in biomass retention and recirculation.

The transport of dissolved substrates is generally governed by diffusion and usually results in a concentration gradient within the aggregate, which favours in turns the formation of a growth rate gradient [START_REF] Nicolella | Wastewater treatment with particulate biofilm reactors[END_REF]. Environmental micro-niches created by diffusion-reaction interactions allow diverse bacterial species to cohabit the same biofilm. Bacterial cells easily adapt to the local surroundings and respond to the concentration gradients, which often results in a stratification where different species occupy specific layers in the biofilm. As a result, biofilms can exhibit various environmental and kinetic characteristics and such diversified microbial activities can be beneficially exploited for treating contemporarily multiple water pollutants. In this case, the overall reactor performance may be closely linked to microscale interactions within the biofilm [START_REF] Morgenroth | Biofilm engineering: linking biofilm development at different length and time scales[END_REF] as the efficiency and robustness of a wastewater treatment plant mainly depend on the composition and activity of its microbial community [START_REF] Wagner | Microbial community composition and function in wastewater treatment plants[END_REF].

Biofilm development

Biofilm formation is a multi-stage process resulting from the balance of several physical (substrate transport), chemical and biological factors (growth yields and substrate conversion rates), and constituted by several stages: adhesion to biofilm support, formation of an attached monolayer and cell proliferation (microcolonies), development of a mature biofilm and lastly dispersion.

The formation of bacterial biofilms starts when a small number of bacteria adhere to a surface. The attachment to the surface can be classified in reversible and irreversible since the initial bound established by the cells with the surface is often weak and bacteria need active mechanisms to reach a more stable association with the surface. The irreversible attachment is characterized by the production of a polysaccharide matrix, which facilitates the bacteria adhesion to the surface. Once attached to the surface, bacteria proliferate, differentiate and form microcolonies which represent the biofilm basic structural units. They are discrete matrix-enclosed communities of bacterial cells that may include cells of one or of many species and their formation depends on clonal growth of attached cells or on active translocation of cells across the surface. Micro colonies grow in size and coalesce to form macro colonies, in many cases constituted by mushroom-like towers separated by fluid-filled voids. This is the point at which the biofilm reaches its maturation, developing a complex three-dimensional dynamic structure, and bacteria profoundly differ from their planktonic state relatively to the number of proteins expressed. At the final stage of biofilm development, macro colonies may be eroded releasing cells from the biofilm into the surrounding medium. Conversely, detachment can be initiated internally leading to the dispersion of individual cells or large clumps. The revert of dispersed biofilm cells to the planktonic state completes this idealistic development cycle [START_REF] El Moustaid | Modeling bacterial attachment to surfaces as an early stage of biofilm development[END_REF][START_REF] Rd Monds | The developmental model of microbial biofilms: ten years of a paradigm up for review[END_REF].

Main processes involved in biofilm systems

Biofilms are both complex ecological and mechanical systems, in which different kinds of processes interact. Biofilm development is determined by "positive" processes (i.e. cell attachment, cell division, and polymer production) which lead to biofilm volume expansion, and "negative" processes, (i.e. cell detachment and cell death) which contribute to biofilm shrinking [START_REF] Picioreanu | Modelling and predicting biofilm structure[END_REF].

Initial attachment of bacterial cells is affected by several factors, including mass transport, surface conditioning, hydrophobicity, surface charge and roughness. Surface colonization is usually performed by some bacterial species which exhibit better capability to attach to the substratum. The attachment and prolification of these bacteria provide conditions that allow additional species to form mixed biofilm. Recent advances in microbial ecology have identified motility as one of the main mediator of multispecies biofilm formation. In particular, it can determine biofilm landscape and contribute to rapid alterations in biofilm populations [START_REF] Houry | Bacterial swimmers that infiltrate and take over the biofilm matrix[END_REF]. Motile bacteria generally have flagella that rotate to propel them and represent one of the features distinguishing the microbes' planktonic form from the sessile lifestyle. The production of flagella and of extracellular matrix seems to be mutually exclusive processes: bacteria give up the ability to move in order to settle down [START_REF] Kolter | Microbial sciences: the superficial life of microbes[END_REF]. Indeed, once a motile bacteria, supplied by the liquid planktonic phase, has successfully infiltrated the biofilm matrix, it can invade a resident community and establish where the environmental conditions are optimal for its growth.

The main biofilm expansion is due to bacterial growth and to extracellular polymer production. The soluble substrates necessary for bacterial growth are dissolved in the liquid flow and to reach the cells, first they pass through a boundary layer (external mass transfer resistance) and then through the biofilm matrix (internal mass transfer resistance) [START_REF] Picioreanu | Modelling and predicting biofilm structure[END_REF].

The transport of nutrients into and through the biofilm matrix is mainly governed by diffusion, although the presence of a porous structure allows advection even within the aggregate. Substrate diffusivity inside the biofilm is less than that in water because of the minimal permeability of biofilm aggregate. In addition, substrate diffusivity has been found to decrease with biofilm depth as a consequence of an increasing density, decreasing porosity, and decreasing permeability with depth [START_REF] Jp Boltz | Systematic evaluation of biofilm models for engineering practice: components and critical assumptions[END_REF]. Besides, the boundary layer is constituted by a thin liquid layer, characterized by a negligible flow over the biofilm/liquid interface. Its thickness mostly depends on the biofilm surface and on the flow regime. The external fluid flow regulates biofilm growth by establishing the concentration of substrates and products at the liquid-solid interface. At the same time the fluid flow shears the biofilm surface eroding the protuberances.

Then, biofilm structure results from the interplay of different interactions, such as mass transfer, conversion rates and detachment forces. An accurate modeling of such a system has to take all of these factors into account.

Biofilm modeling

Biofilm models are considered an active research field, a scientific development which is continuously evolving according to the modern studies in microbiology and with the increasing achievements in their possible applications. According to the IWA Task Group on Biofilm Modelling [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF], existing models have undergone a chronological evolution and can be classified in three main groups:

first generation models, which are based on the simplifying assumptions of steady-state biofilm, uniform distribution of biomass, simplest possible geometry and single limiting substrate. They were aimed at connecting the basic processes occurring in a steady-state biofilm, that is, mass transport from the bulk liquid, biodegradation within the biofilm, biofilm growth and decay;

second generation models, which have been introduced in 1980s and include microbial interactions and non-uniform distribution of biomass still assuming one-dimensional mass transport (by diffusion) and biofilm geometry;

third generation models, which have been developed since 1990s and are based on the idea that some biofilms are characterized by complex and heterogeneous structures that cannot be captured by pioneers models which are based on the assumption of uniform and homogeneous structures. They are aimed at describing spatial non-uniformities that characterize biofilm shape.

Biofilm models can also be categorized in one-dimensional or multidimensional, dynamic or steady-state, single or multispecies etc. More information on the different modeling approaches present in literature, are provided in Chapter 2.

One-dimensional biofilm models have been recognized as sufficiently accurate tools to predict and evaluate biofilm reactor performance. They often allow for mathematical analysis even if they may lack in proper description of biofilm structure and heterogeneity [START_REF] Jp Boltz | Systematic evaluation of biofilm models for engineering practice: components and critical assumptions[END_REF]. Conversely, multidimensional models are characterized by a high level of complexity and require in many cases a demanding numerical treatment and large computing power. A variety of (one-dimensional biofilm) modeling approaches exist in literature. In general, when applied to the case of biofilm reactors, they are able to provide as outputs the following items: biofilm composition in terms of relative mass proportion due to microbial competition and spatial profiles of any number of particulate components in the biofilm; 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; in some cases spatial and dynamic profiles of porosity (consolidation phenomena).

To date, the one-dimensional biofilm model introduced by Wanner and Gujer in 1986 [START_REF] Wanner | A multispecies biofilm model[END_REF] has been implemented in its original and modified versions in several software packages, which have been widely used by engineers for biofilm process design and verification purposes. Such a model model treats biomass as a continuum, is based on the assumption of incompressibility for biofilm matrix and simulates biofilm growth as a convective transport mechanism and substrate transport as a diffusive mechanism. Despite the one dimensionality of the model, ecological aspects of biofilms, such as spatio-temporal population and substrate dynamics, are described relatively in details [START_REF] Abbas | Mathematical contribution to one-dimensional biofilm modeling[END_REF]. However, as most of one-dimensional models, the Wanner-Gujer model includes some elements of uncertainty, mostly related to the fate of particulate (i.e. transport to the biofilm surface, attachment/entrapment, movement inside the biofilm, and hydrolysis and subsequent degradation of particulate substrates), the dynamics and rate of biofilm detachment, external mass transfer boundary layer thickness, and lack of clear and transparent biofilm model calibration protocol [START_REF] Boltz | Mathematical modelling of biofilms and biofilm reactors for engineering design[END_REF]. In addition, the equations constituting the model result in a complex mixed nonlinear hyperbolic-parabolic free boundary value problem, hard to treat with rigorous analytical techniques. Actually, only few studies have addressed the qualitative analysis of the system [START_REF] Abbas | Mathematical contribution to one-dimensional biofilm modeling[END_REF][START_REF] Pritchett | Steady state solutions of a one-dimensional biofilm model[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Szomolay | Analysis of a moving boundary value problem arising in biofilm modelling[END_REF][START_REF] Klapper | Productivity and equilibrium in simple biofilm models[END_REF][START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF].

Scope of the thesis

The overall objective of this dissertation was to study qualitatively and numerically a one-dimensional multispecies multisubstrate biofilm model applied to wastewater treatment and propose a new approach to evaluate the invasion of new species into a constituted biofilm. In particular, the objectives of this study include: 1) qualitative analysis of the free boundary value problem governing the initial phase of biofilm growth; 2) development of numerical simulations to illustrate the model in the case of wastewater treatment and to evaluate interactions within multispecies biofilms; 3) development of an invasion model, able to take into account the colonization phenomenon and providing insights on the movement of particulate components inside the biofilm. The fulfillment of the research objectives pursued in this thesis is described in the following chapters. An extensive literature review of the existing biofilm models is provided in Chapter 2. General guidelines for the selection of the most suitable modeling approach are discussed. A mathematical model based on a continuum approach and able to describe the attachment process during the initial phase of biofilm growth is presented and analyzed qualitatively and numerically in Chapter 3. Chapter 4 and 5 are dedicated to the development of numerical simulations for the evaluation of multispecies biofilm performance. Chapter 6 introduces the "invasion model". In Chapter 7, conclusions and recommendation for future research are presented.

Chapter 2

Continuum and discrete approach in modeling biofilm development and structure: a review

The scientific community has recognized that almost 99% of the microbial life on earth is represented by biofilms. Considering the impacts of this sessile lifestyle on both natural and human activities, extensive experimental activity has been carried out to understand how biofilms grow and interact with the environment. Many mathematical models have also been developed to simulate and elucidate the main processes characterizing the biofilm growth. Two main mathematical approaches for biomass representation can be individuated: continuum and discrete. This review is aimed at evaluating the main features deriving from the application of each approach. Continuum models can simulate the biofilm processes in a quantitative and deterministic way. However, they require a multidimensional formulation to take into account the biofilm spatial heterogeneity, which makes the models quite complicated, significantly increasing the computational efforts. Discrete models are more recent and can represent the typical multidimensional structural heterogeneity of biofilm in good agreement with experimental expectations, but they generate computational results that include elements of randomness, introducing stochastic effects into the solutions. This paper gives general guidelines to select the most suited model, based on the objectives and needs of the model user.

Introduction

Recent advances in quantitative recovery and in direct observation of microbial populations have revealed that biofilms represent the prevailing structures in microbial lifestyle [START_REF] Jw Costerton | Microbial biofilms[END_REF].

In most natural and human environments, biofilms are constituted by highly structured multispecies communities composed of millions of microorganisms that accumulate on surfaces and secrete extracellular polymers which anchor the cells to each other as well as to the surfaces [START_REF] Klapper | Mathematical description of microbial biofilms[END_REF][START_REF] Tolker-Nielsen | Spatial organization of microbial biofilm communities[END_REF][START_REF] Davey | Microbial biofilms: from ecology to molecular genetics[END_REF][START_REF] Watnick | Biofilm, city of microbes[END_REF][START_REF] Jw Costerton | Overview of microbial biofilms[END_REF][START_REF] Klapper | Finger formation in biofilm layers[END_REF][START_REF] Hf Jenkinson | Biofilms adhere to stay[END_REF][START_REF] Stoodley | Biofilms as complex differentiated communities[END_REF][START_REF] Van Loosdrecht | Biofilm structures[END_REF].

The bacteria living in a biofilm are not randomly distributed but they live in distinct niches, benefit from interspecies cooperation [START_REF] Elias | Multi-species biofilms: living with friendly neighbors[END_REF] and show more resistance to toxic substances such as antibiotics, chlorine and detergents thanks to diffusion barriers [START_REF] Jw Costerton | Biofilms, the customized microniche[END_REF]. Biofilms are important components of food chains and are involved in self-purification processes in soil, water and sediments and in the biodegradation of organic compounds including environmental pollutants. They have been used to treat wastewater since the end of the 19th century [START_REF] Nicolella | Wastewater treatment with particulate biofilm reactors[END_REF]. Comparing with suspended cells, the bacteria growing in biofilms show some advantages: i) they cannot be washed away with the water flow but they grow in locations where their food supply remains abundant; ii) they show an increased resistance to antimicrobial agents and allow the achievement of a higher biomass concentration value in bioreactors; iii) their physical structure permits the formation of several bacterial species contributing to the treatment of different organic and inorganic substrates [START_REF] Wang | Review of mathematical models for biofilms[END_REF]. At the same time biofilms can have a significant impact on the surrounding environment, including biofouling, biocorrosion, oil field souring and infections in host tissues or medical implants [START_REF] Klapper | Finger formation in biofilm layers[END_REF].

Biofilm formation is a dynamic process resulting from the balance of several physical (substrate transport, detachment, etc) and biochemical factors (microbial growth, substrate conversion, etc) [START_REF] Rd Monds | The developmental model of microbial biofilms: ten years of a paradigm up for review[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF]. The formation of various biofilm architectures and the related activities are strongly affected by the specific environmental conditions, such as electron donors and acceptors levels, hydrodynamic conditions, carbon source, etc [START_REF] Paul | Effect of shear stress and growth conditions on detachment and physical properties of biofilms[END_REF], in which these sessile communities grow. For example, porous biofilms with channels and voids between the finger-like or mushroom outgrowths are typical of a substrate-transport-limited regime instead compact and smooth biofilms occur when the biomass growth rate is limiting or the shear stress is high [START_REF] Van Loosdrecht | Biofilm structures[END_REF].

Mathematical modeling of biofilm development has been widely performed during the last decades. Biofilm models represent a perfect means to understand the basic principles determining biofilm formation, composition, structure and function [START_REF] Dr Noguera | Biofilm modeling: present status and future directions[END_REF] and therefore they can be used to effectively utilize and control biofilms in industrial and medical settings [START_REF] Picioreanu | Advances in mathematical modeling of biofilm structure[END_REF]. Mathematical models come in many forms that can range from very simple empirical correlations to sophisticated and computationally intensive algorithms that describe three-dimensional (3D) biofilm morphology and activity [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF]. The domain of interest is usually divided in three compartments: the bulk liquid, the boundary layer and the biofilm itself (Figure 2.1). All biofilm models simulate the dynamics of two types of components, particulate (active and inert biomass, EPS) and dissolved (substrates and metabolic products), and generally include three main elements: transport mechanism; consumption and growth mechanism; and loss mechanism [START_REF] Klapper | Finger formation in biofilm layers[END_REF]. The transport of dissolved compounds to the cells within the biofilm matrix is governed by diffusion. It plays a crucial role in biofilm development, since the concentrations of nutrients and products determine the rates of microbial reactions, and all the processes that generate an increase in volume are driven by nutrient availability [START_REF] Picioreanu | Modelling and predicting biofilm structure[END_REF]. Moreover, substrate concentration trends within the biofilm contribute to the formation of different environmental niches [START_REF] Stewart | Diffusion in biofilms[END_REF]. Biomass growth kinetics depend on substrate concentrations; bacteria, by consuming substrate, grow and duplicate, produce exopolymeric substances (EPS) determining an increase in the biofilm volume, usually called biomass spreading. Models proposed different approaches to describe the spreading of the newly formed amount of biomass, but the proper representation of biomass division and spreading still remains one of the most controversial topics [START_REF] Picioreanu | Modelling and predicting biofilm structure[END_REF]. Detachment is a determining factor for biofilm-structure formation [START_REF] Derlon | Growth limiting conditions and denitrification govern extent and frequency of volume detachment of biofilms[END_REF][START_REF] Morgenroth | Detachment -an often over looked phenomenon in biofilm research and modeling[END_REF]. It represents the primary process that balances microbial growth and, thereby, determines the steady state accumulation of the biofilm and the overall biofilm activity [START_REF] Picioreanu | Modelling and predicting biofilm structure[END_REF][START_REF] Stewart | A model of biofilm detachment[END_REF], and it greatly affects the performance and the stability of biofilm reactors [START_REF] Picioreanu | Two-dimensional model of biofilm detachment caused by internal stress from liquid flow[END_REF].

The review presented herein focuses on the description of the different modeling approaches used to represent particulate components and simulate biomass transport mechanism. Biofilm displacement is mainly caused by cell growth and division and EPS production, and can be affected by other processes changing the biofilm volume, such as attachment and detachment [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF]. Based on biomass representation, biofilm models have been classified in two main categories:

continuum models, which do not take directly into account the behavior of an individual microorganism as they treat biomass as a unicuum, based on population-averaged behaviour of different functional groups;

discrete models, which are generally defined bottom-up models, since biofilm structure is not furnished as an input to the model, but the complex morphology of biofilms emerges as a result of the actions and interactions of the biomass units with each other and the environment.

Each category is further subdivided by considering model dimensionality and the way in which diffusion and biomass spreading is treated. In particular, continuum models are classified in one-dimensional (1D) and multidimensional continuum models; discrete models are divided in cellular automaton (CA) models, hybrid differential-discrete CA models and Individual-based Models (IbMs). Continuum models treat the dynamics of biomass spreading by using differential equations, widely used in mechanics and transport phenomena. In discrete models, biomass spreading is assumed to be a stochastic process. This review work focuses on the description of the different modeling approaches present in literature and it is aimed at evaluating the main features of the analyzed biofilm models in order to enable readers to select an appropriate modeling tool based on their own needs.

Continuum models

As their name implies, continuum models consider the domain of interest as a continuum [START_REF] Bolea Albero | A new approach to the simulation of microbial biofilms by a theory of fluid-like pressure-restricted finite growth[END_REF] and biomass spreading as governed by differential equations. All continuum models are based on conservation laws which are formulated as balances of conserved properties (mass, volume, momentum, energy, etc.). For 1D models, these equations come in the form [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF]:

∂D ∂t + ∂J ∂z = R (2.1)
where: z is the space coordinate; t denotes the time variable; D represents a 1D property; R states the net 1D property production rate. Continuum models have undergone an evolution in terms of complexity: from 1D steady-state models, developed during the early 1970s, to multidimensional multispecies dynamic models that have been conceived during the last decade. This evolution has been influenced by the advances in computational and experimental tools and reflects the need of new biofilm models able to provide more complex two or three dimensional descriptions of microbial biofilms, in agreement with experimental observations. Based on their dimensionality, continuum models have been classified in two groups: 1D continuum models and multidimensional continuum models. 1D models consider only the direction perpendicular to the substratum while multidimensional models neglect the concepts of uniform thickness and layering of biomass typical of 1D models, and allow the biofilm matrix to expand in more than one direction.

One dimensional continuum models -Pioneer works

The first continuum biofilm models [START_REF] Atkinson | The overall rate of substrate uptake (reaction) by microbial films. part i-a biological rate equation[END_REF][START_REF] Williamson | A model of substrate utilization by bacterial films[END_REF] have been developed in 1970s in order to evaluate the substrate utilization kinetics in biofilms. These pioneer works were based on the concept that substrate removal from an aqueous phase requires diffusion of reactants into the biofilm, metabolism by microorganisms and diffusion of metabolic products through the biofilm and into the aqueous phase. These models can be considered the first example of continuum models since they were able to catch the essentials of biofilm development, idealizing the processes of substrate utilization, molecular diffusion and mass transport as simultaneous differential equations for a homogeneous layer of bacteria. In [START_REF] Williamson | A model of substrate utilization by bacterial films[END_REF], the authors adopted a schematic representation of the system where the biofilm is assumed to be attached to a flat surface with infinite length and width and characterized by a uniform cell density denoted X f and a locally uniform thickness L f . Substrate concentration within the biofilm changes only in the z direction, assumed perpendicular to the surface, and the rate of reaction is limited by a single substrate named rate-limiting substrate. The decrease in substrate concentration between the bulk liquid and the biofilm surface derives from an incomplete mixing of the liquid phase next to the biofilm surface coupled with mass transfer into biofilm and is modelled by the introduction of a liquid layer adjacent to and permeating the biofilm. In this layer the entire resistance to mass transport from the bulk liquid to the surface is concentrated. The depth of diffusion layer L, is defined as the equivalent depth of liquid through which the actual turbulent mass transport can be described by molecular diffusion alone [START_REF] Chaudhry | A review on the mathematical modeling of biofilm processes: advances in fundamentals of biofilm modeling[END_REF].

The model introduced in [START_REF] Williamson | A model of substrate utilization by bacterial films[END_REF] was able to couple the mass transport from the bulk liquid with the substrate biodegradation within the biofilm. In particular, the following steady-state differential equation has been used to describe substrate utilization within the biofilm:

D f ∂ 2 S f ∂z 2 = kX f S f K S + S f , 0 ≤ z ≤ L f , (2.2)
where:

D f is the molecular diffusivity in biofilms [L 2 T -1 ];
S f is the concentration of rate-limiting substrate at any point in biofilm [ML -3 ]; X f is the bacterial concentration within the biofilm, assumed constant with depth

[ML -3 ]; K S is the Monod half velocity coefficient [ML -3 ]; k is the maximum utilization rate of the rate-limiting substrate [T -1 ].
The rate of substrate utilization within the biofilm has been modelled by a Monod-like bacterial kinetics and the diffusion flux through the diffusion layer and the biofilm by the Fick's law of diffusion. However, this pioneer work did not include any considerations on the growth and decay of the bacteria composing the biofilm.

The model proposed in [START_REF] Williamson | A model of substrate utilization by bacterial films[END_REF] has been improved later by many researchers [START_REF] Be Rittmann | Model of steady-state-biofilm kinetics[END_REF][START_REF] Be Rittmann | Evaluation of steady-state-biofilm kinetics[END_REF][START_REF] Be Rittmann | Substrate flux into biofilms of any thickness[END_REF][START_REF] Be Rittman | The effect of shear stress on biofilm loss rate[END_REF][START_REF] Be Rittmann | Dual limitation of biofilm kinetics[END_REF][START_REF] Be Rittmann | The nonsteady-state-biofilm process for advanced organics removal[END_REF] who amended the basic model of mass transport in a steady state biofilm with additional processes. Rittman and McCarty [START_REF] Be Rittmann | Model of steady-state-biofilm kinetics[END_REF] incorporated the expressions for biofilm growth and decay for a steady-state biofilm, which is defined as a biofilm that for a given bulk liquid substrate concentration has neither net growth nor decay. Later, Rittmann [START_REF] Be Rittman | The effect of shear stress on biofilm loss rate[END_REF] introduced the biofilm loss rate caused by shear stress. This term has been formulated as a first order expression similar to the term used for decay losses.

One dimensional continuum models

With the ongoing progress in experimental methods, more sophisticated multisubstrate-multispecies models have been developed [START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF][START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] B D'acunto | Dynamic modeling of sulfate reducing biofilms[END_REF][START_REF] Kissel | Numerical simulation of mixedculture biofilm[END_REF][START_REF] Lee | One-dimensional mixed-culture biofilm model considering different space occupancies of particulate components[END_REF][START_REF] Rauch | A simplified mixed-culture biofilm model[END_REF][START_REF] Be Rittmann | Development and experimental evaluation of a steady-state, multispecies biofilm model[END_REF][START_REF] Be Rittmann | The transient-state, multiplespecies biofilm model for biofiltration processes[END_REF][START_REF] Tsuno | A simple biofilm model of bacterial competition for attached surface[END_REF][START_REF] Wanner | Competition in biofilms[END_REF][START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF]. These studies neglect the simplifying assumption of single-species biofilms and are mostly centered on the biofilm growth dynamics, including the biofilm thickness, the spatial distribution of microbial species and the substrate concentrations.

The 1D multispecies model of biofilm growth introduced by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Wanner | Competition in biofilms[END_REF] has been successfully applied to many biofilm studies since its development and represents a pioneer work in the understanding of the complex bulk interactions characterizing multispecies biofilms. This model takes into account the following processes: (a) the simultaneous substrate utilization and diffusion within the biofilm; (b) the external mass-transport resistance from the bulk liquid to the biofilm surface; (c) the growth of new biomass proportional to substrate utilization; (d) the biomass loss from endogenous respiration and detachment, and (e) the formation of inert biomass. The following equations have been introduced [START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF]:

∂X i ∂t + ∂ ∂z (uX i ) = ρ i r M,i (z, t, X, S) , i = 1, ... n, 0 ≤ z ≤ L (t) , t > 0, (2.3) ∂u ∂z = n i=1 r M,i = G (z, t, X, S) , 0 < z ≤ L (t) , t > 0, (2.4)
L (t) = u (L (t) , t) + σ a (t) -σ d (t) , t > 0, (2.5
)

∂S j ∂t - ∂ ∂z D j ∂S j ∂z = r S,j (z, t, X, S) , 0 < z < L (t) , t > 0, j = 1, ... , m, (2.6
) where: z is the biofilm growth direction assumed perpendicular to the substratum [L]; ρ i denotes constant density [ML -3 ]; X i (z, t) = ρ i f i denotes the concentration of microorganisms i, X = (X 1 , ..., X n ) [ML -3 ]; f i (z, t) is the volume fraction of microbial species i, n i=1 f i = 1; u(z, t) is the velocity of microbial mass [LT -1 ]; S j (z, t) denotes the concentration of substrate j, S = (S 1 , ..., S m ) [ML -3 ]; r M,i (z, t,X,S) is specific growth rate [ML -3 T -1 ]; L(t) denotes the biofilm thickness, free boundary [L]; σ a (t) is the attachment biomass flux from bulk liquid to biofilm [LT -1 ]; σ d (t) is the detachment biomass flux from biofilm to bulk liquid [LT -1 ];

D j denotes the diffusivity coefficient of substrate j [L 2 T -1 ]; r S,j (z, t, X, S) is the conversion rate of substrate j [ML -3 T -1 ].
Appropriate initial and boundary conditions are required to solve the previous system of nonlinear partial differential equations.

In particular, at the substratum-biofilm interface z = 0 the condition of no concentration gradient is assumed for both the soluble and the particulate components. Equation (2.3) is derived from the mass balance of the ith microbial species set up for a control volume. Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] modeled the spreading of biomass as an advective mass flux of each ith species. In particular, the authors assumed that when the net growth rate is positive in a control volume and the biomass density remains constant, the biomass increases giving rise to a flux of biomass that crosses the control-volume's boundary. Equation (2.4) determines the velocity at which the microbial mass is displaced with respect to the film-support interface. The value of u(z, t) is determined by the mean observed specific growth rate of the biomass and it is assumed identical for all species. Equation (2.5) defines the velocity at which the film-water interface moves; it depends on both the velocity at which the microbial mass is displaced, the velocity at which the biomass is exchanged between the biofilm and the bulk liquid and vice versa, here denoted as σ d (t) and σ a (t). In their work, Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] considered biomass loss due only to: shear stress, modeled setting σ d (t) = λL 2 with λ constant, and sloughing by setting σ d (t) as a δ Dirac function.

The model introduced by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] can be classified as a free boundary value problem which is very complicated to discuss due to the contemporary presence of hyperbolic and parabolic partial differential equations. In addition, when discussed numerically, it needs a special discretization scheme to consider the time-dependent change of the space domain [START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Szomolay | Analysis of a moving boundary value problem arising in biofilm modelling[END_REF][START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF]. The free boundary value problem contains two groups of nonlinear partial differential equations: the first system of n nonlinear hyperbolic partial differential equations describes the growth of microbial species in biofilms (2.3); the second group of m nonlinear parabolic partial differential equations governs the diffusion of substrates (2.6). The two systems are strictly connected as the biological processes they are aimed at modeling. The solution approach used by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] is based on a coordinate transformation that eliminates the moving boundary by introducing the space coordinate ζ(t) = z/L, which describes the distance from the substratum normalized by the biofilm thickness. This mathematical description of an idealized biofilm has been solved by a numerical solution technique based on the method of lines and has been addressed to some case studies. The existence of steady state solutions of this model applied to a single species biofilm has been proved later in [START_REF] Pritchett | Steady state solutions of a one-dimensional biofilm model[END_REF].

Contemporary to Wanner and Gujer, Kissel et al. [START_REF] Kissel | Numerical simulation of mixedculture biofilm[END_REF] formulated a multispecies biofilm model able to describe the competition between microbial species for common substrates within a completely mixed continuous-flow reactor. This model is based on the same continuum approach used by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] but does not include the loss of mass due to detachment. The two models differ in how they describe the net growth of biomass at any position in the biofilm and, consequently, in the numerical treatment adopted for the moving boundary problem. When the biomass increases while the biomass density is kept constant, Kissel et al. [START_REF] Kissel | Numerical simulation of mixedculture biofilm[END_REF] hypothesized that the control volume increases in size (Figure 2.2). Therefore, the numerical modeling of spatial variability in mass fractions and solute concentrations is accomplished by dividing the biofilm into a series of space elements with equal, but variable lengths. After each integration time step, the elements' lengths are recalculated, according to the volume expansion, or contraction obtained for the individual elements. All the equations have been solved numerically by using a fixed-step-size, fourth-order-accurate, Runge-Kutta technique. The model has been addressed only to dynamic state and no attempts have been made to solve equations at steady-state conditions. Figure 2.2: Spreading mechanism adopted in [START_REF] Wanner | A multispecies biofilm model[END_REF] and [START_REF] Kissel | Numerical simulation of mixedculture biofilm[END_REF] (figure adapted from [START_REF] Be Rittmann | Development and experimental evaluation of a steady-state, multispecies biofilm model[END_REF]).

Rittmann and Manem [START_REF] Be Rittmann | Development and experimental evaluation of a steady-state, multispecies biofilm model[END_REF] combined the multispecies biofilm model developed by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] with the steady-state assumption. For a steady-state multispecies biofilm, it is assumed that the growth of all species, deriving from substrate utilization, is equal to all losses. The model is addressed to simulate the competition for space in a multispecies steady-state biofilm and to predict the steady-state substrate fluxes, the biofilm thickness and the species distributions deriving from specific bulk-liquid substrate concentrations. The model contains a set of ordinary differential equations, similar to the mass balance equations derived by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF]. They are converted to be solved in a set of partial differential equations by introducing the pseudotime derivative, a means to perform the iterations required to achieve a correct steady-state solution.

The model introduced by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] was extended later in [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF], in order to simulate the effects of the following biofilm processes neglected in the previous model: advective transport of dissolved components and diffusive movement of particulate components in the biofilm, the development of the biofilm liquid phase volume fraction, the transport of suspended solids within the pore volume of the biofilm and the exchange of cells and particles between the solid matrix and the pore volume, and the simultaneous detachment and attachment to the biofilm surface. In [START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF], two new state variables have been introduced: l and θ. The first one represents the volume fraction of the liquid phase between the particulate components in the biofilm. The second one, also referred to as porosity, is introduced as the ratio of the volume between the biofilm solid matrix and the total biofilm volume. In [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF], the porosity and l represented the same quantity since the transport of suspended solids in the pore volume is neglected. The two variables are related by the following equations [START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF]:

l + nx i=1 P, Si = θ, (2.7)
θ + nx i=1 M, Si = 1, (2.8)
where: θ is the porosity;

l denotes the liquid phase volume fraction; P, Si represents the volume fraction of the solids suspended in the biofilm pore volume;

M, Si is the volume fraction of the biofilm matrix components.

In both [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF], particulate components are assumed to be transported not only by an advective flux, as stated in [START_REF] Wanner | A multispecies biofilm model[END_REF], but an effective diffusive flux of particulate components is introduced. This flux wants to describe the transport of cells and particles in the direction opposite to that of velocity u(z, t). It is independent from microbial growth and accounts for the mixing of cells or particles in the biofilm solid matrix as a result of mechanical deformation of the matrix by hydraulic forces or bioturbation. l is subject to an analogous advective flux since it is assumed that advective transport of particulate components does not change the ratio of liquid to solid phases in the biofilm. Moreover, to compensate the effective diffusive flux of particulate components, a flux of liquid phase in the opposite direction is introduced and a production rate for the liquid phase volume fraction in the biofilm is formulated. The development of the biofilm liquid phase volume fraction is modeled in order to simulate experimental data which show that in some cases porosity decreases from the biofilm surface to the substratum. The dissolved components are assumed to be transported in the liquid phase of the biofilm by a diffusive and advective flux, which is induced by a flux of water that derives from the transport of particulate components. The advective flux assumes a negligibly small value compared to the diffusive one. Moreover, the model [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF] also takes into account the simultaneously attachment and detachment of particulate components at biofilm surface, neglected in the original mixed-culture biofilm model in which only the dominant process was explicitly modeled. More precisely, in [START_REF] Wanner | A multispecies biofilm model[END_REF] the transport of cells and particles happens only towards the biofilm surface and as a consequence, new attaching particulate material can only adsorb at the biofilm surface, but cannot penetrate the biofilm. Modeling simultaneous attachment and detachment is possible only taking into account the diffusive transport of particulate components which reproduces the mixing of cells and particles over the biofilm depth. The partial differential equations introduced in [START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF] have been converted to a system of algebraic and ordinary differential equation and solved by the integration routines and numerical algorithms implemented in AQUASIM, a computer program designed for the identification and simulation of aquatic systems [START_REF] Reichert | Aquasim-a tool for simulation and data analysis of aquatic systems[END_REF]. A very detailed description of this simulation tool has been provided in [START_REF] Wanner | Biofilm modeling with aquasim[END_REF].

Later, Rauch et al. [START_REF] Rauch | A simplified mixed-culture biofilm model[END_REF] introduced a comprehensive simplified model, whose approach consisted in decoupling the modeling of the diffusion process and spatial distribution of bacteria species from the biokinetic reactions. This simplification derives from the need in faster, but sufficiently accurate predictions, avoiding the computational efforts for solving the partial differential equations. Diffusion is modeled as a steady-state phenomenon within each time step since substrate profiles are assumed to reach the equilibrium very fast. The typical concentration boundary layer is neglected. The system is divided in two compartments: bulk liquid and biofilm. The biofilm is constituted by a liquid phase in which dissolved substances are transported by molecular diffusion and a solid matrix which is constituted by several bacterial species, particulate substrate and inert material. The concentration of particulates and density in biofilm are expressed by equation (2.8). Substrates are transported inside the biofilm by molecular diffusion; when they do not penetrate the solid matrix, the reaction is considered as diffusion limited and takes place only over a certain depth of the biofilm. According to Harremoes [START_REF] Harremoes | Biofilm Kinetics[END_REF], the volumetric reaction rate is assumed to be zero-order respect to the concentration of substrate S in the biofilm and the penetration depth is derived from an analytical solution to the diffusion equation. The model is solved by using a two-step procedure: (1) for each conversion process that is influenced by diffusion, the active fraction of the biomass within the biofilm is computed by means of the analytical solution to the diffusion equation; (2) all the conversions within the biofilm are then calculated assuming the biofilm as an ideally mixed reactor but, with only the active fraction of the species contributing. The use of zero-order reaction rates is justified by the need of analytical solution for the substrate penetration depth that represents a basic concept for decoupling the diffusion and biokinetics reaction.

Despite their dimensionality, 1D biofilm models still represent an active topic in biofilm research area as proved by the recent contributions proposed in literature. In [START_REF] Lee | One-dimensional mixed-culture biofilm model considering different space occupancies of particulate components[END_REF], the authors introduced a 1D mixed-culture biofilm model based on the hypothesis that each particulate component has different space occupancy within the biofilm according to its fundamental nature, such as size and density. In this work space occupancy is not defined as the reciprocal of component density, as stated in [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF], but this feature also takes into account the liquid volume that coexists within the biofilm solid matrix. Internal porosity is calculated by the composition of the particulate components, which changes during biofilm growth. The model is based on the same mass balance equations introduced in [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF], but derived for the whole biofilm volume. The concept of effective diffusive transport is introduced and the model has been successfully applied to simulate the consolidation phenomenon. The partial differential equations have been solved by converting them into a system of ordinary differential equations in a dimensionless form, later solved by using the ode15s tool provided in MATLAB software.

Rittmann et al. [START_REF] Be Rittmann | The transient-state, multiplespecies biofilm model for biofiltration processes[END_REF] reported a transient multispecies biofilm model (TMSBM) especially focusing on the kinetics of the growth related microbial products. This model represents a synthesis of the key modeling features used to describe multispecies biofilms and is addressed to biofilms that experience time-varying conditions, particularly including periodic detachment by backwashing. The TMSBM contains non steady-state mass balances for each of the four types of biomass represented (2.3) and for the soluble species in a layer of biofilm (2.6). Similarly to [START_REF] Wanner | A multispecies biofilm model[END_REF], it is assumed that the spreading of the biomass gives rise to a biomass-flux which involves a biomass velocity. This velocity represents the physical rate at which all types of biomass move into or out of a biofilm layer by crossing the layer's boundaries. The partial differential equations constituting the model are solved by using a separate-solution strategy. The calculations for the soluble species are separated from the calculations for the biomass species in order to avoid the accumulation of rounding errors. Substrate and product calculations are performed for fixed biomass distributions, while substrate profiles are kept constant when biomass calculations are performed. This strategy allows to reduce the computing load. To solve the biomass balance equations the TMSBM is based on a hybrid strategy between the approach used in [START_REF] Wanner | A multispecies biofilm model[END_REF] and [START_REF] Kissel | Numerical simulation of mixedculture biofilm[END_REF]. This new method contemporary allows the biomass flux and a change in the layer size. The biofilm is divided in layers having the same size all over the time and the sum of growth and decay for all the biomass species in all layers indicates the overall net growth. Moreover, the model includes a novel accounting scheme that takes into account the net growth of each species in a layer relative to the net growth of all species within that layer and within the adjacent layers.

Recently, D'Acunto and Frunzo [START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF][START_REF] B D'acunto | Dynamic modeling of sulfate reducing biofilms[END_REF] have transformed the partial hyperbolic differential system (2.3) introduced in [START_REF] Wanner | A multispecies biofilm model[END_REF]] into an integral system by using a characteristic-like method, where the characteristics are the lines z = s(z 0 , t) defined by: ∂s ∂t (z 0 , t) = u (s (z 0 , t) , t) , s(z 0 , t) = z 0 0 ≤ z 0 ≤ L 0 .

(2.9)

The same method has been applied later for a qualitative analysis of the attached cell layer in multispecies biofilm formation [START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF]. Compared to the free boundary problem introduced in [START_REF] Wanner | A multispecies biofilm model[END_REF], this biological process is described by a free boundary problem for nonlinear hyperbolic equations where the initial biofilm thickness is set to zero. In this case, the free boundary is represented by a space-like line. An existence and uniqueness theorem of solution to the systems (2.3) -(2.6) has been proved by the fixed point theorem [START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF]. The method of characteristics has been used also for numerical purposes to simulate the dynamics of multispecies biofilms [START_REF] B D'acunto | Dynamic modeling of sulfate reducing biofilms[END_REF][START_REF] Mr Mattei | Mathematical modeling of competition and coexistence of sulfate-reducing bacteria, acetogens, and methanogens in multispecies biofilms[END_REF].

Klapper and Szomolay [START_REF] Klapper | An exclusion principle and the importance of mobility for a class of biofilm models[END_REF] have demonstrated by an exclusion principle that Wanner-Gujer model [START_REF] Wanner | A multispecies biofilm model[END_REF], evaluated under steady-state conditions, leads to restrictions on ecological structure since it neglects downward microbial motility. The introduction of a diffusion flux for motile species [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Reichert | Movement of solids in biofilms: significance of liquid phase transport[END_REF] may be able to negate conditions leading to the exclusion principle. In a recent contribution, D'Acunto et al. [START_REF] B D'acunto | Modeling multispecies biofilms including new bacterial species invasion[END_REF] have been able to take into account the biological process of colonization of new species and transport from bulk liquid to biofilm (or vice-versa) keeping the equations hyperbolic.

Parallel to the development of 1D dynamic models, more complex steady-state models have been developed during the last decade. Perez et al. [START_REF] Pérez | Modeling biofilm and floc diffusion processes based on analytical solution of reaction-diffusion equations[END_REF] developed a biofilm model based on the assumption of zero and first order kinetics in biofilms to compute fluxes of substrate into the biofilm. More precisely, the approach used is based on the weighted average of the analytical solutions for first and zero-order reaction kinetics. Compared to numerical models, the use of analytical solutions, in addition to simplicity, allows the analysis of the effects that each term, or parameter can have on the overall flux. Beyenal and Lewandowski [START_REF] Lewandowski | Modeling mass transport and microbial activity in stratified biofilms[END_REF] introduced a model able to reproduce biofilm heterogeneity by using a 1D continuum approach. The biofilm is modeled as composed by a finite number of layers characterized by different nutrient concentration, effective diffusivity and density. Each of these layers is modeled as a uniform biofilm and the effective diffusivity is recognized as the control parameter for space discretization. The model is aimed at quantifying mass transfer in stratified biofilm and comparing the results with a homogeneous biofilm model. The effective diffusivity is expressed as a linear function of space coordinate z and biofilm density; its gradient is used to append the equation quantifying mass transfer in homogeneous biofilms by a factor representing biofilm heterogeneity. More recently, Gonzo et al. [START_REF] Ee Gonzo | Continuum heterogeneous biofilm model-a simple and accurate method for effectiveness factor determination[END_REF][START_REF] Ee Gonzo | The continuum heterogeneous biofilm model with multiple limiting substrate monod kinetics[END_REF] developed a new approach to model steady-state activity of heterogeneous biofilm. The main difference with the work of Beyenal and Lewandowski [START_REF] Lewandowski | Modeling mass transport and microbial activity in stratified biofilms[END_REF] is that the new approach does not require numerical simulations.

Multidimensional continuum models

The development of multidimensional continuum models is relatively recent and reflects the need of reproducing the complex morphology of biofilms, which arises from the interaction with the surrounding liquid and the dynamics of transport and consumption of substrates. These models are amenable to mathematical analysis and do not rely on ad hoc rules to simulate growth processes.

As stated in [START_REF] Hj Eberl | A finite difference scheme for a degenerated diffusion equation arising in microbial ecology[END_REF], in the 1D case, dynamic biofilm models have been mostly formulated as free boundary value problems and were based on the assumption that newly produced biomass is converted into new biofilm volume which moves according to a convective transport mechanism. The increasing biofilm thickness and the speed of propagation of the biofilm/liquid interface normal to the substratum can be calculated from the production terms by integration over the biofilm thickness. In the multidimensional case, this approach requires the introduction of an evolution equation for the convective biomass transport velocity. This equation can be derived by introducing the idea that biofilm growth generates a pressure field within the biofilm, which is responsible for the establishment of a spreading velocity. Therefore, a further unknown variable (pressure) has to be modeled. To solve this issue different modeling approaches have been introduced: the description of the biofilm as a rigid/elastic/viscoelastic Figure 2.3: 2D continuum models: a) spreading mechanism adopted by Eberl et al. [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF]; b) spreading mechanism adopted by Dockery and Klapper [START_REF] Klapper | Finger formation in biofilm layers[END_REF]. solid or highly viscous fluid has been investigated and the idea of estimating material properties of the biomass and incorporating them in model equations has been taken into account in order to better understand biofilm structural stability.

The first attempt to model biofilm growth as a convective transport mechanism has been introduced in [START_REF] Wood | Diffusion and reaction in biofilms[END_REF][START_REF] Wood | Cellular growth in biofilms[END_REF] where the authors developed a macroscopic description of microbial growth by using the sub-cell-scale information of mass transport and intracellular reactions. However, the mechanistic problem arising from the calculation of the convective field has been solved only by introducing an empiricism which requires the experimental determination of a growth coefficient.

An alternative approach to the convective transport mechanism has been introduced in [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF]. In this work, the authors developed a spatio-temporal continuum model in which the biomass spreading is described by a nonlinear density-dependent diffusion mechanism (Figure 2.3).

The model is aimed at describing hydrodynamics, transport and consumption of nutrients and biomass production for a single species biofilm. Biomass spreading occurs only when the biomass density m(t, x) reaches a known a priori maximum value (waiting time behavior); as a consequence a density-dependent diffusion coefficient is introduced. The system is divided in two regions separated by an interface Γ : Ω 1 represents the liquid region and Ω 2 is the solid biofilm region. The distinction between the liquid region Ω 1 and Ω 2 is made by the biomass density m(x, t) = 0 or m(x, t) > 0, respectively. The model is governed by the following equations [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF]:

∇u = 0, ∂u ∂t + u • ∇u = - 1 ρ ∇p + ∇ 2 u, in Ω 1 = {x Ω|m (t, x) = 0} , (2.10)
u = 0, in Ω 2 = {x Ω|m (t, x) > 0} , (2.11
)

∂c ∂t +u•∇c = ∇• (d 1 (m) ∇c)-f (c, m) in Ω 2 = {x Ω|m (t, x) > 0} , (2.12) ∂m ∂t = ∇ • (d 2 (m) ∇m) + g (c, m) in Ω 2 = {x Ω|m (t, x) > 0} , (2.13) f (c, m) = k 1 cm k 2 + c , g(c, m) = k 3 (f (c, m) -k 4 m) , (2.14
d 2 (m) = m max -m a m b (2.15)
The equations (2.10)-(2.15) solve both the hydrodynamics and the biofilm evolution; they are strictly connected since the regions Ω 1 and Ω 2 both depend on m(t, x). Biomass production is assumed to be established only by reaction kinetics and the biomass diffusivity is assumed to vanish as m(t, x) becomes small, but it increases as m(t, x) grows thanks to biochemical reactions (2.12)-(2.15). The model (2.10)-(2.15) is mathematically rather complicated and difficult to handle analytically. To solve it, the authors assumed hydrostatic state, since the major difficulties of the model derive from the Navier-Stokes equations (2.10) and introduced dimensionless dependent variables. In this way the model reduces to a spatio-temporal predator prey model for biomass and nutrients. The model behavior has been validated only by numerical simulations carried out in 1D and 3D by using a finite difference scheme which is solved explicitly for the slower biomass spreading process and implicitly for the faster nutrient transport process. The numerical analysis is aimed at showing the sensitivity of the biofilm behavior to crucial parameters and confirms that model results are in good agreement with previous experimental and modeling experience.

Further analysis and application of the model introduced in [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF] have been performed in [START_REF] Hj Eberl | A finite difference scheme for a degenerated diffusion equation arising in microbial ecology[END_REF][START_REF] Duvnjak | Time-discretisation of a degenerate reactiondiffusion equation arising in biofilm modeling[END_REF][START_REF] Hj Eberl | A transient density dependent diffusion-reaction model for the limitation of antibiotic penetration in biofilms[END_REF][START_REF] Ma Efendiev | Existence and longtime behavior of solutions of a nonlinear reaction-diffusion system arising in the modeling of biofilms[END_REF]. In these studies, the authors mainly focus on the evaluation of different numerical schemes able to handle the diffusion singularity effects arising from equation (2.15). Recently, Eberl and Sudarsan [START_REF] Hj Eberl And R Sudarsan | Exposure of biofilms to slow flow fields: The convective contribution to growth and disinfection[END_REF] have extended the degenerated diffusion-reaction model for biofilm growth and disinfection introduced in [START_REF] Hj Eberl | A transient density dependent diffusion-reaction model for the limitation of antibiotic penetration in biofilms[END_REF] to account for the convective transport of substrates in the bulk liquid. A thin-film approximation to the Navier-Stokes equations was considered so that the fluid flow could be solved analytically. However this simplification can be used only for slow flows in narrow channels. The concept of biofilm as a homogeneous, viscous, and incompressible fluid of constant density, satisfying Darcy's law has been firstly introduced in [START_REF] Klapper | Finger formation in biofilm layers[END_REF]. In this pioneer work, Dockery and Klapper introduced an equation that regulates the state variable p (pressure) in the biofilm phase:

λ∇ 2 p + g (u(S)) = 0, (2.16)
where:

p is the osmotic pressure [ML -1 T -2 ]; λ is the Darcy constant [T L 3 M ];
S is the concentration of the rate-limiting substrate

[ML -3 ]; g is a prescribed growth function [T -1 ]; u is the substrate uptake rate [ML -3 T -1 ].
The pressure equation is solved in the biofilm region after setting specific boundary conditions: the aqueous region is supposed to be static near the biofilm surface and so a constant pressure is assumed for the bulk liquid (Figure 2.3). Equation (2.16) is coupled with the solution of a nutrient diffusion-reaction mass balance, which provides the field of concentration S. The substrate is assumed to diffuse through the bulk region into the biofilm, where it also spreads and is consumed. The model does not take into account the internal chemical signaling for biofilm growth and behavior and the influence of fluid dynamics. The equations have been solved numerically on a uniform two-dimensional (2D) rectangular grid. The biofilm-bulk liquid interface evolution is tracked by using the level set method.

In a further work, Klapper et al. [START_REF] Klapper | Viscoelastic fluid description of bacterial biofilm material properties[END_REF] investigated biofilm material properties by using a system of viscoelastic fluid equations, coupled to a linear Jeffrey's viscoelastic stress-strain law. More precisely, the standard fluid equations of momentum (2.17) and continuity (2.18) have been simplified under the hypotheses of quasi-static evolution of phenomena (slow growth in the biofilm matrix), high viscosity and divergence-free velocity field in the biofilm compartment: ∇σ = 0, (2.17)

∇(ρu) = 0, (2.18)
where:

σ is the stress tensor for the biofilm matrix [ML -1 T -2 ]; ρ is the biofilm matrix density [ML -3 ; u is the biofilm matrix velocity field [MT -1 ].

The system of equations is based on elements of associated polymer physics and results consistent with experimental results of biofilm deformation, failure and detachment in response to mechanical forces.

The work of [START_REF] Klapper | Effect of heterogeneous structure in mechanically unstressed biofilms on overall growth[END_REF] represents the sequel of [START_REF] Klapper | Finger formation in biofilm layers[END_REF] and examines the formation of biofilm fingers and mushrooms. According to [START_REF] Klapper | Viscoelastic fluid description of bacterial biofilm material properties[END_REF], the hypothesis of biofilm as a viscoelastic fluid is adopted, but the analysis is restricted to the case of static or nearly static bulk fluid. Therefore, the substrate is assumed only to diffuse from bulk liquid into the biofilm and the shear stress and the associated viscoelastic response is not considered. Under the hypothesis of incompressibility of the biofilm matrix and the assumption of Darcy's law for the biofilm interface velocity, the continuity equation reduces to equation (2.16). According to [START_REF] Klapper | Finger formation in biofilm layers[END_REF], the biofilm-liquid interface z = h(x, y, t) evolution follows equation (2.19) [START_REF] Klapper | Effect of heterogeneous structure in mechanically unstressed biofilms on overall growth[END_REF]:

dh dt = -(∇p • n) (ẑ • n) (2.19)
The performed nonlinear analysis suggests that in the case of biofilms free of external mechanical stress, overall growth is inhibited by the presence of growing perturbations in the linear stage. A generalization of the previous 2D model [START_REF] Klapper | Finger formation in biofilm layers[END_REF] and of the earlier 1D model [START_REF] Wanner | A multispecies biofilm model[END_REF] has been proposed by Alpkvist and Klapper [START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF] who developed a continuum model for the heterogeneous growth in biofilm systems with multiple species and multiple substrates. This model represents the early work of the rigorous mathematical treatment of continuum multidimensional multispecies biofilm models and it is based upon a combination of the approach introduced in [START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Klapper | Finger formation in biofilm layers[END_REF]. The domain is subdivided in two regions: the biomass region B t and the liquid region Ω Bt where Ω is defined as an open sunset of R 3 . The domain has two moving boundaries: the biomass-liquid interface defined by the curve Γ t and the bulk liquid interface at a fixed height Γ Hb above B t , defined by the curve Γ Hb . In the region above the curve Γ Hb , fluid mixing is able to replenish or remove diffusive components faster than they are used or produced.

The model takes into account N b different components or phases for the biomass region and N c different substrates. The model consists of a series of partial differential equations derived on the basis of conservation laws and reaction kinetics. As in [START_REF] Wanner | A multispecies biofilm model[END_REF], the transport of biomass is governed by an advective process characterized by a volumetric flow u(t, x) equal for all species. According to [START_REF] Klapper | Finger formation in biofilm layers[END_REF], the biofilm is modeled as a homogeneous, viscous, incompressible fluid with a velocity given by Darcy's law [START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF]:

u = -λ∇p, (2.20)
where:

p = p(t, x) is the pressure [ML -1 T -2 ]; λ is the Darcy constant [T L 3 M ].
The model is based on semilinear Poisson equations for substrate concentrations (2.21), linear Poisson equation for pressure (2.22), and advective equations for the biomass volume fractions (2.23) [START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF]:

-D j ∇ 2 C j = r j j = 1, ..., N c , (2.21)
-∇ 2 p = N b i=1 g i ρ * i i = 1, ..., N b , (2.22
)

∂ϑ i ∂t -∇p • ∇ϑ i = g i ρ * i - N b i=1 g i ρ * i i = 1, ..., N b , (2.23)
where:

C j is the substrate concentration [ML -3 ]; D j is the assumed constant substrate diffusivity [L 2 T -1 ]; r j is the substrate uptake rates [ML -3 T -1 ]; g i is the biomass growth or loss rate [ML -3 T -1 ];
ρ * i is the individual density for biomass components, assumed to be constant in time and space [ML -3 ]; θ i = θ i (t, x) denotes the volume fraction of the ith species.

Applied to a planar biofilm system the model reduces to a 1D model equivalent to Wanner and Gujer system. Model simulations have been based upon accepted numerical methods with an existing error analysis. In particular, the time evolution of the biomass region is calculated by using a level set function as in [START_REF] Klapper | Finger formation in biofilm layers[END_REF]. The model has been used to simulate in 2D and 3D biofilm growth in growth-limited and transport-limited regimes.

A continuum model for a single species bacterial colony growth has been introduced in [START_REF] Alpkvist | A new mathematical model for chemotactic bacterial colony growth[END_REF]. The biomass density is supposed to assume zero outside and a distinct value inside the biomass region. Therefore, the biofilm growth problem is solved by tracking the boundary of the biomass domain. In this case, the equation for the time-varying velocity field has been derived by considering a chemotactic growth, i.e. the flow of biomass has been assumed to move in the direction of the increasing nutrient concentration. This equation is formulated such that the equation of continuity for biomass density holds and it is coupled to the reaction-diffusion equation for substrates. According to [START_REF] Klapper | Finger formation in biofilm layers[END_REF], the equations have been solved by using the level set method for what concerns the moving boundary of biofilm region and a second-order finite difference scheme for the substrate field.

Another deterministic approach to model biofilm growth has been derived from material mechanics [START_REF] Hj Dupin | Pore-scale modeling of biological clogging due to aggregate expansion: A material mechanics approach[END_REF]. The biofilm is modeled as a continuous, uniform, isotropic, and hyper-elastic material, whose expansion and deformation are governed by material stress-strain relations. The density is kept constant by deforming the biofilm matrix; this means that the pressure generated by cell division has to meet the resistance of the EPS matrix surrounding microbial cells.

The work of Alpkvist and Klapper [START_REF] Alpkvist | Description of mechanical response including detachment using a novel particle model of biofilm/flow interaction[END_REF] represents one of the first attempts made to include the liquid flow hydrodynamics in a biofilm model. The authors focus mainly on the interaction between fluid flow and biomass and therefore neglect biofilm growth. Biofilm is characterized as a viscous body embedded in a viscous fluid in a 3D domain and the whole system is governed by the Navier-Stokes equations of momentum and incompressibility-induced mass balance:

∂u ∂t + u • ∇u = - 1 ρ ∇p + ν∇ 2 u + 1 ρ F, (2.24)
∇ • u = 0, (2.25)
where:

p = p(t, x) is the osmotic pressure [ML -1 T -2 ]; u = u(t, x) is the flow field [LT -1 ]; F = F(t, x)
is the internal viscoelastic force density within the fluid-biofilm system

[ML -2 T -2 ]; ρ is the density of water [ML -3 ]; ν is the kinematic viscosity of water [M 2 T ].
The macroscopic behaviour of biofilms is coupled to the flow of biomass and liquid through the function F = F(t, x) in equation (2.24). The equations have been solved by using the Immersed Boundary Method which allows the authors to evaluate the behavior of the whole system by directly solving the Navier-Stokes equations.

Further works have tried to reproduce both biofilm growth and liquid hydrodynamics; they are based on simplifying assumptions for what concerns Navier-Stokes equations and make use of specific methods to solve the partial differential equations on irregular shape domains. One example is furnished by the work of Duddu et al. [START_REF] Duddu | A combined extended finite element and level set method for biofilm growth[END_REF] who proposed a continuum model to estimate substrate concentration, biomass advection velocity and biomass volume fraction and extended later the model to fluid flow velocity field calculation [START_REF] Duddu | A two-dimensional continuum model of biofilm growth incorporating fluid flow and shear stress based detachment[END_REF]. The biofilm is characterized as a homogeneous isotropic elastic material constituted by two components, the active and the inactive biomass, while the fluid is assumed to behave as Newtonian with constant viscosity and in laminar flow. The fluid flow and the stress deformation problems are uncoupled under the hypothesis of small stress induced deformation [START_REF] Picioreanu | Two-dimensional model of biofilm detachment caused by internal stress from liquid flow[END_REF].

The biofilm growth is supposed to be irrotational; therefore the velocity field is derived from a potential function. The system of partial differential equations governing the fluid hydrodynamics, substrate transport at steady-state, the mass balance for total biomass written in terms of the growth velocity potential and the mass balance equation for active biomass have been solved by using the extended finite element method while the location of the biofilm/fluid interface is evaluated by the level set method.

In [START_REF] Cumsille | A novel model for biofilm growth and its resolution by using the hybrid immersed interface-level set method[END_REF] a Hele-shaw type-like modeling is introduced: the modeled system is assumed to be composed by two fluids (biofilm and liquid) characterized by different viscosities and separated by a moving interface. The velocity field in the liquid compartment is assumed to be divergence-free due to incompressibility while the velocity field in the biofilm compartment is not divergence-free and equation (2.16) is assumed to hold. Compared to [START_REF] Klapper | Finger formation in biofilm layers[END_REF], this work solves the pressure equation in the entire domain by imposing transmission conditions on the biofilm/bulk liquid interface, including the effects of the fluid motion induced by the evolution of the biofilm/liquid interface and accounts for advective substrate transport in and out of the biofilm. The mathematical problem has been solved by coupling the immersed interface method with the Level-Set method. Similarly, Cogan [START_REF] Ng Cogan | Two-fluid model of biofilm disinfection[END_REF] developed a model in which the biofilm is treated as a viscous fluid immersed in a fluid of lower viscosity but the equations governing the system are assumed to be Stokes type-like due to the low values of Re. Other works deal with the biofilm treated as a biological gel composed of EPS and water and in which the bacteria are immersed [START_REF] Ng Cogan | The role of the biofilm matrix in structural development[END_REF] or investigate the biofilm by using a phase-field model in 1D and 2D respectively [START_REF] Zhang | Phase field models for biofilms. i. theory and one-dimensional simulations[END_REF][START_REF] Zhang | Phase field models for biofilms. ii. 2-d numerical simulations of biofilm-flow interaction[END_REF]. Clarelli et al. [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF] have recently introduced a fluid dynamics model based on the mixture theory which considers the biofilm as a multiphase fluid. In contrast with most of the existing models, this work considers a finite speed of propagation for the hyperbolic equations.

Discrete Models

Discrete models for biofilm research started to be developed in 1990s. Biofilms are assumed to be living systems inherently stochastic and researchers have been devising ways to express this stochasticity [START_REF] Cs Laspidou | Cellular-automata and individual-based approaches for the modeling of biofilm structures: Pros and cons[END_REF]. They use approaches where the large-scale dynamics are emergent from the processes occurring at a small-scale and are generally defined bottom-up models since biofilm structure is not furnished as an input to the model but the complex morphology of biofilms emerges as a result of the actions and interactions of the biomass units with each other and the environment. The rules used to model interactions at a local level can be motivated purely from biological principles, instead of analysis from a mathematical and physical framework [START_REF] Alpkvist | Threedimensional biofilm model with individual cells and continuum eps matrix[END_REF]. The basic idea consists in splitting the biomass accumulation and transport in two separate processes: the biomass growth kinetics are still governed by ordinary differential equations as for the continuum models, while the biomass transport mechanism is realized in a discrete way. Therefore, discrete models have been classified in three groups based on biomass representation and spreading mechanism adopted:

-CA models; -Hybrid differential-discrete CA models, in which the mass transport is described by using differential equations while the biofilm structure development is treated by using a CA approach; -IbMs.

In CA models, the biomass is represented in an array of small compartments (usually rectangular), as opposed to the agent-based representation of the IbM that uses particles located anywhere in space and characterized by essential state variables like cell mass and cell volume. As well as continuum models, CA models use volume averaging properties (density or concentration) as state variable for the biomass and are so called biomass-based models. Wood and Whitaker [START_REF] Wood | Diffusion and reaction in biofilms[END_REF][START_REF] Wood | Cellular growth in biofilms[END_REF] have provided an analysis of conditions under which biomass averaging is a valid computational tool. The two groups differ also on the biomass spreading rules used (see Sections 2.3.1, 2.3.2, 2.3.3).

Cellular Automaton Models

Modeling biomass growth and spreading has been widely performed by using a CA approach [START_REF] Barker | A cellular automaton model of microbial growth[END_REF][START_REF] Chang | A three-dimensional, stochastic simulation of biofilm growth and transport-related factors that affect structure[END_REF][START_REF] Rl Colasanti | Cellular automata models of microbial colonies[END_REF][START_REF] Sw Hermanowicz | A model of two-dimensional biofilm morphology[END_REF][START_REF] Sw Hermanowicz | Two-dimensional simulations of biofilm development: effects of external environmental conditions[END_REF][START_REF] Sw Hermanowicz | A simple 2d biofilm model yields a variety of morphological features[END_REF][START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF][START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF][START_REF] Pizarro | Quantitative cellular automaton model for biofilms[END_REF][START_REF] Ge Pizarro | Two-dimensional cellular automaton model for mixed-culture biofilm[END_REF][START_REF] Tang | An improved cellular automaton method to model multispecies biofilms[END_REF][START_REF] Wimpenny | A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models[END_REF]. CA models were originally developed for the Game of Life (conceived by the mathematician John Horton Conway in 1970) and were based on simple rules for building complex structures from simple and repetitive elements [START_REF] Wimpenny | A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models[END_REF]. In particular, the basic idea consists in miming the physical laws by a series of simple rules that are easy to compute quickly and in parallel. More properly, a CA model consists of a simulation which is discrete in time, space and state [START_REF] Gb Ermentrout And L Edelstein-Keshet | Cellular automata approaches to biological modeling[END_REF]. Usually the model space is discretized in a grid of rectangular elements (often squares in 2D or cubes in 3D). Each grid element has four first-order neighbors and another four second-order neighbors in the 2D rectangular space discretization [START_REF] Picioreanu | Modelling and predicting biofilm structure[END_REF]. The grid cell is allowed to fill up to a predetermined maximum and a simple rule-based system is employed to locate the extra biomass in a new compartment [START_REF] Cs Laspidou | Cellular-automata and individual-based approaches for the modeling of biofilm structures: Pros and cons[END_REF]. Substrate diffusion is usually simulated by random walks of individual substrate particles; while biofilm growth is described as the multiplication of individual microbial cells when they consume substrate particles.

CA models can be divided in three classes: (1) deterministic or Eulerian automata; (2) lattice gas models; and (3) solidification models [START_REF] Gb Ermentrout And L Edelstein-Keshet | Cellular automata approaches to biological modeling[END_REF]. For the first class of models, the spatial domain is divided into a fixed lattice and each lattice point has a state associated with it. The state at the next step is determined by earlier states of the cell and its neighbor. This type of CA model reproduces evolution equation with a partial differential equation or an integral equation. Lattice gas models are called particle systems and consist of a discrete spatial grid on which particles move and interact in some prescribed fashion. Solidification models resemble lattice gas models except for the concept of bound state.

The Diffusion-Limited Aggregation (DLA) models represent the first attempt to model bacteria colonies structures using a discrete approach [START_REF] Fujikawa | Diversity of the growth patterns of bacillus subtilis colonies on agar plates[END_REF][START_REF] Fujikawa | Fractal growth of bacillus subtilis on agar plates[END_REF][START_REF] Matsushita | Diffusion-limited growth in bacterial colony formation[END_REF][START_REF] Tolman | Cluster-size distribution in the incremental growth of dla clusters[END_REF][START_REF] Ta Witten | Diffusion-limited aggregation, a kinetic critical phenomenon[END_REF]. These models are based on the same grid system as standard CAs, but the array contains particles that can move between the squares in a prescribed fashion. They are based on an analogy between crystal growth and biofilm accumulation. In particular, DLA models assume that both crystallization and biofilm formation are driven by the mass transfer of some essential dissolved compounds from bulk liquid to a solid surface. These models are mostly focused on the important role played by concentration gradients in the growth mechanism of bacterial colonies. The biofilm growth is assumed to be determined by the deposition of new layers of material on an existing surface. Dissolved matter has to diffuse through boundary layers; when it reaches a reactive surface, a surface reaction transforms it in solid phase. The basic idea of these models consists in choosing a seed particle as the origin of a square lattice on a plane. Biofilm growth occurs when another particle, released far from the origin and allowed to move randomly, arrives at the nearest neighboring site to the origin and sticks to the site. Later, these two particles are frozen in this position and another particle is released. Repeating this procedure, the cluster grows assuming in many cases an open and branched structure. DLA models are based on the simplifying assumption that nutrients diffuse only across a liquid boundary layer; actually nutrients diffuse also into the biofilm, leading to the appearance of a reaction zone in the bulk biofilm. This means that biofilm does not grow only at the surface but also in volume and the expansion of the solid-liquid biofilm interface is caused by internal pressure generated by the growing biomass. DLA models have been applied to simulate the growth of bacterial colonies both at very low nutrient level on an agar plate and under higher nutrient concentrations. Although the shapes of DLA patterns may resemble those of certain bacterial colonies , the biological mechanism is clearly distinct since cells are added through division of nearby cells.

Later, Wimpenny and Colasanti [START_REF] Wimpenny | A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models[END_REF] have developed a model that adds biological rules to DLA models. The stationary particle used in DLA models is replaced by a microbial cell. This microbial cell can occupy a single square and can produce copies of itself that will occupy neighboring squares. The cells consume resource units that can randomly diffuse over a predetermined range of neighboring compartments. In this model, growth occurs only if there is available free space in the neighborhood of the cell. This mechanism generates growth only in the outermost cell layer, just like in crystal formation and neglects any growth occurring inside the biofilm matrix. Moreover, the model does not take into account the conservation laws of the substrate amount converted into biomass. Despite these shortcomings, the model was able to demonstrate how changings in the concentration of a rate-limiting substrate can cause different morphology varying from dense structure to biofilms penetrated by water channels.

A quantitative CA model for homogeneous biofilms has been introduced in [START_REF] Pizarro | Quantitative cellular automaton model for biofilms[END_REF]. The main objective of this work was to link the stochastic CA parameters with the physical and kinetic parameters used in biofilm modeling in order to obtain quantitative predictions of macroscale activity. The CA model is presented as constituted by six different elements: lattice, cell, states, time, rules and neighborhood. The biofilm system is divided into two lattices: the first describing the spatial location of food particles, the second describing the spatial location of the microbial particles that constitute the biofilm. In the substrate lattice, cells are subdivided into layers that represent the possible directions of displacement of substrate particles during diffusion. Each layer in the substrate lattice can have one of two states, describing absence (0) or presence (1) of food. The number of food particles in a local neighborhood defines the concentration of substrate at that location. In the microbial lattice, the cells can assume three different states, absence (0), presence of one microbial particle (1), or presence of two microbial particles [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF]. The latter state describes the situation right after a reproduction event. The information on each lattice is updated at discrete time intervals. The dynamics of this update are governed by the CA rules, which represent the interaction of each cell with its neighborhood of cells with the corresponding cells in the superimposed lattice. Each rule represents the application of most important processes occurring in biofilm, namely diffusion, substrate utilization, bacterial growth, bacterial decay, and microbial distribution. The rules are applied to the two lattices in a sequential manner. Substrate diffusion is modeled by a random movement of food particles in the lattice. This movement is simulated in two steps, mixing and transport. Substrate utilization is modeled by introducing the probability that during a time interval, a microbial particle will consume a food particle. Microbial growth is modeled according to [START_REF] Wimpenny | A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models[END_REF]. The probability of a microbial particle disappearing from the lattice at a given time step is evaluated by a first order coefficient which takes into account microbial decay and detachment. After the growth and decay steps, the microbial is updated according to the biomass distribution rules: 1) conversion of a cell with two microbial particles into two cells with one microbial particle each; 2) elimination of the empty cells. The CA approach introduced by Pizarro et al. [START_REF] Pizarro | Quantitative cellular automaton model for biofilms[END_REF] has been applied later to incorporate the formation of inert biomass within a structurally heterogeneous multi-species biofilm, the mechanisms of inert biomass decay and to include a self-organizing development of the biofilm structure [START_REF] Ge Pizarro | Two-dimensional cellular automaton model for mixed-culture biofilm[END_REF].

Hybrid differential-discrete cellular automaton models

Hybrid differential-discrete CA are a class of models in which nutrient diffusion is modeled by using a differential equation usually assumed at a steady-state with respect to the bacterial growth, while the biomass spreading is treated by CA rules [START_REF] Picioreanu | Two-dimensional model of biofilm detachment caused by internal stress from liquid flow[END_REF][START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF][START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF][START_REF] Hunt | A computer investigation of chemically mediated detachment in bacterial biofilms[END_REF][START_REF] Laspidou | Modeling the development of biofilm density including active bacteria, inert biomass, and extracellular polymeric substances[END_REF][START_REF] Laspidou | Evaluating trends in biofilm density using the umcca model[END_REF][START_REF] Laspidou | Finite element modeling to expand the umcca model to describe biofilm mechanical behavior[END_REF][START_REF] Dr Noguera | Simulation of multispecies biofilm development in three dimensions[END_REF][START_REF] Picioreanu | Discrete-differential modelling of biofilm structure[END_REF][START_REF] Picioreanu | A theoretical study on the effect of surface roughness on mass transport and transformation in biofilms[END_REF][START_REF] Picioreanu | Effect of diffusive and convective substrate transport on biofilm structure formation: a twodimensional modeling study[END_REF]. This type of model presents the same features of CA models, thus it is characterized by the same drawbacks. However, the use of finite difference methods for solving the nutrient field can lead to a faster and more realistic model solution [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF].

A first tentative of combining continuous models with discrete ones to simulate complex biological structures, has been introduce by Ben-Jacob et al. [START_REF] Ben-Jacob | Generic modelling of cooperative growth patterns in bacterial colonies[END_REF], who developed a detailed model of bacterial colony growth using CA systems. The model includes the following generic features: diffusion of nutrients; movement of the bacteria; reproduction and sporulation; local communication. Nutrient diffusion has been modeled by solving a diffusion equation on a triangular lattice. Bacterial cells are divided in groups called "walkers" which can move on a triangular lattice within an envelope. Each walker is described by its location and an internal energy which affects its activity. The walker can lose or gain energy; when this energy drops to zero the walker becomes stationary while when this amount increases thanks to the consumption of nutrients and reaches a threshold, the walker duplicates. The model involves elements of cell-cell communication and chemotaxis and begins to reflect some of the complexity of a microbial community.

Hermanowicz [START_REF] Sw Hermanowicz | A model of two-dimensional biofilm morphology[END_REF][START_REF] Sw Hermanowicz | Two-dimensional simulations of biofilm development: effects of external environmental conditions[END_REF][START_REF] Sw Hermanowicz | A simple 2d biofilm model yields a variety of morphological features[END_REF]] developed a 2D model in which biofilm is represented by a 2D array of "cells". Each model cell can be "occupied", i.e. occupied by the biomass, or "empty", i.e. filled with water; mathematically it is represented by a dynamic variable that changes according to prescribed rules. The work is aimed at demonstrating how the CA approach is able to model the formation of self-organized structures based on simple development rules on a small scale and to evaluate the effect of the external environmental conditions. Model cells occupied by biomass can growth, divide or detach themselves according to a set of rules. Cell division depends on the probability of division, evaluated as a function of the environmental condition, such as nutrient concentration [START_REF] Sw Hermanowicz | A simple 2d biofilm model yields a variety of morphological features[END_REF]:

P = c c + K = (c/K) (c/K) + 1 , (2.26)
where: P is the probability of division; c is the local concentration of the limiting substrate

[ML -3 ]; K is the Monod half-saturation constant [ML -3 ].
A dividing grid cell spawns a daughter cell which will occupy one of the eight neighboring grid units with the following rules: if the grid cell is empty, the daughter grid cell will occupy it; if there is more than one empty cells, the choice is random; if all the neighboring cells are occupied the movement occurs in the direction of least resistance. This direction is evaluated calculating the shortest distance from each occupied grid cell to the biomass interface. In this case the daughter cell will push a whole line of cells in the direction of the nearest biofilm surface, to make place for itself. This mechanism of displacement resembles the concept of biomass advective flux formulated by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF], but in this case the biomass does not move with a uniform velocity, but jumps in a stochastic manner. Grid cells consume substrates that diffuse inside the boundary layer, whose thickness is a model input, and the biomass, while substrate concentrations remain constant outside the boundary layer. A matrix representing nutrient concentrations is superimposed on the working space containing water and biomass. The concentration field can be described by a Poisson equation. This equation is not solved numerically, but a modified version of the analytical solution obtained for 1D biofilm and zero-order nutrient uptake kinetics is introduced [START_REF] Sw Hermanowicz | A simple 2d biofilm model yields a variety of morphological features[END_REF]:

c = ⎛ ⎜ ⎝cS 1/2 - ⎛ ⎝ k 2D 1 8 8 i=1 1 d 2 i -1 ⎞ ⎠ 1/2 ⎞ ⎟ ⎠ 2 , (2.27)
where:

C S is nutrient concentration maintained constant outside the boundary layer

[ML -3 ];
k is the uptake rate for zero-order kinetics

[ML -3 T -1 ]; d i is a penetration distances [L].
Based on the idea that the resistance to mass transport in this case is a function of the penetration distance, equation (2.27) is derived modeling the overall resistance as a harmonic means of the resistances evaluated in the eight directions of nutrient supply considered for each point inside the biofilm. Detachment occurs randomly at a fluid/biomass interface with a probability increasing proportional with the biofilm thickness. More precisely, the probability of cell erosion is evaluated as a function of the hydrodynamic shear stress and biofilm cohesion. This approach allows the consideration of the detachment of larger clusters.

Picioreanu et al. [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF] introduced the first, so defined, hybrid-differential approach suitable for modeling immobilized cells, growing in a gel matrix. The model is still based on a CA approach for biomass spreading, but it evaluates the substrate field by solving a common reaction-diffusion equation. This model tries to overcome a common drawback of CA models deriving from the use of abstract parameters such as units of resource or random-walk distance, and relies on physical/chemical/biological parameters commonly used to describe biofilm systems (yields, concentrations, rates, fluxes of nutrients). Moreover, the combination of differential with discrete models allows the authors to predict the correct time evolution of biofilm growth, concentrations, fluxes and conversion rates, despite typical CA algorithms which work in a completely abstract time and space. In this pioneer work, the state of the system is represented by using two variables: the soluble limiting substrate concentration and the biomass density, coupled to a matrix which stores information about the grid element occupation. The model takes into account the three main processes characterizing biofilm development in hydrostatic conditions (i.e. diffusion-reaction-growth) and is aimed at demonstrating the validity of the new combined differential-discrete approach in studying biofilm development. Substrate transport occurs only by diffusion through a concentration boundary layer and further into the biofilm matrix and it is expressed in dimensionless form by the following equation [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF]:

∂S ∂t = D d 2 ∂ 2 S ∂X 2 + ∂ 2 S ∂Y 2 + ∂ 2 S ∂Z 2 -ρ S (C, S) , (2.28)
where: S = c S /c S0 is the dimensionless substrate concentration and c S0 is the substrate concentration in the bulk liquid; C = c X /c Xm is the dimensionless biomass concentration and c Xm is the maximum biomass density in a colony;

D is the diffusion coefficient [L 2 T -1 ]; d is the characteristic length (in this case the bead diameter) [L]; X = x/d, = y/d, Z = z/d are the space coordinates in dimensionless form; ρ S (C, S) is the normalized rate of substrate consumption [T -1 ].
The solution of diffusion-reaction processes is uncoupled from the calculations regarding the slower process of biomass spreading. In particular, equation (2.28) is solved by using relaxation algorithms and maintaining the matrices of biomass density and occupation state at a frozen level.

The biomass density is evaluated by solving the following equation [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF]:

∂C ∂t = ρ X (C, S) , (2.29)
where:

ρ X (C, S) is the normalized rate of biomass accumulation [T -1 ].
The occupation matrix is updated after solving the biomass balance. In particular, the biomass is redistributed when the maximum density is achieved in an elemental volume (x, y, z). The biomass is divided in two equal parts, redistributed in the neighboring space with no preferential direction according to simple CA rules (Figure 2.4). The pressure exerted by the biomass growing in the biofilm depth generates displacement of cells towards the biofilm-liquid interface. A single-cell release mechanism for detachment is considered only for the biomass located outside the carrier sphere.

The model introduced in [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF] has been applied later to simulate the biofilm growth on solid flat surfaces [START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF]. Despite continuum models, biofilm structure properties such as shape, porosity and density do not need to be initiated as input data but are generated by the model itself. Simulations at different substrate conversion/ transport rate ratios have been performed to evaluate their effect on biofilm structure. The biofilm surface shape has been characterized by using statistical quantities, such as biofilm surface enlargement, roughness, fractal dimension of biofilm surface while biofilm structure complexity has been evaluated as solids hold-up and biofilm compactness.

An extension of the previous model has been presented in [START_REF] Picioreanu | Discrete-differential modelling of biofilm structure[END_REF][START_REF] Picioreanu | A theoretical study on the effect of surface roughness on mass transport and transformation in biofilms[END_REF][START_REF] Picioreanu | Effect of diffusive and convective substrate transport on biofilm structure formation: a twodimensional modeling study[END_REF] and accounts for biomass growth and spreading, diffusive and convective transport and transformation of substrates, flow around the biofilm structure. The 2D model is fully quantitative, being based on first principles as Navier-Stokes equations, Figure 2.4: Schematic representation of the spreading rules adopted by a) Picioreanu et al. [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF][START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF]; b) Laspidou and Rittmann [START_REF] Laspidou | Modeling the development of biofilm density including active bacteria, inert biomass, and extracellular polymeric substances[END_REF][START_REF] Laspidou | Evaluating trends in biofilm density using the umcca model[END_REF]. Figure adapted from [START_REF] Cs Laspidou | Cellular-automata and individual-based approaches for the modeling of biofilm structures: Pros and cons[END_REF][START_REF] Tang | An improved cellular automaton method to model multispecies biofilms[END_REF].

substrates mass balances and kinetic laws for biomass growth. The biomass growth and spreading is modeled following the approach introduced in [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF]. The mass balance of the substrate is modeled by a convection-diffusion-equation and the flow field is governed by the incompressible Navier-Stokes equations in laminar regimes (2.10). The flow field and biofilm shape are interdependent since flow field shears the biofilm surface, erodes the protuberances and regulates the substrate concentrations at the biofilm-liquid interface. At the same time, changes in the biofilm shape determine a new boundary condition for the flow field and thus, different flow and substrate concentration. In this work a new strategy to follow biofilm development in time is presented. It is based on the idea that there is a clear separation of time scales in the biofilm growth. In this contest each process is solved assuming all the other processes occurring at different time scales at steady state.

Detachment has been incorporated later in the hybrid discrete-differential approach previously described [START_REF] Picioreanu | Two-dimensional model of biofilm detachment caused by internal stress from liquid flow[END_REF]. In this work, two known biofilm detachment mechanisms, i.e. erosion (loss of small biofilm parts, eventually only cells, mainly from the biofilm surface) and sloughing (loss of massive biofilm chunks, often broken from the substratum surface), have been modeled in a unitary way assuming that detachment is caused by stress developed in the biofilm structure [START_REF] Van Loosdrecht | Mathematical modelling of biofilm structures[END_REF]. The authors assumed that the biofilm detachment results from the combined effect of liquid shear and biofilm strength. The liquid flow above the biofilm exerts forces on the biofilm structure, both in normal and tangential to its surface, so the biofilm structure is subjected to a state of stress. The biofilm is assumed to be a homogeneous, isotropic, elastic material in the state of plane strain and the criterion of maximum distortion energy is applied to evaluate where the biofilm will break. In particular biofilm breakage is supposed to occur when the equivalent stress, expressed as a function of the normal and shear stresses, exceeds the cohesion strength.

Similarly to [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF], Laspidou and Rittman [START_REF] Laspidou | Modeling the development of biofilm density including active bacteria, inert biomass, and extracellular polymeric substances[END_REF][START_REF] Laspidou | Evaluating trends in biofilm density using the umcca model[END_REF] developed a multi-component cellular automaton model which combines the discrete representation of the solid phase by CA with classical continuous methods for soluble components. The model is based on the theory developed and quantified in [START_REF] Laspidou | A unified theory for extracellular polymeric substances, soluble microbial products, and active and inert biomass[END_REF][START_REF] Laspidou | Non-steady state modeling of extracellular polymeric substances, soluble microbial products, and active and inert biomass[END_REF] and considers three solid species including bacteria, EPS and inert residual biomass, two soluble microbial products, one limiting-growth substrate and an electron acceptor. They are all quantified in dimensionless form according to [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF]. The model reproduces in a 2D domain the growth of active biomass, the EPS and utilization-associated products formation, the EPS hydrolysis to biomass-associated products and their utilization as electron donors and the endogenous decay of active biomass to residual dead cells. The same solution strategy of Picioreanu et al. [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF] is used for solving the substrate field, but two new concepts are introduced for the cellular automaton algorithm: the composite density CompDen i,j and the biofilm consolidation which is aimed at describing the increases in biofilm density that occurs over time and deeper in the biofilm. The composite density varies with time and space and it is calculated for each CA cell according to equation (2.30):

CompDen i,j = X i,j a χ a,max + EP S i,j eps max + X i,j res χ res,max ,

where: X i,j a is the dimensionless density of the active biomass in i,jth cell; X r es i,j is the dimensionless density of the true residual inert biomass in i,jth cell; EP S i,j is the dimensionless concentration of EPS i,jth cell; χ a,max is the maximum active biomass packing density [ML-3]; χ res,max is the maximum of the true residual inert biomass packing density [ML-3]; eps max is the maximum EPS packing density [ML -3 ].

Each of the solid-phase components is computed from mass-balance equations and redistributed according to the CA algorithm except for the residual inert biomass, which is expected to accumulate only at the bottom of the biofilm.

The spreading of active biomass and EPS is simulated as the division of mother cells in daughter cells. In particular, the excess biomass is redistributed from one cell or compartment when the composite density exceeds a maximum value. The maximum composite density is specific for each compartment and increases over time with bioage (age of each biofilm department) in order to simulate the consolidation phenomenon.

A consolidation ratio is calculated for each compartment as an exponential function of biofilm age. It represents the degree of maximum packing density. The excess biomass is redistributed when the sum of the dimensionless density of the solid-phase components exceeds the consolidation ratio. The model distributes the excess biomass identifying the shortest path or path of least resistance (Figure 2.4). Moreover, a first-order detachment is included for the outmost layer of the biofilm. The outputs of the model have been used later to perform the biofilm stress analysis aimed at evaluating the biofilm's strength and resistance to detachment [START_REF] Laspidou | Finite element modeling to expand the umcca model to describe biofilm mechanical behavior[END_REF].

Individual Based Models

The term IbMs is addressed to a class of multidimensional models whose objective is to describe the actions and properties of the individuals constituting the bacterial population or community [START_REF] Ferrer | Individual-based modelling: an essential tool for microbiology[END_REF][START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF][START_REF] Picioreanu | Particle-based multidimensional multispecies biofilm model[END_REF][START_REF] Xavier | A framework for multidimensional modelling of activity and structure of multispecies biofilms[END_REF]. IbMs use a bottom-up approach and can also be classified as spatially structured population methods. This type of model was born with the aim of surmounting all the drawbacks deriving from the application of discrete rules, typical of the CA approach, to biofilm spreading [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF]. IbMs allow cells movement only on a continuous set of directions and distances while CA models typically require that biomass moves only in the finite number of lattice directions. Usually transport and reaction of a solute species, local microbial growth rates are modeled using differential approaches [START_REF] Xavier | A modelling study of the activity and structure of biofilms in biological reactors[END_REF].

In all IbMs, bacterial cells represent the fundamental entities and they are modeled as hard spheres in continuous 3D space, each of them having a variable volume, mass and a set of mutable growth parameters. These spherical agents act independently, analogously to how individual bacterial cells behave within biofilms. IbM models do not specify any global (population level) laws such as exponential population growth. The behavior of the agents is defined explicitly with a set of rules that mimic the behavior of individual bacterial cells, i.e. growth through the consumption of substrates, reproduction through cell division, production of metabolites etc.

The first attempt to model bacterial colony growth by using this approach has been proposed in [START_REF] Ju Kreft | Bacsim, a simulator for individualbased modelling of bacterial colony growth[END_REF][START_REF] Ju Kreft | Applications of individual-based modelling in microbial ecology[END_REF] and it has been used later to simulate a multispecies biofilm [START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF]. The use of IbM for biofilms can be classified as a more realistic approach that quantitatively incorporates the physiology of individual cells [START_REF] Hellweger | A bunch of tiny individuals-individual-based modeling for microbes[END_REF]. In [START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF], the authors introduced a fully quantitative IbM based on BacSim, which consists of two main parts: one deals with the simulation of the growth and behavior of individual bacteria as autonomous agents; the other one deals with the simulation of substrate and product reaction and diffusion. According to [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF], the bacterial growth has been simulated assuming the diffusion process at a pseudo-steady state (between each biomass spreading iteration) since biofilm growth is usually a much slower process than diffusion of substrates into the biofilm. Each cell performs "actions" as a result of the environmental conditions and its internal state. It grows by consuming the substrates and divides when a certain volume is reached; it moves as a consequence of being pushed by its neighbors [START_REF] Xavier | Assessment of three-dimensional biofilm models through direct comparison with confocal microscopy imaging[END_REF]. The model considers the random variation of cell parameters, the maximal uptake rate and the volume-at-division, by using a Gaussian distribution with a coefficient variation (CV) of 10%. The pressure buildup due to the growth of biomass is released by maintenance of a minimum distance between the neighboring cells. For each cell, the vector sum of all positive overlap radii with the neighboring cells is calculated and then the position of the cell is shifted in the direction opposite to this vector. Therefore, the biomass packing in the biofilm is defined by the shoving parameter K shov which represents the spacing among cells. When a cell reaches a critical volume, it divides resulting in the creation of another cell, the "daughter", and the mass of the original cell is distributed slightly unevenly between these two spheres. The random choice of the direction for the placement of daughter cells and the uneven division of mass between cells makes the model stochastic.

The substrate concentration is governed by a reaction-diffusion equation which is solved by using a relaxation method. The uptake rates calculated for each bacterium occupying a certain grid element of the substrate field are averaged on an area percentage basis. Despite the unilateral shoving mechanism adopted in [START_REF] Ju Kreft | Bacsim, a simulator for individualbased modelling of bacterial colony growth[END_REF], Kreft et al. [START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF] introduced a mutual mechanism which minimizes the effect of a bias making the shoving independent of the sequence with which the bacteria are accessed. This result is achieved by inverting the sequence of who shoves whom every 10 steps. The performed simulations revealed that the mutual scheme is better in avoiding overlaps and relocating cells and helps in reducing sequential execution biases. The BacSim framework has been compared with the population-level model introduced by Picioreanu et al. [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF][START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF]. The simulations show similar results in principle, due to modeling the same physical processes, but differ in details of biofilm shape and growth of minority species. In particular, the IbM method results better suited to the description of multispecies systems since CA models determine the production of too much internal mixing of species within a colony or the generation of anisotropic colonies.

The IbM approach introduced by Kreft et al. [START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF] has been used later to model the mechanism of production and spreading of extracellular polymeric substances [START_REF] Kreft | Effect of eps on biofilm structure and function as revealed by an individual-based model of biofilm growth[END_REF][START_REF] Kreft | Modeling biofilms with extra-cellular polymeric substances[END_REF] and to test several evolutionary and ecological hypotheses [START_REF] Ju Kreft | Biofilms promote altruism[END_REF]. In [START_REF] Kreft | Effect of eps on biofilm structure and function as revealed by an individual-based model of biofilm growth[END_REF], the EPS formation is stoichiometrically coupled to growth; the EPS produced is first bound to the bacterial agent forming a protective layer and then excreted as a separate agent that will participate in the shoving mechanism along with the bacterial agents. In [START_REF] Ju Kreft | Biofilms promote altruism[END_REF], the IbM approach has been applied to study the development of altruistic behavior by bacteria in biofilms.

The very detailed level of biofilm description characterizing IbMs can represent a disadvantage in modeling systems with large-scale heterogeneity. In order to extend the spatial scale of the previous IbMs, Picioreanu et al. [START_REF] Picioreanu | Particle-based multidimensional multispecies biofilm model[END_REF] introduced a multidimensional particle-based modeling approach which considers the presence of larger biomass particles, but keeps the rules for biomass redistribution and shoving introduced in [START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF]. In this model, the biofilm biomass is divided in spherical particles containing only one type of active biomass and a fraction of inert biomass resulting from the decay of all active biomass types. The size of these spherical agents is chosen to represent a cell cluster of similar cells and no variability in metabolic parameters is included for all biomass particles of the same type. Biomass division and spreading is based on the same mechanism introduced in [START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF]. A simplified biomass detachment model is introduced: it consists of removing every particle, which is shifted above an imposed biofilm thickness limit, due to a shoving step. The substrate field is governed by a dynamic-state diffusion reaction equation, which is uncoupled from the solution of biomass evolutions, as stated in [START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF]. A steady-state solution of the partial differential equation for mass balances of soluble substrates in the biofilm is found by a nonlinear multigrid algorithm. The numerical simulations reveal that the IbM framework for cell transport describes its continuum counterpart at least in a 1D case.

The modeling approaches introduced in [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF][START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF][START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF][START_REF] Picioreanu | Particle-based multidimensional multispecies biofilm model[END_REF], have been integrated to provide a framework that defines the structure for multidimensional, multispecies dynamic modeling of biofilm systems [START_REF] Xavier | A framework for multidimensional modelling of activity and structure of multispecies biofilms[END_REF]. This IbM takes into account the concept of structured biomass which is constituted by multiple bacterial species, inert biomass and EPS. Three spatial scales are considered: a) individual scale, which deals with the behavior of biomass agents; b) biofilm scale, which works on a community level; c) system scale, which takes into account the interactions between the bulk liquid compartment and the biofilm. According to Picioreanu et al. [START_REF] Picioreanu | Particle-based multidimensional multispecies biofilm model[END_REF], biomass particles can represent either a single cell or a cluster of cells of the same species. Each biomass particle is correlated to a pDocument which defines the number and type of particulate species constituting the agent. The mass of each particulate species varies according to a bioconversion equation and determines changes in the agent volume. The agent duplicates when the maximum particle radius is reached. The masses of all particulate species contained in the dividing agent are then redistributed between the two resulting agents. The EPS production and excretion are modeled according to [START_REF] Kreft | Effect of eps on biofilm structure and function as revealed by an individual-based model of biofilm growth[END_REF] (Figure 2.5). In the case of EPS decay or inactivation of the bacterial biomass, the framework includes net reduction of the biofilm volume. A multidimensional extension of the method used in [START_REF] Wanner | A multispecies biofilm model[END_REF] is introduced to study detachment and other biomass losses. In particular, a continuous detachment speed function is used to model both erosion and sloughing, as extensively described in [START_REF] Xavier | A general description of detachment for multidimensional modelling of biofilms[END_REF]. The computation of solute concentration fields is decoupled from biomass dynamics, as adopted in [START_REF] Picioreanu | A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads[END_REF][START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF][START_REF] Ju Kreft | Individual-based modelling of biofilms[END_REF][START_REF] Picioreanu | Particle-based multidimensional multispecies biofilm model[END_REF]. The solute concentration field is computed by a multigrid solver, as previously applied in [START_REF] Picioreanu | Particle-based multidimensional multispecies biofilm model[END_REF]. Bulk concentration of solute species can be: constant, in which case the bulk liquid is assumed to be an infinite and constant supply of the solute; intermittent, by alternating feast and famine cycles, or can be computed from a mass balance equation applied to the whole system. An alternative method to treat EPS has been introduced by Alpkvist et al. [START_REF] Alpkvist | Threedimensional biofilm model with individual cells and continuum eps matrix[END_REF]. In particular, the authors used a continuous representation of EPS combined with an IbM of individual bacteria. According to [START_REF] Klapper | Finger formation in biofilm layers[END_REF][START_REF] Alpkvist | Threedimensional biofilm model with individual cells and continuum eps matrix[END_REF], EPS is modeled as a viscous fluid, which is well justified by both experimental facts and physical grounds. On the other hand, the IbM approach is used to model the behavior of each bacterial cell, the local interactions between different microbial species and individual variation of microbial cells. The movement of EPS and cells on a global level in the biofilm is governed by an advection speed which is assumed to follow the Darcy's law. At the same time, the individual cells (biomass spheres) undergo a shoving mechanism when they get too close to each other; the cell shoving introduces small local deviations in the flow field. The model has been applied to study the consolidation process in mature biofilms. This process seems to derive from the presence of a negative pressure in the lower region of the biofilm which is generated by EPS and cell degradation processes and results in cell transport towards the substratum.

A new modeling platform dedicated to IbM of microbial communities has been recently introduced by Lardon et al. [START_REF] La Lardon | idynomics: next-generation individual-based modelling of biofilms[END_REF]. In this work, the authors have tried to combine most features of previous models and incorporate various improvements in order to provide a common basis for further developments. In order to address this aim, an open-source software called individual-based Dynamics of Microbial Communities Simulator (iDynoMiCS) has been developed. The iDynoMiCS structure emerges from the combination of the previous modeling approaches but presents inherent differences. Primarily iDynoMiCS allows for the introduction of non-bacterial agents (archea, protozoa, algae or fungi). Microbial agents are structured in compartments including all intracellular components (active or inactive biomass, storage compounds, etc.), bounded by an outer layer of capsular EPS. All the agents are updated in a random order, which is changed for each time step in order to remove any bias. An individual agent can carry out a different suite of reactions than other individuals of the same species. The EPS excretion is represented by a particulate method: the EPS produced is continuously released into the environment and distributed to the EPS particles of the same type, present in the neighborhood of the EPS-producing agent; if no such agents are found, a new EPS particle is created. To recreate a continuum representation of EPS, smaller radii of EPS particles are adopted (Figure 2.5). According to [START_REF] Klapper | Finger formation in biofilm layers[END_REF][START_REF] Alpkvist | Threedimensional biofilm model with individual cells and continuum eps matrix[END_REF], a pressure field to model biomass spreading or consolidation, is introduced. As all IbMs, iDynoMiCS is affected by stochasticity in the choice of the initial agent locations and masses, the cell division threshold volume, the cell death threshold volume, the daughter cell orientation and size, the excretion direction of new EPS particles and the order of agents update. During the last decade, the IbM approach has been widely used to predict several structural features of microbial biofilms and the results match experimental observations. The IbM approach has been used to evaluate the biomass production/consumption and transport of biofilm for microbial fuel cells [START_REF] Picioreanu | A computational model for biofilm-based microbial fuel cells[END_REF]. The effect of microbial motility on biofilm morphology has been analyzed in [START_REF] Picioreanu | Microbial motility involvement in biofilm structure formation-a 3d modelling study[END_REF] and the concepts of IbM have been applied to describe and optimize a biofilm and granular reactor [START_REF] Xavier | Multiscale individual-based model of microbial and bioconversion dynamics in aerobic granular sludge[END_REF].

Discussion

Due to biofilm involvement in a large range of human activities and natural processes, developing an effective mathematical modelling approach may be essential for elucidating the processes involved in biofilm formation and maturation, as well as for developing a strategy to minimize the biofilm related risks and exploit their technical possibilities. The heterogeneity of biofilm structure and the interdependence of physical/chemical/biological processes, occurring at different time and space scales, make mathematical modeling of biofilm growth and structure a special challenge for researchers.

A detailed overview of the wide range of modeling approaches developed during the last decades has been presented above. Selecting an appropriate model may represent a challenging issue for both researchers and practitioners [START_REF] Morgenroth | Biofilm models for the practitioner[END_REF]. The scope and output of the model constitute a discriminator factor: practitioners are interested in developing models able to predict quantitatively the performance and responses of biofilm reactors; researchers consider modeling as a powerful tool to understand the fundamental mechanisms regulating the formation and performance of biofilms. Therefore, practitioners aim at simpler biofilm models, which can be easily calibrated by using the data provided by experimental activity, while in research the degree of complexity is increasing over time.

On the basis of the model classification proposed in this work, general guidelines for the selection of the most suited modeling tool, based on the specific needs of the model user, are provided. In particular, two general questions are answered in the following sections, which should cover the doubts arising in choosing a modeling approach: i) when to use 1D, 2D or 3D models? ii) Should Continuum or Discrete approach be used? However, there is no general answer to these questions because the use of a specific modeling tool is objective depending. For instance, different environmental conditions produce a variety of structurally very different biofilms; some mathematical models are able to capture this heterogeneity; others, based on simplifying assumptions, are not. The extent to which simplifications and idealizations must or can be introduced depends on the particular purpose of the mathematical model [START_REF] Hermann | What do biofilm models, mechanical ducks, and artificial life have in common?[END_REF].

When to use 1D, 2D or 3D models?

As widely described above, 1D models consider only the direction perpendicular to the substratum: this represents a valid simplification when vertical gradients of variables and parameters are orders of magnitude higher than those in the directions parallel to the carrier surface [START_REF] Xavier | A general description of detachment for multidimensional modelling of biofilms[END_REF]. This hypothesis verifies in the case of uniform bulk liquid conditions over the whole substratum area, when the substratum area is regular and large enough compared to biofilm thickness or for smooth biofilm surfaces. Since this applies to many (not all) of engineering biofilm systems, 1D models have been widely used to predict the whole process dynamics of biofilm reactors and are increasingly used as educational material in engineering curricula [START_REF] Boltz | Mathematical modelling of biofilms and biofilm reactors for engineering design[END_REF]. In addition, the choice of 1D models may reflect the need of keeping the computational effort at a low level. An inherent limitation of 1D models rely on the simplified modeling of bulk liquid as a completely mixed compartment. The calculations regarding the flow field are neglected as well as the interactions between liquid flow and biofilm surface. However, recent improvements in microscopy and imaging techniques have revealed that numerous biofilms are not uniform and spatial irregularities in real biofilms cannot be interpreted using a conceptual model of biofilms in which microorganisms are uniformly distributed in a continuous matrix of extracellular polymers [START_REF] Lewandowski | Modeling mass transport and microbial activity in stratified biofilms[END_REF]. Biofilms have been recognized as complex 3D heterogeneous entities characterized by a highly porous structure filled with fluid which supplies nutrients to microorganisms and erodes biofilm surface resulting in removal of biomass. This spatially heterogeneous architecture can induce complex flow patterns and affect mass transfer [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF]. 1D models result inappropriate to describe the dynamics of biofilm activity when structural dependent factors such as external mass transfer coefficient or porosity vary significantly with time. Therefore, 2D and 3D models have been developed to capture this heterogeneity. Multidimensional models are able to evaluate the substrate removal and biomass production rates of dynamic biofilm systems, but they can be used also to evaluate the interaction among biofilm shape, fluid flow, biomass decay and detachment, by taking into account the fluid dynamics modeling of the liquid phase and the effect of biofilm geometry on the external mass transfer rates. They are suitable to study the effect of different environmental conditions on biofilm structure and to evaluate the multidimensional interactions between microbial communities, such as microbial segregation. 2D models are not representative of a 3D domain where flow can by-pass dense biofilm structures [START_REF] Wanner | Mathematical Modeling of Biofilms[END_REF]. However, simulations performed by Eberl et al. [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF] highlighted that in most cases, 2D and 3D models lead to equivalent external mass transfer coefficient. On the other hand, 3D models can be useful for biofilms consisting in isolated colonies where advective transport becomes not negligible. The use of highly accurate 2D and 3D models requires a detailed description of the biofilm structure at a meso-scale, which can be accomplished only by using the modern investigation techniques, and the solution of a nonlinear system of partial differential equations in a complex domain. In particular, multidimensional biofilm models have been singularly used as research tools, where an accurate resolution of the processes occurring inside the biofilm is required, and their application as engineering tools is limited by the high spatial resolution and the level of detail required for model calibration.

Should continuum or discrete approach be used?

Continuum models represent a valid alternative to the discrete approach since all the drawbacks that characterize discrete models seem to arise from the discreteness of the spreading mechanism adopted [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF]. As above described, a continuum biofilm model: i) is characterized by a continuous representation of biomass; ii) is based on differential equations widely used in physics to model the dynamics of biomass spreading and iii) generates deterministic solutions. The main advantages of continuum biofilm descriptions derive from the use of the powerful framework of partial differential equations. Indeed, the use of differential calculus allows achieving quantitative results for substrate transport that can be compared with data measured in real systems.

The 1D continuum models persist today as widely used methods to describe macroscopic conversions and for interpreting and predicting biofilm reactors performance [START_REF] Xavier | Assessment of three-dimensional biofilm models through direct comparison with confocal microscopy imaging[END_REF]. The 1D dynamic multispecies model of Wanner and Reichert [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF] implemented in the software AQUASIM, is up to date the most widely used biofilm model applied to engineering design, as it is sufficiently accurate for predicting global mass conversion rates for a full bioreactor. However, 1D models do not give any knowledge on the local spatial architecture of the biofilm and multidimensional continuum models have been developed with the aim of covering this gap.

The main challenges in developing multidimensional continuum models rely on the presence of moving boundaries, i.e. the biofilm-fluid interface, fluid flow, non-linear growth kinetics and discontinuous gradients across the boundary biofilm-fluid interface [START_REF] Duddu | A combined extended finite element and level set method for biofilm growth[END_REF]. In particular, the use of multidimensional continuum models implies high computational efforts and sometimes requires simplifying assumptions to solve the differential equations describing biofilm evolution on irregular domains. Flow field calculations are usually much more computationally expensive than simulations of biofilm growth; therefore bulk flow hydrodynamics has been usually neglected in many multidimensional continuum models. Moreover, a variety of resolution methods, all of them characterized by significant computational efforts, for solving elliptic equations on irregularly shaped domains in conjunction with moving interfaces have been investigated. The formulation and derivation of continuum models require a comprehensive mathematical skill, a higher computational effort compared to discrete methods and the computational algorithms could be sometimes not trivial [START_REF] Alpkvist | Threedimensional biofilm model with individual cells and continuum eps matrix[END_REF].

Despite the high computational efforts, multidimensional continuum models are more convenient than discrete approaches when applied to mechanical problems. Indeed, in this case discrete models can hardly describe the global relationships between the individuals constituting the biofilm.

On the other hand, discrete models are able to represent the typical multidimensional structural heterogeneity of biofilm in good agreement with experimental expectations, but they generate computational results that include elements of randomness. Their output depends on the sequence of execution of methods on the discrete objects and introduces stochastic effects into the solutions. In the case of IbMs, the stochasticity manifests in two occasions: (a) for the random choice of direction for the placement of "daughter" particles and (b) for the uneven division of mass between in-cell division [START_REF] Xavier | Assessment of three-dimensional biofilm models through direct comparison with confocal microscopy imaging[END_REF]. For CA, the rules used for biomass spreading are sometimes formulated arbitrary and might lead to aesthetically driven, rather than to physically motivated, model formulation [START_REF] Hj Eberl | A new deterministic spatiotemporal continuum model for biofilm development[END_REF]. Generally they are lattice dependent and not invariant to changes of coordinate system. Moreover the same initial conditions can lead to different model outputs and error analyses are non-trivial [START_REF] Alpkvist | Threedimensional biofilm model with individual cells and continuum eps matrix[END_REF]. Therefore, for a discrete model, several runs with the same initial state are needed to average the stochastic effect before conclusions are drawn.

Despite the aforementioned disadvantages, both discrete approaches, CA and IbM represent powerful modeling tools which have been applied not only in ecology, but in many other disciplines such as social, economical, demographical and political sciences. Cellular automaton models avoid the mathematical and computational burden of continuous models, and are more straightforward [START_REF] Tang | An improved cellular automaton method to model multispecies biofilms[END_REF]. Similarly to multidimensional continuum models, they work on a larger scale than IbMs, which are usually used for studies at the scale of micrometers to centimeters and therefore are computationally more intensive. When CA are applied to the case of single species biofilms, they produce similar results of continuum models; while in the case of multispecies biofilms, these algorithms generate a lot of internal mixing within the colony. However, over the years, several types of spreading rules have been adopted in order to minimize the excess of mixing [START_REF] Dr Noguera | Biofilm modeling: present status and future directions[END_REF][START_REF] Tang | An improved cellular automaton method to model multispecies biofilms[END_REF][START_REF] Laspidou | Modeling the development of biofilm density including active bacteria, inert biomass, and extracellular polymeric substances[END_REF]. According to [START_REF] Cs Laspidou | Cellular-automata and individual-based approaches for the modeling of biofilm structures: Pros and cons[END_REF], CA are especially suitable to model old and aged biofilms, characterized by the presence of cavities and experiencing the phenomenon of consolidation.

IbMs started to be developed with the aim of simulating individual or localized behaviors which are usually not adequately described by population-average models. These models are mainly addressed to capture the micro-scale level and to describe how individual processes, interactions and local variability affect the macroscopic structure of biofilms. One of the main drawbacks deriving from the application of this approach, rely on the assumption of individual microorganisms as hard spheres and on the use of a predetermined shoving parameter to model the direction of the biomass movement and porosity. Furthermore information on the individual heterogeneity of growth parameters, the volume fraction occupied by cells in colonies and the biomass spreading mechanism adopted by different microorganisms are sometimes missed. In addition, despite their computational demand, the IbM approach can incorporate rare species or rare events, it can make a distinction between spreading mechanisms adopted by different bacteria and operates at the highest spatial resolution level relevant in a biofilm [START_REF] Picioreanu | Use of mathematical modelling to study biofilm development and morphology[END_REF]. IbMs represent a big promise in modelling multispecies biofilms and in incorporating concepts such as cell-to-cell signaling, quorum sensing and cellular motility.

Conclusions

In this work, a comparison of the two types of approach, namely the continuum and discrete, used to simulate the development of biofilm structure has been conducted in order to elucidate the main advantages/disadvantages deriving from the application of each approach. Continuum models benefit from the framework of differential calculus and represent a valuable tool to understanding biofilm processes in a quantitative and deterministic way. However, 1D continuum models assume a planar geometry and therefore, cannot take into account the biofilm spatial heterogeneity. On the other hand, the formulation of multidimensional continuum models requires a comprehensive mathematical skill and sometimes, very costly or demanding numerical solutions. Discrete models are more recent in time in biofilm research and are based on the idea that biofilms can be characterized as stochastic living systems. These models have shown their capability of representing biofilm structure heterogeneity in good agreement with experimental results. However, they introduce elements of randomness, mostly in modeling the spreading of biomass and leading to shapes resembling biofilm structure but that may not simulate reality exactly.

All biofilm models constitute valuable tools in predicting biofilm growth and structure and their choice depends mostly on the situation that needs to be modeled. To date, extensive experimental activities are being carried out to understand how biofilm grows and interacts with the environment. Evaluating the complexity of intracellular interactions and communications represents one of the future challenges in biofilm research. Only the collaboration among researchers with different expertise will lead to the definition and development of new modeling approaches, able to take into account the advances in biofilm ecology.

Chapter 3

Analysis and simulations of the initial phase in multispecies biofilm formation

The work presents a mathematical modelling approach to study dynamic competition during the attachment phenomena in the initial phase of biofilm growth. Biofilm development is described by a set of nonlinear hyperbolic partial differential equations. Diffusion of substrates through biofilm is modeled by a set of semilinear parabolic partial differential equations. The two sets of equations are mutually connected. The resulting mathematical problem is a free boundary value problem, which is essentially hyperpolic. A characteristic-like method is introduced to convert differential equations into integral equations. Fixed-point theorem is used to obtain existence, uniqueness and properties of solutions. The model has been applied to the competition of heterotrophic-autotrophic bacteria in a multispecies biofilm. The effects of different attachment rates on the biofilm dynamics including biofilm thickness, volume fractions of bacterial species and substrate concentration trends have been investigated. The simulations show that the different attachment rates influence biofilm thickness, of course. However, the volume fractions of bacterial species mainly depend on biofilm internal dynamics and substrate concentration trends. The bulk concentrations of microbial species play a relative important role only in the outermost layers of biofilm.

Introduction

Biofilms are commonly defined layer like aggregations of microorganisms and are involved in a variety of scenarios including pollution, corrosion, biofouling, biomedical applications, bacterial growth in water distribution systems, attached growth systems for wastewater treatment. Biofilm metabolism is characterized by some inherent features that provide several advantages and some challenges for applications. Indeed, a single biofilm can exibit varying environmental and kinetic characteristics, that is, it can include a variety of microbial groups contributing to the conversion of different organic and inorganic substrates; biofilm cells are at least 500 times more resistant to antimicrobial agents [START_REF] Jw Costerton | Microbial biofilms[END_REF] and benefit from interspecies cooperation.

Biofilm development is determined by "positive" processes, like cell attachment, cell division, and polymer production, which leads to biofilm volume expansion, and "negative" processes, like cell detachment and cell death, which contribute to biofilm shrinking. The main biofilm expansion is due to bacterial growth and extracellular polymer production. The nutrients necessary for biofilm growth are dissolved in the bulk liquid and are transported by molecular diffusion first through the boundary layer, where the external mass transfer resistance is concentrated, and then through the biofilm matrix. The external fluid flow regulates biofilm growth by establishing the concentration of substrates and products at the liquid-solid interface. At the same time the fluid flow shears the biofilm surface, eroding the protuberances. So biofilm structure results from the interplay of different interactions, such as mass transfer, conversion rates and detachment forces. An accurate modeling of such a system have to take into account all of these factors, since these factors strongly affect the overall performance of biofilm-based systems [START_REF] Picioreanu | Modelling and predicting biofilm structure[END_REF].

The transport of microorganims to and from the biofilm (attachment and detachment) is particularly important since it defines the microbial ecology of the biofilm and plays a crucial role in the start up of biofilm reactors. In particular, biofilm formation occurs when specified bacterial species, able to make the first colonization of substratum surface, attach to the surface and start producing an extrapolymeric matrix that will allow the attachment of other microbial species [START_REF] Jw Costerton | Overview of microbial biofilms[END_REF]. This bacterial adhesion can be performed only by few microbial species and it is especially rapid and specific if the surface in question is itself a nutrient. After the first colonization, a thin layer of biofilm constitutes and the flux of micro organisms from bulk liquid to biofilm continues playing a crucial role in both biofilm development and in competition of microbial species for substrates and place.

The attachment process is influenced by both physical and biological factors, such as bulk liquid characteristics in terms of kind and concentration of microbial species, flow velocity and turbulence, geometry of the substratum. Therefore, the attachment flux can assume different values determining different biofilm development pathways. At the same time, the growth of the existing thin layer of biofilm is influenced by substrate concentrations and coupled with attachment flux determines the overall biomass volume increase.

Mathematical modelling of the attachment process represents a useful tool to evaluate this phenomenon and to predict the structural biofilm development. Many of biofilm models developed during the last decades are able to reproduce the complex interactions existing between the main processes, including cell attachment, involved in the formation of biofilm structure. Biofilm models have undergone a temporary evolution and mostly differ on the way biomass spreading is treated. Continuum models [START_REF] Klapper | Mathematical description of microbial biofilms[END_REF][START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Wanner | Competition in biofilms[END_REF] consider biomass as an unicuum and are based on conservation principles. These studies are mostly centered on the biofilm growth dynamics including the biofilm thickness and spatial distribution of microbial species and substrate concentration. These continuum models can be related to the underlying description offered by models at the microscopic scale as documented in [START_REF] Bellomo | Multiscale biological tissue models and flux-limited chemotaxis from binary mixtures of multicellular growing systems[END_REF]. Later, discrete models have been developed to reproduce biofilm spatial heterogeneities by using simple rules [START_REF] Picioreanu | Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach[END_REF][START_REF] Wimpenny | A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models[END_REF]. These models can capture the various biofilm growth patterns observed in experiments and strongly suggest that the biofilm structure is largely determined by the surrounding substrate concentration. In this work a mathematical model based on a continuum approach and able to describe the attachment process during biofilm growth is presented. In particular, the objectives of this study include:

-to propose a mathematical modelling approach to study population dynamics during the attachment phenomena in the initial phase of biofilm growth;

-to provide a qualitative analysis to the solutions of the corresponding free boundary value problem;

-to develop numerical simulations to illustrate the model.

Initial phase of biofilm formation

Consider a multispecies biofilm formed by n microbial species. 1D model of biofilm growth, based on the continuum description, is governed by the following equations [START_REF] Klapper | Mathematical description of microbial biofilms[END_REF][START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Wanner | Competition in biofilms[END_REF],

∂X i ∂t + ∂ ∂z (uX i ) = ρ i r M,i (z, t, X, S), i = 1, ..., n, 0 ≤ z ≤ L(t), t > 0, (3.1) ∂u ∂z = n i=1 r M,i = G(z, t, X, S), 0 < z ≤ L(t), t > 0, (3.2)
L(t) = u(L(t), t) + σ(t), t > 0, (3.3)
∂S j ∂t - ∂ ∂z (D j ∂S j ∂z ) = r S,j (z, t, X, S), 0 < z < L(t), t > 0, j = 1, ..., m, (3.4)
where: z is the biofilm growth direction, assumed perpendicular to substratum; f i (z, t) is the volume fraction of microbial species i, n i=1 f i = 1; ρ i denotes constant density; X i (z, t) = ρ i f i denotes the concentration of microorganism i, X = (X 1 , ..., X n ); u(z, t) is the velocity of the microbial mass; S j (z, t) denotes the concentration of substrate j, j = 1, ..., m, S = (S 1 , ..., S m ); r M,i (z, t,X,S) is the specific growth rate; L(t) denotes biofilm thickness; σ(t) is the biomass flux from bulk liquid to biofilm; D j denotes the diffusivity coefficient of substrate j; r S,j (z, t, X, S) is the conversion rate of substrate j.

Free boundary value problem

As outlined in Section 3.1, we want to discuss the free boundary value problem for the initial phase of biofilm development. We consider the situation where a thin layer of biofilm has been already formed. This leads to assume a strictly positive initial thickness for the biofilm: L(0) > 0. The case where L(0) = 0 was discussed in [START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF]. In addition, it is assumed that there is no biomass flux at the support and this implies u(0, t) = 0. Therefore, the following initial conditions will be associated to equations (3.1)-(3.3)

X i (z, 0) = ϕ i (z), i = 1, ..., n, 0 ≤ z ≤ L(t), (3.5)
u(0, t) = 0, t > 0,

L(0) = L 0 > 0, (3.6)
where L 0 denotes the initial thickness of biofilm and ϕ i (z), i = 1, ..., n, the initial concentrations of microbial species.

The initial phase of biofilm growth is strongly influenced by the attachment. This process can be defined as the immobilization of cells suspended in the bulk liquid to biofilm or substratum. Mathematically, the attachment is modelled as a flux, usually denoted by σ. Moreover, the biofilm development in this phase also depends on substrate availability, since a thin layer of biofilm has been already constituted. For mature biofilms, the most significative biomass flux occurs from biofilm to bulk liquid. This biological process, known as detachment, mostly depends on hydrodynamic conditions and biofilm thickness. However, in the initial phase, attachment is the prevailing process. So, σ is assumed to be a strictly positive function of time in this work. The attachment affects microbial species distribution into biofilms and can even introduce new species existing in the bulk liquid in the already developed biofilm structure. So, the following boundary condition is needed [START_REF] Wanner | A multispecies biofilm model[END_REF],

X i (L(t), t) = ψ i (t), i = 1, ..., n, t > 0, (3.8)
where ψ i (t) denotes the concentration of the microbial species i in the bulk liquid. Finally, we introduce the initial-boundary conditions for equations (3.4). The initial conditions are quite general

S j (z, 0) = S j0 (z), 0 ≤ z ≤ L 0 , j = 1, ..., m. (3.9)
Moreover, no flux boundary conditions are assumed at the biofilm support and prescribed boundary condition on the free boundary

∂S j ∂z (0, t) = 0, S j (L(t), t) = S jL (t), t > 0, j = 1, ..., m, (3.10)
S j0 (L 0 ) = S jL (0), S j0 (0) = 0. (3.11)
The value assumed by substrate concentrations at biofilm/bulk-liquid interface take into account the effect of mass-transport resistance. Equations (3.10) 2 determine the substrate concentration trends into biofilm and reflect the bulk liquid substrate availability. Microbial growth rate depends on substrate concentration so the boundary conditions strongly affect bacterial species distribution, determining the predominance of some bacterial species over others.

The free boundary problem considered in this work is summarized by equations (3.1)-(3.4) with initial-boundary conditions (3.5)-(3.11).

Characteristic-like method

When 0 ≤ z 0 ≤ L 0 , we consider the characteristic-like lines for system (3.1), fig.

3.1, z = z(z 0 , t), 0 ≤ z 0 ≤ L 0 , t > 0, (3.12)
defined by

∂z ∂t (z 0 , t) = u(z(z 0 , t), t), z(z 0 , 0) = z 0 , 0 ≤ z 0 ≤ L 0 , t > 0. (3.13)
When z 0 = L(t 0 )(> L 0 ), 0 < t 0 ≤ t, we consider the characteristic-like lines, fig.

3.1,

z = z(L(t 0 ), t) = c(t 0 , t), 0 < t 0 ≤ t, t > 0, (3.14)
defined by ∂c ∂t (t 0 , t) = u(c(t 0 , t), t), c(t 0 , t 0 ) = L(t 0 ), 0 < t 0 ≤ t, t > 0.

(3.15)

d dt X i (z(z 0 , t), t) = F i (z(z 0 , t), t, X, S), i = 1, ..., n, 0 ≤ z 0 ≤ L 0 , t > 0, (3.16
) with initial conditions

X i (z(z 0 , 0), 0) = X i (z 0 , 0) = ϕ i (z 0 ), 0 ≤ z 0 ≤ L 0 , i = 1, ..., n, (3.17)
where

F i = ρ i r M,i -X i n h=1
r M,h = F i (z, t, X, S), i = 1, ..., n.

(3.18)

Using (3.15) in (3.1) gives the following system of nonlinear differential equations

d dt X i (c(t 0 , t), t) = F i (c(t 0 , t), t, X, S), i = 1, ..., n, 0 < t 0 ≤ t, t > 0, (3.19)
with initial conditions

X i (c(t 0 , t 0 ), 0) = X i (L(t 0 ), t 0 ) = ψ i (t 0 ), 0 < t 0 ≤ t, t > 0, i = 1, ..., n. (3.20) Differential system (3.16)-(3.17) is equivalent to the integral system ⎧ ⎨ ⎩ X i (z(z 0 , t), t) = ϕ i (z 0 ) + t 0 F i (z(z 0 , τ), τ, X(z(z 0 , τ), τ), S(z(z 0 , τ), τ)) dτ, i = 1, ..., n, 0 ≤ z 0 ≤ L 0 , t > 0, (3.21
) which incorporates the initial conditions (3.17).

Differential system (3.19)-(3.20) is equivalent to the integral system ⎧ ⎨ ⎩ X i (c(t 0 , t), t) = ψ i (t 0 ) + t t0 F i (c(t 0 , τ), τ, X(c(t 0 , τ), τ), S(c(t 0 , τ), τ)) dτ, i = 1, ..., n, 0 < t 0 ≤ t, t > 0, (3.22) which incorporates the initial conditions (3.20).

Solving systems (3.21) and (3.22) gives X i along the lines (3.12) and (3.14), respectively. Therefore, the complete resolution requires the knowledge of those lines. Now, the integral equations for lines (3.12) and (3.14) are derived. From (3.13)

z(z 0 , t) = z 0 + t 0 u(z(z 0 , τ), τ) dτ, 0 ≤ z 0 ≤ L 0 , t > 0, (3.23)
From (3.2)

u(z(z 0 , t), t) = z(z 0 ,t) 0 G(ζ, t, X(ζ, t), S(ζ, t)) dζ, (3.24)
where equation (3.6) has been employed. By considering the change of variable

ζ = z(ζ 0 , t), equation (3.24) becomes u(z(z 0 , t), t) = z 0 0 G(z(ζ 0 , t), t, x(ζ 0 , t), s(ζ 0 , t)) ∂z ∂ζ 0 (ζ 0 , t) dζ 0 , (3.25)
where the following notations have been used

x(ζ 0 , t) = X(z(ζ 0 , t), t), s(ζ 0 , t) = S(z(ζ 0 , t), t).

(3.26)

Inserting equation (3.25) into (3.23) gives the integral equation for z(z 0 , t)

⎧ ⎨ ⎩ z(z 0 , t) = z 0 + t 0 dτ z0 0 G(z(ζ 0 , τ), τ, x(ζ 0 , τ), s(ζ 0 , τ)) ∂z ∂ζ 0 (ζ 0 , τ)dζ 0 , 0 ≤ z 0 ≤ L 0 , t > 0, (3.27
) which incorporates the initial condition z(z 0 , 0) = z 0 . Moreover, since ∂z/∂z 0 is involved, we also need

∂z ∂z 0 (z 0 , t) = 1 + t 0 G(z(z 0 , τ), τ, x(z 0 , τ), s(z 0 , τ)) ∂z ∂z 0 (z 0 , τ) dτ, (3.28)
which follows easily from (3.27).

Consider the characteristic-like lines when z 0 > L 0 . From (3.15)

c(t 0 , t) = L(t 0 ) + t t0 u(c(t 0 , τ), τ) dτ, 0 < t 0 ≤ t, t > 0. (3.29) From (3.2) u(c(t 0 , t), t) = z(L 0 ,t) 0 G(ζ,
+ t 0 0 G(c(τ 0 , t), t, x(τ 0 , t), s(τ 0 , t)) ∂c ∂τ 0 (τ 0 , t) dτ 0 ,
where notations (3.26) have been used, and in addition

x(τ 0 , t) = X(c(τ 0 , t), t), s(τ 0 , t) = S(c(τ 0 , t), t). G(c(τ 0 , τ), τ, x(τ 0 , τ), s(τ 0 , τ)) ∂c ∂τ 0 (τ 0 , τ) dτ 0 , 0 < t 0 ≤ t, t > 0, which incorporates the initial condition c(t 0 , t 0 ) = L(t 0 ). Another version of this equation will be provided in Section 3.4, where also ∂c/∂ t 0 will be derived. Finally, by using notations (3.26) and (3.32), equations (3.21) and (3.22) are rewritten as

⎧ ⎨ ⎩ x i (z 0 , t) = ϕ i (z 0 ) + t 0 F i (z(z 0 , τ), τ, x(z 0 , τ), s(z 0 , τ)) dτ, i = 1, ..., n, 0 ≤ z 0 ≤ L 0 , t > 0, (3.34)
⎧ ⎨ ⎩ x i (t 0 , t) = ψ i (t 0 ) + t t0 F i (c(t 0 , τ), τ, x(t 0 , τ), s(t 0 , τ)) dτ, i = 1, ..., n, 0 < t 0 ≤ t, t > 0.
(3.35)

Free boundary

Consider free boundary equation (3.3), rewritten as

L(t 0 ) = u(L(t 0 ), t 0 ) + σ(t 0 ) = u(c(t 0 , t 0 ), t 0 ) + σ(t 0 ). (3.36)
Use equation (3.31)

L(t 0 ) = σ(t 0 ) + L 0 0 G(z(ζ 0 , t 0 ), t 0 , x(ζ 0 , t 0 ), s(ζ 0 , t 0 )) ∂z ∂ζ 0 (ζ 0 , t 0 ) dζ 0 (3.37) + t 0 0
G(c(τ 0 , t 0 ), t 0 , x(τ 0 , t 0 ), s(τ 0 , t 0 )) ∂c ∂τ 0 (τ 0 , t 0 ) dτ 0 .

Hence,

L(t 0 ) = L 0 + t 0 0 σ(θ)dθ (3.38) + t 0 0 dθ L 0 0 G(z(ζ 0 , θ), θ, x(ζ 0 , θ), s(ζ 0 , θ)) ∂z ∂ζ 0 (ζ 0 , θ) dζ 0 + t 0 0 dθ θ 0 G(c(τ 0 , θ), θ, x(τ 0 , θ), s(τ 0 , θ)) ∂c ∂τ 0 (τ 0 , θ) dτ 0 , 0 < t 0 ≤ t.
Insert the expression above into equation (3.33) and obtain

c(t 0 , t) = L 0 + t 0 dθ L 0 0 G(z(ζ 0 , θ), θ, x(ζ 0 , θ), s(ζ 0 , θ)) ∂z ∂ζ 0 (ζ 0 , θ) dζ 0 (3.39) + t 0 0 σ(θ)dθ + t t 0 dθ t 0 0 G(c(τ 0 , θ), θ, x(τ 0 , θ), s(τ 0 , θ)) ∂c ∂τ 0 (τ 0 , θ) dτ 0 + t 0 0 dθ θ 0
G(c(τ 0 , θ), θ, x(τ 0 , θ), s(τ 0 , θ)) ∂c ∂τ 0 (τ 0 , θ) dτ 0 , 0 < t 0 ≤ t, t > 0.

Finally, we easily derive the equation for ∂c/∂t 0 from (3.39)

∂c ∂t 0 (t 0 , t) = σ(t 0 ) + t t0
G(c(t 0 , τ), τ, x(t 0 , τ), s(t 0 , τ)) ∂c ∂t 0 (t 0 , τ) dτ. (3.40)

Special problem

Equations (3.34), (3.35), (3.27), (3.28), (3.38), (3.39) and (3.40), which describe the free boundary problem, are mutually connected. In addition, they depend on substrates S j , j = 1, ..., m. So, also diffusion equations (3.4) must be involved. In this section we discuss a special mathematical problem by neglecting the dependence on substrates. The general situation will be considered in Section 3.6. Precisely, in this section we analyze the free boundary problem governed by the following system of integral equations

x i (z 0 , t) = ϕ i (z 0 ) + t 0 F i (z(z 0 , τ), τ, x(z 0 , τ)) dτ, i = 1, ..., n, (3.41)
z(z 0 , t) = z 0 + t 0 dτ z0 0 G(z(ζ 0 , τ), τ, x(ζ 0 , τ)) ∂z ∂ζ 0 (ζ 0 , τ) dζ 0 , (3.42)
∂z ∂z 0 (z 0 , t) = 1 + t 0 G(z(z 0 , τ), τ, x(z 0 , τ)) ∂z ∂z 0 (z 0 , τ) dτ, (3.43)
when 0 ≤ z 0 ≤ L 0 , 0 < t ≤ T, T > 0, (3.44)
and

x i (t 0 , t) = ψ i (t 0 ) + t t0 F i (c(t 0 , τ), τ, x(t 0 , τ)) dτ, i = 1, ..., n, (3.45)
c(t 0 , t) = t 0 dθ L 0 0 G(z(ζ 0 , θ), θ, x(ζ 0 , θ)) ∂z ∂ζ 0 (ζ 0 , θ) dζ 0 (3.46) +L 0 + t 0 0 σ(θ)dθ + t t 0 dτ t 0 0 G(c(τ 0 , τ), τ, x(τ 0 , τ)) ∂c ∂τ 0 (τ 0 , τ) dτ 0 + t 0 0 dτ τ 0 G(c(τ 0 , τ), τ, x(τ 0 , τ)) ∂c ∂τ 0 (τ 0 , τ) dτ 0 , ∂c ∂t 0 (t 0 , t) = σ(t 0 ) + t t0
G(c(t 0 , τ), τ, x(t 0 , τ)) ∂c ∂t 0 (t 0 , τ) dτ, (3.47) where x n+1 (z 0 , t) = z(z 0 , t), x n+2 (z 0 , t) = ∂z ∂z 0 (z 0 , t).

L(t 0 ) = t 0 0 dθ L 0 0 G(z(ζ 0 , θ), θ, x(ζ 0 , θ)) ∂z ∂ζ 0 (ζ 0 , θ) dζ 0 (3.48) +L 0 + t 0 0 σ(θ)dθ + t 0 0 dθ θ 0 G(c(τ 0 , θ), θ, x(τ 0 , θ)) ∂c ∂τ 0 (τ 0 , θ) dτ 0 , when z 0 ≥ L 0 and 0 < t 0 ≤ t, 0 < t ≤ T, T > 0. ( 3
(3.51)

In addition, let

ϕ n+1 (z 0 ) = z 0 , F n+1 (τ, x(ζ 0 , τ)) = G(z(ζ 0 , τ), τ, x 1 , ..., x n ) ∂z ∂ζ 0 (ζ 0 , τ), (3.52)
ϕ n+2 (z 0 ) = 1, F n+2 (τ, x(z 0 , τ)) = G(z(z 0 , τ), τ, x 1 , ..., x n ) ∂z ∂z 0 (z 0 , τ). (3.53)
Assume F i continuous and bounded

M i = max |F i |, 0 < τ ≤ t, 0 ≤ z 0 < L 0 , |x i -ϕ i | < ρ i , i = 1, ..., n +2, (3.54)
where ρ i are positive constants. Setting

T = min{ρ 1 /M 1 , ..., ρ n /M n , ρ n+1 /(M n+1 L 0 ), ρ n+2 /M n+2 }, (3.55)
suppose that F i satisfy the Lipschitz condition

|F i (τ, x) -F i (τ, x)| < λ i n+2 h=1 |x h -xh |, λ i > 0, i = 1, ..., n + 2, (3.56)
on Proof: Denote by Σ 1 the space of continuous vectors x which satisfy the inequalities

D 1 = {0 < τ ≤ t < T, 0 < ζ 0 ≤ z 0 < L 0 , |x i -ϕ i | < ρ i , i = 1, ..., n + 2}. ( 3
|x i -ϕ i | < ρ i , i = 1, ..., n + 2, 0 ≤ z 0 ≤ L 0 , 0 ≤ t ≤ T, (3.58)
and consider the norm

||x|| = n+2 i=1 max exp(-γ 1 z 0 -γ 2 t)|x i (z 0 , t)|, 0 ≤ z 0 ≤ L 0 , 0 ≤ t ≤ T, (3.59)
where γ 1 and γ 2 are positive constants that will be specified later on. Consider the map y = Ax on Σ 1 defined by

y i (z 0 , t) = ϕ i (z 0 ) + t 0 F i (τ, x(z 0 , τ)) dτ, i = 1, ..., n, (3.60)
y n+1 (z 0 , t) = ϕ n+1 (z 0 ) + t 0 dτ z0 0 F n+1 (τ, x(ζ 0 , τ)) dζ 0 , (3.61)
y n+2 (z 0 , t) = ϕ n+2 (z 0 ) + t 0 F n+2 (τ, x(z 0 , τ)) dτ. ( 3
|y i -ỹi | exp(-γ 1 z 0 -γ 2 t) ≤ (3.63) λ i t 0 exp(-γ 2 (t -τ )) n+2 h=1 exp(-γ 1 z 0 -γ 2 τ )|x h (z 0 , τ), τ) -xh (z 0 , τ), τ)| dτ ≤ λ i ||x -x|| t 0 exp(-γ 2 (t -τ )) dτ ≤ λ i γ 2 ||x -x||, i = 1, ..., n, |y n+1 -ỹn+1 | exp(-γ 1 z 0 -γ 2 t) ≤ (3.64) λ n+1 t 0 dτ z0 0 e -γ1(z0-ζ0)-γ2(t-τ ) n+2 h=1 e -γ1ζ0-γ2τ |x h (ζ 0 , τ), τ)-x h (ζ 0 , τ), τ)| dζ 0 ≤ λ n+1 ||x -x|| t 0 dτ z0 0 e -γ1(z0-ζ0)-γ2(t-τ ) dζ 0 ≤ λ n+1 γ 1 γ 2 ||x -x||, |y n+2 -ỹn+2 | exp(-γ 1 z 0 -γ 2 t) ≤ (3.65) λ n+2 t 0 e -γ2(t-τ ) n+2 h=1 e -γ1z0-γ2τ |x h (z 0 , τ), τ) -xh (z 0 , τ), τ)| dτ ≤ λ n+2 ||x -x|| t 0 e -γ2(t-τ ) dτ ≤ λ n+2 γ 2 ||x -x||.
λ i γ 2 + λ n+1 γ 1 γ 2 + λ n+2 γ 2
can be made < 1 if γ 1 and γ 2 are selected large enough. Therefore, y = Ax is a contractive map and the theorem is proved. Similar reasonings are now applied to system (3.45)-(3.47). In addition, we assume that the functions x i (z 0 , t), i = 1, ...n, and z(z 0 , t) are known, as determined in the previous discussion. In particular, it is assumed to be known the following function

ψ n+1 (t 0 , t) = L 0 + t 0 0 σ(θ)dθ+ t 0 dθ L 0 0 G(z(ζ 0 , θ), θ, x(ζ 0 , θ)) ∂z ∂ζ 0 (ζ 0 , θ) dζ 0 .
(3.66) Define the functions x n+1 (t 0 , t) and x n+2 (t 0 , t) as follows

x n+1 (t 0 , t) = c(t 0 , t), x n+2 (t 0 , t) = ∂c ∂t 0 (t 0 , t), (3.67
F n+1 (τ, x(τ 0 , τ)) = G(c(τ 0 , τ), τ, x 1 , ..., x n ) ∂c ∂τ 0 (τ 0 , τ), (3.69)
ψ n+2 (t 0 , t) = σ(t 0 ), F n+2 (τ, x(t 0 , τ)) = G(c(t 0 , τ), τ, x 1 , ..., x n ) ∂c ∂t 0 (t 0 , τ).

(3.70) Assume F i continuous and bounded

K i = max |F i |, 0 < τ < t, 0 < τ 0 < t, |x i -ψ i | < μ i , i = 1, ..., n + 2, (3.71)
where μ i are positive constants. Furthermore,

T = min{μ 1 /K 1 , ..., μ n /K n , μ n+1 /(2K n+1 ), μ n+2 /K n+2 }, (3.72)
and F i satisfy the Lipschitz condition

|F i (τ, x) -F i (τ, x)| < λ i n+2 h=1 |x h -xh |, λ i > 0, i = 1, ..., n + 2, (3.73)
on Proof: Denote by Σ 2 the space of continuous vectors x(t 0 , t) which satisfy the condition

D 2 = {0 < τ 0 < t 0 , 0 < τ < t 0 ≤ t < T, |x i -ψ i | < μ i , i = 1, ..., n + 2}. ( 3
|x i -ψ i | < μ i , i = 1, ..., n + 2, 0 ≤ t 0 ≤ T, 0 ≤ t ≤ T, (3.75)
and consider the norm

||x|| = n+2 i=1 max exp(-γ 3 t 0 -γ 4 t)|x i (t 0 , t)|, 0 ≤ t 0 ≤ T, 0 ≤ t ≤ T, (3.76)
where γ 3 and γ 4 are positive constants that will be specified later on.

Consider the map y = Bx on Σ 2 defined by [START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF]. Precisely, n i=1 f i = 1 at any time if it is so initially. In addition, f i ≥ 0 at any time.

y i (t 0 , t) = ψ i (t 0 ) + t t0 F i (τ, x(t 0 , τ)) dτ, i = 1, ..., n, (3.77)
y n+1 (t 0 , t) = ψ n+1 (t 0 , t) + t t 0 dτ t0 0 F n+1 (τ, x(τ 0 , τ)) dτ 0 , (3.78)
+ t 0 0 dτ τ 0 F n+1 (τ, x(τ 0 , τ)) dτ 0 , y n+2 (t 0 , t) = ψ n+2 (t 0 ) + t t0 F n+2 (τ, x(t 0 , τ)) dτ. ( 3
|y i (t 0 , t)) -ỹi (t 0 , t))|e -γ3t0-γ4t (3.80) ≤ λ i t t0 e -γ4(t-τ ) n+2 h=1 e -γ3t0-γ4τ |x h (t 0 , τ) -xh (t 0 , τ)|dτ ≤ λ i ||x -x||/γ 4 , i = 1, ..., n, |y n+1 (t 0 , t) -ỹn+1 (t 0 , t)|e -γ3t0-γ4t (3.81) ≤ λ n+1 t t0 dτ t0 0 e -γ3(t 0 -τ 0 ) e -γ4(t-τ ) n+2 h=1 e -γ3τ 0 -γ4τ |x h (τ 0 , τ) -xh (τ 0 , τ)|dτ 0 +λ n+1 t0 0 dτ τ 0 e -γ3(t 0 -τ 0 ) e -γ4(t-τ ) n+2 h=1 e -γ3τ 0 -γ4τ |x h (τ 0 , τ) -xh (τ 0 , τ)|dτ 0 ≤ 2λ n+1 ||x -x||/(γ 3 γ 4 ), |y n+2 (t 0 , t) -ỹn+2 (t 0 , t)|e -γ3t0-γ4t (3.82) ≤ λ n+2 t t0 e -γ4(t-τ ) n+2 h=1 e -γ3t0-

Effect of substrates

Biofilm development is a complex process strongly influenced by substrates availability. Substrate concentration trends inside biofilm derive from the combination of two processes that take place in the same time within biofilm: microbial conversion of substrates and transport substrate by molecular diffusion. Molecular diffusion of substrates into biofilm is influenced by several factors, such as flux of substrates received from bulk, diffusivity of substrate species, structural characteristics of biofilm.

Substrate trends determine spatially heterogeneous growth of microbial species that contributes to the formation of environmental microniches allowing the coexistence of different microbial groups. Distinct chemical niches exist at different depths in biofilms, there are several studies in which chemical gradients have been related to the distribution of specific bacterial species. For example, the stratified distributions of the bacteria that constitute methanogenic consortia have also been described and can be understood in terms of the diffusive exchange of metabolites among species.

In this section the influence of substrates on biofilm growth is considered and it is assumed r M,i = r M,i (z, t, X, S). (3.83) Note that equations (3.41)-(3.49) still hold, but F i and G must be modified according to (3.83)

F i = F i (z, τ, x, s), G = G(z, τ, x, s), (3.84)
where notations (3.26) and (3.32) have been used. The diffusion of substrates is governed by equations (3.4) with initial-boundary conditions (3.9)- (3.11). The solution can be expressed in terms of integral equations by using known results on the heat equation in general regions, e.g. [START_REF] Li Rubinstein | The Stefan Problem[END_REF]. So, we obtain

S j (z, t) = L0 0 S j0 (ζ 0 )N j (z, ζ 0 , t)dζ 0 + t 0 D j w j (τ )N j (z, L(τ ), t -τ )dτ (3.85) + t 0 S jL (τ )[N j (z, L(τ ), t -τ ) L(τ ) -D j N jζ (z, L(τ ), t -τ )] dτ + t 0 dτ L(τ ) 0 r S,j (ζ, τ, X(ζ, τ ), S(ζ, τ ))N j (z, ζ, t -τ )dζ, j = 1, ..., m, w j (t) = 2 L0 0 S j0 (ζ 0 )G j (L(t), ζ 0 , t) dζ 0 + 2 t 0 ṠLj (τ )G j (L(t), L(τ ), t-τ ) dτ (3.86) +2 t 0 dτ L(τ ) 0 r Sj (ζ, τ, X(ζ, τ ), S(ζ, τ ))N jz (L(t), ζ, t -τ )dζ +2 t 0 D j w j (τ )N jz (L(t), L(τ ), t -τ ) dτ, j = 1, ..., m,
where the following notations have been used

w j (t) = ∂S j ∂z (L(t), t), K j (z, t) = exp(-z 2 )/4D j t 4πD j t , N j (z, ζ, t -τ ) = K j (z -ζ, t -τ ) + K j (z + ζ, t -τ ), G j (z, ζ, t -τ ) = K j (z -ζ, t -τ ) -K j (z + ζ, t -τ ).
s j (z 0 , t) = L0 0 S j0 (ζ 0 )N j (z(z 0 , t), ζ 0 , t)dζ 0 + t 0 D j w j (τ )N (z(z 0 , t), L(τ ), t-τ )dτ (3.87) + t 0 S jL (τ )[Nj(z(z 0 , t), L(τ ), t -τ ) L(τ ) -D j N jζ (z(z 0 , t), L(τ ), t -τ )] dτ + t 0 dτ L0 0 r S,j (ζ(ζ 0 , τ), τ, x(ζ 0 , τ), s(ζ 0 , τ))N j (z(z 0 , t), ζ(ζ 0 , τ), t-τ ) ∂ζ ∂ζ 0 dζ 0 + t 0 dτ τ 0
r S,j (c(τ 0 , τ), τ, x(τ 0 , τ), s(τ 0 , τ))N j (z(z 0 , t), c(τ 0 , τ), t-τ ) ∂c ∂τ 0 dτ 0 .

Moreover, for z 0 > L 0 and 0 < t 0 ≤ t, system (3.85) reduces to r S,j (c(τ 0 , τ), τ, x(τ 0 , τ), s(τ 0 , τ))N j (c(t 0 , t), c(τ 0 , τ), t-τ ) ∂c ∂τ 0 dτ 0 .

s j (t 0 , t) = L0 0 S j0 (ζ 0 )N j (c(t 0 , t), ζ 0 , t)dζ 0 + t 0 D j w j (τ )N (c(t 0 , t), L(τ ), t-τ )dτ (3.88)
+ t 0 S jL (τ )[Nj(c(t 0 , t), L(τ ), t -τ ) L(τ ) -D j N jζ (c(t 0 , t), L(τ ), t -τ )] dτ
System (3.86) can be treated similarly. Notice that integral system (3.85)-(3.86) has been widely discussed since it is involved in the well-known Stefan problem. This system can be associated to equations (3.41)-(3.49), modified according to (3.84), and a result of uniqueness and existence to solutions can be deduced.

Numerical simulations

Heterotrophic-autotrophic competition for space with oxygen as common substrate proposed in [START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF] has been used to provide numerical simulations of the free boundary problem introduced in Section 3.2. The proposed numerical example, described in details by the equations in Table 3.1, is based on mass balance equations for substrates, products, and bacterial groups and includes the bio-chemical reactions of heterotrophic-autotrophic competition. The model considers the kinetics of microbial growth and decay and takes into account two groups of bacteria Heterotrophic Bacteria (X 1 ) and Autotrophic Bacteria (X 2 ), and three components (substrates), Ammonia (S 1 ), Organic Carbon (S 2 ) and Oxygen (S 3 ). Inert is modelled as another microbial species, whose growth derives from the heterotrophic and autotrophic biomass decay. Oxygen is used for organic carbon oxidation, nitrification and endogenous respiration. Oxidation of ammonia nitrogen to nitrate by the autotrophs provides energy for autotrophic growth.

Process X 1 X 2 X 3 S 1 S 2 S 3 Process rate r j HG 1 - --1 Y1 - -α1-Y1 Y1 μ m,1 X 1 S1 K1,1+S1 S3 K3,1+S3 AG - 1 - - -1 Y2 -α2-Y2 Y2 μ m,2 X 2 S2 K2,1+S2 S3 K3,2+S3 HER -1 - - - - - b res,1 X 1 S3 K3,1+S3 AER --1 - - - - b res,2 X 2 S3 K3,2+S3 AD -1 - 1 - - - b m,1 X 1 HD --1 1 - - - b m,2 X 2
Table 3.1: Stoichiometry and rate laws for microbial processes. HG = heterotroph growth; AG = autotroph growth; HER = heterotroph endogenous respiration; AER = autotroph endogenous respiration; HD = heterotroph decay; AG = autotroph decay.

The differential system (3.1)-(3.4) has been integrated numerically [START_REF] Acunto | Computational Methods for PDE in Mechanics[END_REF][START_REF] Acunto | Computational Partial Differential Equations for Engineering Science[END_REF]. An initial biofilm thickness L 0 = 0.3 mm has been assumed and Dirichlet-Neumann boundary conditions have been adopted.

The oxygen concentration at the interface biofilm/bulk liquid, the values of ammonia and acetate fluxes from bulk liquid to biofilm and biological parameters are reported in Table 3 The values of kinetics and stoichiometric parameters reported in Table 3.1 and used in numerical simulations are the following:

μ m,1 = 25; μ m,2 = 5; b m,1 = 1; b m,2 = 1; b res,1 = 1; b res,2 = 0.5; Y 1 = 0.8; Y 2 = 0.6; K 1,1 = 5; K 2,1 = 1; K 3,1 = 0.1; K 3,2 = 0.1.
Numerical simulations have been developed in order to predict the biomass distribution and substrate concentration trends over biofilm depth. The results are shown in Figure 3.3 and Figure 3.2, respectively. In particular, three sets of simulations at different attachment rates but at the same simulation time have been performed. The objective is the evaluation of the effects of attachment rate on biofilm growth in terms of biofilm thickness, bacterial species distribution and substrate concentration trends. When the attachment rates is equal to 5 mmd -1 (Figure 3.3 (A)) biofilm thickness is more than two time the biofilm thickness when the attachment rates is less than 1 mmd -1 (Figure 3. 3 (A,B)). This difference determines different substrate concentration trends into biofilm and, consequently, different biofilm structure occurs. In this case heterotroph bacteria are found to be predominant at the outmost layer of biofilm.

It is interesting to note the sharp variation on biofilm volume fraction in the superficial layer of the biofilm. This was expected, since, after the initial phase of attachment, the biofilm volume fraction is mainly determined by internal bacteria metabolism and not by the external biomass flux concentration. The sum equal to one of the three different biofilm volume fractions emphasizes the quality of numerical integration. The diffused substrate concentration trends in the biofilm, for three different attachment rates, are shown in Figure 3.2. The different thickness of biofilm and different kinds of bacterial species growing into biofilm determine different substrate concentration trends. When the attachment rates is equal to 5 mmd -1 (Figure 3.2 (A)) there is a sharper decrease of oxygen concentration than attachment rates is less than 1 mmd -1 (Figure 3. 2 (A,B)). This occurs since a great concentration of heterotrophs at the outmost layer of biofilm implies a greater consumption of oxygen.

Chapter 4

Mathematical modeling of competition and coexistence of sulfate-reducing bacteria, acetogens and methanogens in multispecies biofilms

This work presents an integrated mathematical model able to simulate the physical, chemical and biological processes prevailing in a sulfate reducing biofilm under dynamic conditions. The model includes sulfate reduction by complete and incomplete sulfate reducing bacteria; lactate removal by sulfate reduction and by acetogenic bacteria and acetate consumption via methanogenesis. Numerical integration based on the method of characteristics has been developed. The major problem of sulfate-reducing fixed-growth reactors

is the formation of undesired bacterial species which compete for space and substrate within the biofilm with sulfate reducing bacteria. The effect of COD/SO 2- 4 ratio on the reactor performances in terms of bacterial species distribution and substrate diffusion trends in the biofilm has been assessed. The simulation results reveal a stratification of microbial activities in biofilm reflecting the different ecological niches created by substrate gradients.

Introduction

High Sulfate Containing Wastewaters (HSCW) are generated from various industrial activities, such as pulp and paper industries, mining and mineral processing, production of explosives, scrubbing of flue gases, food processing and petrochemical industries [START_REF] Pnl Lens | Biotechnological treatment of sulfate-rich wastewaters[END_REF]. These anthropogenic activities have contributed to local imbalances in the natural sulfur cycle, resulting in acidification, leaching of toxic metals, elevated sulfate levels in natural waters, potential production of corrosive and toxic sulfide, emissions of SO 2 , H 2 S and odorous volatile sulfur compounds, cat clays and heavy metal release upon oxygen exposure of sediments after dredging [START_REF] Lens | Perspectives of sulfate reducing bioreactors in environmental biotechnology[END_REF]. HSCWs often contain elevated concentrations of metals (iron, aluminium and manganese and other heavy metals) and metalloids, deriving from the mining and processing of metal ores and coals, which increase the complexity of the degradation routes [START_REF] Db | Acid mine drainage remediation options: a review[END_REF][START_REF] Papirio | Acid mine drainage treatment in fluidized-bed bioreactors by sulfate-reducing bacteria: a critical review[END_REF].

During the last years numerous psychochemical and biological techniques have been investigated for the neutralization and removal of metals and sulfate from wastewaters. Two main categories can be individuated: passive and active processes. Passive treatment processes commonly replace the conventional neutralisation techniques involving the addition of a chemical-neutralising agent. Passive treatment processes require less energy and chemicals and relatively low maintenance costs. Among these it is possible to enumerate natural wetlands, aerobic and anaerobic wetlands, open limestone channels. Although the passive processes are considered low-cost treatment technologies, their efficiency is not very high compared to very expensive surface requests in terms of land. Active treatment processes are, instead, much more efficient. The treatment efficiency is improved through the application of energy, chemical and biological agents. Active treatment processes require certainly higher maintenance cost and manpower when compared to passive processes, but these costs are offset by the high treatment efficiency. Technologies such as reverse osmosis, ion exchange, limestone and chemical neutralisation and active biological treatment represent typical examples of active treatment processes. In particular, active biological sulfate removal from HSCWs represents a valid and cost-effective alternative to the costly and sometimes complex physico-chemical sulfate removal methods.

Biological sulfate removal can be accomplished in two steps: a dissimilatory sulfate reduction to sulfide performed by Sulfate Reducing Bacteria (SRB), followed by sulfide removal through partial oxidation to sulfur or precipitation of heavy metals sulfide. The dissimilatory sulfate reduction can take place in methanogenic or sulfidogenic bioreactors. The production of sulfide has been shown to be inhibitory for anaerobic digestion. As a consequence many studies have been carried out to assess the sulfide toxicity, individuate the most-suitable strategies to prevent it and steer the competition between SRB, acetogenic and methanogenic microorganisms in the direction of methanogenesis. On the other hand, sulfidogenesis can be seen as an ideally suited process to remove both sulfate and heavy metals from HSCWs and the interest in the application of this process as the main step for the biological treatment of specific wastestreams from chemical, mining and galvanic industries as well as scrubbing water for flue-gas desulfurization, has been growing [START_REF] Hulshoff | New developments in reactor and process technology for sulfate reduction[END_REF].

Biological sulfate reduction is mediated by heterotrophic or autotrophic SRB, able to reduce sulfate to sulfide in the presence of a carbon source (CO 2 , acetate, lactate, propionate, etc.). Different types of carbon can be used as energy sources; most of the substrates are typical fermentation products or intermediate breakdown products of larger molecules [START_REF] Ah Kaksonen | Sulfate reduction based bioprocesses for the treatment of acid mine drainage and the recovery of metals[END_REF]. A minimum chemical oxygen demand (COD) to sulfate mole ratio of 0.67 is required for achieving theoretically possible removal of sulfate [START_REF] Liamleam | Electron donors for biological sulfate reduction[END_REF]. Lens et al [START_REF] Pnl Lens | Biotechnological treatment of sulfate-rich wastewaters[END_REF] reported that SRB are very diverse in their carbon source utilization and metabolic activities. The availability of carbon and energy source provides the energy for the growth and maintenance of SRB. SRB carry out sulfate reduction basing on the following reaction [START_REF] Herlihy | The importance of sediment sulfate reduction to the sulfate budget of an impoundment receiving acid mine drainage[END_REF]:

SO 2- 4 + 8e -+ 4H 2 O → S 2-+ 8OH - (4.1)
In most cases the electron donor and carbon source are the same compound. However when hydrogen is used as the electron donor, CO 2 can be used as carbon source by SRB. The selection of the electron donor depends on the ability of SRB to utilize the substrate, its costs per unit of reduced sulfate, the availability in sufficient quantities and the remaining pollution load of the additive in the wastestream [START_REF] Pnl Lens | Biotechnological treatment of sulfate-rich wastewaters[END_REF][START_REF] Ah Kaksonen | Sulfate reduction based bioprocesses for the treatment of acid mine drainage and the recovery of metals[END_REF]. The choice of a suitable carbon source and electron donor for this process is still a challenge.

SRB can be classified into two groups based on their functional ability to oxidise the organic compounds completely to CO 2 -SRB completely oxidizers (SRB (C) ) or incompletely to acetate and CO 2 -SRB incompletely oxidizers (SRB (I) ). Postgate [START_REF] Jr Postgate | The sulphate-reducing bacteria[END_REF] indicated that lactate offers potential advantages as carbon source and electron donor in the sulfate reduction process. Lactate can be used by many SRB species; its oxidation results in high biomass yield and high alkalinity production. However, the potential accumulation of acetate in the effluent due to the incomplete oxidation of lactate to acetate and CO 2 represents the main disadvantage of using lactate as carbon source. For this inconvenient a large amount of lactate is needed to achieve complete reduction of sulfate, contributing to increase the costs of bioreactors performance. In addition, due to the release of acetate, the COD of the effluent stream increases. The incomplete oxidation of carbon sources to acetate can be attributed to the lower value of free energy for the oxidation of acetate to carbon dioxide which prevents further oxidation of acetate to carbon dioxide [START_REF] Jr Postgate | The sulphate-reducing bacteria[END_REF]. Furthermore, the presence of acetate and lactate can allow the development of both methanogenic archae and cetogenic bacteria that can ferment lactate, resulting in the production of acetate. Due to their kinetic properties, high levels of lactate encourage the growth of acetogenic bacteria. On the other hand lactate oxidation becomes dominant under conditions of lactate limitation and excess sulfate [START_REF] Hj Laanbroek | Alcohol conversion by desulfobulbus propionicus lindhorst in the presence and absence of sulfate and hydrogen[END_REF]. Indeed, in an investigation based on a full scale anaerobic digester [START_REF] Zellner | Degradation of lactate by an anaerobic mixed culture in a fluidized-bed reactor[END_REF], lactate oxidizers were shown to have lower K s and μ max values than lactate fermenters.

Numerous reactor designs dedicated to biological sulfate reduction have been reported [START_REF] Ah Kaksonen | Sulfate reduction based bioprocesses for the treatment of acid mine drainage and the recovery of metals[END_REF]. They can be classified into two main groups: i) suspended growth reactors, that involve the growth of planktonic bacteria, such as batch reactors, baffled reactors, up-flow anaerobic sludge bed reactors and gas-lift reactors; and ii) attached growth reactors, that involve a bacterial biomass attached to media (biofilm), i.e. fixed bed reactors or fluidized bed reactors. Various immobilized biomass reactors have gained increasing attention due to the advantages of displacing biomass in biofilms. Bacteria growing in biofilms cannot be washed out with the water flow. This allows to retain the biomass within the reactor and therefore to operate at shorter hydraulic retention time (HRT). Maximal biomass retention is desirable for process stability and minimal sludge production. Moreover the high biomass retention and concentration characterizing biofilm reactors strongly affects the achievable loading rates, with the possibility of obtaining high treatment efficiencies [START_REF] Ah Kaksonen | Sulfate reduction based bioprocesses for the treatment of acid mine drainage and the recovery of metals[END_REF]. In addition, biofilms show good tolerance for shocks of hydraulic and organic loading and can allow treating contemporary different pollutants thanks to niche differentiation.

Biological sulfate reduction in anaerobic fixed growth reactors has been investigated extensively at lab-scale. In designing these biofilm reactors, in predicting their behavior under different operating conditions and in understanding the complex microbial relations existing in anaerobic environments in the presence of sulfate, mathematical modelling seems to be essential. Indeed, mathematical models can be used to estimate parameters that cannot be observed directly in experiments and develop an online control strategy. Therefore the use of mathematical modelling clearly benefits engineers, designers and operators [START_REF] Frunzo | Dynamic mathematical modeling of sulfate reducing gas-lift reactors[END_REF][START_REF] Garuti | Full-scale ananox (r) system performance[END_REF].

The scope of this work is to evaluate the SRB growth in multispecies biofilms by modelling the competition between the different bacterial groups involved in the lactate metabolism under biosulfidogenic conditions. In particular this work is aimed at evaluating the dynamical response of the model under established boundary conditions assessing the effect of different COD/SO 2- 4 ratios on microbial population shifts. The numerical simulations have been obtained with great accuracy by the method of characteristics. Simulation results show that the model can predict the short-term responses of biofilm performance to substrate variations in the bulk liquid as well as the long-term development of film thickness and microbial species.

Statement of the problem

Sulfate reducing applications usually utilize mixed cultures comprising of SRB and anaerobic fermentative microorganisms, such as methanogens and acetogens [START_REF] Ah Kaksonen | Sulfate reduction based bioprocesses for the treatment of acid mine drainage and the recovery of metals[END_REF]. To perform a complete reduction of sulfate to sulfide, SRB have to effectively compete with the other anaerobic bacteria for the available organic substrate. The presence of sulfate seems to be crucial in this competition. As stated in [START_REF] Muyzer | The ecology and biotechnology of sulphatereducing bacteria[END_REF], the degradation of organic matter in sulfate-reducing environments is different from the degradation in methanogenic environments. Macromolecules, such as proteins, polysaccharides and lipids are hydrolysed by hydrolytic bacteria. Subsequently, the monomers amino acids, sugars and fatty acids are fermented by fermentative bacteria into a range of fermentation products, such as acetate, propionate, butyrate, lactate and hydrogen. In the presence of sulfate, SRB consume these fermentation products. However, in the absence of sulfate, hydrogen and acetate, the acetate having been produced directly by fermentation or indirectly by acetogenesis, are consumed by the methanogens. Among simple organic substrates, SRB have been demonstrated to use lactate, ethanol, methanol, acetate, proprionate and butyrate [START_REF] Papirio | Biological inverse fluidized-bed reactors for the treatment of low ph-and sulphate-containing wastewaters under different cod conditions[END_REF]. Lactate can support the growth of a wide spectrum of SRB, encouraging microbial diversity and consequent treatment system resilience [START_REF] Oo Oyekola | Kinetic analysis of biological sulphate reduction using lactate as carbon source and electron donor: effect of sulphate concentration[END_REF]. Lactate can be metabolized via fermentation or sulfate-reducing oxidation or both by a wide range of microorganisms. Lactate fermentation is the anaerobic degradation of lactate, independent of sulfate reduction [START_REF] Oo Oyekola | Effect of culture conditions on the competitive interaction between lactate oxidizers and fermenters in a biological sulfate reduction system[END_REF]. Lactate is oxidized either incompletely or completely in the presence of sulfate by a diverse range of SRB strains [START_REF] Oo Oyekola | Study of anaerobic lactate metabolism under biosulfidogenic conditions[END_REF]. According to [START_REF] Papirio | Biological inverse fluidized-bed reactors for the treatment of low ph-and sulphate-containing wastewaters under different cod conditions[END_REF], only particular species of SRB are able to oxidize lactate to CO 2 whereas others oxidize lactate to acetate and very few can use acetate as carbon source. Besides the limited capability of SRB to degrade it, acetate can accumulate in solution even if other microorganisms, such as methanogens, are present. Competition between the different microbial groups depends on the kinetic properties of the interacting microorganisms, such as the maximum specific growth rate (μ max ) and substrate affinity (K s ) [START_REF] Oo Oyekola | Kinetic analysis of biological sulphate reduction using lactate as carbon source and electron donor: effect of sulphate concentration[END_REF]. Extensive experimental efforts have been devoted to the kinetic study of lactate metabolic pathway under biosulfidogenic conditions in chemostat cultures [START_REF] Oo Oyekola | Kinetic analysis of biological sulphate reduction using lactate as carbon source and electron donor: effect of sulphate concentration[END_REF][START_REF] Oo Oyekola | Study of anaerobic lactate metabolism under biosulfidogenic conditions[END_REF]. In these studies, the effects of different sulfate concentrations, lactate concentrations and volumetric loading rates on the kinetics of lactate utilization and the stoichiometry of biological sulfate reduction have been investigated. In the case of immobilized biomass reactors, the competition between the different microbial groups is regulated not only by kinetic properties and dilution rates, but substrate diffusion and niche differentiation have been found to have a crucial role in dictating lactate utilization pathway. In this work the dynamics of the anaerobic sulfate reduction in a multispecies biofilm are discussed. Chemical, physical and biological transient processes are analyzed. In particular the model takes into account the bioprocess pathways reported in Figure 4.1.

The model can simulate the activities of microorganisms living in a sulfate reducing multispecies biofilm and evaluate the interactions between the related processes: lactate and acetate consumption, sulfate reduction, and bacterial growth and decay. Three reacting components are simultaneously considered: lactate, sulfate and acetate. The proposed model takes into account the growth of two types of sulfate reducing bacteria classified into two groups based on their functional ability to oxidise the lactate completely to carbon dioxide (LDSRB (C) ) or incompletely to acetate and carbon dioxide (LDSRB (I) ). The presence of lactate in an anaerobic environment allows the development of acetogenic bacteria (AB) with the production of Acetate and Hydrogen. The undesired acetate production by both incomplete sulfate reducing bacteria and acetogenic bacteria allows the growth of methanogenic archea (MA) that produce methane as a final metabolic product. Inert residues (Inert), deriving from microbial biomass decay, are also taken into account. AB compete for space and lactate with sulfate reducing bacteria, while MA compete only for space. The growth of these microorganisms is favoured by the formation of zones in biofilm characterized by different substrate concentration levels.

The mathematical model

The biofilm growth is governed by the following equations [START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF]:

∂X i ∂t + ∂ ∂z (uX i ) = ρ i r M,i (z, t, X, S) , 0 ≤ z ≤ L (t) , t > 0, i = 1, 2, 3, 4, 5, (4.2
)

∂u ∂z = 5 i=1 r M,i (z, t, X, S) , 0 < z ≤ L (t) , t > 0, (4.3)
where X i = ρ i f i (z, t) denotes the concentration of the microbial species and inert residues; f i denotes the volume fraction of microbial species i = 1, 2, 3, 4, 5; ρ i is the density assumed constant; u(z, t) is the velocity of the microbial mass displacement with respect to the biofilm support interface; the term r M,i (z, t, X, S) represent the biomass growth rate, X = (X 1 , X 2 , X 3 , X 4 , X 5 ) and S = (S 1 , S 2 , S 3 , S 4 , S 5 ). The net biomass growth rates are given by:

r M,1 = (μ 1 -K d,1 )X 1 , (4.4) r M,2 = (μ 2 -K d,2 )X 2 , (4.5) r M,3 = (μ 3 -K d,3 )X 3 , (4.6) r M,4 = (μ 4 -K d,4 )X 4 , (4.7)
while for inert residues

r M,5 = K d,1 X 1 + K d,2 X 2 + K d,3 X 3 + K d,4 X 4 , (4.8)
where μ 1 , μ 2 , μ 3 and μ 4 are the biomass growth rates for biomass X 1 , X 2 , X 3 and X 4 . K d,1 , K d,2 , K d,3 and K d,4 are the decay-inactivation rates for the single microbial species. The biomass growth rates are given by:

μ 1 = μ max,1 S 1 K 1,1 + S 1 S 2 K 1,2 + S 2
, (4.9)

μ 2 = μ max,2 S 1 K 2,1 + S 1 S 2 K 2,2 + S 2
, (4.10)

μ 3 = μ max,3 S 2 K 3,2 + S 2 , (4.11
)

μ 4 = μ max,4 S 3 K 4,3 + S 3 , (4.12)
where: μ max,i is the maximum growth rate for biomass i; K i,j is the half saturation constant for substrate j(S 1 , S 2 , S 3 ) of biomass i. The diffusion of substrates is governed by the equations:

∂S j ∂t -D j ∂ 2 S j ∂z 2 = r S,j (z, t, X, S) , 0 < z < L (t) , 0 < t ≤ T, j = 1, 2, 3, (4.13
) where D j denotes the diffusivity coefficient and r S,j (z, t, X, S) the net conversion rate of substrate j, expressed by:

r S,1 = -1.5 (1 -Y 1,2 ) Y 1,2 μ 1 -1.5 (1 -Y 2,2 ) Y 2,2 μ 2 , (4.14) r S,2 = - 1 Y 1,2 μ 1 - 1 Y 2,2 μ 2 - 1 Y 3,2 μ 3 , (4.15) r S,3 = -0.8 (1 -Y 2,3 ) Y 2,3 μ 2 -0.8 (1 -Y 3,3 ) Y 3,3 μ 3 - 1 Y 4,3 μ 4 , (4.16)
where Y i,j denotes the yield for biomass i and substrates j. The following initial-boundary conditions will be considered for equations (4.2, 4.3, 4.13):

X i (z, 0) = ϕ i (z) , u(0, t) = 0, 0 ≤ z ≤ L 0 , t ≥ 0, i = 1, 2, 3, 4, 5, (4.17)
S j (z, 0) = S 0j (z) , 0 ≤ z ≤ L 0 , j = 1, 2, 3, (4.18
)

∂S j ∂z (0, t) = 0, S j (L (t) , t) = G j (t), t > 0, i = 1, 2, 3, (4.19)
The functions ϕ i (z) represent the initial concentrations of biomass i, the functions S 0j (z) represent the initial substrate concentrations into biofilm, G j (t)represent the values assumed by substrates S j at the biofilm-bulk liquid interface. The free boundary evolution is governed by the following ordinary differential equation:

L (t) = u (L (t) , t) -σ (L (t) , t) , t > 0, (4.20)
with the following initial condition:

L (0) = L 0 , (4.21)
where L 0 denotes the initial biofilm thickness and σ(L(t), t) represents the velocity at which biomass is exchanged between biofilm and bulk liquid [START_REF] Wanner | A multispecies biofilm model[END_REF], the expression used in this work is:

σ(L(t), t) = λL 2 (t) (4.22)
The qualitative analysis of system (4.2, 4.3, 4.13) developed in [START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF] was based on the characteristic method. The numerical method proposed in [START_REF] B D'acunto | Dynamic modeling of sulfate reducing biofilms[END_REF] has been applied. The procedure can be briefly summarized as follows: from the initialboundary conditions u is computed; then L, X, S, are computed in this order; next the computational process is repeated and the solution at the final time is obtained. Numeric integration of the system (4.2, 4.3, 4.13, 4.17, 4.18, 4.19, 4.20, 4.21) has been performed using original software. The schematic representation of the microbial process is reported in Table 4.1.

Results and discussion

The mathematical model proposed in this paper has been applied to simulate the sulfate reduction process in a multispecies biofilm with an initial thickness of 300 μm. The initial conditions and biological parameters, used in the model, are reported in Table 4 The model assumes appropriate boundary conditions for the biological process modelled. An initial arbitrary biomass distribution has been adopted in order to evaluate the dynamical response of the model as the system tends to reach an equilibrium that is not affected by the initial conditions. The microbial equilibrium is only governed by the boundary conditions, as for each dynamic model. Kinetics and stoichiometric parameters, and diffusion coefficients reported in Table 4.3 have been adopted. 4 ratios in the range 0.5 -2 have been investigated. The simulations have been performed to evaluate the dynamical response of the biofilm in terms of volume fractions of bacteria and concentration trends of substrates. In particular the results show the model capability to reveal the microbial stratification in the biofilm, evaluating the effect of substrate diffusion on biomass growth.

d -1 ] 1 Sulfate reduction by X 1 -1.5 (1-Y1,2) Y1,2 -1 Y1,2 1 μ 1 X 1 2 Sulfate reduction by X 2 -1.5 (1-Y2,2) Y2,2 -1 Y2,2 0.8 (1-Y2,2) Y2,2 1 μ 2 X 2 3 Sulfate reduction by X 3 -1 Y3,2 0.8 (1-Y3,2) Y3,2 1
μ 3 X 3 4 Sulfate reduction by X 4 -1 Y4,3 1 μ 4 X 4 5 Decay of X 1 -1 1 Kd 1 X 1 6 Decay of X 2 -1 1 Kd 2 X 2 7
Decay of X LDSRB (C) and LDSRB (I) prevail in the outer layer of biofilm where sulfate and lactate remain abundant. In the deepest zone of the biofilm, characterized by a low level of sulfate and lactate, due to substrate diffusion coupled with microbial consumption, the MA compete for space with other microbial species. Indeed the MA are mostly present in the deepest zone of the biofilm where the optimal conditions for their growth are established. With a COD/SO 2- 4 ratio of 0.5 (Figure 4.2) the AB are present at the inner layer of the biofilm, where the concentration of sulfate is lower while both LDSRB (C) and LDSRB (I) are found to be predominant at the outmost layer of the biofilm. In particular, LDSRB (C) and LDSRB (I) represent the most abundant species in the biofilm showing that in the presence of excess sulfate, the quantitative oxidation of lactate to acetate or CO 2 coupled to sulfate reduction is the dominant reaction.

A similar result has been achieved in [START_REF] Lb Celis-García | Performance of a downflow fluidized bed reactor under sulfate reduction conditions using volatile fatty acids as electron donors[END_REF], where the authors experienced high participation of SRB, and in particular of LDSRB (I) , on COD removal in a down-flow fluidized-bed reactor. Concerning sulfate reduction, it is possible to note that sulfate is not completely depleted within biofilm, probably due to the presence of incomplete oxidizers, while lactate-COD concentration drops to zero in 200 μm (Figure 4.2 (A)). As shown in Figure 4.2 (B), methanogenesis and acetogenesis are not completely suppressed; however the volume fraction of AB is sensitively reduced respect to the initial condition and this trend is expected to exacerbate with time. On the other hand, the formation of a zone in the inner part of biofilm characterized by abundance of acetate and lack of lactate-COD could support the methanogenic metabolism allowing the methanogens to remain present in the biofilm. According to the experience of [START_REF] Papirio | Biological inverse fluidized-bed reactors for the treatment of low ph-and sulphate-containing wastewaters under different cod conditions[END_REF], acetate production can be recognized as the rate-limiting step in such a sulfate-reducing process.

In Figure 4.3 and 4.4 is shown the response of the multispecies biofilm to the increasing COD/SO 2- 4 ratio. As experienced in [START_REF] Papirio | Biological inverse fluidized-bed reactors for the treatment of low ph-and sulphate-containing wastewaters under different cod conditions[END_REF], the excess of lactate over sulfate continuously guaranteed the required carbon source for SRB to reduce sulfate to sulfide. The exposure to higher COD/SO 2- 4 ratios was enough for the development of substantial sulfidogenesis leading to sulfate depletion (Figures 4.3 (A), 4.4 (A)). In this condition, acetogens do not experience competition for the remaining COD; therefore the area of acetogenic within the biofilm becomes broader at increasing COD/SO

Conclusions

A mathematical model able to simulate the physical, chemical and biological processes prevailing in a multispecies sulfate reducing biofilm under dynamic conditions has been presented. Special attention has been given to the competition between sulfate reduction, acetogenesis and methanogenesis. The effects of the variations of the operational conditions in terms of COD/SO

Modelling microbial population dynamics in multispecies biofilms including Anammox bacteria A 1-D mathematical model for analysis and prediction of microbial interactions within multispecies biofilms including Anammox pathway is presented. The model combines the related processes of organic carbon oxidation, denitrification, nitrification and Anammox and phenomena of substrate reaction and diffusion, biomass growth and advection, detachment. The biofilm growth process is governed by nonlinear hyperbolic PDEs and substrate dynamics are dominated by semilinear parabolic PDEs. It follows a complex system of PDEs on a free boundary domain. Equations are integrated numerically by using the method of characteristics as strongly suggested by the qualitative analysis of the free boundary value problem. Mass conservation equation plays an important role in checking the accuracy of simulations. The model has been applied to simulate Anammox competition and to evaluate the influence of substrate diffusion on microbial stratification. Specific scenarios are analyzed. The results reveal that in a thick multispecies biofilm, including heterotrophic, aerobic autotrophic nitrifying and Anammox bacteria, oxygen diffusion limitation determines the formation of both aerobic and anoxic microenvironments favouring interspecies competition. In contrast, oxygen excess causes a disturbance on microbial interactions leading to Anammox bacteria loss. The model predictions may help engineers or operators to have a better insight into biofilm dynamics in order to optimize process design or practical operation.

A modified version of this chapter was published as: Mattei, M.R., Frunzo, L., D'Acunto, B., Esposito, G. and Pirozzi, F. (2015). Modelling microbial population dynamics in multispecies biofilms including Anammox bacteria. Ecological Modelling, 2015, vol. 304, pp. 44-58, DOI: 10.1016/j.ecolmodel.2015.02.007.

Introduction

Biofilm research had been neglected for a long time until microbiologists rediscovered these fascinating communities almost 40 years ago [START_REF] Tj Battin | Microbial landscapes: new paths to biofilm research[END_REF]. In the past decades the number of studies performed on surface-associated microbes has increased considerably and today, we recognize that most, if not all, microbial species can form biofilms [START_REF] Kolter | Microbial sciences: the superficial life of microbes[END_REF][START_REF] Deygout | Homogenised model linking microscopic and macroscopic dynamics of a biofilm: Application to growth in a plug flow reactor[END_REF]. These sessile communities form anywhere there is a surface with a little moisture and some nutrients and nearly always harbour a multitude of microbial species, which compete or coexist thanks to niche differentiation.

The elaborate biofilm structure constitutes the ideal environment for the development of different microbial groups which may be confronted with dynamic changes in nutrient profile, either due to environmental changes, or due to the metabolism and migration of other populations [START_REF] Moons | Bacterial interactions in biofilms[END_REF]. In multispecies biofilms, different species may be separated into discrete layers, according to their metabolic activities, affinity for substrates, growth rates and sensitivity towards inhibiting substances [START_REF] Cd Nadell | The sociobiology of biofilms[END_REF]. A striking example is provided by the stratification which manifests itself between aerobic and anaerobic species in the oxygenated and anoxic regions of biofilms of wastewater treatment reactors. In these systems, the relative abundance of different microbial populations strongly depends on the metabolic potential of the inhabiting strains which contribute themselves to local nutrient composition through metabolism and growth. As a consequence, while bacteria maximize their growth on the available nutrients, they might change the ratio of components therein, thereby creating conditions that favor the development of other microbial groups [START_REF] Moons | Bacterial interactions in biofilms[END_REF].

The inherent synergical interactions characterizing natural biofilm communities have shown to facilitate the simultaneous removal of various pollutants in wastewater treatment reactors. In several biofilm-based processes, redox stratification is experienced due to the formation of strong concentration gradients of both electron donors and acceptors, and the accumulation of metabolic waste products, that can be used as growth substrates by other microorganisms. This has been highlighted, for instance, in the case of the ANaerobic AMMonia OXidation (Anammox) process, which has become one of the most promising innovative techniques for the biological removal of nitrogen from wastewater. In the absence of molecular oxygen, Anammox bacteria catalyze this novel process, where ammonia is anaerobically oxidized to nitrogen gas, with nitrite as the electron acceptor. In their environments, Anammox bacteria receive the key substrates (ammonium and nitrite) in cooperation and competition with other N-cycle microorganisms. Environmental studies have provided strong evidence for the close cooperation between aerobic ammonia-oxidizing bacteria, which inhabit the aerobic regions of biofilms, and Anammox bacteria, the first producing one of the substrates (nitrite) for the second [START_REF] Terada | Redox-stratification controlled biofilm (rescobi) for completely autotrophic nitrogen removal: The effect of co-versus counter-diffusion on reactor performance[END_REF]. This cooperative action has laid the foundation for a new wastewater treatment technology, whose performance is closely linked to complex and delicate metabolic interactions within biofilms. Careful microenvironment control is required due to the high sensitivity of Anammox bacteria to oxygen and nitrite, and the slow specific growth rate of this microbial strain. Consequently, evaluating the competition in population dynamics represents a crucial step in optimizing the control criteria of Anammox process.

Mathematical models represent useful tools to explore the microbial competition and coexistence in multispecies biofilms. Many biofilm models including Anammox metabolism have been presented over the last decades. In 2002, Hao et al. [START_REF] Hao | Sensitivity analysis of a biofilm model describing a one-stage completely autotrophic nitrogen removal (canon) process[END_REF] introduced the first biofilm model aimed at evaluating the dynamics of a multispecies biofilm including Anammox metabolism and performing a completely autotrophic nitrogen removal. In this model, the dynamics of heterotrophs were not included. This gap has been filled in by the same authors in 2004 [START_REF] Hao | Model-based evaluation of cod influence on a partial nitrification-anammox biofilm (canon) process[END_REF], by investigating the effect of heterotrophic growth on influent organic substrate on the performance of a partial nitritration-anammox biofilm reactor. An extension of the previous model was presented by Lackner et al. [START_REF] Lackner | Heterotrophic activity compromises autotrophic nitrogen removal in membrane aerated biofilms: results of a modelling study[END_REF], who evaluated the growth of heterotrophic bacteria on influent organic carbon and microbial decay products in both co-diffusion and counter diffusion systems. Recently, Mozumder et al. [START_REF] Msi Mozumder | Effect of heterotrophic growth on autotrophic nitrogen removal in a granular sludge reactor[END_REF] have proposed a different version of the model in which they evaluated the influence of heterotrophic growth on autotrophic nitrogen removal in a granular sludge reactor.

In this contribution the same approach introduced in [57] was followed, and a model able to study the Anammox competition in a multispecies biofilm performing nitrogen removal through partial-nitritation Anammox and COD degradation has been developed. The effect of heterotrophic growth on the performance of a nitrifying biofilm has been investigated in previous works. In most of the cases, heterotrophic denitrification has been modelled via both nitrite and nitrate.

However, as shown by experimental results heterotrophic denitrification is a sequential process [START_REF] Khin | Novel microbial nitrogen removal processes[END_REF]. Indeed only in a recent contribution [START_REF] Msi Mozumder | Effect of heterotrophic growth on autotrophic nitrogen removal in a granular sludge reactor[END_REF], a sequential denitrification mechanism has been adopted. The production of nitrite due to the reduction of nitrate by denitrifying microorganisms has been found to play a crucial role in the development of Anammox bacteria. Therefore, we applied the kinetics of denitrification proposed by Kaelin et al. [START_REF] Kaelin | Extension of asm3 for twot-step nitrification and denitrification and its calibration and validation with bacth tests and pilot scale data[END_REF] in the case of planktonic growth, to a multispecies biofilm model. For biomass decay, the concept of endogeneous respiration has been adopted instead of the death-regeneration model [START_REF] Hao | Model-based evaluation of cod influence on a partial nitrification-anammox biofilm (canon) process[END_REF].

The performed experimental results on Anammox biofilm reactors have shown that the growth of this microbial species is strongly affected by the oxygen concentration within biofilms resulting from the preceeding interaction between microbial metabolisms and diffusion. A contributory effect is provided by the shear stress which erodes the biofilm surfaces exposing the Anammox bacteria to higher oxygen concentrations internally. In this work, we focused on the dynamic model behavior of the system under different bulk oxygen concentrations and shear stress conditions.

The proposed model is based on the widely used 1D dynamic multi-species biofilm model introduced by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] which is suited best to model situations in which competition between the microbial constituents of the biofilm matrix is significant. In [START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF], the authors presented an analysis of solutions to a free boundary value problem related to the multispecies biofilm model introduced in [START_REF] Wanner | A multispecies biofilm model[END_REF]. The system of partial differential equations characterizing the model was discussed by using the method of characteristics. The latter represents an efficient resolution method for finding analytical and numerical solutions to hyperbolic partial differential equations. If this method is applied to the free boundary value problem, arising from biofilm growth, it does not require a coordinate transformation to fix the size of the domain in the 1D set-up, as proposed by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF]. Consequently, if the initial phase of biofilm development, where the biofilm initial length is zero, is considered, no coordinate transformation is initially possible, but the method of characteristics can still be applied. Different modelling scenarios have been experimented in order to evaluate microbial population shifts upon changes in operating conditions. Changes in the bacterial populations, substrate profiles within biofilm, and hence in, the nitrogen removal characteristics, have been monitored over time and at different DO levels and shear stress conditions. The simulations show the ability of such a model to examine such scenarios and their dynamic effects on biofilm behavior and performance.

Biological problem

Nutrient removal has become a big concern for wastewater treatment since the end of the 20 th century due to the discharge of nitrogen and phosphorous compounds-laden wastewater to natural water bodies, that results in many cases in eutrophication, emissions of nitrous oxide to atmosphere during oxidation of ammonia and toxicity to aquatic invertebrate and vertebrate species [START_REF] Okabe | Spatial microbial distributions of nitrifiers and heterotrophs in mixed-population biofilms[END_REF][START_REF] Philips | Origin, causes and effects of increased nitrite concentrations in aquatic environments[END_REF]. The eutrophication may be managed and controlled primarily by restricting the nutrient inputs to natural bodies, by means of appropriate wastewater treatment plants (WWTPs) and this nutrient loading restriction can be accomplished by a wide variety of external and internal controls [START_REF] Vh Smith | Eutrophication: impacts of excess nutrient inputs on freshwater, marine, and terrestrial ecosystems[END_REF].

Nutrient removal efficiency in WWTPs depends on several factors, including treatment technologies used, influent wastewater characteristics, mechanical and operational failures and facility design limitations. It represents a special challenge for WWTPs due to the additional costs associated with the complex treatment technology and technical designs required to produce effluent containing low nutrient concentrations [START_REF] Ro | Contribution of wastewater treatment plant effluents to nutrient dynamics in aquatic systems: A review[END_REF]. Consequently, new approaches and techniques have been studied and tested in order to meet the increasingly stringent discharge standards and to achieve a more sustainable and cost efficient nutrient removal for a variety of so-called hard-to-treat concentrated nutrient streams, including sludge digester supernatant, manure, piggery wastewaters and several industrial wastewater [START_REF] Garuti | Full-scale ananox (r) system performance[END_REF][START_REF] Wyffels | Modeling and simulation of oxygen-limited partial nitritation in a membrane-assisted bioreactor (mbr)[END_REF]. Ammonium nitrogen is conventionally removed from wastewater by the sequential biochemical processes of aerobic autotrophic nitrification, using molecular oxygen as electron acceptor (conversion of NH + 4 to NO - 2 and further to NO - 3 ), and anoxic heterotrophic denitrification (conversion of NO - 2 and NO - 3 to gaseous nitrogen) using organic carbon as electron donor [START_REF] Khin | Novel microbial nitrogen removal processes[END_REF][START_REF] Hjuler | A model of the simultaneous nitrogen and carbon removal in a biofilm[END_REF]. In general, the combined system nitrification/denitrification requires considerable amounts of resources (4.57 kg of O 2 and 2-4 kg COD per kg of ammonium nitrogen), results in a high production of sludge (1 kgV SS/kgN ), depends on external addition of carbon source to achieve complete denitrification and leads to high costs of construction, operation and maintenance, since the two steps are usually accomplished in separate oxic and anoxic units [START_REF] Paredes | New aspects of microbial nitrogen transformations in the context of wastewater treatment -a review[END_REF][START_REF] Ao Sliekers | Completely autotrophic nitrogen removal over nitrite in one single reactor[END_REF].

An innovative and more sustainable biological nitrogen-removal technology, alternative to the traditional nitrification/denitrification system, the Anammox process, has been recently discovered in a denitrifying fluidized bed reactor [START_REF] Mulder | Anaerobic ammonium oxidation discovered in a denitrifying fluidized bed reactor[END_REF] and ever since then, Anammox has been extensively researched as a promising method able to overcome the shortcomings of conventional treatments [START_REF] Msm Jetten | -2004: 10 years of research on the anaerobic oxidation of ammonium[END_REF]. In this process ammonium is oxidized with nitrite serving as the electron acceptor under anaerobic conditions, producing nitrogen gas and nitrate [START_REF] Msm Jetten | The anaerobic oxidation of ammonium[END_REF]. No addition of external carbon is required due to the autotrophic nature of these bacteria, negligible sludge is produced thanks to the low biomass yield, oxygen requirements are reduced two-fold since only half of ammonium is oxidized to nitrite and low energy is required as compared to the conventional nitrification/denitrification process [START_REF] Rc Jin | The importance of the substrate ration in the operation of the anammox process in upflow biofilter[END_REF][START_REF] Strous | Effects of aerobic and microaerobic conditions on anaerobic ammonium-oxidizing (anammox) sludge[END_REF].

The growth of Anammox bacteria is strongly influenced by the interaction with other microbial strains since they need a nearby nitrite source. Indeed, nitrite plays a crucial dual role in Anammox reaction: it acts as electron acceptor for the ammonium oxidation and as electron donor for the CO 2 reduction to biomass. Therefore, the application of the Anammox process requires a combination of aerobic and anaerobic conditions that might be performed in a system consisting of either two reactors or a single reactor. In the first case the process takes places in two reactors in series: a partial nitrification reactor where 50% of the wastewater ammonium content is oxidized to nitrite, and a separate unit for the anaerobic oxidation of ammonia. In this way the two biological reactors can be controlled separately, saving 50% on oxygen requirement, 100% on organic carbon source requirement for denitrification, and producing less sludge compared to the conventional nitrification-denitrification process [START_REF] Swh Van Hulle | Construction start-up and operation of a continuously aerated laboratory-scale sharon reactor in view of coupling with an anammox reactor[END_REF]. In the second case both processes occur in the same reactor, like sequencing batch reactors or biofilm reactors, and the oxygen concentration becomes a key control for these types of applications, usually named oxygen-limited nitrogen removal processes. Two processes have been identified to operate in oxygen-limited conditions in one single reactor and in absence of organic carbon source: i) completely autotrophic nitrogen removal over nitrite (CANON) [START_REF] Ao Sliekers | Completely autotrophic nitrogen removal over nitrite in one single reactor[END_REF], which has been extensively evaluated by simulations [START_REF] Hao | Sensitivity analysis of a biofilm model describing a one-stage completely autotrophic nitrogen removal (canon) process[END_REF][START_REF] Hao | Model-based evaluation of temperature and inflow variations on a partial nitrification-anammox biofilm process[END_REF], and ii) oxygen limited autotrophic nitrification denitrification (OLAND) [START_REF] Kuai | Ammonium removal by oxygen-limited autotrophic nitrification-denitrification system[END_REF]. CANON and OLAND processes show only minimal difference in stoichiometry, suggesting that they are based on the same removal pathway. In biofilm-based systems, many reactions are combined in the same reactor and complex physical interactions exist between the microorganisms living in the biofilm. All these interactions are linked together by nitrogen conversions and diffusion gradients. In general, a biofilm-based reactor operates in continuous conditions since the oxygen levels are governed by gradients in biofilm systems; the classical nitrification, performed by the ammonium oxidizers, takes place in the outer aerobic layers while the anaerobic oxidation occurs in the deeper zones of the biofilm. Anoxic conditions required for anammox metabolism are established by the oxygen respiration operated by aerobic ammonium oxidizers. Anammox bacteria, in turn, remove the toxic nitrite and convert the remaining ammonium into nitrogen gas. These biofilm systems are characterized by a long biomass retention time which perfectly fits with the long doubling time of Anammox bacteria and allows the formation of different substrate concentration gradients in the same reactor [START_REF] Gong | Feasibility of a membrane-aerated biofilm reactor to achieve single-stage autotrophic nitrogen removal based on anammox[END_REF][START_REF] Ka Third | The canon system (completely autotrophic nitrogen-removal over nitrite) under ammonium limitation: interaction and competition between three groups of bacteria[END_REF]. The application and industrialization of the Anammox process have been restricted by the slow growth rate of the Anammox bacteria and the widespread inhibition factors existing in nitrogen-rich wastewater. Indeed, this type of bacteria show high sensitivity to changing environmental conditions and to the composition of wastewater making the process more difficult to initiate and recover from inhibition. Moreover Anammox activity is based on the harmonious and balanced interaction with other bacteria which can be disturbed and interfere with nitrogen removal. The behaviour of such a microbial community is complex and depends on multiple parameters. It is evident that substrate diffusion in the biofilm plays a crucial role in defining the composition, diversity and dynamics of such biofilm bacterial communities since it facilitates a variety of microhabitats. In particular, Dissolved Oxygen (DO) has been recognized as a critical operational parameter for the Anammox process as it can strongly affect the coexistence and co-performance of microbial populations with complementary and/or opposed environmental requirements. Anammox inhibitions can be controlled by proper measures. Substrate concentration and loading rate control, pH adjustment, sludge acclimatization, DO and ORP control are all effective in preventing and relieving Anammox inhibition [START_REF] Rc Jin | The inhibition of the anammox process: A review[END_REF]. However evaluating Anammox process through experiments would take a long time due to bacteria slow growing rates. Thus, the use of mathematical modelling is crucially helpful in testing a large variation of environmental and operational conditions that might influence the process [START_REF] Hao | Model-based evaluation of cod influence on a partial nitrification-anammox biofilm (canon) process[END_REF].

Model construction and numerical approach

The proposed model simulates the dynamics of the biological ammonium removal in a multispecies biofilm. Chemical, physical and biological transient processes are analyzed. In particular, the model takes into account the oxygen-limited nitrogen removal process coupled with soluble organic carbon removal by aerobic and denitrifying microrganisms, as shown in Figure 5.1. The competition between heterotrophs and autotrophs in completely autotrophic ammonium removal processes has been usually neglected [START_REF] Hao | Sensitivity analysis of a biofilm model describing a one-stage completely autotrophic nitrogen removal (canon) process[END_REF][START_REF] Hao | Model-based evaluation of temperature and inflow variations on a partial nitrification-anammox biofilm process[END_REF]; however recent studies have shown that in autotrophic biofilms up to 50% of the biomass can be heterotrophic [START_REF] Okabe | Spatial microbial distributions of nitrifiers and heterotrophs in mixed-population biofilms[END_REF] and Anammox bacteria can coexist with denitrifying bacteria [START_REF] Hao | Model-based evaluation of cod influence on a partial nitrification-anammox biofilm (canon) process[END_REF][START_REF] Lackner | Heterotrophic activity compromises autotrophic nitrogen removal in membrane aerated biofilms: results of a modelling study[END_REF][START_REF] Hao | Model-based evaluation of oxygen consumption in a partial nitrification-anammox biofilm process[END_REF]. The inhibition of Anammox bacteria by organic matter can be defined concentration dependent: at high concentrations of organic carbon heterotrophic bacteria are able to outcompete Anammox bacteria due to their faster growth; at low concentrations of organic matter heterotrophs can not dominate and outcompete Anammox bacteria and different biological reactions are promoted. However, the microbial interactions in such a multispecies biofilm are very complex and the coexistence of Anammox bacteria with other processes is related to many parameters. In this work, the influence of organic carbon and nitrogen loading rates and DO concentration on biofilm dynamics have been considered.

Figure 5.1: Main microbial interactions of the simulated biological process

The proposed model takes into account the activities of heterotrophs and autotrophs living in a biofilm, evaluating the interactions between the related processes: organic carbon oxidation, denitrification, nitrification and Anammox. Five reacting components are simultaneously considered: ammonium (S 1 ), nitrite (S 2 ), nitrate (S 3 ), soluble organic carbon (S 4 ) and oxygen (S 5 ). Organic carbon concentration is expressed in terms of Chemical Oxygen Demand (COD) which represents a measurement conventionally used in environmental chemistry to characterize indirectly the amount of organic compounds [START_REF] Tchobanoglous | Wastewater Engineering: Treatment and Reuse[END_REF].

The oxygen-limited nitrogen removal process involves the growth of three autotrophic organisms: aerobic ammonia-oxidizing bacteria AeAOB (X 1 ), aerobic nitrite-oxidizing bacteria NOB (X 2 ) and anaerobic ammonia-oxidizing bacteria AnAOB (X 3 ), taking also into account the related decay and endogenous respiration processes.

COD is removed by both aerobic oxidation and denitrification, performed by facultative heterotrophic bacteria HB (X 4 ). Inert residues Inerts (X 5 ) are also taken into account. According to [START_REF] Kaelin | Extension of asm3 for twot-step nitrification and denitrification and its calibration and validation with bacth tests and pilot scale data[END_REF], anoxic denitrification involves two reactions including a first step from nitrate to nitrite and second from nitrite to molecular nitrogen. Nitrite production as an intermediate of the denitrification process becomes of interest since Anammox competes for nitrite with denitrification. Only anoxic endogenous respiration over nitrate is considered for the autotrophic microorganisms, while all anoxic processes for heterotrophic microorganisms are doubled in a first step from nitrate to nitrite and second from nitrite to molecular nitrogen. The biofilm growth is governed by the following equations [START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF]:

∂X i ∂t + ∂ ∂z (uX i ) = ρ i r M,i (z, t, X, S), 0 ≤ z ≤ L(t), t > 0, i = 1, 2 , 3, 4, 5, (5.1)
∂u ∂z = 5 i=1 r M,i (z, t, X, S), 0 < z ≤ L(t), t > 0, (5.2)
where X i = ρ i f i (z, t) denotes the concentration of the microbial species and inert residues i = 1, 2, 3, 4, 5, f i is the volume fraction of microbial species i

5 i=1 f i = 1, L(t)
is the thickness of biofilm, ρ i the biofilm constant density, S j (z, t) is the concentration of substrate j = 1, 2, 3, 4, 5, u(z, t) is the velocity of the microbial mass displacement with respect to the biofilm support interface, r M,i (z, t, X, S) represent the biomass growth rates, X = (X 1 , X 2 , X 3 , X 4 , X 5 ) and S = (S 1 , S 2 , S 3 , S 4 , S 5 ).

The biomass growth rates are given by:

r M,1 = (μ 1 -b 1 -c 1 )X 1 , (5.3)
r M,2 = (μ 2 -b 2 -c 2 )X 2 , (5.4)
r M,3 = (μ 3 -b 3 -c 3 )X 3 , (5.5)
r M,4 = (μ 4,1 + μ 4,2 + μ 4,3 -b 4 -c 4,1 -c 4,2 )X 4 , (5.6)
while for inert residues

r M,5 = f I (b 1 + c 1 )X 1 + f I (b 2 + c 2 )X 2 + f I (b 3 + c 3 )X 3 + f I (b 4 + c 4,1 + c 4,2 )X 4 , (5.7
) where μ 1 , μ 2 , and μ 3 are the net biomass growth rates for biomass X 1 , X 2 and X 3 ; μ 4,1 , μ 4,2 and μ 4,3 are the net biomass growth rates for biomass X 4 in aerobic (μ 4,1 ) and anoxic conditions over nitrate (μ 4,2 ) and nitrite (μ 4,3 ); b 1 , b 2 , b 3 and b 4 are the Aerobic Endogenous respiration rates for the single microbial species; c 1 , c 2 and c 3 are the Anoxic Endogenous respiration rates over nitrate for the autotrophic microorganisms; c 4,1 and c 4,2 are the Anoxic Endogenous respiration rates for X 4 , over nitrate and nitrite respectively. They are given by:

μ 1 = μ max,1 S 1 K 1,1 + S 1 S 5 K 1,5 + S 5 , (5.8)
μ 2 = μ max,2 S 2 K 2,2 + S 2 S 5 K 2,5 + S 5 , (5.9)
μ 3 = μ max,3 K 3,5 K 3,5 + S 5 S 1 K 3,1 + S 1 S 2 K 3,2 + S 2 , (5.10)
μ 4,1 = μ max,4 S 4 K 4,4 + S 4 S 5 K 4,5 + S 5 , (5.11)
μ 4,2 = β 1 • μ max,4 K 4,5 K 4,5 + S 5 S 4 K 4,4 + S 4 S 3 K 4,3 + S 3 , (5.12
)

μ 4,3 = β 2 • μ max,4 K 4,5 K 4,5 + S 5 S 4 K 4,4 + S 4 S 2 K 4,2 + S 2 , (5.13)
b 1 = b m,1 S 5 K 1,5 + S 5 , (5.14)
b 2 = b m,2 S 5 K 2,5 + S 5 , (5.15)
b 3 = b m,3 S 5 K 3,5 + S 5 , (5.16)
b 4 = b m,4 S 5 K 4,5 + S 5 , (5.17)
c 1 = η • b m,1 K 1,5 K 1,5 + S 5 S 3 K 4,3 + S 3 , (5.18)
c 2 = η • b m,2 K 2,5 K 2,5 + S 5 S 3 K 4,3 + S 3 , (5.19)
c 3 = η • b m,3 K 3,5 K 3,5 + S 5 S 3 K 4,3 + S 3 , (5.20)
c 4,1 = η 1 • b m,4 K 4,5 K 4,5 + S 5 S 3 K 4,3 + S 3 , (5.21)
c 4,2 = η 2 • b m,4 K 4,5 K 4,5 + S 5 S 2 K 4,2 + S 2 , (5.22)
where μ max,i denotes the maximum net growth rate for biomass i, K i,j the affinity constant of substrate j for biomass i, β 1 and β 2 the reduction factor for denitrification NO 3 -NO 2 and NO 2 -N 2 respectively, b m,i the decay-inactivation rate for biomass i, η the anoxic reduction factor for b m,i for autotrophic microorganisms, η 1 and η 2 the reduction factors for b m,4 in anoxic conditions NO 3 -NO 2 and NO 2 -N 2 respectively, f I the inert content in lysis of biomass. The diffusion of substrates is governed by the equations:

∂S j ∂t -D j ∂ 2 S j ∂z 2 = r S,j (z, t, X, S), 0 < z < L(t), 0 < t ≤ T, j = 1, 2 , 3, 4, 5 (5.23
) where D j denotes the diffusivity coefficient and r S,j (z, t, X, S) the net conversion rate of substrate j.

These are expressed by:

r S,1 = (- 1 Y 1 -i N,B )μ 1 X 1 +(- 1 Y 3 -i N,B )μ 3 X 3 -i N,B (μ 2 X 2 +μ 4,1 X 4 +μ 4,2 X 4 +μ 4,3 X 4 ) +(i N,B -i N,I f I )[(b 1 + c 1 )X 1 +(b 2 + c 2 )X 2 +(b 3 + c 3 )X 3 +(b 4 + c 4,1 + c 4,2 )X 4 ], (5.24)
r S,2 = 1 Y 1 μ 1 X 1 - 1 Y 2 μ 2 X 2 + (- 1 Y 3 - 1 1.14 )μ 3 X 3 -(1 - 1 Y 4 ) 1 1.14 μ 4,2 X 4 + (1 - 1 Y 4 ) 1 1.72 μ 4,3 X 4 + 1 -f I 1.14 c 4,1 X 4 - 1 -f I 1.72 c 4,2 X 4 , (5.25)
r S,3 = ( 1 1.14 )μ 3 X 3 + (1 - 1 Y 4 ) 1 1.14 μ 4,2 X 4 + 1 Y 2 μ 2 X 2 - 1 -f I 2.86 (c 1 X 1 + c 2 X 2 + c 3 X 3 ) - 1 -f I 1.14 c 4,1 X 4 , (5.26)
r S,4 = - 1 Y 4 (μ 4,1 X 4 + μ 4,2 X 4 + μ 4,3 X 4 ) (5.27) r S,5 = (1 - 3.43 Y 1 )μ 1 X 1 - 1 Y 4 μ 4,1 X 4 + (1 - 1.14 Y 2 )μ 2 X 2 -(1 -f I )(b 1 X 1 + b 2 X 2 + b 3 X 3 + b 4 X 4 ) (5.28)
where Y i denotes the yield for biomass i, f I denotes the inert content in lysis biomass i, i N,B is Nitrogen content in biomass, i N,I is Nitrogen content in inert biomass.

The following initial-boundary conditions will be considered for equations (5.1), (5.2) and (5.21)

X i (z, 0) = ϕ i (z), u(0, t) = 0, 0 ≤ z ≤ L 0 , t ≥ 0, i = 1, 2, 3, 4, 5, (5.29)
S j (z, 0) = S j0 (z), 0 ≤ z ≤ L 0 , j = 1, 2, 3, 4, 5, (5.30)

∂S j ∂z (0, t) = 0, S j (L(t), t) = S jL (t), 0 < t ≤ T, j = 5, (5.31)

∂S j ∂z (0, t) = 0, ∂S j ∂z (L(t), t) = G j , 0 < t ≤ T, j = 1, 2, 3, 4 (5.32)
The functions ϕ i (z), i = 1, ..., 5, represent the initial concentrations. Condition (5.29) 2 follows from the relationship g i (0, t) = u(0, t)X i (0, t) of the biomass flux at z = 0. The functions S j0 (z) represent the initial values of substrates. The function S jL (t) in (5.31) 2 is the oxygen value in the bulk liquid. G j in (5.32) 2 represent the substrate fluxes and reproduce the operational conditions of such a biofilm-based system. The free boundary evolution is governed by the following ordinary differential equation:

dL dt (t) = u(L(t), t) -σ(L(t), t), (5.33)
with initial condition:

L(0) = L 0 , (5.34)
where L 0 denotes the initial biofilm thickness and σ(L(t), t) represent the velocity at witch biomass is exchanged between biofilm and bulk liquid and is assumed as a known function of L and t. in this work it takes the expression of σ(L(t), t) = λL 2 (t), (5.35) where λ is the shear constant. As proved in [START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF], by using the method of characteristics the equations (5.1) can be written as:

d dt X i (z(z 0 , t), t) = F i (z, t, X, S), 0 ≤ z 0 ≤ L 0 , t > 0, i = 1, ..., 5 (5.36)
where the characteristic lines z = z(z 0 , t) are defined as

∂z ∂t = u(z(z 0 , t), t), z(z 0 , 0) = z 0 , 0 ≤ z 0 ≤ L 0 , t > 0.
(5.37)

From this definition and from the no-flux condition at the substratum (5.29) 2 , it follows that the characteristic starting at z 0 = 0, is by necessity z = 0. Consequently, since the method of characteristic is able to provide the solution along any characteristic, this method of course gives the solution of the system for z = 0 at any t. So, the resolution of the biological problem does not require any boundary condition for z = 0, since the solution is derived uniquely by the initial condition. Mathematical models including prescribed boundary conditions on f i for z = 0 and z = L(t) could be considered, although it is not required by the biological process and it could lead to cases without solutions. The numerical simulations have been performed by applying the method of characteristics to the system of non-linear equations (5.1), and by using the Euler explicit method for the system of semi-linear equations (5.23) [START_REF] Acunto | Computational Methods for PDE in Mechanics[END_REF]. This method was first introduced in [START_REF] B D'acunto | Dynamic modeling of sulfate reducing biofilms[END_REF], showing great efficiency and accuracy, and could be applied to evaluate the dynamics of biofilms with zero initial thickness, overcoming the lack of the current models [START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF]. The great accuracy of the method is strictly connected to the existence of an invariant in the numerical integration. In particular, let consider equation (5.2) and note that it is equivalent to 5 i=1 f i = 1, which states the conservation of mass for the whole system. The evaluation of each bacterial species volume fraction is performed separately along the characteristics at each time step integration, allowing us to use equation 5 i=1 f i = 1 as a check control and to evaluate the error. The overall model stoichiometry and kinetics for the one-dimensional biofilm model and the respective parameter values used for numerical simulations, derive from published literature and are reported in Table 5.1 and Table 5.2.

For all the simulations reported in this paper, the model assumes an initial biofilm thickness of 300 μm and appropriate boundary conditions for the biological process modeled. Ammonium and COD were supplied from the bulk liquid at defined loading rates; DO at the biofilm-bulk liquid interface was fixed as a model input parameter so the oxygen load was not adjusted to reach optimal nitrogen removal. The influent contained no nitrite or nitrate. Ammonium and nitrate do not inhibit the Anammox process but nitrite concentration exceeding 100 mgN/L can inactivate it. However, in reality, nitrite concentration rarely reaches this level, therefore the nitrite inhibition on the Anammox process was not included in the model. An initial arbitrary biomass distribution has been adopted in order to prove the dynamic response of the model which tends to reach an equilibrium that is not affected by the initial conditions. The microbial equilibrium is only governed by the boundary conditions, as for any dynamic model. The adopted initial biofilm volume fractions are AeAOB 20%, AnAOB 17%, NOB 33%, HB 30%, Inerts 0%. In all simulations, inert is modeled as another microbial species whose growth derives from the heterotrophic and autotrophic biomass decay.

The model capability of covering different ranges of input values of different parameters represents a sort of priority and the model providing more accurate forecast of the system's behavior due to inputs of different concentrations with lowest errors is considered to be the most desirable one [START_REF] Vafajoo | Model-based evaluation of operating parameters on canon process in a membrane-aerated biofilm reactor[END_REF]. Moreover, an accurate model should also be able to simulate the effect of variation of operational parameters, which can strongly affect microbial population dynamics. For this reason, four different simulation scenarios have been considered. The first and second scenarios examine the effect of different DO concentrations on the biofilm dynamics and in particular on microbial interactions and stratification over time. The following boundary conditions for equation (5.23) have been adopted: ammonium surface load of 2 gm -2 d -1 , DO concentration of 3 mg/L (Scenario 1) and 5 mg/L (Scenario 2) in the bulk liquid, COD surface load of 0.2

gm -2 d -1 .
The third scenario (Scenario 3) examines the influence of a DO step-change on biofilm dynamics. The case study is achieved by changing the concentration of DO in the bulk liquid keeping constant ammonium and COD surface loads. In particular, during the course of a simulation the DO in the bulk liquid has been shifted from 3 to 5 mg/L and kept constant in order to study the long-term exposure effects on population dynamics. The effect of fluctuating conditions occurs on a much smaller scale than that of biofilm growth. Thus, these conditions can be observed only as short-term effects. However in a biofilm-based process including the Anammox pathway, it is important to know the limitations of the system and whether or not severe disturbances to the system would lead to the loss of this microbial strain.

The fourth scenario (Scenario 4) evaluates the influence of a shear stress change due to a difference in hydrodynamics conditions, in the composition, diversity and dynamics of biofilm bacterial communities. The approach used is similar to Scenario 3: during the course of a simulation the shear stress constant has been changed from 50 m -1 d -1 to 150 m -1 d -1 and the microbial population dynamics have been analyzed. An overview of the different modeling scenarios considered in this study is reported in Table 5.3

Results

Numerical solutions to the free boundary problem stated in Section 5.3 have been obtained by using the method of characteristics, e.g. [START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF]. The accuracy of the solutions was checked by considering the mass conservation equation n i=1 f i = 1. All the simulations in this section have been performed by using an original software developed in MatLab platform.

The simulations have been performed to evaluate the dynamic response of the system in terms of volume fractions of bacteria and substrate concentration trends within biofilm. In particular the results show the model capability to reveal microbial stratification in biofilm, evaluating the effect of substrate diffusion on biomass growth. DO in the bulk liquid has been postulated as the determining factor for bacteria survival or out-competition. So, different oxygen levels have been examined in order to explain the selection of different type of bacteria over time as well as the population shift deriving from a perturbation in the operational conditions. The sum equal to one of all different biofilm volume fractions in all the simulations validates the quality of numerical integration. In all the Figures, the concentration of ammonium, nitrite and nitrate is expressed as mg N/L.

Scenario 1

The simulation results for the multispecies biofilm performance over time at DO = 3 mg/L are shown in Figures 5. Substrate concentration trends derive from the combination between diffusion and microbial metabolism and are characterized by: high ammonia levels in the entire biofilm since it represents the most abundant substrate in the bulk liquid (highest loading rate) and diffuses rapidly; it is consumed by AeAOB where there is abundance of oxygen and by AnAOB where oxygen concentration drops to zero; penetrated oxygen profile with high levels in the outer part of biofilm but the formation of an anoxic zone starts to be visible; high nitrate and nitrite levels in all the biofilm since they are produced by AeAOB and NOB and diffuse within the biofilm without be consumed by AnAOB and denitrifying HB, whose metabolism is still inhibited by oxygen concentration. In Figure 5.2 (C,D), the biomass stratification appears visible but the biofilm has not reached the equilibrium condition since the substrate concentrations, inside biofilm, are still high and could be used for the growth of new biomass. AeAOB prevail in the outer layer of biofilm where oxygen and ammonium remain abundant. In the deepest zone of biofilm, characterized by a low level of oxygen, the AnAOB compete for space with other microbial species and their volume fraction increases with time. At 100 days (Figure 5.2 (E,F)) the stratification becomes significant and oxygen concentration drops to zero in the bottom of the biofilm. Nitrite and nitrate concentrations are still high; nitrate represents the outcome of the biological process and is consumed only by denitrifying bacteria, whose growth is limited by the low COD concentration. Nitrite is produced in the outer part of the biofilm by AeAOB; it can be consumed by NOB, which are out-competed by AOB and aerobic HB for oxygen, denitrifying HB and AnAOB which experience respectively limitation for COD and ammonium. At 150 days (Figure 5.2 (G,H)) the stratification is more defined than after 100 days but not completely developed. In Figure 5.3 (I,M) it is already evident that in the inner layer of biofilm NOB are outcompeted by AnAOB and Inerts volume fraction is increasing over time. NOB tend to occupy the central zone where oxygen is not zero and where there is abundance of nitrite which are produced by AeAOB. Figure 5.3 (O,P) shows the typical stratification of a mature biofilm with a simulation time of 300 days: the AeAOB occupy the outer layer of biofilm out-competing the other microbial species; the HB are mostly present in the central zone of biofilm where the optimal conditions for their growth are established; in the deepest zone of biofilm the AnAOB out-compete NOB and HB since their supply (nitrite and ammonium) is more abundant. Since most of the influent COD is oxidized at the outmost layer of the biofilm, denitrification by heterotrophs plays only a minor role in this biofilm system.

Scenario 2

Figures 5.4 and 5.5 display variations of biomass distribution and substrate concentration trends within biofilm during 300 days of system operation at DO = 5 mg/L. Similarly to Figure 5.2, after few days biofilm population structure is still influenced by the initial condition (Figure 5.4 (A)), but as showed in Figure 5.4 (B), the biofilm experiences a fully penetrated oxygen profile. This diversity strongly reflects on microbial interactions and in particular on AnAOB metabolism. In fact, the AnAOB volume fraction decreases over time until reaching the completely loss from the reactor (Figure 5.4 (A,C,E,G), Figure 5.5 (I,M,O)). Figure 5.4 (B,D,F,H) and Figure 5.5 (L,N,P) provide the variations of substrate concentrations over time. Ammonium and COD concentrations are low within biofilm in all the simulations. Nitrate concentration keeps higher compared to the results achieved for DO =3 mg/L; this trend is strictly connected with the new biofilm microbial composition: nitrate is produced by NOB which experience a perfect survival condition thanks to the simultaneous presence of excess nitrite, produced by AeAOB, and oxygen. The NOB competition for nitrite with AnAOB and for oxygen with AeAOB is relieved, and the cells can start to produce nitrate. Nitrite concentration keeps lower compared to the previous results as it is oxidized by NOB whose metabolism is characterized by higher rate than AnAOB consumption rate. The biofilm stratification starts to be evident after 50 days and reaches the following expected configuration ( indicates that there is still anaerobic ammonium activity after 250 days of oxygen excess, albeit at a low rate, and the interactions between AeAOB and AnAOB could be restablished as soon as the oxygen limitation had been relieved. The ability of the system to tolerate oxygen excess for up to 250 days without irreversible damage shows that such a biofilm-based system represents a robust process for ammonium removal. However, control of aeration plays a crucial role in the biofilm population dynamics and a good DO control is essential to obtain a good N-removal efficiency in this type of biofilm-based systems.

Scenario 3

Scenario 4

The main objective of this scenario is to study the population dynamics of a biofilm subjected to variable shear stress. The simulation has been run considering as initial condition the results of the 100 days time simulation with DO = 3 mg/L ( The effect of shear stress on biofilm physical and microbial properties may be explained by physiological adaptation mechanisms of the same microbial species in the biofilm and/or change in the microbial composition of biofilms. In particular at low shear (Figures 5.2 and 5.3) biofilm displays a high diversity level suggesting that the biofilm maturation stage has already begun; at high shear stress (Figures 5.8 and 5.9) diversity decreases so that one species can dominate the bacterial community. Thereby, decrease in diversity in response to shear stress alteration may be explained by the fact that high detachment forces lead to the formation of a thinner biofilm characterized by a fully penetrated oxygen profile ( The higher concentration of oxygen inhibits Anammox metabolism and determines the loss of this bacterial group after 300 days (Figure 5.9 (O,P)). At the same time, the oxygen availability leads to the formation of a thicker aerobic zone, ideal for the growth of AeAOB and NOB (Figure 5.8 (A,C,E,G), Figure 5.9 (I,M,O)). The concentration of nitrite keeps lower than Figure 5.2 (B,D,F,H) and Figure 5.3 (L,N,P) while nitrate still represents the main outcome of the overall process.

Discussion

The potential of the Anammox technology to perform autotrophic nitrogen removal has been widely demonstrated by the extensive experimental activity carried out during the last decades. A summary of the experimental studies described in literature can be found in [START_REF] Van Hulle | Engineering aspects and practical application of autotrophic nitrogen removal from nitrogen rich streams[END_REF]. Nevertheless, the application and industrialization of the Anammox process have been restricted by the long start-up periods due to the slow growth rate of the Anammox bacteria and the widespread inhibition factors existing in nitrogen-rich wastewater. Indeed, this type of bacteria show high sensitivity to changing environmental conditions and to the composition of wastewater making the process more difficult to initiate and recover from inhibition [START_REF] Dapena-Mora | Evaluation of activity and inhibition effects on anammox process by batch tests based on the nitrogen gas production[END_REF]. Moreover, Anammox activity is based on the harmonious and balanced interaction with other bacteria which can be disturbed and interfere with nitrogen removal. The behavior of such a microbial community is complex and depends on multiple parameters [START_REF] Van Hulle | Engineering aspects and practical application of autotrophic nitrogen removal from nitrogen rich streams[END_REF].

Therefore, the mathematical model could serve as a support tool to gain essential information in the identification of the key factors affecting the efficiency and stability of the process. The mathematical model proposed in this paper has been applied to simulate the Anammox bacteria interactions in a multispecies biofilm including heterotrophic and autotrophic nitrifying bacteria. The main objective was to provide a better understanding of the ecophysiological interactions establishing between nitrifiers, Anammox bacteria and heterotrophs in a bofilm treating wastewaters characterized by low organic matter content and high nitrogen concentration.

The simulation results have shown that substrate diffusion in the biofilm plays a crucial role in defining the composition, diversity and dynamics of such biofilm bacterial communities since it facilitates the formation of various microenvironments [START_REF] Van Hulle | Engineering aspects and practical application of autotrophic nitrogen removal from nitrogen rich streams[END_REF]. In particular, Dissolved Oxygen (DO) has been effectively recognized as a critical operational parameter for the Anammox process [START_REF] Jung | Factors affecting the activity of anammox bacteria during start up in the continuous culture reactor[END_REF][START_REF] Guo | Long-term effect of dissolved oxygen on partial nitrification performance and microbial community structure[END_REF][START_REF] Okabe | Development of long-term stable partial nitrification and subsequent anammox process[END_REF] as it can strongly affect the coexistence and co-performance of microbial populations with complementary and/or opposed environmental requirements [START_REF] Egli | Microbial composition and structure of a rotating biological contactor biofilm treating ammonium-rich wastewater without organic carbon[END_REF][START_REF] Szatkowska | A one-stage system with partial nitritation and anammox processes in the moving-bed biofilm reactor[END_REF]. More precisely, the growth of Anammox bacteria has been found to be favored by the formation of zones in the biofilm characterized by different substrate concentration levels, as shown in Figures 5.2 and 5.3. This is confirmed by extensive in-situ experimental activity carried out on biofilm based reactors, which clearly indicates a distinct microdistribution of microorganisms within the biofilm [START_REF] Okabe | Development of long-term stable partial nitrification and subsequent anammox process[END_REF][START_REF] Egli | Microbial composition and structure of a rotating biological contactor biofilm treating ammonium-rich wastewater without organic carbon[END_REF][START_REF] Nielsen | Kinetics, diffusional limitation and microscale distribution of chemistry and organisms in a canon reactor[END_REF][START_REF] Vázquez-Padín | Microbial community distribution and activity dynamics of granular biomass in a canon reactor[END_REF][START_REF] Vlaeminck | Aggregate size and architecture determine microbial activity balance for one-stage partial nitritation and anammox[END_REF].

The effects of DO on microbial stratification are different depending on the conditions applied. Under oxygen-limited conditions (Figures 5.2, 5.3), the oxygen is consumed in the outer layer of the biofilm by AeAOB, NOB and heterotrophs and thus does not penetrate the biofilm completely. Therefore, the Anammox process can be performed in the anoxic layers making use of the produced nitrite that diffuses further into the biofilm. Similar results have been found in [START_REF] Nielsen | Kinetics, diffusional limitation and microscale distribution of chemistry and organisms in a canon reactor[END_REF] where microbiological analysis demonstrated that the aerobic ammonia oxidizers are located at the surface layer of the aggregate, while anaerobic ammonia oxidizers occupy most the anoxic interior parts. Competitions between AeAOB, NOB and AnAOB occur for oxygen, ammonium as well as nitrite [START_REF] Pynaert | Characterization of an autotrophic nitrogen-removing biofilm from a highly loaded lab-scale rotating biological contactor[END_REF]. The regions where oxygen and nitrite levels are low, while ammonium remains not limiting, constitute the ideal competitive environment between these microbial species, since under these conditions NOB are not favored as they have to compete for both electron donor (nitrite) and electron acceptor (oxygen). As shown in Figures 5.2 and 5.3, NOB do not compete very effectively with Anammox bacteria for nitrite and are most likely limited by the available oxygen [START_REF] Egli | Microbial composition and structure of a rotating biological contactor biofilm treating ammonium-rich wastewater without organic carbon[END_REF]. This implies that only a part of the produced nitrite is further oxidized to nitrate by NOB. Therefore, both ammonium and nitrite can diffuse to the lower anoxic part where the Anammox process takes place [START_REF] Egli | Microbial composition and structure of a rotating biological contactor biofilm treating ammonium-rich wastewater without organic carbon[END_REF]. A possible explanation for the limited activity of NOB could be related to the difference in oxygen half saturation constant for AeAOB and NOB [START_REF] Hanaki | Nitrification at low levels of dissolved oxygen with and without organic loading in a suspended-growth reactor[END_REF] or to the inhibition that hydroxylamine severely exerts on NOB [START_REF] Yang | Investigation of batchwise nitrite build-up by an enriched nitrification culture[END_REF]. In the deepest anoxic zones of the biofilm, AnAOB compete with denitrifying microorganisms and depend on AeAOB for nitrite, if ammonium is not limiting for the process. Heterotrophic bacteria experience a double limitation: they compete for oxygen with AeAOB and NOB in the aerobic layer of the biofilm, and for nitrite in the anoxic zone. As shown in Figures 5.2 and 5.3, heterotrophic activity is negatively influenced by the assumed low influent organic matter concentration: heterotrophs are not able to dominate in this system and outcompete Anammox organisms [START_REF] Ruscalleda | Heterotrophic denitrification on granular anammox sbr treating urban landfill leachate[END_REF]. However, when a higher concentration of organic carbon is present in combination with ammonium and nitrite, we expect a faster growth of heterotrophic denitrifiers which would eliminate Anammox bacteria [START_REF] Sabumon | Anaerobic ammonia removal in presence of organic matter: a novel route[END_REF][START_REF] Tang | Suppression of anaerobic ammonium oxidizers under high organic content in high-rate anammox uasb reactor[END_REF]. The threshold concentration for organic carbon in which denitrifiers outcompete Anammox bacteria differs from report to report [START_REF] Van Hulle | Engineering aspects and practical application of autotrophic nitrogen removal from nitrogen rich streams[END_REF]. Under excess of oxygen, the AeAOB and heterotrophic bacteria are not able to consume all the available oxygen, which can inhibit Anammox bacteria determining temporary nitrite accumulation. The simultaneous presence of excess nitrite and oxygen favors the growth of NOB, and the co-performance of AnAOB and AeAOB is drastically reduced (Figures 5.4, 5.5). Vàsquez-Padìn et al. [START_REF] Vázquez-Padín | Microbial community distribution and activity dynamics of granular biomass in a canon reactor[END_REF] used microsensors to measure the oxygen and nitrite concentrations inside the biofilm at different DO levels. For nitrite, they measured a clearly visible nitrite peak in the nitrification zone, corresponding to the external layers of the biofilm, and a reduced concentration in the Anammox zone where nitrite were consumed together with ammonium under anoxic conditions. For oxygen, they observed in all the cases a curve profile due to the combination of the internal diffusion with the biological reaction. Similar results are visible in Figures 5.4 -5.9. Under oxygen-limited conditions, they found out that the high oxygen demand in the surface layer, due to the activity of AeAOB, allowed the formation of an anoxic zone where the Anammox process could take place. They also observed an increased oxygen penetration depth with a DO level in the bulk liquid ranging from 1.5 to 35.2 mgO 2 L -1 . However in all the cases, the biofilm never experienced a fully penetrated oxygen profile and the Anammox bacteria could always survive in the inner part of the biofilm. Conversely, show that the prolonged exposure to a DO level of 5 mgO 2 L -1 would lead to the loss of Anammox activity, in agreement with Egli et al. [START_REF] Egli | Enrichment and characterization of an anammox bacterium from a rotating biological contactor treating ammonium-rich leachate[END_REF], who found out that the oxygen inhibits Anammox metabolism reversibly at low oxygen levels (air saturation level of 0.25-2 %) but probably irreversibly at higher levels. This can be explained by considering that the optimal oxygen concentration for performing a completely autotrophic nitrogen removal and allow the formation of an anoxic zone, depends on multiple parameters, such as the biofilm thickness and density, the COD content, the ammonium surface load, the temperature etc. On this basis, Vlaeminck et al. [START_REF] Vlaeminck | Aggregate size and architecture determine microbial activity balance for one-stage partial nitritation and anammox[END_REF] studied the effect of aggregate size on the proportion of microbial nitrite production and consumption. According to our results (Figures 5.8 and 5.9), they observed a higher Anammox activity with increasing aggregate size.

Conclusions

A mathematical model able to simulate the physical, chemical and biological processes prevailing in a multispecies biofilm for ammonia removal is presented. The model has been focused on the competition between AeAOB, NOB, AnAOB, and HB into biofilm and is able to evaluate the significant control that diffusion and DO exert on microbial stratification. Moreover, the proposed model adequately considers the bioconversion processes in the biological biofilm under dynamic conditions and it is able to predict the biomass stratification and substrate concentration trends for different DO levels and operational conditions. High DO level is responsible for loss of AnAOB from the reactor.

The model reproduces the different timescales at which changes in the concentrations of soluble components, total biomass and interspecies variations occur. Equations governing the free boundary value problem have been integrated numerically and the simulations reveal that the model is able to evaluate properly the effects that boundary conditions exert on bacterial competition.

In conclusion, the model can be used to optimize the anaerobic ammonia removal and control the bacterial stability in the system as it provides a qualitative understanding of microbial species and their activity within the biofilm.

-3.43 Y1 1 μ1X1 2 Growth of X3 -1 Y3 -i N,B -1 Y3 -1 1.14 1 1.14 1 μ3X3 3 Aerobic Growth of X4 -i N,B -1 Y4 -1 Y 4 1 μ4,1X4 1
Anoxic Growth of X4(NO 3 -NO 2 ) -i N,B -1 Y4 -(1 -1 Y4 ) 1 1.14 (1 -1 Y4 ) 1 1.14 1 μ4,2X4 4
Anoxic Growth of X4(NO 2 -N 2 ) -i N,B -1 Y4 (1 -1 Y4 ) 1 1.72 1 μ4,3X4 5
Growth of X2 -i N,B -1 Y2 1 Y2 1 -1.14 Y2 1 μ2X2 7 Aerobic Endogenous respiration of X1 i N,B -i N,I f I -(1 -f I ) -1 f I b1X1 8 Aerobic Endogenous respiration of X3 i N,B -i N,I f I -(1 -f I ) -1 f I b3X3 9 Aerobic Endogenous respiration of X4 i N,B -i N,I f I -(1 -f I ) -1 f I b4X4 10 Aerobic Endogenous respiration of X2 i N,B -i N,I f I -(1 -f I ) -1 f I b2X2 11 Anoxic Endogenous respiration of X1 i N,B -i N,I f I - 1-f I 2.86 -1 f I c1X1 12 Anoxic Endogenous respiration of X3 i N,B -i N,I f I - 1-f I 2.86 -1 f I c3X3 13 Anoxic Endogenous respiration of X4(NO 3 -NO 2 ) i N,B -i N,I f I 1-f I 1.14 - 1-f I 1.14 -1 f I c4,1X4 6
Anoxic Endogenous respiration of X4(NO 2 -N 2 ) i N,B -i N,I f I - 1-f I 1.72 -1 f I c4,2X4 14

Modelling multispecies biofilms including new bacterial species invasion

A mathematical model for multispecies biofilm evolution based on continuum approach and mass conservation principles is presented. The model can describe biofilm growth dynamics including spatial distribution of microbial species, substrate concentrations, attachment, and detachment, and, in particular, is able to predict the biological process of colonization of new species and transport from bulk liquid to biofilm (or vice-versa). From a mathematical point of view, a significant feature is the boundary condition related to biofilm species concentrations on the biofilm free boundary. These data, either for new or for already existing species, are not required by this model, but rather can be predicted as results. Numerical solutions for representative examples are obtained by the method of characteristics. Results indicate that colonizing bacteria diffuse into biofilm and grow only where favorable environmental conditions exist for their development.

Introduction

In both natural and artificial environments microorganisms often exist in an organized form known as a biofilm. Microbial biofilms are highly structured habitats consisting of surface-associated microorganisms enclosed in an exopolysaccharide matrix and organized into microcolonies. In some cases cell clusters are separated by interstitial voids and channels, which create a characteristic porous structure. An advantage of attaching to a surface is the ability to anchor to a preferred environment for bacterial growth. Moreover, the whole community benefits from the close spatial arrangement of different bacterial species and the potential for interaction and co-metabolism [START_REF] Wimpenny | Heterogeneity in biofilms[END_REF]. This structure might even be considered as an immobilized enzyme system in which the milieu and the enzyme activities are constantly changing and evolving to an approximately steady state [START_REF] Iw Isutherland | The biofilm matrix-an immobilized but dynamic microbial environment[END_REF]. On the way to development of mature biofilm, substrate concentrations become heterogeneous, allowing formation of microniches characterized by particular environmental conditions. These microniches provide growth conditions suitable for new species. The presence of relatively large channels and pores within the matrix structure might allow the entry of colonizing cells, present in the bulk liquid, and their establishment within the biofilm [START_REF] Iw Isutherland | The biofilm matrix-an immobilized but dynamic microbial environment[END_REF]. The newly colonizing cells can swim in the channels within the biofilm matrix and possibly even within the biofilm itself [START_REF] Wilson | Statistical quantification of detachment rates and size distributions of cell clumps from wild-type (pa01) and cell signalling mutant (jp1) pseudomonas aeruginosa biofilms[END_REF] and find favorable environmental conditions for growth. Successive colonization processes influence biofilm structure starting from the initial colonization phase of the substratum and, over time, affect microbial species distribution, both on the surface and within the biofilm structure.

The whole process of biofilm formation can be viewed in the context of primary ecological succession. According to the theory [START_REF] Cr | Changes in community properties during microbial succession[END_REF] and supported by experimental evidence, biofilm formation is initiated by pioneering microbial species which attach stochastically to the surface forming a monolayer. This first stage in biofilm formation process is characterized by high level of diversity and can be followed by a secondary colonization of bacteria that benefit from a protective environment in the biofilm and/or feed on the remnants of other bacteria [START_REF] Ac Martiny | Longterm succession of structure and diversity of a biofilm formed in a model drinking water distribution system[END_REF]. However, species richness is usually found to decline in the early-middle stages of biofilm development due to microbial competition for space and other needed resources.

Then, subsequently, in more mature communities, the number of niches again increases due to the formation of substrate gradients, accumulation of metabolic waste products that can be used as growth substrates by other microrganisms, and the development of a complex heterogeneous architecture with greater area for attachment and growth, all of which promote a second increase of species diversity connected to the appearance of more specialized populations.

In parallel to experimental investigations, complex mathematical models and numerical simulations have been developed to investigate development, structures, and ecological interactions of biofilms. However, little attention has been directed towards successional invasion. Here a mathematical model for multispecies biofilm formation and development is presented that does so. The model is based on a continuum approach for one space dimension [START_REF] Wanner | A multispecies biofilm model[END_REF][START_REF] Acunto | Qualitative analysis and simulations of a free boundary problem for multispecies biofilm models[END_REF][START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Atkinson | The overall rate of substrate uptake (reaction) by microbial films. part i-a biological rate equation[END_REF][START_REF] Kissel | Numerical simulation of mixedculture biofilm[END_REF] and then generalized to three-dimensions in [START_REF] Klapper | Mathematical description of microbial biofilms[END_REF][START_REF] Klapper | Finger formation in biofilm layers[END_REF][START_REF] Alpkvist | A multidimensional multispecies continuum model for heterogeneous biofilm development[END_REF] with the intention of predicting biofilm growth, space distribution of species, substrate trends, attachment and detachment. Important biological processes are added, specifically colonization of new species diffusing from bulk liquid to biofilm and developing of latent microbial species within the biofilm. The first was already included in preceding models, e.g. [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF][START_REF] Klapper | An exclusion principle and the importance of mobility for a class of biofilm models[END_REF], but only at a high price since boundary conditions for the invading species were needed on the free boundary. The present model does not need such data. In fact, rather, boundary values of all species, including the new ones, are determined self-consistently by the model. The diffusion of colonizing species from bulk liquid into the biofilm has been described by using a diffusion-reaction equation.

Spread through diffusion supposes a random character of mobility. Future prospectives include the possibility of characterizing the movement of bacteria by using chemotactically-driven mobility.

As an example of a target problem, consider latent anoxic bacteria, present in the bulk, diffusing into the biofilm and growing in the inner layer of biofilm where the oxygen concentration is equal to zero. In this particular case the anoxic bacteria concentration in the bulk is non-zero, while the value of this species on the free boundary is equal to zero. Notably, the growth process is hyperbolic and the diffusion process of new or latent species parabolic. The two are mutually connected but governed through different equations that are coupled by introducing a growth rate term arising from mobile species concentrations in the hyperbolic equations governing biofilm development and connected with parabolic equations governing concentration diffusive processes. In this way the model can take into account all colonizing cells of all bacterial species present in the bulk liquid that diffuse within the biofilm and grow when environmental conditions allow, thus selecting species that are best able to grow in the particular conditions of any specific microenvironment in the biofilm. The model is quite general and can handle any number of species, included the colonizing and latent species, and any number of substrates. Growth processes are governed by a system of nonlinear hyperbolic partial differential equations and concentration diffusive processes by parabolic partial differential equations. Substrate trends are governed by a system of semi-linear parabolic partial differential equations. All equations are mutually connected and lead to a free boundary value problem, presented in Sec. 6.2. Also, the generalization to 3D is briefly described. Some qualitative properties of solutions are proved, mainly to show the consistency of the model, in Sec. 6.3. Numerical solutions are based on the method of characteristics and the accuracy was checked by comparison to the equation n i=1 f i (z, t) = 1, which states that the sum of volume fractions gives one. In Sec. 6.4 the model is applied to a well-known biological process of heterotrophic and autotrophic bacteria competition. The results confirm the capability of the model to predict biomass distribution, substrate concentration trends over biofilm depth, and formation of new bacterial species.

Invasion model

Equations for biofilms

Consider multispecies biofilm growth in one space dimension. Denote by z the biofilm growth direction, assumed perpendicular to a substratum located at z = 0. The dynamics are governed by the equations

∂X i ∂t + ∂ ∂z (uX i ) = ρ i r M,i (z, t, X, S) + ρ i r i (z, t, ψ, S), (6.1
) i = 1, ..., n, 0 ≤ z ≤ L(t), t > 0, ∂u ∂z = n i=1 (r M,i + r i ), 0 < z ≤ L(t), t > 0, (6.2
)

∂ψ i ∂t - ∂ ∂z D M,i ∂ψ i ∂z = r ψ,i (z, t, ψ, X, S), (6.3)
i = 1, ..., n, 0 < z < L(t), t > 0, L(t) = u(L(t), t) + σ a (t) -σ d (L(t), t), t > 0, (6.4)
where: X i (z, t) = ρ i f i denotes the concentration of microorganism i, X = (X 1 , ..., X n ); f i (z, t) is the volume fraction of microbial species i; hence n i=1 f i = 1; ρ i denotes constant density; S j (z, t) denotes the concentration of substrate j, j = 1, ..., m, S = (S 1 , ..., S m ); u(z, t) is the velocity of the microbial mass; r M,i (z, t,X,S) is the specific growth rate; ψ i (z, t) concentration of planktonic species diffusing from bulk liquid to biofilm; ψ = (ψ 1 , ..., ψ n ); r i (z, t, ψ, S) is the specific growth rate due to planktonic species; r ψ,i (z, t, ψ, X, S) is the conversion rate of motile species; L(t) denotes biofilm thickness, free boundary; D M,i denotes the diffusivity coefficient of planktonic species; σ a (t) is the attachment biomass flux from bulk liquid to biofilm; σ d (L(t), t) is the detachment biomass flux from biofilm to bulk liquid.

Equations (6.1) follow from the mass balance for a generic control volume; though the term r i is new, these equations have formed the core of many biofilm models. The functions r i (z, t, ψ, S) denote specific growth rates due to planktonic species and derives from the diffusion of colonizing bacterial species from the bulk liquid. Equations (6.2) follow from n i=1 f i = 1. Note that it is assumed that the volume fractions occupied by planktonic species are negligible. Equations (6.3) are newly introduced here. They govern diffusion of bacterial cells through channels and pores and within biofilms. The functions r ψ,i (z, t, ψ, S) represent the loss due to the conversion of motile bacteria to non-motile state as well as gain from biofilm cells becoming planktonic. We do not include the latter effect here as it is likely not important for invasion of new species, and so will suppress the X dependence below. Note that if there is in fact conversion of biofilm to planktonic cells, there would be an accompanying sink in (6.1) though that term would likely be small. Finally, free boundary equation (6.4) follows from global mass balance and includes attachment and detachment terms.

Suitable initial-boundary conditions are associated to equations (6.1)-(6.4). Some of them are immediate; for example, no flux conditions on substratum imply u(0, t) = 0, ∂ψ i ∂z (0, t) = 0, t > 0. (6.5) Some others are constrained, since the initial concentration of colonizing species must be assigned to be identically zero. Designating the new bacterial species by indexes i = n 1 + 1, ..., n, the initial conditions for X i are

X i (z, 0) = ϕ i (z), i = 1, ..., n 1 , ϕ i (z) = 0, i = n 1 + 1, ..., n. 0 ≤ z ≤ L(t), (6.6)
For equations (6.3), we assume that the invading biological process starts at t = 0 and so the initial-boundary conditions

ψ i (z, 0) = ψ i0 (z) = 0, 0 ≤ z ≤ L(0), (6.7
)

ψ i (L(t), t) = ψ iL (t), t > 0, (6.8)
are added to (6.5) 2 , where the quite general functions ψ iL (t) represent the concentrations of the new species in the bulk liquid. Also, about the possible connection between dispersal and water channels, we note that this is a layer of complexity beyond the current model. Finally, the initial biofilm thickness is prescribed for equation (6.4) L(0) = L 0 . (6.9)

Remark 2.1 Equation (6.3) could be considered with different initial-boundary conditions depending on context. For example,

ψ i (z, 0) = 0, 0 ≤ z ≤ L 0 , ∂ψ i ∂z (0, t) = 0, ∂ψ i ∂z (L(t), t) = ψ iL (t), t > 0,
describe inoculation flux of a new species.

Remark 2.2 Explicitly note that no boundary conditions for X and its derivatives have been assigned on the free boundary z = L(t).

Equations for substrates

Substrate profiles are governed by the following equations

∂S j ∂t - ∂ ∂z D j ∂S j ∂z = r S,j (z, t, X, S), 0 < z < L(t), t > 0, (6.10)
where j = 1, ..., m, r S,j (z, t, X, S) is the conversion rate of substrate j, and D S,j denotes the diffusivity coefficient of substrate j, with initial conditions

S j (z, 0) = S j0 (z), 0 ≤ z ≤ L 0 , j = 1, ..., m. (6.11)
No substrate flux is assumed at the substratum boundary z = 0, i.e., ∂S j ∂z (0, t) = 0, t > 0, j = 1, ..., m. (6.12)

On the free boundary z = L(t), Dirichlet conditions

S j (L(t), t) = S jL (t), t > 0, j = 1, ..., m, (6.13)
or Neumann conditions

∂S j ∂z (L(t), t) = S jL (t), t > 0, j = 1, ..., m, (6.14)
or mixed condition, depending on the problem, can be prescribed.

3D Model

The 1D model presented in the previous sections can be generalized to 3D by starting from the model described in [START_REF] Klapper | Finger formation in biofilm layers[END_REF]. Denote by B t the 3D region occupied by the biofilm and let x = (x 1 , x 2 , x 3 ) be a generic point. Then,

X i = X i (x, t), f i = f i (x, t), Sj = S j (x, t) , u = u(x, t), r M,i = r M,i (x, t,X,S), ψ i = ψ i (x, t), r i = r i (x, t, ψ, S).
If u = ∇p, where p denotes the pressure within the biofilm, the equations governing biofilm and substrate evolution are written as

∂X i ∂t + ∇ • (X i ∇p) = ρ i r M,i (x, t, X, S) + ρ i r i (x, t, ψ, S), x ∈ B t , (6.15) ∇ 2 p = n i=1 (r M,i + r i ), x ∈ B t , (6.16
)

∂ψ i ∂t -D M,i ∇ 2 ψ i = r ψ,i (x, t, ψ, S), x ∈ B t , (6.17)
∂S j ∂t -D j ∇ 2 S j = r S,j (z, t, X, S), x ∈ B t . (6.18)

Qualitative properties of solutions

Consider the free boundary value problem described in Sections 6.2.1 and 6.2.2 and suppose that a uniqueness and existence theorem exists. Note, under suitable assumptions and some modifications, such a result could be proved as in [START_REF] B D'acunto | Analysis and simulations of the initial phase in multispecies biofilm formation[END_REF].

Let us show some simple properties of solutions to the new model in order to emphasize its consistency and capability to predict the formation of new bacterial species.

Property 3.1 Consider equations (6.1), rewritten for convenience,

∂X i ∂t + ∂ ∂z (uX i ) = ρ i r M,i (z, t, X, S) + ρ i r i (z, t, ψ, S), 0 ≤ z ≤ L(t), t > 0, (6.19
) where, now, i is a prefixed index. Suppose that

ϕ i (z) = 0, r M,i | X i =0 = 0, r i | ψ i =0 = 0. (6.20)
Then equation (6.19) admits the unique solution X i = 0. Therefore the new species, indexed i, does not develop (ϕ i = 0). Note that hypothesis (3.13) 2 is usually satisfied for microbial species. Thus the term r i is essential for colonization to occur in the model.

Property 3.2 Consider equation (6.3), rewritten for convenience,

∂ψ i ∂t - ∂ ∂z D M,i ∂ψ i ∂z = r ψ,i (z, t, ψ, S), 0 < z < L(t), t > 0, (6.21)
where i is the same index as in equation (6.19). Assume homogeneous initialboundary conditions

ψ i (z, 0) = 0, 0 ≤ z ≤ L(0), ∂ψ i ∂z (0, t) = 0, ψ i (L(t), t) = 0, t > 0. (6.22)
Then, equation (6.21) admits the unique solution ψ i = 0. If a species is not present in the bulk liquid, then it cannot diffuse into the biofilm and develop, as outlined in Property 3.1. (Also recall that the model does not allow biofilm cells to become planktonic, that is, r ψ,i is independent of X.) Note further that hypothesis (3.13) 3 provides an important constraint on term r i ; as a consequence, even if the function ψ j = 0 for another species j, j = i, still species i is not affected.

Consider now system (6.1) rewritten in terms of volume fraction

∂f i ∂t + u ∂f i ∂z = r M,i + r i -f i ∂u ∂z , i = 1, ..., n, 0 ≤ z ≤ L(t), t > 0. (6.23)
Introduce the characteristic-like lines z = z(z 0 , t) defined by

∂z ∂t (z 0 , t) = u(z(z 0 , t), t), z(z 0 , 0) = z 0 , 0 ≤ z 0 ≤ L 0 . (6.24)
Insert (6.24) in system (6.1) to obtain

ḟi = r M,i + r i -f i ∂u ∂z , i = 1, ..., n, z = z(z 0 , t), 0 ≤ z 0 ≤ L 0 , t > 0. (6.25)
Summing on i and using (6.2) results in

ḟ = (1 -f ) n i=1 (r M,i + r i ), z = z(z 0 , t), 0 ≤ z 0 ≤ L 0 , t > 0, (6.26)
where f = n i=1 f i . Since it is apparent that equation (6.26) admits the solution f = 1 we can state the following property. Property 3.3. Consider system (6.1) and assume that (6.2) holds. Then, under the hypothesis f (z 0 , 0) = 1 we have f (z(z 0 , t), t) = 1, 0 ≤ z 0 ≤ L 0 , t > 0.

Numerical solutions and applications

Numerical solutions to the free boundary problem stated in Section 3.2 have been obtained by using the method of characteristics, e.g. [START_REF] Acunto | Free boundary problem for an initial cell layer in multispecies biofilm formation[END_REF][START_REF] B D'acunto | Dynamic modeling of sulfate reducing biofilms[END_REF]. Accuracy was checked by comparison to the equation n i=1 f i (z, t) = 1. The simulations in this section have been performed by original software developed for the model presented in this work. Heterotrophic-autotrophic competition for space with oxygen as common substrate proposed in [START_REF] Wanner | A multispecies biofilm model[END_REF] is considered. This example is based on mass balance equations for substrates, products, and bacterial groups and includes the bio-chemical reactions of heterotrophic-autotrophic competition. The model considers the kinetics of microbial growth and decay and takes into account two groups of bacteria: heterotrophic bacteria (X 1 ) and autotrophic bacteria (X 2 ), and three components (substrates): ammonia (S 1 ), organic carbon (S 2 ) and oxygen (S 3 ). Oxygen is used for both organic carbon oxidation and nitrification. Oxidation of ammonia to nitrate by the autotrophs provides energy for autotrophic growth. Inert material is modelled as another microbial species, whose growth derives from the decay of heterotrophic and autotrophic biomass. Organic carbon concentration is expressed in terms of Chemical Oxygen Demand (COD) which represents a measure conventionally used in environmental chemistry to characterize indirectly the amount of organic compounds. In this specific context COD measures the amount of molecular oxygen per unit of material that would be needed to oxidize it to CO 2 [START_REF] Tchobanoglous | Wastewater Engineering: Treatment and Reuse[END_REF]. The stoichiometry and process rates of the model including expressions for r M,i and r S,j and the relevant parameter values come from [START_REF] Wanner | A multispecies biofilm model[END_REF]. The specific growth rate via planktonic species is expressed as

Parameter

r i = k col,i ψ i k ψ,i + ψ i S i k S i + S i S 3 k S3 + S 3 , i = 1, 2, (6.27)
where k col,i is the maximum colonization rate of motile species, k ψ,i is the kinetic constant for motile bacteria, and k S i is the half-saturation constant for substrate S i .

The planktonic species loss term is expressed as

r ψ,i = - 1 Y ψ,i r i = - k col,i Y ψ,i ψ i k ψ,i + ψ i S i k S i + S i S 3 k S3 + S 3 , i = 1, 2, (6.28)
where Y ψ,i is the yield of non-motile bacteria on motile species.

Numerical simulations demonstrate the capability of the model to predict biomass distribution, substrate concentration profiles over biofilm depth, and invasion of new bacterial species. The results are shown in Figures 6.1-6.4 (note that biofilm is growing from left to right). In particular, two sets of simulations using different colonizing bacterial species have been performed: simulation set 1 illustrates the dynamics of heterotrophic colonization and simulation set 2 analyzes the time evolution of autotrophic colonization. The objective is the evaluation of the effects of colonization of different bacterial species on biofilm growth, bacterial species distribution and substrate concentration profiles.

The values of ammonia, COD, and oxygen concentrations at the biofilm/bulk liquid interface used for model simulations are reported in Table 6.1.

Simulation Set 1: Heterotrophic colonization

The first simulation monitors over time the heterotrophic colonization of a biofilm initially constituted only by autotrophic bacteria. The initial biofilm thickness is 0.3 mm and the concentration of colonizing heterotrophic bacteria in the bulk liquid has been set to ψ 1 (L(t), t) = 1mgCODl -1 . Figure 6.1: Effect of heterotrophic colonization on bacterial volume fractions, substrate concentration profiles, and ψ 1 profile within biofilm after 1 (A1,A2,A3) , 2 (B1,B2,B3) , 3 (C1,C2,C3), 5 (D1,D2,D3) days.

In addition, ψ 2 (L(t), t) = 0. Simulation results are reported in Figures 6.1 A3, B3, C3, D3) and 6.1 (E3, F3, G3, H3), ψ 1 never reaches zero within biofilm, indicating that merely the contemporary presence of substrates and colonizing motile species can lead to the growth of non-motile bacteria. This fact is consistent with observations that substrate concentrations have a regulatory effect on the dynamics of biofilm structure since the colony size can be directly correlated with the substrate concentration profiles in the biofilm. Figures 6.1 (A1, B1, C1, D1) and 6.2 (E1, F1, G1, H1) show that heterotrophic bacteria develop only in the outmost part of biofilm where oxygen and COD are present. The higher growth rate of heterotrophic bacteria there allows the suppression of autotrophs. When the external invasion takes place, (Figures 6.2 (E1, F1)), it is possible to note that the heterotrophic bacteria penetrate into biofilm (Figures 6.2 (G1, H1)). The simulations were stopped at 30 days. By going on with the simulations, it would be possible to see that heterotrophs continue to penetrate into biofilm. As it is reasonable to expect, the COD shows a constant profile when heterotrophic bacteria are not present in biofilm, Figures 6.1 (A2, B2, C2). As soon as heterotrophic bacteria start to growth within biofilm, the COD concentration decreases, Figures 6.1 (D2) and 6.2 (E2, F2, G2, H2). In the outmost part of biofilm where oxygen concentration is high (Figures 6.2 (E2, F2, G2, H2)), according to experimental results [START_REF] Horn | Growth and decay in an auto-/heterotrophic biofilm[END_REF][START_REF] Wäsche | Influence of growth conditions on biofilm development and mass transfer at the bulk/biofilm interface[END_REF], ammonia and COD show decreasing profiles.

Simulation Set 2: Autotrophic colonization

Complementing the previous set, simulation set 2 investigates the autotrophic colonization of a biofilm constituted only by heterotrophic bacteria with an initial thickness of 0.3 mm, Figures 6. A concentration of colonizing autotrophic bacteria in the bulk liquid equal to ψ 2 (L(t), t) = 1mgCODl -1 has been considered. In addition, ψ 1 (L(t), t) = 0. The new species develops only in the inner part of the biofilm where the environmental conditions are favorable for autotrophs. See Figures 6.3 (A2, B2, C2, D2) and 6.4 (E2, F2, G2, H2) for substrate concentration profiles within biofilm at 1, 2, 3, 5 days simulation time and 7.5, 10, 20, 30 days simulation time respectively.

It is evident that autotrophs develop in the inner part of the biofilm where COD is not present (and thus the heterotrophs cannot grow) but where high levels of ammonia are available. Figures 6.3 (A1, B1, C1, D1) and 6.4 (E1, F1, G1, H1) show the biofilm structure in terms of bacterial volume fractions at different simulation times: 1, 2, 3, 5, 7.5, 10, 20, 30 days. After 7.5, 10 and 20 days, Figures 6.4 (E1, F1, G1), the biofilm is characterized by a typical heterotrophs-autotrophs stratification: autotrophs are dominant in the innermost part of biofilm and heterotrophs develop at the outmost layers [START_REF] Horn | Growth and decay in an auto-/heterotrophic biofilm[END_REF]. The selected boundary conditions, in particular the high ammonia concentration determine the strong autotrophic development, thus the biomass distribution in Figures 6.4 (H1), in confirmation of the consistency of the proposed model.

Conclusions

A mathematical model describing the invasion of new bacterial species into an already constituted biofilm has been presented, motivated by the importance of invasion of mobile, diffusing organisms. The biological process of microbial invasion is likely a frequent one in biofilms and can be viewed in the context of ecological succession or simply as a consequence of life in a fluctuating environment. Invasion may be an important event; even beyond the ability of a new species simply to establish itself in an already present community, synergistic interactions could arise between colonizing and resident species resulting in cooperative fitness benefits. In other cases the invasion of a new species may lead to the loss of resident species, contributing to a decrease in ecological richness. Depending on context, either of these might be desirable or undesirable -consider for example introduction of pathogens or probiotics to the gut microbiome.

Understanding the mechanisms that govern invasion is important in order to clarify community functions. This work represents a new approach for modelling the colonization process while using a continuum model based on hyperbolic equations. Qualitative properties of solutions showed the consistency of the model. A simple well-known example has been studied using numerical simulations developed by the method of characteristics. The results show that the model is able to predict the colonizing process in a reasonable way. A next step might be the application of the model to instances of multispecies biofilms in which resident species benefit from metabolic commensalism and an invading species might interfere, either positively or negatively.

Chapter 7 Conclusion and Future Work 7.1 Conclusion

Biofilm models have been recognized as useful tools for studying and exploring fundamental processes of such a complex system on a wide range of temporal and spatial scales as well as elucidating microbial competition and coexistence. Biofilm models are characterized by different levels of complexity: one-dimensional models assume the biofilm as planar; multidimensional models are able to describe the spatial heterogeneity of biofilm structure in multiple directions allowing us to capture most of the features of a biofilm system.

In Chapter 2, the different modeling approaches present in literature have been extensively analyzed with the aim of highlighting the main advantages and disadvantages characterizing each approach. In particular, biofilm models have been classified in two main categories, namely the continuum and discrete models; for each class of models the required inputs, the outputs and the approach used to simulate the main biofilm processes/features have been summarized. Based on these elements and the results achieved by comparative studies [START_REF] Morgenroth | Comparing biofilm models for a single species biofilm system[END_REF][START_REF] Hj Eberl | Modelling a spatially heterogeneous biofilm and the bulk fluid: Selected results from benchmark problem 2(bm 2)[END_REF][START_REF] Be Rittmann | Results from the multi-species benchmark problem(bm 3) using one-dimensional models[END_REF], general guidelines have been outlined in order to provide support to the selection of the most suitable modeling approach. Multidimensional models are characterized by high requirements in terms of input data, computational resource and mathematical skills of the user; their use has been exclusively confined to research while their application as engineering tools has been mostly limited by the high spatial resolution and the level of detail required for model calibration. Besides, multidimensional models usually describe local biofilm development and result not feasible for large scale applications, such as the reactor scale in the case of wastewater treatment. Despite the limitation for their dimensionality, one-dimensional models are characterized by a good results-to-effort-ratio and are often recommended for practical use [START_REF] Mašić | Investigation of a biofilm reactor model with suspended biomass[END_REF]. In addition, one-dimensional models, as well as all continuum models, are suitable for mathematical analysis and generate deterministic solutions. However, on the basis of the elucidations furnished in Chapter 2, they should not be pointed out 123 merely as simplified mathematical descriptions of biofilm growth, since modeling of bacterial biofilms, even in the one-dimensional case, lead to complex mixed hyperbolic-parabolic free boundary value problems, not easily accessible to qualitative analysis techniques. For this reason, they have been normally solved numerically by using custom made softwares.

In Chapter 3, a qualitative analysis of the solutions to the free boundary value problem governing bacterial population dynamics during the attachment phenomenon in the initial phase of biofilm formation has been provided. The exchange of particulate species between biofilm and bulk liquid, namely the attachment and detachment process, plays a crucial role in biofilm technology as it contributes to defining the microbial ecology of biofilm ecosystems. The discussion has been addressed to the case of an already constituted thin layer of biofilm, whose initial thickness is assumed to be strictly positive. The partial differential equations constituting the free boundary value problem have been converted into relatively simple ordinary differential equations. The key of this transformation relies on the introduction of the characteristic-like lines. The obtained ODE-type biofilm model has been converted into the equivalent integral form, incorporating the initial conditions. Properties of the solutions have been investigated: by using the fixed-point theorem, their existence and uniqueness have been proved. These results demonstrate that the method of characteristics may be used to obtain a qualitative insight into dynamic biofilm processes which so far have been studied primarily numerically. Besides, such a method has been used to integrate numerically the differential system with the aim of investigating the effects of attachment rate on biofilm performance.

The method of characteristics has been recognized as a valid alternative to the commonly known coordinate transformation method introduced by Wanner and Gujer [START_REF] Wanner | A multispecies biofilm model[END_REF] and applied to the free boundary value problem in order to fix the size of the domain in the 1D set-up. Actually, the method of characteristics may be applied even to the case of initial biofilm formation, while the coordinate transformation is not possible as the biofilm initial length needs to be set to zero.

In the further chapters, the transformed simpler to handle ODE-type biofilm model has been used to study with numerical techniques population dynamics in specific multispecies biofilms which experience microbial coexistence and/or competition. For what concerns the exchange with the bulk liquid, only detachment has been taken into account as it represents the prevailing transfer process for a mature biofilm. The accuracy of the numerical integration has been checked by testing the invariance of the mass conservation evaluated for the whole system. Microbial interactions within multispecies biofilms can be classified in antagonistic, such as competition over nutrients and growth inhibition, or synergistic. Synergistic interactions among different species usually predominate over antagonistic ones, revealing in many cases in metabolic cooperation with one species utilizing the metabolite produced by a neighboring species [START_REF] Elias | Multi-species biofilms: living with friendly neighbors[END_REF]. This phenomenon is usually defined as metabolic commensalism and allows the coexistence of different microbial groups through niche differentiation.

The inherent synergical interactions characterizing natural biofilm communities have shown to facilitate the simultaneous removal of various pollutants in wastewater treatment reactors.

In several biofilm-based processes, redox stratification is experienced due to the formation of strong concentration gradients of both electron donors and acceptors, and the accumulation of metabolic waste products, that can be used as growth substrates by other microorganisms. On the other hand, microbial population may experience antagonistic interactions, reflecting on the overall efficiency of the wastewater treatment process. A striking example is provided by the reduction of sulphate removal efficiencies in sulfidogenic reactors, treating acidic sulfate-containing wastewater, mainly due to the accumulation of acetate and the presence of different species competing for a common substrate. Based on these considerations, in Chapter 4 the dynamical response of a multispecies biofilm performing a sulfidogenic process at different COD/SO 2- 4 ratios has been investigated. The bacterial community has been supposed to be constituted by sulfate-reducing bacteria, acetogens and methanogens; the reactor is simulated to be fed with a mixture of sulfate, acting as the electron acceptor and lactate as electron donor. The simulation results have found confirmation in experimental research. COD/sulfate ratio has been recognized as a crucial variable in the optimization of lactate utilization via oxidation in preference to fermentation and in the maximization of the efficiency of biological sulfate reduction. The model could be further used to understand how the operational conditions can affect biofilm growth and microbial competition, recommend start-up and feeding strategies as well as provide insights regarding biofilm reactor configurations. Microbial interactions in multispecies biofilms also play a crucial role in the case of the ANaerobic AMMonia OXidation (Anammox) process, which represents one of the most promising innovative techniques for the biological removal of nitrogen from wastewater. In Chapter 5, the one-dimensional mathematical model has been applied to analyse and predict microbial interactions within multispecies biofilms including Anammox pathway. Based on the results achieved in previous models, the related processes of organic carbon oxidation, denitrification, nitrification and Anammox have been combined in order to evaluate the influence of bulk substrate concentrations and diffusion on microbial stratification. Specific scenarios have been analyzed: first the effect of different bulk liquid oxygen concentration on microbial population stratification has been assessed; second, the model response to a variation of operational parameters (dissolved oxygen level and shear stress condition) has been investigated. The results reveal that in a thick multispecies biofilm, including heterotrophic, aerobic autotrophic nitrifying and Anammox bacteria, oxygen diffusion limitation determines the formation of both aerobic and anoxic microenvironments favouring interspecies competition. In contrast, oxygen excess causes a disturbance on microbial interactions leading to Anammox bacteria loss. Moreover, it is evident that biofilm erosion, which is mainly dictated by the specific reactor hydrodynamics, provides a contributory effect to the loss of Anammox bacteria as it exposes them to higher oxygen concentrations internally. The combined effect of lower shear stress and oxygen may be seen as a possible operating strategy to facilitate the coexistence of the different microbial groups. Nevertheless, some assumptions made in modeling this complex biological system may lead to an oversimplified representation of inert mass production and exchange of microbial products.

Indeed, the incorporation of microbial by-products and EPS formation, kinetics expression of heterotrophic growth on hydrolysis by-products may pave for a more accurate modeling of Anammox process in order to effectively control the performance of Anammox systems regarding effluent quality and process dynamics.

In Chapter 6, a new modeling approach to simulate bacterial species colonization of a constituted multispecies biofilm has been presented. The model includes all the main biofilm processes, such as substrate diffusion, attachment, detachment, microbial growth, biomass spreading, but it is mainly focused on the description of new species invasion and transport from bulk liquid to biofilm (or vice-versa). The model has been conceived in the framework of continuum models with the intention of providing a modeling alternative to the widely used Wanner-Recheirt model [START_REF] Wanner | Mathematical modeling of mixed-culture biofilms[END_REF] for the specifically colonization of new species diffusing from bulk liquid to biofilm. The basic idea constists in coupling the diffusion of colonizing species with the main biofilm development by introducing a growth rate term which depends on the concentration of motile species. This expedient has allowed us to keep the growth process hyperbolic with no need of special boundary conditions for the growth of microbial species within the biofilm. The resulting model is constituted by a system of non linear hyperbolic partial differential equations, a system of semi-linear parabolic equations governing substrate diffusive processes and a system of diffusion-reaction equations for simulating the spread of motile bacteria within biofilm matrix. The models is amenable to qualitative analysis and some properties on the consistency of solutions have been shown. By using the method of characteristics, the model has been solved in numerical form and addressed to study the colonization of heterotrophic-autotrophic biofilm. The latter represents only one of the possible applications as microbial invasion represents a likely frequently biological process in biofilms.

Future directions

One-dimensional biofilm models are considered to be a valuable tool in biofilm process research. They can address the qualitative and/or quantitative analysis and constitute a foundation and a framework within which further modifications and developments can be made [START_REF] Mašić | Investigation of a biofilm reactor model with suspended biomass[END_REF]. In this context, mathematical modeling can be used to solidify existing or derive new models, and to analyze and explain their behaviour. Based on the results achieved in this study, biofilm research areas which may benefit from further progress can be individuated. They are briefly presented in the following sections.

Invasion of new species

The invasion of colonizing species into a constituted biofilm represents a frequent biological process in multispecies biofilms.

The consequences of the establishment of new species in a pre-existing microbial community are varied: the resident microbial species could benefit from synergistic interactions with the invading species, the first producing metabolites which inhibit their metabolism at high concentration, the second growing on such products. This results in an increased ecological richness and in a cooperative fitness. Not all the interactions established between resident and colonizing species have a positive feedback: the invasion may lead to the loss of resident species contributing to a decrease in microbial diversity. The mechanisms governing this biological process are poorly understood: the invading species may be driven by a random movement within the biofilm or a more complex chemotactic approach may be used. An extension of the study to the case of chemotactic movement of motile bacteria within the porous structure of biofilms would be desirable. To the best of our knowledge, the movement of motile species driven by chemotaxis has never been included in a biofilm model but most of the literature studies address the case of planktonic bacterial state. Furthermore, the invasion model could be combined with the modelling of the initial phase of biofilm formation.

Effect of sloughing events and volumetric detachment on microbial composition

In Chapter 5, the effect of continuous detachment, modeled as a function of biofilm thickness, on microbial population dynamics in multispecies biofilms including Anammox pathway has been investigated. However, not only erosion but also single sloughing events and volumetric detachment may influence the results.

Morgenroth and Wilderer [START_REF] Morgenroth | Influence of detachment mechanisms on competition in biofilms[END_REF] have shown in the case of autotrophic/heterotrophic biofilm that detachment has a significant impact on competition and structure of the microbial population. Therefore, a natural extension of the developed model for simulating Anammox competition would be to implement different detachment mechanisms instead of simply considering the erosion performed by liquid flow.

Calibration and sensitivity analysis

For adequately representing experimental data with a mathematical model, calibration of the model to collected data is inevitable and required to effectively use biofilm models in engineering design [START_REF] Brockmann | Biofilm model calibration and microbial diversity study using monte carlo simulations[END_REF]. The biofilm reactor model calibration process represents a "hot" topic in biofilm research: efforts to define the type of information that is required to calibrate a biofilm reactor model, sampling and tests procedures are limited [START_REF] Boltz | Mathematical modelling of biofilms and biofilm reactors for engineering design[END_REF]. Actually, only few studies are dedicated to the definition of a specific calibration protocol. This lack may be related to the availability of widely range of significantly different biofilm modeling approaches, various empirical mass transfer boundary layer formulations, and host of system-specific model parameters. In future, more efforts are planned to be dedicated to the development of a calibration algorithm and a specific sensitivity analysis approach.

4 . 4 6 . 3 121 List of Tables 3 . 1

 446312131 Figure adapted from these two papers. . . . . 3.1 Characteristic-like lines . . . . . . . . . . . . . . . . . . . . . . . 3.2 Effect of attachment rate (σ) on the volumetric fraction of the bacterial species in biofilm. A: σ = 5 mmd -1 ; B: σ = 1 mmd -1 ; C: σ = 0.5 mmd -1 . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Effect of attachment rate (σ) on the substrate trends in biofilm. A: σ = 5 mmd -1 ; B: σ = 1 mmd -1 ; C: σ = 0.5 mmd -1 . . . . . . . 4.1 Main pathways of the biological process. . . . . . . . . . . . . . 4.2 Substrate trends in the biofilm (A) and bacterial volumetric fractions (B) in the biofilm for a COD/SO4 2-ratio = 0.5. Dotted line: sulfate concentration; dashdot line: COD; continuous line: acetate concentration. . . . . . . . . . . . . . . . . . . . . . 4.3 Substrate trends in the biofilm (A) and bacterial volumetric fractions (B) in the biofilm for a COD/SO 2- 4 ratio = 1.5. Dotted line: sulfate concentration; dashdot line: COD; continuous line: acetate concentration. . . . . . . . . . . . . . . . . . . . . . . . ix x Substrate trends in the biofilm (A) and bacterial volumetric fractions (B) in the biofilm for a COD/SO 2- 4 ratio = 2. Dotted line: sulfate concentration; dashdot line: COD; continuous line: acetate concentration. . . . . . . . . . . . . . . . . . . . . . . . 5.1 Main microbial interactions of the simulated biological process . . 5.2 Effects of applied DO (3 mg/L) on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 10 (A,B), 50 (C,D), 100 (E,F), 150 (G,H) days time simulation. . . 5.3 Effects of applied DO (3 mg/L) on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 200 (I,L), 250 (M,N), 300 (O,P) days time simulation. . . . . . . . 5.4 Effects of applied DO (5 mg/L) on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 10 (A,B), 50 (C,D), 100 (E,F), 150 (G,H) days time simulation. . . 5.5 Effects of applied DO (5 mg/L) on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 200 (I,L), 250 (M,N), 300 (O,P) days time simulation. . . . . . . . 5.6 Effects of a DO change from 3 mg/L (maintained over 100 days) to 5 mg/L on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 10 (A,B), 50 (C,D), 100 (E,F), 150 (G,H) days time simulation. . . . . . . . . . 5.7 Effects of a DO change from 3 mg/L (maintained over 100 days) to 5 mg/L on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 200 (I,L), 250 (M,N), 300 (O,P) days time simulation. . . . . . . . . . . . . . . 5.8 Effects of a change in applied shear stress constant from 50m -1 d -1 (maintained over 100 days) to 150m -1 d -1 on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 10 (A,B), 50 (C,D), 100 (E,F), 150 (G,H) days time simulation. . . . . . . . . . . . . . . . . . . 5.9 Effects of a change in applied shear stress constant from 50m -1 d -1 (maintained over 100 days) to 150m -1 d -1 on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 200 (I,L), 250 (M,N), 300 (O,P) days time simulation. . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Effect of heterotrophic colonization on bacterial volume fractions, substrate concentration profiles, and ψ 1 profile within biofilm after 1 (A1,A2,A3) , 2 (B1,B2,B3) , 3 (C1,C2,C3), 5 (D1,D2,D3) days. 6.2 Effect of heterotrophic colonization on bacterial volume fractions, substrate concentration profiles, and ψ 1 profile within biofilm after 7.5 (E1,E2,E3) , 10 (F1,F2,F3) , 20 (G1,G2,G3), 30 (H1,H2,H3) days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Effect of autotrophic colonization on bacterial volume fractions, substrate concentration profiles, and ψ 2 profile within biofilm after 1 (A1,A2,A3) , 2 (B1,B2,B3) , 3 (C1,C2,C3), 5 (D1,D2,D3) days. 120 6.4 Effect of autotrophic colonization on bacterial volume fractions, substrate concentration profiles, and ψ 2 profile within biofilm after 7.5 (E1,E2,E3) , 10 (F1,F2,F3) , 20 (G1,G2,G3), 30 (H1,H2,H3) days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Stoichiometry and rate laws for microbial processes. HG = heterotroph growth; AG = autotroph growth; HER = heterotroph endogenous respiration; AER = autotroph endogenous respiration; HD = heterotroph decay; AG = autotroph decay. . . . . . . . . . . 3.2 Operational parameters used for model simulations. . . . . . . . . 4.2 Operational parameters used for model simulation . . . . . . . . . 4.1 Petersen Matrix of the proposed model . . . . . . . . . . . . . . . 4.3 Kinetic, stoichiometric and diffusion coefficients used in the model 5.1 Kinetic and Stoichiometric Parameters used for Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Peterson matrix of the proposed model . . . . . . . . . . . . . . . 5.3 Overview of the different modeling scenarios executed in this study 6.1 Operational parameters used for model simulations . . . . . . . .

Figure 2 . 1 :

 21 Figure 2.1: Schematic representation of biofilm model classification adopted in this review and based on biomass representation and dimensionality: a) 1D continuum models; b) multidimensional continuum models; c) Cellular Automata; d) Individual based Models (for multidimensional models only the 2D representation has been reported).

Figure 2 . 5 :

 25 Figure 2.5: Spreading mechanisms adopted by: a), c) Xavier et al. [131]; b) Lardon et al. [141]. Figure adapted from these two papers.

Figure 3 . 1 :

 31 Figure 3.1: Characteristic-like lines

(3. 32 )

 32 Inserting equation(3.31) into(3.29) gives the integral equation for c(t 0 , t)c(t 0 , t) = L(t 0 ) + ζ 0 , τ), τ, x(ζ 0 , τ), s(ζ 0 , τ)

. 49 )

 49 In this situation systems (3.41)-(3.43) and (3.45)-(3.47) and equation (3.48) can be solved in series as shown in the following theorems. Consider system (3.41)-(3.43) and introduce some new notations. Firstly, let us redefine the vector x x = (x 1 , ..., x n , x n+1 , x n+2 ) (3.50)

. 57 ) 1

 571 Theorem If hypotheses(3.54)-(3.57) hold and ϕ i ∈ C(0, L 0 ), then there exists a unique continuous solution to system (3.41)-(3.43), x i ∈ C((0, L 0 ) × (0, T )).

. 62 )

 62 Map (3.60)-(3.62) maps Σ 1 into itself because of hypotheses (3.54)-(3.55). In addition, it is a contractive map. Indeed, setting ỹ = Ax we derive

Summing ( 3 .

 3 63)-(3.65) gives ||y -ỹ|| ≤ λ||x -x||,

. 74 ) 2

 742 Theorem If hypotheses(3.71)-(3.74) hold and ψ i , σ ∈ C(0, T ), i = 1, ..., n, then there exists a unique continuous solution x i ∈ C((0, T ) × (0, T )), i = 1, ..., n + 2 to system (3.46)-(3.48).

Now, system ( 3 .

 3 85)-(3.86) is suitably transformed by using positions (3.26), (3.32) and the change of variables ζ = z(ζ 0 , t), ζ = c(τ 0 , t) introduced in Section 3.3. Therefore, for 0 ≤ z 0 ≤ L 0 we get

  j (ζ(ζ 0 , τ), τ, x(ζ 0 , τ), s(ζ 0 , τ))N j (c(t 0 , t), ζ(ζ 0 , τ), t-τ )

Figure 3 . 2 :

 32 Figure 3.2: Effect of attachment rate (σ) on the volumetric fraction of the bacterial species in biofilm. A: σ = 5 mmd -1 ; B: σ = 1 mmd -1 ; C: σ = 0.5 mmd -1 .

Figure 3 . 3 :

 33 Figure 3.3: Effect of attachment rate (σ) on the substrate trends in biofilm. A: σ = 5 mmd -1 ; B: σ = 1 mmd -1 ; C: σ = 0.5 mmd -1 .

Figure 4 . 1 :

 41 Figure 4.1: Main pathways of the biological process.

Figures 4. 2

 2

  , 4.3 and 4.4 show the results of model simulations, named respectively set A, set B and set C, performed to assess the COD/SO 2- 4 ratio effect on the reactor performances in terms of bacterial volume fractions (Figures 4.2 (B), 4.3 (B), 4.4 (B)) and concentration trends of substrates (Figures 4.2 (A), 4.3 (A), 4.4 (A)) within biofilm for a 10 day simulation. COD/SO 2-

1 ]Figure 4 . 2 : 1 ]Figure 4 . 3 :

 142143 Figure 4.2: Substrate trends in the biofilm (A) and bacterial volumetric fractions (B) in the biofilm for a COD/SO4 2-ratio = 0.5. Dotted line: sulfate concentration; dashdot line: COD; continuous line: acetate concentration.

Figures 4. 2

 2 Figures 4.2, 4.3 and 4.4 show biomass distribution and substrate concentration trends at 0.5, 1.5 and 2 COD/SO 2- 4 ratio respectively. As shown in Figures 4.2 (B), 4.3 (B) and 4.4 (B), after 10 days, the biomass stratification appears visible:LDSRB (C) and LDSRB (I) prevail in the outer layer of biofilm where sulfate and lactate remain abundant. In the deepest zone of the biofilm, characterized by a low level of sulfate and lactate, due to substrate diffusion coupled with microbial consumption, the MA compete for space with other microbial species. Indeed the MA are mostly present in the deepest zone of the biofilm where the optimal conditions for their growth are established. With a COD/SO 2- 4 ratio of 0.5 (Figure4.2) the AB are present at the inner layer of the biofilm, where the concentration of sulfate is lower while both LDSRB (C) and LDSRB (I) are found to be predominant at the outmost layer of the biofilm. In particular, LDSRB (C) and LDSRB (I) represent the most abundant species in the biofilm showing that in the presence of excess sulfate, the quantitative oxidation of lactate to acetate or CO 2 coupled to sulfate reduction is the dominant reaction.A similar result has been achieved in[START_REF] Lb Celis-García | Performance of a downflow fluidized bed reactor under sulfate reduction conditions using volatile fatty acids as electron donors[END_REF], where the authors experienced high participation of SRB, and in particular of LDSRB (I) , on COD removal in a down-flow fluidized-bed reactor. Concerning sulfate reduction, it is possible to

  Figures 4.2, 4.3 and 4.4 show biomass distribution and substrate concentration trends at 0.5, 1.5 and 2 COD/SO 2- 4 ratio respectively. As shown in Figures 4.2 (B), 4.3 (B) and 4.4 (B), after 10 days, the biomass stratification appears visible:LDSRB (C) and LDSRB (I) prevail in the outer layer of biofilm where sulfate and lactate remain abundant. In the deepest zone of the biofilm, characterized by a low level of sulfate and lactate, due to substrate diffusion coupled with microbial consumption, the MA compete for space with other microbial species. Indeed the MA are mostly present in the deepest zone of the biofilm where the optimal conditions for their growth are established. With a COD/SO 2- 4 ratio of 0.5 (Figure4.2) the AB are present at the inner layer of the biofilm, where the concentration of sulfate is lower while both LDSRB (C) and LDSRB (I) are found to be predominant at the outmost layer of the biofilm. In particular, LDSRB (C) and LDSRB (I) represent the most abundant species in the biofilm showing that in the presence of excess sulfate, the quantitative oxidation of lactate to acetate or CO 2 coupled to sulfate reduction is the dominant reaction.A similar result has been achieved in[START_REF] Lb Celis-García | Performance of a downflow fluidized bed reactor under sulfate reduction conditions using volatile fatty acids as electron donors[END_REF], where the authors experienced high participation of SRB, and in particular of LDSRB (I) , on COD removal in a down-flow fluidized-bed reactor. Concerning sulfate reduction, it is possible to

1 ]Figure 4 . 4 :

 144 Figure 4.4: Substrate trends in the biofilm (A) and bacterial volumetric fractions (B) in the biofilm for a COD/SO 2- 4 ratio = 2. Dotted line: sulfate concentration; dashdot line: COD; continuous line: acetate concentration.

2

 2

  and 5.3. After 10 days (Figure5.2 (A,B)) biomass stratification is not visible in the inner layers of biofilm structure; in effect microbial species distribution is still affected by initial conditions.

Figure 5 . 2 :

 52 Figure 5.2: Effects of applied DO (3 mg/L) on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 10 (A,B), 50 (C,D), 100 (E,F), 150 (G,H) days time simulation.

Figure 5 . 3 :

 53 Figure 5.3: Effects of applied DO (3 mg/L) on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 200 (I,L), 250 (M,N), 300 (O,P) days time simulation.

Figure 5 . 4 (

 54 Figures 5.4 and 5.5 display variations of biomass distribution and substrate concentration trends within biofilm during 300 days of system operation at DO = 5 mg/L. Similarly to Figure5.2, after few days biofilm population structure is still influenced by the initial condition (Figure5.4 (A)), but as showed in Figure5.4 (B), the biofilm experiences a fully penetrated oxygen profile. This diversity strongly reflects on microbial interactions and in particular on AnAOB metabolism. In fact, the AnAOB volume fraction decreases over time until reaching the completely loss from the reactor (Figure5.4 (A,C,E,G), Figure5.5 (I,M,O)). Figure5.4 (B,D,F,H) and Figure5.5 (L,N,P) provide the variations of substrate concentrations over time. Ammonium and COD concentrations are low within biofilm in all the simulations. Nitrate concentration keeps higher compared to the results achieved for DO =3 mg/L; this trend is strictly connected with the new biofilm microbial composition: nitrate is produced by NOB which experience a perfect survival condition thanks to the simultaneous presence of excess nitrite, produced by AeAOB, and oxygen. The NOB competition for nitrite with AnAOB and for oxygen with AeAOB is relieved, and the cells can start to produce nitrate. Nitrite concentration keeps lower compared to the previous results as it is oxidized by NOB whose metabolism is characterized by higher rate than AnAOB consumption rate. The biofilm stratification starts to be evident after 50 days and reaches the following expected configuration (Figure 5.4 (E,G), Figure 5.5 (I,M,O)): the AeAOB continue to occupy the outer layer of biofilm out-competing

Figure 5 . 4 :

 54 Figure 5.4: Effects of applied DO (5 mg/L) on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 10 (A,B), 50 (C,D), 100 (E,F), 150 (G,H) days time simulation.

Figure 5 . 5 :

 55 Figure 5.5: Effects of applied DO (5 mg/L) on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 200 (I,L), 250 (M,N), 300 (O,P) days time simulation.

Figures 5 .

 5 Figures 5.6 and 5.7 show the results of a biofilm which experienced an increase in the DO value. The simulation has been run considering as initial condition the results of the 100 days time simulation with DO = 3 mg/L (Figure 5.2 (E,F)) and assuming as boundary conditions DO = 5 mg/L, ammonium surface load of 2 gm -2 d -1 , COD surface load of 0.2 gm -2 d -1 . The program run for 300 days and the results clearly indicate that there is no significant difference in the biomass population between Figure 5.5 (O) and Figure 5.7 (O). This is caused by the fact that the dynamics of the problem are influenced only by the boundary conditions, which determine the achievement of the same final microbial distribution. This result marks the great accuracy of the numerical approach used.Figure 5.6 (A,C,E,G) and Figure 5.7 (I,M,O) show that prolonged exposure to high oxygen level can lead to the loss of AnAOB from the reactor. In particular, Figure 5.7 (I)

  Figures 5.6 and 5.7 show the results of a biofilm which experienced an increase in the DO value. The simulation has been run considering as initial condition the results of the 100 days time simulation with DO = 3 mg/L (Figure 5.2 (E,F)) and assuming as boundary conditions DO = 5 mg/L, ammonium surface load of 2 gm -2 d -1 , COD surface load of 0.2 gm -2 d -1 . The program run for 300 days and the results clearly indicate that there is no significant difference in the biomass population between Figure 5.5 (O) and Figure 5.7 (O). This is caused by the fact that the dynamics of the problem are influenced only by the boundary conditions, which determine the achievement of the same final microbial distribution. This result marks the great accuracy of the numerical approach used.Figure 5.6 (A,C,E,G) and Figure 5.7 (I,M,O) show that prolonged exposure to high oxygen level can lead to the loss of AnAOB from the reactor. In particular, Figure 5.7 (I)

Figure 5 . 6 :

 56 Figure 5.6: Effects of a DO change from 3 mg/L (maintained over 100 days) to 5 mg/L on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 10 (A,B), 50 (C,D), 100 (E,F), 150 (G,H) days time simulation.

Figure 5 . 7 :

 57 Figure 5.7: Effects of a DO change from 3 mg/L (maintained over 100 days) to 5 mg/L on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 200 (I,L), 250 (M,N), 300 (O,P) days time simulation.

Figure 5 . 2 (

 52 E,F)). The shear stress constant has been increased from 50 m -1 d -1 to 150 m -1 d -1 and kept constant for all the course of the simulation. The results in terms of biofilm composition and substrate profiles as a function of time are shown in Figures 5.8 and 5.9. Results clearly demonstrate a link between shear stress, diffusion and composition of the microbial communities.

Figure 5 . 8 :

 58 Figure 5.8: Effects of a change in applied shear stress constant from 50m -1 d -1 (maintained over 100 days) to 150m -1 d -1 on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 10 (A,B), 50 (C,D), 100 (E,F), 150 (G,H) days time simulation.

Figure 5 . 8 (

 58 B,D,F,H), Figure5.8 (L,N,P)).

Figure 5 . 9 :

 59 Figure 5.9: Effects of a change in applied shear stress constant from 50m -1 d -1 (maintained over 100 days) to 150m -1 d -1 on bacterial population distribution (left) and substrate concentration trends (right) within biofilm after 200 (I,L), 250 (M,N), 300 (O,P) days time simulation.

and 6 . 2 .

 62 In particular, the results are in terms of bacterial volume fractions Figures 6.1 (A1, B1, C1, D1), substrate concentration profiles Figures 6.1 (A2, B2, C2, D2), and ψ 1 profiles Figures 6.1 (A3, B3, C3, D3) within biofilm at 1, 2, 3, 5 days simulation time and bacterial volume fractions Figures 6.1 (E1, F1, G1, H1), substrate concentration profiles Figures 6.1 (E2, F2, G2, H2), and ψ 1 profiles

Figures 6 . 1 (

 61 Figures 6.1 (E3, F3, G3, H3) within biofilm at 7.5, 10, 20, 30 days simulation time respectively. The simulations show that colonizing bacteria diffuse into biofilm and grow only where there are favorable environmental conditions for their development, (Figures 6.1 (C1, D1) and Figures 6.2 (E1, F1, G1, H1)), as determined by substrates trends (Figures 6.1 (A2 B2 C2 D2), Figures 6.2 (E2 F2 G2 H2)).

Figure 6 . 2 :

 62 Figure 6.2: Effect of heterotrophic colonization on bacterial volume fractions, substrate concentration profiles, and ψ 1 profile within biofilm after 7.5 (E1,E2,E3) , 10 (F1,F2,F3) , 20 (G1,G2,G3), 30 (H1,H2,H3) days.

  (

Figure 6 . 3 :

 63 Figure 6.3: Effect of autotrophic colonization on bacterial volume fractions, substrate concentration profiles, and ψ 2 profile within biofilm after 1 (A1,A2,A3) , 2 (B1,B2,B3) , 3 (C1,C2,C3), 5 (D1,D2,D3) days. The diffused substrate concentration profiles, for eight different time simulations, are shown in Figures 6.1 (A2, B2, C2, D2) and 6.2 (E2, F2, G2, H2).

3

 3

  and 6.4.

Figure 6 . 4 :

 64 Figure 6.4: Effect of autotrophic colonization on bacterial volume fractions, substrate concentration profiles, and ψ 2 profile within biofilm after 7.5 (E1,E2,E3) , 10 (F1,F2,F3) , 20 (G1,G2,G3), 30 (H1,H2,H3) days.

  By considering the change of variable ζ = z(ζ 0 , t) in the first integral and ζ = c(τ 0 , t) in the second one

					t, X(ζ, t), S(ζ, t)) dζ	(3.30)
		+	c(t 0 ,t)	G(ζ, t, X(ζ, t), S(ζ, t)) dζ,
				z(L 0 ,t)
	u(c(t 0 , t), t) =	0	L 0	G(z(ζ 0 , t), t, x(ζ 0 , t), s(ζ 0 , t))	∂z ∂ζ 0	(ζ 0 , t) dζ 0	(3.31)

  γ4τ |x h (t 0 , τ) -xh (t 0 , τ)|dτ. ≤ λ n+2 ||x -x||/γ 4 .

	Summing (3.80)-(3.82) gives
		||y -ỹ|| ≤ λ||x -x||,
	where	λ =	n i=1	λ i γ 4	+	2λ n+1 γ 3 γ 4	+	λ n+2 γ 4
	Finally, consider equation (3.48). By using Th. 1 and Th. 2, we immediately
	obtain the solution for this equation.
	Interesting properties of solutions to systems (3.41)-(3.43) and (3.46)-(3.48)
	can be proved as in

can be made < 1 if γ 3 and γ 4 are selected large enough. Therefore, y = Bx is a contractive map and the theorem is proved.

  .2.

	Parameter	Unit	Set A Set B Set C
	COD Flux	gm -2 d -1	0.4	0.4	0.4
	Ammonia Flux	gm -2 d -1	0.8	0.8	0.8
	Oxygen Concentration	mgl -1	8	8	8
	Time Simulation	h	24	24	24
	Initial Biofilm thickness	cm	0.03	0.03	0.03
	Attachment Rates	mmd -1	5	1	0.5
	Initial Volume Fraction of HB	-	0.65	0.65	0.65
	Initial Volume Fraction of AB	-	0.34	0.34	0.34
	Initial Volume Fraction of Inert	-	0.01	0.01	0.01
	Table 3.2: Operational parameters used for model simulations.

Table 4 . 2

 42 .2. The model has been addressed to evaluate the microbial structure of a sulfate-reducing biofilm as affected by changing COD to sulfate ratios.

	Parameter	Unit	Set A Set B Set C
	COD Concentration	mgCODl -1	0.2	0.3	0.1
	Sulfate concentration	mgl -1	0.1	0.2	0.2
	Time Simulation	d	10	10	10
	Initial Biofilm thickness	μm	300	300	300
	Initial Volume Fraction of LDSRB (C)	-	0.3	0.3	0.3
	Initial Volume Fraction of LDSRB (I)	-	0.3	0.3	0.3
	Initial volume fraction of AB	-	0.2	0.2	0.2
	Initial volume fraction of MA	-	0.2	0.2	0.2
	Initial volume fraction of Inert	-	0	0	0
	Shear constant	m -1 d -1	2000 2000 2000

: Operational parameters used for model simulation

Table 4 .

 4 3: Kinetic, stoichiometric and diffusion coefficients used in the model

	2-4	ratio on the

Table 5 .

 5 1: Kinetic and Stoichiometric Parameters used for Numerical Simulations

	1	Y1

Table 5 . 3 :

 53 Overview of the different modeling scenarios executed in this study

	Chapter 6

Table 6 .

 6 1: Operational parameters used for model simulations

		Unit	Set 1 Set 2
	COD Concentration	mgl -1	3	3
	Ammonia Concentration	mgl -1	13	13
	Oxygen Concentration	mgl -1	8	8
	Time Simulation	d	10	10
	Initial Biofilm thickness	mm	0.3	0.3
	k ψ,1 k ψ,2	d -1 d -1	0.002 0.002 0.002 0.002
	Initial Volume Fraction of X 1	-	0.0	1.0
	Initial Volume Fraction of X 2	-	1.0	0.0

This chapter was published as: D'Acunto, B., Esposito, G., Frunzo, L., Mattei, M.R. and Pirozzi, F. (2013). Analysis and simulations of the initial phase in multispecies biofilm formation. Communications in Applied and Industrial Mechanics, 2013, pp. 1-23, DOI: 10.1685/journal.caim.448.

ratio = 1.5. Dotted line: sulfate concentration; dashdot line: COD; continuous line: acetate concentration.

Acknowledgments

This chapter was published as: Mattei, M.R., D'Acunto, B., Esposito, G., Frunzo, L. and Pirozzi, F. (2014). Mathematical modeling of competition and coexistence of sulfate-reducing bacteria, acetogens and methanogens in multispecies biofilms. Desalination and Water Treatment, 2014, pp. 1-9, DOI: 10.1080/19443994.2014.937764.

Components