Particle-Based Multidimensional Multispecies Biofilm Model

In this paper we describe a spatially multidimensional (two-dimensional[2-D] and three-dimensional [3-D]) particle-based approach formodeling the dynamics of multispecies biofilms growing on multiplesubstrates. The model is based on diffusion-reaction mass balancesfor chemical species coupled with microbial growth and spreadingof biomass represented by hard spherical particles. Effectively,this is a scaled-up version of a previously proposed individual-basedbiofilm model. Predictions of this new particle-based modelwere quantitatively compared with those obtained with an establishedone-dimensional (1-D) multispecies model for equivalent problems.A nitrifying biofilm containing aerobic ammonium and nitriteoxidizers, anaerobic ammonium oxidizers, and inert biomass waschosen as an example. The 2-D and 3-D models generally gavethe same results. If only the average flux of nutrients needsto be known, 2-D and 1-D models are very similar. However, thebehavior of intermediates, which are produced and consumed indifferent locations within the biofilm, is better describedin 2-D and 3-D models because of the multidirectional concentrationgradients. The predictions of 2-D or 3-D models are also differentfrom those of 1-D models for slowly growing or minority speciesin the biofilm. This aspect is related to the mechanism of biomassspreading or advection implemented in the models and shouldreceive more attention in future experimental studies.

INTRODUCTION

Top

Introduction

Biofilm formation is determined by a variety of biological,chemical, and physical processes. During biofilm development,the large number of chemical and biological species that occurand interact simultaneously over a broad range of length andtime scales can easily confound human intuition. By generatingquantitative predictions from descriptions of biofilm characteristicsthat have been turned into the rational form of a complete setof equations, mathematical models can lead to deeper and broaderunderstanding. The appropriate level of mathematical sophisticationthat is chosen depends on the required accuracy and the availableresources. For practical engineering purposes, simple modelsare often sufficient to make safe design decisions. The firstbiofilm models were conceived with this goal. They describedbiofilms as uniform steady-state films containing a single typeof organism, in which substrate concentrations were governedexclusively by one-dimensional (1-D) mass transport and biochemicalreactions (30). Later, stratified models able to represent thedynamics of multisubstrate and multispecies biofilms were developed(33). This generation of 1-D models has helped enormously inunderstanding complex interactions in multispecies biofilms;the application of such models has even reached the stage wherethey are increasingly used as educational material in engineeringcurricula (e.g., AQUASIM software [29]). Yet 1-D models areunable to generate the characteristically multidimensional biofilmmorphology. Mathematical description of the more accurate pictureof biofilm structural heterogeneity that is emerging from numerousexperimental observations requires new biofilm models that providemore complex two-dimensional (2-D) and three-dimensional (3-D)descriptions of a microbial biofilm. Although different mathematicalmodels have recently been proposed for explaining formationof the geometric biofilm structure, the effects of physicalfactors, such as substrate gradients, on 3-D microbial speciesdistribution in biofilms has almost been neglected. Advancesin experimental techniques (confocal laser scanning microscopy,microelectrodes, etc.) have been paralleled by advances in computationalpower. However, we owe the rise of multidimensional biofilmmodels not only to more powerful hardware but also to newlydeveloped, highly efficient numerical methods.

In general, a quantitative representation of the system studiedis superior to a merely qualitative picture. This is one ofthe reasons for expressing hypotheses in mathematical form.A mathematical model consists of the full set of equations abstractingthe information required to simulate a system. A rigorous andmechanistic representation of the real system must be basedon the fundamental laws of physics, chemistry, and biology,which are also called first principles. The development of modelsbased on first principles enforces systematic and comprehensiveclarification of concepts and assumptions and imposes rationalmethods for approaching a problem. For example, biofilm modelsbased on reaction and transport principles (33) have provento be useful not only for testing the soundness of differentscientific concepts but also for establishing rational strategiesfor designing biofilm systems. Models based on first principlesprovide a unified view of microbial growth systems, includingbiofilms. They promote lateral transfer of insight between variousscientific domains. Diffusion-reaction models routinely usedin chemical engineering are now widely used to simulate biofilmsystems (33, 34). Fluid mechanics methods have been used tostudy biofilm rheology (5, 31) and hydrodynamic conditions inthe liquid environment surrounding a biofilm matrix (3, 4, 6,7, 23, 24, 26). Laws of structural mechanics and finite elementanalysis methods of civil engineering, have been used to studybiofilm growth and detachment (6, 26).

When a quantitative mathematical model of biofilm structureand function is constructed, it is advantageous to constructthe model from submodels, each of which describes one of thevarious ongoing processes in a biofilm, including (i) biomassgrowth and decay, (ii) biomass division and spreading, (iii)substrate transport and reactions, (iv) biomass detachment,(v) liquid flow past the biofilm, and (vi) biomass attachment.Attachment is an important process because it determines theinitial pattern of colonization of the substratum and the possibleimmigration of any type of cell from the liquid phase to variouslocations in the existing biofilm. The advantages of a modularbiofilm model are manifold and include a better understandingof specific phenomena, better validation of individual modelcomponents, the possibility of exchanging routines or submodelswith other biofilm models, more flexibility in solving decoupledmodel equations, and reflection of the modular structure ofbiofilm communities and processes.

For the processes described above, proper representation ofbiomass division and spreading is one of the most controversialtopics. The difficulty in modeling the spreading of microbialcells inside colonies is that there must be a mechanism to releasethe pressure generated by the growing bacteria. Different solutionshave been proposed, but given the lack of experimental evidence,none can claim to be correct. Current models for biofilm structuredeal with bacteria in two different ways, depending mostly onthe intended spatial scale of the model. One approach, individual-basedmodeling (IbM), attempts to model the biofilm community by describingthe actions and properties of individual bacteria (13, 14).IbM allows individual variability and treats bacterial cellsas the fundamental entities. Essential state variables are,for example, the cell mass (m), the cell volume (V), etc. Theother approach treats biofilms as multiphase systems and usesvolume averaging to develop macroscopic equations for biomassdynamics. These models can be called biomass-based models becausethey use the mass of cells per unit of volume (density or concentration[C_X]) as the state variable for biomass. A comprehensive analysisof conditions under which biomass averaging is a valid computationaltool has been performed previously (36, 37). The biomass-basedmodels can be further divided into two classes on the basisof the mechanism used for biomass spreading. One subclass includesdiscrete biofilm models (i.e., cellular automata), in whichbiomass can be shifted only stepwise along a finite number ofdirections according to a set of discrete rules (11, 12, 19,21, 22, 27, 35). The other subclass of biomass-based modelstreats biomass as a continuum, and biomass spreading is generallymodeled by differential equations widely used in physics (5,6, 8).

Perhaps the easiest way to model biofilm spreading is by implementinga discrete dynamic model, often called a cellular automaton.Here, discrete means that space, time, and properties of thesystem can have only a finite number of states. The space isfirst discretized into boxes (grid elements), forming a lattice.Typical for cellular automaton models is the requirement thatbiomass can move only in the finite number of lattice directions.A variety of discrete rules, although easy to implement, mayproduce qualitatively different and rather arbitrary structures.They might easily mislead the researcher to aesthetically drivenrather than physically motivated model formulation. Eberl etal. (8) pointed out a number of inherent drawbacks that discreteor stochastic mathematical models for biofilm spreading possess,as follows: (i) due to the small set of possible directionsfor movement, the models are not invariant with changes in thecoordinate system; (ii) a colony of bacteria growing in a uniformenvironment is not always round, as it should be; (ii) ad hocrules for priority must be specified for simultaneous shiftsinto the same grid cell; and (iv) many different and somewhatarbitrary rules can be formulated for spreading. In addition,we have seen from numerous computational experiments (14) that,when applied to multispecies biofilms, the cellular automaton-basedmodel described by Picioreanu et al. (21, 22) produces too muchinternal mixing (dispersion) of species within a colony. Thisis a direct consequence of the stochastic pushing rules usedfor biomass redistribution. Multispecies biofilm simulationswith a variant of the discrete algorithm that eliminates thisproblem tend to generate anisotropic colonies instead (25).

Many of the drawbacks of discrete spreading directions and discretedisplacement distances described above can be surmounted byallowing cell movement with a continuous set of directions anddistances. A realistic model of this kind was proposed previously(13) for bacterial colony growth, and it has been used recentlyfor simulation of a multispecies biofilm (14). Bacterial cellsare represented as hard spheres, with each cell having, besidesvariable volume and mass, a set of variable growth parameters.Each cell grows by consuming nutrients and divides when a certainvolume is reached. In this model, the pressure buildup due tothe increase in biomass is relaxed by minimizing the overlapof cells. The value of this IbM approach is strengthened bythe relative ease with which a mechanism for production andspreading of extracellular polymeric substances can be implemented(15). In an extension of the model described above (15, 16),how biofilm structure can be influenced by extracellular polymericsubstance production was demonstrated.

The IbM approach seems to be very appealing to microbiologistsbecause it allows individual variability and, by treating bacterialcells as the fundamental units, it can incorporate rare speciesor rare events, such as mutations. However, the very detailedlevel of biofilm description can also be a disadvantage whensystems with large-scale heterogeneity are modeled. It is evidentthat IbMs require more computer resources than biomass-basedmodels at the same spatial resolution require. However, a reconciliationbetween modeling scope and required resources is possible. Byallowing the presence of larger biomass particles (e.g., particlesthat are 10 to 20 µm in diameter), one can give up thecell level biomass representation while still keeping the shovingor pushing principle for biomass redistribution. Below, we callthis approach a particle-based modeling approach, which wasthe focus of this theoretical study.

The aims of this study were (i) to describe the particle-basedmodeling approach and solution of model equations and (ii) tocompare quantitatively and critically results of the new 2-Dand 3-D particle-based models with results for equivalent 1-Dproblems solved with established computational tools, such asAQUASIM (29).

MATERIALS AND METHODS

Top

Materials and Methods

Biofilm model description.
The new biofilm model, which is capable of describing microbialand solute distributions in two and three spatial dimensions,was constructed by combining ideas described previously (14,21, 22) and extending the spatial scale of the model by implementinga very efficient multigrid solver. In order to facilitate thecomparison with 1-D multispecies biofilm models (33, 34), abiofilm reactor system similar to that described by Reichert(29) was considered. The biofilm reactor contains two compartments:(i) bulk liquid and (ii) biofilm (Fig. 1A). The bulk liquidcompartment having a volume V (expressed in cubic meters) iscontinuously fed at a constant flow rate Q (in cubic metersper day) with a solution containing a number N_S of soluble substrates(n) at concentrations

(in grams per cubic meter). The bulk liquid is assumed to be completely mixed andcontinuously withdrawn from the reactor at the same flow rate(Q). The substrate concentration in the bulk liquid is

. The biofilm, with average thickness L_F, grows ona planar support with area A_F. The biofilm and bulk liquid compartmentsare in contact and exchange solutes only by diffusion. The bulkliquid volume (V) is very large compared with the biofilm volume(V_F), and thus V can be considered constant. For simplicity,biological activity and other reactions are not considered inthe bulk liquid. Biomass is not present in the inflow.

View larger version (62K):
[in this window]
[in a new window]
FIG. 1. Model biofilm system. (A) Continuously stirred tank reactor containing a liquid phase in contact with a biofilm growing on a planar surface. (B) Rectangular computational domain (2-D or 3-D) enclosing a small part of the whole biofilm. (C) Rectangular uniform grid (mesh size, h = ${Delta}$ x = ${Delta}$ y = ${Delta}$ z) used for solution of the partial differential equations for substrate diffusion and reaction. (D) Biofilm biomass contained in spherical particles holding one type of active biomass, as well as inactive biomass. All biomass particles situated in one grid element (i,j,k) experience the same substrate concentrations.

The computational domain is only a small part of the whole systemand includes both a small fraction of the bulk liquid and thebiofilm (Fig. 1B). The biofilm grows in a rectangular domain(box) with dimensions L_X by L_Y by L_Z on a planar, inert substratumwith the surface at z = 0. In 2-D, the size of the domain isL_X x L_Z. The processes that contribute to biofilm developmentthat are included in the model are as follows.

(i) Biomass growth and decay.
It is assumed that the biofilm consists of N_B types of activebiomass (with b = 1,2,...,N_B) and only one type of inert biomass(with b = 0) that results from the decay of the active biomass.The separation of active biomass into different categories isusually based on differences in metabolism. For example, thecase studies presented here include aerobic ammonium oxidizers,aerobic nitrite oxidizers, and anaerobic ammonium oxidizers.The 3-D biofilm structure is represented by a collection ofN_P nonoverlapping hard biomass spheres, also called biomassparticles. Each spherical particle (p) contains only one typeof active biomass and a fraction of inert biomass. The totalmass [] (in grams per particle) of particle p is therefore the sum of its active biomass [] and the inert biomass []:

(1)
Itis assumed that all particles have a constant density ( ${rho}$ _X) (ingrams of biomass per cubic meter of particle). When the biomassof the particle changes over time, the volume and radius [R^(p)]change accordingly. In 3-D this is:

(2)
Inthe 2-D model reduction the biomass particles are in fact cylindersspanning the whole domain (i.e., length [ ${Delta}$ x] and radius):

The growth of each biomass type (b) (active or inert)with time is described by an ordinary differential equationrepresenting the mass balance for each biomass particle:

(3)
The net reaction rates for generation of activebiomass (r_X,b) and inert biomass (r_X,0) are typically functionsof the active biomass of the particle [] and concentrations of various substrates present at the center(x,y,z) of the biomass sphere []. If a set of N_S relevant substrates was chosen, then n = 1,...,N_S.

An initial distribution of biomass particles at t = 0 must bedefined. It is assumed that all N_P,0 initial biomass particleshave mass m₀ and, correspondingly, radius R₀. They are randomlydistributed on the planar substratum surface at z = 0, and theircenters are placed at z = R₀. Therefore, the initial coordinatesof biomass particles are (x,y,R₀) with x [0,L_X] and y [0,L_Y].There are two main differences from the model described by Kreftet al. (14). First, the size of the biomass spheres can be chosento represent a cluster of similar cells (for example, a microcolonywith a diameter of $~$ 10 to 20 µm) rather than a single bacterialcell (typical size, $~$ 1 µm). Second, for the sake of simplicity,the metabolic (kinetic and stoichiometric) parameters are thesame for all biomass particles of the same type. Consequently,the biomass representation as hard spherical particles usedin this study can be considered a simplification and scale-upof the IbM approach (13, 14). This description of biomass asparticles can also be seen as an unconventional biomass discretizationscheme suitable for differential spreading of the various speciesin a multispecies biofilm. Biomass spreading in this model doesnot result in undue dispersal of clonal clusters; the extentof such mixing may affect competition and cooperation withinand between species in a biofilm. One can find analogies inother multiparticle models in theoretical physics and biology(2, 10).

(ii) Biomass division.
If there are sufficient nutrients in the environment, the bacterialmass in each particle grows according to equation 3. However,the total biomass (i.e., active biomass plus inert biomass)in a spherical biomass particle is assumed to be limited toa maximum value (M_X) (in grams of biomass per cubic meter ofparticle) (m_X < M_X) that is independent of the growth ratefor the sake of simplicity (13, 14). M_X is conveniently chosento obtain the desired total biomass density in the biofilm (C_X,max)(in grams of biomass per cubic meter of biofilm). When the M_Xin a sphere is reached, a new daughter sphere is created, whichtouches the mother sphere in a randomly chosen direction. Partof the biomass contained in the mother sphere is redistributedto the daughter sphere. There are two significant problems thatshould be mentioned here, which are related to mass redistributionand to the direction of division. First, if particles are splitinto precisely equal parts, then synchronized divisions occuruntil, due to spatial gradients of the substrate, biomass particlesgrow at different rates and the synchrony is gradually broken(13). Synchrony can be broken from the beginning by prescribingan unequal splitting ratio (e.g., by uniformly randomly redistributingbetween 40 and 60% of the mother particle). Second, if isotropicspreading of the bacterial colony is desired (producing a roundcluster of biomass particles growing at the same rate in alldirections), then the new particle must be placed at randomaround the dividing particle. Although here we consider onlythis coccoid type of growth, there are cases in which a preferentialgrowth direction prevails (for example, in the development ofthe filamentous bacteria that form activated sludge flocs).

(iii) Biomass spreading.
Biofilm spreading as a result of biomass growth occurs by shovingof the spheres when they get too close to each other. In thismodel, the pressure that builds up due to biomass growth isrelaxed by minimizing the overlap of spheres. For each particle,the vector sum of all overlap radii with neighboring particlesis calculated, and then the position of the particle is shiftedin the direction opposite the direction of this vector (similarto the algorithm described by Kreft et al. [13]): x_j,new = x_j– ${delta}$ x_j, where x_j is the coordinate of the current particlep (i.e., x₁ = x, x₂ = y, and x₃ = z). The three components ofthe displacement vector are ${delta}$ x_j (or explicitly, ${delta}$ x₁ = ${delta}$ x, ${delta}$ x₂ = ${delta}$ y, and ${delta}$ x₃ = ${delta}$ z):

(4)
where R^(p) is theradius of the current particle, R⁽ⁿ⁾ is the radius of a neighboringcell, is the coordinate of an overlapping neighboring particle, and d⁽ⁿ⁾ is the distance between particles:

(5)
The shove radius (k) is a multiplier on the particleradius that allows adjustment of the minimal spacing betweenparticles.

Treatments of domain boundaries constitute special cases ofbiomass spreading. If the particle's surface extends into thesubstratum (the particle's center is at z < R), then theparticle is shifted to z = R. Furthermore, undesired edge effectscan simply be avoided by not having edges, that is, by usingperiodic or wrapped-around boundaries in the directions parallelto the substratum (directions x and y). This means, for instance,that in the x direction, particles close to the domain boundaryat x = L_X can also have particles as neighbors at the boundaryopposite, x = 0. Consequently, particles shoved out at x = L_X+ ${delta}$ x reenter the computational domain at x = ${delta}$ x, and those pushedout to x = – ${delta}$ x are resited at x = L_X – ${delta}$ x.

(iv) Biomass detachment.
Detachment is implemented in this model only in a simplisticway, with the sole purpose of keeping the biofilm thicknessconstant in some of the example cases. The detachment modelsimply consists of removing every particle which is shiftedabove an imposed biofilm thickness limit (center z > L_F)due to a shoving step. The removed particle is not relocatedbut merely considered lost. Examples of biofilm systems growingin similar detachment conditions are the constant-depth filmfermentor described by Peters and Wimpenny (20) and, to someextent, the biofilm particles with a surface smoothed due tohigh abrasion rates in the biofilm airlift reactor during steadystate (17, 32).

(v) Mass balances of soluble substrates in the biofilm.
The spatial distribution (the field) of the concentrations ofa given set of substrates, products, etc. influences the growthrate of a particular biomass type. Conversely, the spatial distributionof bacterial activity affects the substrate concentration fields.Due to growth, division, spreading, and detachment, the spatialdistribution of biomass varies over time. This causes temporalvariation of the substrate fields as well. For this reason,equations of dynamic state diffusion-reaction mass balancesmust be written for each of the n substrates in the biofilmvolume:

(6)
where r_S,n is the net reactionrate and D_n is the effective diffusivity of substrate n in thebiofilm. The temporal derivative on the left is the rate ofsubstrate accumulation. The right side contains the isotropic3-D diffusion rate, with uniform diffusivity (a spatial variationof the diffusion coefficient, depending on, for instance, thelocal bacterial density, can easily be taken into account).The net reaction rate is the algebraic sum of the rates of allindividual processes that produce or consume the specific substrate.For example, if we consider a nitrifying biofilm with threemicrobial populations (see Results and Discussion), the netreaction rate for nitrite is the rate of nitrate productionby ammonia oxidizers minus the rates of consumption by nitriteoxidizers and by anammox bacteria. Reaction rates (r_S) frequentlydepend on multiple solute concentrations (C_S,n), and they regularlyinclude nonlinear saturation and inhibition factors.

The boundary conditions for equation 6 were as follows: (i)constant substrate concentrations at the top boundary (the biofilm-bulkliquid interface) (in order to simplify the analysis of modelresults, the external resistance to mass transfer was neglected;that is, the concentrations of solutes were assumed to be identicalto those in the bulk liquid at any point above the plane z =z_F, where z_F is the maximum biofilm height [Fig. 1B], whichincreases as the biofilm thickness increases over time)

(7)
(ii) zero flux of substrate through the impermeablesupport material (the biofilm-support interface)

(8)
and (iii) connected lateral boundaries (alsocalled periodic or cyclic boundary type)

(9)
Forthe initial state at t = 0, a uniform distribution of concentrationsthroughout the whole domain was assumed:

(10)

(vi) Mass balances of soluble substrates in the bulk liquid.
Mass balances for solutes in the bulk liquid are written fordynamic conditions:

(11)
What equation11 simply accounts for is that the accumulation of dissolvedmatter in the bulk volume (left side) is the result of (i) flowinto and out of the bioreactor (first term on the right side)and (ii) the global conversion rates in the biofilm volume (thetriple integral on the right side). The initial condition forequation 11 at t = 0 is:

(12)
Thesolution of equation 11 for the bulk concentrations [] is used as a boundary condition for equation 6.As an exception, the oxygen concentration in the bulk liquidcompartment was kept constant, simulating controlled aerationof the bioreactor.

Solution methods. (i) 1-D biofilm model.
The 1-D model for multispecies biofilms (33, 34) was used toevaluate the 2-D and 3-D particle-based biofilm model describedin this study. In particular, we used the 1-D model implementationknown as the computer program AQUASIM (29).

(ii) 2-D and 3-D biofilm model.
The dynamics of soluble substrate concentrations due to transportand reaction are very fast compared with the characteristictimes of biomass growth, division, and spreading. This observationallowed us to decouple the solution of biomass evolution fromthat of substrate diffusion-reaction mass balances, as explainedin more detail elsewhere (24). Due to this time separation assumption,the algorithm for the six model processes can be executed sequentially,as explained in the pseudocode below.

Set initial condition for biomass [N_P,0 particles with massm₀ and radius R₀]

Set initial condition for solutes in bulk liquid and biofilm(equations 10 and 12)

Biomass dynamics

Biomass growth and decay (equations 3and 2)
Biomass division
Biomass spreading (equation4)
Biomass detachment

Solute dynamics

Mass balances of solutes in the bulk liquid(equation11)
Mass balances of solutes in the biofilm(equations 6 to9 solved for steady state)

Time step, t $<-$ t + ${Delta}$ t

while t < t_end

First, a number of inoculum biomass particles are attached tothe planar support surface at randomly chosen positions. Valuesfor inlet concentrations of solutes are set in the whole computationaldomain, discretized in N_X x N_Y x N_Z cubic elements (or N_X xN_Z squares in 2-D) of size ${Delta}$ x = L_X/(N_X – 1) (Fig. 1C).Then, for biomass growth and decay, the biomass balances foreach particle are integrated directly for one time step ( ${Delta}$ t)to find . Division of biomass is executed for the particles whose new mass exceeds the threshold, i.e.,when . Following the growth and division steps, many particles overlap. Therefore, spreading accordingto the shoving mechanism described above is necessary. Aftersome of the particles are removed in the detachment step, thebiomass dynamics in the time interval ${Delta}$ t have been completed;that is, all new positions, masses, and particle radii at t+ ${Delta}$ t are known. With this new biomass distribution, we now needto compute the new fields of solute concentrations. Becausethe biomass particles are assumed to be circular and the gridelements of substrate fields are rectangular, the average ofthe bacterial masses (m_Xi) of all cells occupying a certaingrid element must be calculated prior to solution of the substratemass balances. The biomass concentration [] is used in the reaction rate expressions in the mass balances(equations 6 and 11). First, solute mass balances in bulk liquid(equation 11) are integrated with a simple forward Euler step.We use the dynamic form of equation 11 because attempts to couplea steady-state solute balance for the bulk liquid with the timeiterations for biomass and solute balances in the biofilm produceda numerically unstable scheme. The overall solute transformationrates in the biofilm (the triple integral in equation 11) areapproximated by summation. The same time step was used for bothbiomass and bulk solute balances. Second, a steady-state ( ${partial}$ C_S,n/ ${partial}$ t= 0) solution is found for the elliptic form of the partialdifferential equations (PDEs) for mass balances of soluble substratesin the biofilm (equation 6) with boundary conditions (equations7 to 9). The PDEs are solved by a nonlinear multigrid method(28). The multigrid algorithm has the advantage of simple, flexible,and very efficient solution of the systems of elliptic PDEsthat usually occur in spatial diffusion-reaction problems. Sincemultigrid algorithms scale linearly with system size (N), theyallow a dramatic increase in computational speed for 3-D applicationscompared with the alternating direction implicit algorithm usedin the earlier versions of the model, which scales with N² (13,14). The steady-state solution provides concentration fields[] that are used in the next time step for the solution of biomass balances (equation 3). The rate of biomassdetachment (per substratum area) at any moment t was calculatedby averaging over the previous 100 time steps (0.5 day).

Case studies. (i) Reaction model.
To assess the new biofilm modeling approach against a classic1-D biofilm model, a multispecies nitrifying biofilm communitywas chosen as a test case (9). This biofilm contains three activemicrobial types: aerobic ammonium oxidizers (NH), aerobic nitriteoxidizers (NO), and anaerobic ammonium oxidizers (AN; anammoxbacteria). In addition, the inert biomass is a fourth type ofparticulate material resulting from the decay of all activebiomass types. The reaction model includes the substrate dependenciesof the various species, which lead to competition among aerobicammonium oxidizers, aerobic nitrite oxidizers, and anammox bacteriafor oxygen, ammonium, and nitrite in the biofilm and productionof nitrate and N₂. Each of the three microbial types is capableof growth, as well as maintenance (endogenous respiration underaerobic and anaerobic conditions). All model parameters areshown in Table 1. The reaction stoichiometry and kinetics ofmodel processes are summarized in Table 2. The biofilm systemanalyzed in this study was chosen only to illustrate the model'scapability, and it does not represent any real system. Othermore realistic but more complex systems can be easily constructed,for example, by including heterotrophic bacteria, as shown byNoguera and Picioreanu (18).

View this table:
[in this window]
[in a new window]
TABLE 1. Model parameters^a

View this table:
[in this window]
[in a new window]
TABLE 2. Matrix of the stoichiometries and kinetics for the biological processes modeled^a

(ii) Test conditions.
Three growth regimens were used to evaluate the microbial distributionin biofilms as predicted by multidimensional models. Case Iis an oxygen-limited regimen achieved by a low oxygen concentrationin the reactor influent (2 g of O₂ m^–3) and a sufficientelectron donor concentration (30 g of N-NH₄ m^–3 and noother N compounds). In addition, case I* starts with a microbialdensity eight times lower than that in case I. Case II simulatesthe biofilm formation at high oxygen levels (10 g of O₂ m^–3)while maintaining the same nonlimiting ammonia concentration.In case III, the influent concentrations are approximately thosecalculated for the effluent in case II, simulating the behaviorof a second biofilm reactor connected in series with the firstreactor. In this case, the biofilm forms under nitrogen-limitingconditions [= 4 g of N-NH₄ m^–3 and= 3 g of N-NO₂ m^–3]. In all cases,the oxygen concentration in the bulk liquid compartment wasassumed to be controlled by aeration, maintaining the constantconcentration of the reactor influent.

RESULTS AND DISCUSSION

Top

Results and Discussion

The multidimensional biofilm model was evaluated with respectto 1-D models by means of several indicators, such as the spatialdistribution of substrate concentrations and fluxes, biofilmstructure, total biofilm biomass, dynamics of biomass distribution,dynamics of biomass detachment, dynamics of solute concentrationsin the bulk liquid, and steady-state concentrations of biomassand solutes in the biofilm. An example of the spatial distributionof substrate concentrations in 3-D is presented for case IIIin Fig. 2 and 3. The biofilm geometric structure and total biomassaccumulated in cases I and I* are shown in Fig. 4 and 5, respectively.The dynamics of biomass distribution, of biomass detachment,and of solute concentrations in the bulk are shown in Fig. 6(case I), Fig. 7 (case II), and Fig. 8 (case III). Figure 9shows the biomass spatial distributions for cases I to III in2-D simulations. The steady-state concentrations of biomassand solutes in the biofilm are shown in Fig. 10 (case I) andFig. 11 (case II). Animations of biofilm development can befound at the web page http://www.bt.tudelft.nl/content/ebt/researchgroup_ebt.html.

View larger version (97K):
[in this window]
[in a new window]
FIG. 2. Biomass distribution and multidimensional gradients of the intermediate nitrite, produced by ammonium oxidizers (blue) and consumed by nitrite oxidizers (red) and anammox bacteria (green), during biofilm development on day 5 (a and d), day 10 (b and e), and day 25 (cand f) under N-limiting conditions. For case III,

,

= 10 kg/m³,

= 4 kg/m³,

= 3 kg/m³, and

= 23 kg/m³. Panels d to f show nitrite concentrations (white, maximum; black, minimum). The top view maps show the concentrations at the bottom of the biofilm (z = 0); the lateral view maps are vertical sections along the lines indicated. The corresponding biomass distributions are shown in panels a to c.

View larger version (124K):
[in this window]
[in a new window]
FIG. 3. Oxygen (a) and nitrite (b) fluxes (indicated by the lengths and directions of the arrows) and 2-D concentration fields (darker contour areas indicate lower concentrations). The biomass particles are ammonium oxidizers (blue circles) and nitrite oxidizers (red circles).

View larger version (92K):
[in this window]
[in a new window]
FIG. 4. Spatial biomass distribution just before the onset of detachment when the thickness cutoff of 500 µm was reached. The biofilms grew under oxygen-limiting conditions (case I) from a high (a) or low (b and c) initial density. (a) 2-D model with high inoculum density

, day 10; (b) 2-D model with low inoculum density

, day 12 (arrows indicate the locations of ammonium-oxidizing inoculum particles); (c) 3-D model with low inoculum density

, day 14.

View larger version (28K):
[in this window]
[in a new window]
FIG. 5. Comparison of the development of areal biomass and bulk NH₄⁺ concentrations in simulations of the 1-D AQUASIM model (lines), as well as the new pseudo-1-D (*), 2-D ( ${circ}$ ), and 3-D (•) biofilm models. The biofilms grew under oxygen-limiting conditions[

= 10 kg/m³,

= 30 kg/m³] from a high initial density (case I;

) (a) or a low initial density (case I*;

) (b).

View larger version (29K):
[in this window]
[in a new window]
FIG. 6. Biofilm development under oxygen-limiting conditions[case I;

,

= 2 kg/m³,

= 30 kg/m³]. (a) Concentrations of dissolved chemical species in the bulk liquid; (b) areal biomass concentrations; (c) biomass detachment rates. In panels a and b, the lines indicate the 1-D AQUASIM results, while symbols indicate the 2-D biofilm model results. NH, aerobic ammonium oxidizers; NO, aerobic nitrite oxidizers; AN, anaerobic ammonium oxidizers; I, inert biomass; T, total biomass.

View larger version (31K):
[in this window]
[in a new window]
FIG. 7. Biofilm development with abundant oxygen and nitrogen[case II;

,

= 10 kg/m³,

= 30 kg/m³]. (a) Concentrations of dissolved chemical species in the bulk liquid; (b) areal biomass concentrations; (c) biomass detachment rates. In panels a and b, the lines indicate the 1-D AQUASIM results, while the symbols indicate the 2-D biofilm model results. NH, aerobic ammonium oxidizers; NO, aerobic nitrite oxidizers; AN, anaerobic ammonium oxidizers; I, inert biomass; T, total biomass.

View larger version (39K):
[in this window]
[in a new window]
FIG. 8. Biofilm development under N-limiting conditions [case III;

,

= 10 kg/m³,

= 4 kg/m³,

= 3 kg/m³,

= 23 kg/m³]. (a) Vertical profiles of biomass concentration: 1-D AQUASIM model (lines) and horizontally averaged values for the 2-D biofilm model (symbols). (b) Growth of the biofilm (areal biomass concentration) computed with the 1-D AQUASIM model (lines) and with the 2-D biofilm model (symbols). NH, aerobic ammonium oxidizers; NO, aerobic nitrite oxidizers; AN, anaerobic ammonium oxidizers; I, inert biomass; T, total biomass.

View larger version (87K):
[in this window]
[in a new window]
FIG. 9. Biomass spatial distributions for cases I to III in 2-D simulations: ammonium oxidizers (blue), nitrite oxidizers (red), and anammoxbacteria (green). The higher the fraction of inert biomass in a biomass particle, the lighter the color. (a) Case I.

= 24,

= 2 kg/m³,

= 30 kg/m³. (b) Case II.

= 24,

= 10 kg/m³,

= 30 kg/m³. (c) Case III.

= 3,

= 10 kg/m³,

= 4 kg/m³,

= 3 kg/m³,

= 23 kg/m³. Biomass distributions are shown for three phases of development, the onset of detachment (a), an intermediate state (b), and during the steady state (c).

View larger version (36K):
[in this window]
[in a new window]
FIG. 10. Vertical concentration profiles for biofilm growth underoxygen-limiting conditions [case I;

= 24,

= 2 kg/m³,

= 30 kg/m³]. (a and b) biomass; (c) dissolved chemical species. The lines indicate values computed with the 1-D AQUASIM model, while the symbols indicate values obtained with the 2-D biofilm model. NH, aerobic ammonium oxidizers; NO, aerobic nitrite oxidizers; AN, anaerobic ammonium oxidizers; I, inert biomass; T, total biomass.

View larger version (34K):
[in this window]
[in a new window]
FIG. 11. Vertical concentration profiles for biofilm growth withabundant oxygen and ammonia [case II;

= 24,

= 10kg/m³,

= 30 kg/m³]. (a and b) Biomass; (c) dissolved chemical-species. The lines indicate values computed with the 1-D AQUASIM model, while the symbols indicate values obtained with the 2-D biofilm model. NH, aerobic ammonium oxidizers; NO, aerobic nitrite oxidizers; AN, anaerobic ammonium oxidizers; I, inert biomass; T, total biomass.

Spatial distribution of substrate concentrations and fluxes.
An essential question is, why do we need spatial (e.g., 2-Dor 3-D) biofilm models? What is the relative importance of concentrationgradients in directions other than normal to the substratumsurface? This computational study and other studies of thiskind (7, 23) show that the answer is also a relative one. Dependingon the environmental conditions and biofilm structure, for somechemical species there are 3-D gradients, whereas for othersthe 1-D model description may be sufficient. A 3-D simulationwas carried out under the nitrogen-limiting conditions of caseIII. The spatial distribution of nitrifying bacteria at threemoments in time and the corresponding spatial distribution ofthe intermediate nitrite are shown in Fig. 2. There are obviouslyregions where NO₂^– accumulates due to production by ammoniumoxidizers and places where NO₂^– is depleted due to consumptionby nitrite oxidizers or anammox bacteria. Magnified representationsof the 2-D oxygen and nitrite fluxes and concentration fieldsin bacterial colonies are shown in Fig. 3. As a general observation,chemical species that are only consumed (e.g., oxygen and ammonium)or only produced (e.g., nitrogen) in the biofilm are more likelyto form only vertical, one-directional concentration gradients(Fig. 3a), unless the biofilm geometry is very porous and channeled.Conversely, in the concentration fields of intermediates (e.g.,nitrite), which may be produced and consumed in different places(e.g., by different microbial populations), strong multidirectionalgradients are likely to be found. This is clearly shown in Fig.3b, in which the nitrite flux is oriented along the biofilmsurface from the ammonium-oxidizing colonies to the nitrite-oxidizingcolonies.

Biofilm geometric structure.
Results of multidimensional models can be understood by lookingat the differences in geometric biofilm structure predictedfor different conditions (Fig. 4). As discussed previously (22),2-D and 3-D models based on the processes of diffusion, reaction,and growth suggest that formation of porous biofilm structuresis favored under nutrient-limiting conditions. This can be explainedas an unavoidable consequence of the existence of strong substrategradients, which lead to a much greater microbial growth rateat the top of the biofilm than at the bottom of the biofilm.When the initial microbial distribution on the surface is nonuniform,a wavy biofilm surface develops due to a self-enhancing process;this biofilm surface instability was also studied from a mathematicalperspective by Dockery and Klapper (5). Intrinsically, planar1-D biofilm models assume a planar biofilm surface and concentrationboundary layers, as well as 1-D substrate gradients and fluxes.Therefore, an unstable biofilm surface is impossible, and acompact layer with greater total biomass per area develops dueto the absence of substrate-depleted valleys.

A sparse initial biomass distribution, as shown previously (24),leads to rugged biofilms with deep vertical channels. This iswhy the biofilm shown in Fig. 4b, corresponding to case I*,which was inoculated with three spheres of each biomass type,was more porous than the biofilm shown in Fig. 4a, correspondingto case I, which was inoculated with 24 spheres of each type.

Biofilm total biomass.
The total biomass accumulation predicted by the 1-D, 2-D, and3-D models is an important indicator of biofilm structure. Themost direct comparison between the growth curves for the 1-DAQUASIM models and the 2-D and 3-D models can be made underoxygen-limiting conditions because the oxygen gradients areessentially 1-D (case I) (Fig. 5a). If the biofilm thicknessis not leveled off by detachment, the multidimensional modelspredict less biomass than the 1-D models predict. Reduced biomassformation is shown in Fig. 5, which shows that the total biomasswas present only until the biofilm thickness reached 500 µm,the critical thickness for the onset of detachment. When theexperiment was started with less biomass on the carrier (caseI*), not only was the amount of biomass that accumulated lowerthan the amount that accumulated in case I (and lower than 1-Dpredictions), but also there was a significant difference betweenthe 2-D and 3-D model results (Fig. 5b). The total biomassespredicted by 1-D, 2-D, and 3-D models do match when compactstructures with a flat, horizontal surface form, as they doin the absence of nutrient limitation (case II) (Fig. 7b, inset).

For narrow 2-D systems, in which the length (L_X) is so smallthat surface instability cannot develop, a pseudo-1-D systemis obtained as a degenerated case of the multidimensional system.The total biomass and bulk concentrations obtained with sucha pseudo-1-D model (Fig. 5) are approximately the same as thosecalculated with 1-D AQUASIM (Fig. 5).

Dynamics of biomass distribution.
The development or ecological succession of the total amountof each biomass type attached to the substratum was recordedfor all three test cases. The areal biomass concentrations (ingrams of chemical oxygen demand [COD] biomass per square meter)were calculated by integration of biomass concentration profilesfor the 1-D biofilm:

(13)
and by summationof all particulate biomasses for the 2-D and 3-D models foreach b = 0, 1, ..., N_B

(14)
Initially,the areal biomass concentrations computed with 1-D and 2-D modelsmatch very well (Fig. 6b, 7b, and 8b). However, due to porousbiofilm formation in nutrient-limited case I (Fig. 6b), lessbiomass accumulates in the 2-D model until the detachment begins.During the detachment, the total biomass concentration staysclose to 5 g of COD m^–2, but the ratio of the populationschanges, progressing to a steady state. The fast-growing population(in this case the ammonium consumers) dominates initially. Thethick layer of aerobic bacteria then creates conditions favorablefor development of the anaerobic population of anammox bacteria(namely, very low O₂ concentrations and very high NH₄⁺ and NO₂^–concentrations). Gradually, the anammox bacteria displace theammonium oxidizers, as shown in Fig. 9a. After 57 days, theanaerobic ammonium oxidizers dominate. Note that this is a typicalexample of ecological succession, where the fast-growing pioneeror opportunist population is replaced by a more slowly growinggleaner population that can sustain a higher standing stockdue to, for example, more efficient use of the same resourcesor changes in resource use (1). Here, the ammonia oxidizers,which use the increasingly limiting oxygen, are replaced bythe anammox bacteria, which use nitrite as the alternative resource.The aerobic NO₂^– oxidizers are completely displaced fromthe biofilm (Fig. 9) because they have to compete on two sides(for oxygen with the aerobic ammonium oxidizers and for nitritewith the anammox bacteria). Good agreement between the biomassratios of 1-D and 2-D models is shown in Fig. 6b. However, moreinert biomass accumulated in the 2-D particle-based model simulationsthan in the 1-D continuum-based model simulations.

At a high O₂ concentration (case II; 10 g m^–3), the anaerobiclayer is much thinner and deeper in the biofilm. Consequently,the anaerobic ammonium oxidizers develop slowly, and the aerobicnitrite oxidizers occupy an important fraction of the biofilm.The initial fast growth of NH₄⁺ oxidizers and the catching upby the NO₂^– oxidizers agree well in 1-D and 2-D modelsuntil about 100 days (Fig. 7b). Again, at later stages, formationof inert biomass is much more pronounced in the 2-D model. The1-D AQUASIM model predicts a steady-state biofilm with verysimilar proportions of the three populations (Table 3, totalbiomass in case II), and there is a trend towards a diffusethree-layer distribution with NH₄⁺ oxidizers on the top, NO₂^–oxidizers in the middle, and anammox bacteria on the bottom.However, in the 2-D particle-based model, much of the anammoxbacterial layer is replaced by inert biomass (Fig. 9b), becausethe inert biomass is less mobile in the particle-based model(i.e., it has a lower advective velocity) than in the 1-D continuummodel, which reduces the flow of inert biomass to the sheared-offbiofilm top. Because of the particle spreading mechanism, theobserved trend in all simulations was that the slowest and thereforedeepest microbial species were not washed out of the biofilmin the 2-D simulations as quickly as they were in the 1-D simulations.

View this table:
[in this window]
[in a new window]
TABLE 3. Comparison of the 1-D multispecies AQUASIM biofilm model and the new 2-D multispecies biofilm model in terms of steady-state predictions for the biofilm development indicators areal biomass concentration [X_b] and biomass detachment rate [J_X,b^(det)]

In N-limited conditions (case III), the whole biofilm developmentis slowed down. As a result, the maximum thickness (500 µm)and maximum biomass (5 g of COD m^–2) are reached muchlater, only after 200 days. More importantly, the largest biofilmfraction is composed of inert material (Fig. 8b). Active nitrifierssurvive only in the top biofilm layers, within the penetrationdepth of NH₄⁺ and NO₂^– (Fig. 9c), and the bottom biofilmlayer is inactive and frozen. This lack of growth and upwardpushing of biomass eliminates the reason for the differencesbetween the 1-D and 2-D model results in this case (Fig. 8b).

Dynamics of solute concentrations.
Generally, the substrate concentrations in the reactor bulk[] relax faster to a steady state than the total biomass concentrations relax (Fig. 7a). This is a consequenceof the fact that the are determined by fluxes to or from the biofilm, which in turn depend mainly on the biomassdistribution in the top biofilm layers. Because of the highersubstrate concentrations near the biofilm surface, the biomassgrowth rates are also higher; thus, a stable biomass profilenear the surface is established sooner.

Dynamics of biomass detachment.
By imposing a maximum biofilm thickness, detachment of biomassparticles from the biofilm surface occurs continuously afterz_F = 500 µm has been reached. This results in a gradualfilling of biofilm pores with time and, eventually, developmentof a compact biofilm. The biomass detachment rate during a timestep ( ${Delta}$ t) is, in the 2-D model, related to the carrier area (L_X ${Delta}$ x)and results from the removal of N_det,t biomass particles:

(15)
The detachment rate increases during the porousbiofilm stage, as the surface area in contact with the z_F thresholdalso increases and more biomass particles are pushed out ofthe upper biofilm boundary (Fig. 6c). The rate of total biomassdetachment indicates the rate of overall biomass growth in the biofilm. As expected, under fully aerobicconditions (case II), the flux of detached biomass is higher(0.8 g of COD m^–2day^–1) (Fig. 7c) than it is inthe low-oxygen setting (case I; 0.6 g of COD m^–2day^–1)(Fig. 6c). The fast-growing population of aerobic NH₄⁺ oxidizersgrows preferentially in the top biofilm layer and, consequently,has a high detachment rate. The slowly growing bacteria accumulatingin deeper layers must push through the layer of aerobic NH₄⁺oxidizers first and, therefore, begin to detach later (anaerobicNH₄⁺ oxidizers begin to detach 60 days after aerobic NH₄⁺ oxidizers[Fig. 6c], and aerobic nitrite oxidizers begin to detach 5 daysafter aerobic NH₄⁺ oxidizers [Fig. 7c]). Due to the small sizeof the model biofilm system, noisy oscillations in detachmentrates occur due to the synchronized removal of discrete particles.

Steady-state concentrations and fluxes.
A perfect steady state with respect to the spatial distributionof biomass and solutes is not reached in the particle-basedmultidimensional biofilm model. This is due to the stochasticcharacteristics of the model, as the algorithm uses random numberson several occasions, including (i) for the initial distributionof biomass particles on the surface, (ii) for the random choiceof direction for the placement of daughter particles, and (iii)for uneven cell division, which is a desynchronizing divisionevent in the biofilm.

Four variables that have practical importance, substrate concentration, substrate fluxes (J_S,i), total biomassin the biofilm (X_b), and detachment rates , were selected for comparison of the 1-D multispecies model withthe 2-D particle-based model at steady state (Tables 3 and 4).Solute concentrations in the bulk liquid were in general verysimilar in the two models, while a greater difference was observedfor nitrite and nitrate in oxygen-limited case I. Profiles forsubstrate concentrations in the biofilm [C_S,i(z)] are shownin Fig. 10c (case I) and Fig. 11c (case II). The biomass distributionin the biofilm, C_X,b(z), displays some differences (Fig. 8a,10, and 11). In general, the agreement among biomass profilesis better for the upper biofilm layers, where, due to the presenceof higher substrate concentrations, the biomass growth ratesand consequently the advective velocities are higher. The substrateprofiles computed for the two models agree very well, indicatingthat they are not very sensitive to the biomass distributionin the biofilm depth. As a consequence, the substrate fluxesto the biofilm (Table 4) are very similar. Finally, the goodmatch in biomass profiles near the biofilm surface is reflectedin good agreement for biomass detachment rates at a steady state(Table 3).

View this table:
[in this window]
[in a new window]
TABLE 4. Comparison of the 1-D multispecies AQUASIM biofilm model and the new 2-D multispecies biolfilm model in terms of steady-state predictions for the biofilm development indicators bulk substrate concentration [C_S,i^(b)] and solute flux [J_S,i]

These observations suggest that complete knowledge of the 2-Dor 3-D distribution of microbial species is probably not neededto calculate overall conversion rates (where the interest isusually from an engineering perspective), but species localizationmight have a great effect on inherently local microbial competition,cooperation, and selection (where the interest of microbiologistsoften is).

Conclusions.
In this paper we describe a spatially multidimensional (2-Dand 3-D) particle-based approach for dynamic modeling of multispeciesbiofilms growing on multiple substrates. Predictions of thenew 2-D and 3-D particle-based models were quantitatively andcritically compared with results for equivalent 1-D problemsobtained with the classic, 1-D multispecies model (33, 34) implementedin the AQUASIM software (29). For comparison, biofilms containingthree active microbial populations (aerobic ammonium oxidizers,aerobic nitrite oxidizers, and anaerobic ammonium oxidizers),one type of inert biomass, and five soluble chemical species(oxygen, ammonium, nitrite, nitrate, and N₂) were studied underthree growth regimens.

The multidimensional models predicted less biomass per carrierarea (i.e., more porous biofilms), especially under nutrient-limitedconditions and when biofilm development started from a spatiallysparse inoculum.

Chemical species that are either consumed or produced in geometricallyhomogeneous biofilms mainly form only one-directional gradientsof concentration. For intermediates, which are produced andconsumed in different locations in the biofilm, strong multidirectionalgradients are likely to be found. These gradients can be predictedaccurately only by 2-D and 3-D models.

Comparison with the 1-D predictions showed that full knowledgeof the 2-D or 3-D distribution of microbial species may notbe required for calculation of solute fluxes and overall conversionrates. However, from an ecological and evolutionary point ofview, the distribution of species might have great effects onlocal microbial competition and cooperation, as well as selectionprocesses.

Especially for the slowest growing species and for the inertbiomaterial the 1-D and 2-D models predict different dynamicsand depth profiles. This is mainly caused by differences inbiomass spreading mechanisms (overall advective velocity throughthe entire 1-D slice versus local pushing of hard spheres ina 2-D or 3-D space), combined with different coupling of inertbiomass to active biomass. Experimental observations of thespreading of color-tagged bacteria through growing biofilmsshould help improve this aspect of biofilm models.

It is believed that IbM of 3-D biofilm structure and the particle-basedapproach as its scaled-up version, although computationallymore demanding than 1-D models, go one step further towardsmodeling biofilm systems from first principles. Not only physicaland chemical processes, such as substrate transport and reactions,but also biomass behavior can be described in a mechanisticallymeaningful way as the result of small-scale actions and interactions,such as individual microbial growth, division, and movement.By its extendible nature, this new generation of biofilm modelscould explain a large variety of phenomena from microbiologyand microbial ecology, retaining at the same time the quantitativefeatures required in engineering applications.

FOOTNOTES

* Corresponding author. Mailing address: Department of Biotechnology, Delft University of Technology, Julianalaan 67, 2628 BC Delft, The Netherlands. Phone: 31 15 2781551. Fax: 31 15 2782355. E-mail: C.Picioreanu{at}tnw.tudelft.nl.

REFERENCES

Top