**DOI:**10.1128/AEM.70.5.3024-3040.2004

## ABSTRACT

In this paper we describe a spatially multidimensional (two-dimensional [2-D] and three-dimensional [3-D]) particle-based approach for modeling the dynamics of multispecies biofilms growing on multiple substrates. The model is based on diffusion-reaction mass balances for chemical species coupled with microbial growth and spreading of biomass represented by hard spherical particles. Effectively, this is a scaled-up version of a previously proposed individual-based biofilm model. Predictions of this new particle-based model were quantitatively compared with those obtained with an established one-dimensional (1-D) multispecies model for equivalent problems. A nitrifying biofilm containing aerobic ammonium and nitrite oxidizers, anaerobic ammonium oxidizers, and inert biomass was chosen as an example. The 2-D and 3-D models generally gave the same results. If only the average flux of nutrients needs to be known, 2-D and 1-D models are very similar. However, the behavior of intermediates, which are produced and consumed in different locations within the biofilm, is better described in 2-D and 3-D models because of the multidirectional concentration gradients. The predictions of 2-D or 3-D models are also different from those of 1-D models for slowly growing or minority species in the biofilm. This aspect is related to the mechanism of biomass spreading or advection implemented in the models and should receive more attention in future experimental studies.

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 occur and interact simultaneously over a broad range of length and time scales can easily confound human intuition. By generating quantitative predictions from descriptions of biofilm characteristics that have been turned into the rational form of a complete set of equations, mathematical models can lead to deeper and broader understanding. The appropriate level of mathematical sophistication that is chosen depends on the required accuracy and the available resources. For practical engineering purposes, simple models are often sufficient to make safe design decisions. The first biofilm models were conceived with this goal. They described biofilms as uniform steady-state films containing a single type of organism, in which substrate concentrations were governed exclusively by one-dimensional (1-D) mass transport and biochemical reactions (30). Later, stratified models able to represent the dynamics of multisubstrate and multispecies biofilms were developed (33). This generation of 1-D models has helped enormously in understanding complex interactions in multispecies biofilms; the application of such models has even reached the stage where they are increasingly used as educational material in engineering curricula (e.g., AQUASIM software [29]). Yet 1-D models are unable to generate the characteristically multidimensional biofilm morphology. Mathematical description of the more accurate picture of biofilm structural heterogeneity that is emerging from numerous experimental observations requires new biofilm models that provide more complex two-dimensional (2-D) and three-dimensional (3-D) descriptions of a microbial biofilm. Although different mathematical models have recently been proposed for explaining formation of the geometric biofilm structure, the effects of physical factors, such as substrate gradients, on 3-D microbial species distribution in biofilms has almost been neglected. Advances in experimental techniques (confocal laser scanning microscopy, microelectrodes, etc.) have been paralleled by advances in computational power. However, we owe the rise of multidimensional biofilm models not only to more powerful hardware but also to newly developed, highly efficient numerical methods.

In general, a quantitative representation of the system studied is superior to a merely qualitative picture. This is one of the reasons for expressing hypotheses in mathematical form. A mathematical model consists of the full set of equations abstracting the information required to simulate a system. A rigorous and mechanistic representation of the real system must be based on the fundamental laws of physics, chemistry, and biology, which are also called first principles. The development of models based on first principles enforces systematic and comprehensive clarification of concepts and assumptions and imposes rational methods for approaching a problem. For example, biofilm models based on reaction and transport principles (33) have proven to be useful not only for testing the soundness of different scientific concepts but also for establishing rational strategies for designing biofilm systems. Models based on first principles provide a unified view of microbial growth systems, including biofilms. They promote lateral transfer of insight between various scientific domains. Diffusion-reaction models routinely used in chemical engineering are now widely used to simulate biofilm systems (33, 34). Fluid mechanics methods have been used to study biofilm rheology (5, 31) and hydrodynamic conditions in the liquid environment surrounding a biofilm matrix (3, 4, 6, 7, 23, 24, 26). Laws of structural mechanics and finite element analysis methods of civil engineering, have been used to study biofilm growth and detachment (6, 26).

When a quantitative mathematical model of biofilm structure and function is constructed, it is advantageous to construct the model from submodels, each of which describes one of the various ongoing processes in a biofilm, including (i) biomass growth 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 the initial pattern of colonization of the substratum and the possible immigration of any type of cell from the liquid phase to various locations in the existing biofilm. The advantages of a modular biofilm model are manifold and include a better understanding of specific phenomena, better validation of individual model components, the possibility of exchanging routines or submodels with other biofilm models, more flexibility in solving decoupled model equations, and reflection of the modular structure of biofilm communities and processes.

For the processes described above, proper representation of biomass division and spreading is one of the most controversial topics. The difficulty in modeling the spreading of microbial cells inside colonies is that there must be a mechanism to release the pressure generated by the growing bacteria. Different solutions have been proposed, but given the lack of experimental evidence, none can claim to be correct. Current models for biofilm structure deal with bacteria in two different ways, depending mostly on the intended spatial scale of the model. One approach, individual-based modeling (IbM), attempts to model the biofilm community by describing the actions and properties of individual bacteria (13, 14). IbM allows individual variability and treats bacterial cells as the fundamental entities. Essential state variables are, for example, the cell mass (*m*), the cell volume (*V*), etc. The other approach treats biofilms as multiphase systems and uses volume averaging to develop macroscopic equations for biomass dynamics. These models can be called biomass-based models because they use the mass of cells per unit of volume (density or concentration [*C _{X}*]) as the state variable for biomass. A comprehensive analysis of conditions under which biomass averaging is a valid computational tool has been performed previously (36, 37). The biomass-based models can be further divided into two classes on the basis of the mechanism used for biomass spreading. One subclass includes discrete biofilm models (i.e., cellular automata), in which biomass can be shifted only stepwise along a finite number of directions according to a set of discrete rules (11, 12, 19, 21, 22, 27, 35). The other subclass of biomass-based models treats biomass as a continuum, and biomass spreading is generally modeled by differential equations widely used in physics (5, 6, 8).

Perhaps the easiest way to model biofilm spreading is by implementing a discrete dynamic model, often called a cellular automaton. Here, discrete means that space, time, and properties of the system can have only a finite number of states. The space is first discretized into boxes (grid elements), forming a lattice. Typical for cellular automaton models is the requirement that biomass can move only in the finite number of lattice directions. A variety of discrete rules, although easy to implement, may produce qualitatively different and rather arbitrary structures. They might easily mislead the researcher to aesthetically driven rather than physically motivated model formulation. Eberl et al. (8) pointed out a number of inherent drawbacks that discrete or stochastic mathematical models for biofilm spreading possess, as follows: (i) due to the small set of possible directions for movement, the models are not invariant with changes in the coordinate system; (ii) a colony of bacteria growing in a uniform environment is not always round, as it should be; (ii) ad hoc rules for priority must be specified for simultaneous shifts into the same grid cell; and (iv) many different and somewhat arbitrary 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-based model described by Picioreanu et al. (21, 22) produces too much internal mixing (dispersion) of species within a colony. This is a direct consequence of the stochastic pushing rules used for biomass redistribution. Multispecies biofilm simulations with a variant of the discrete algorithm that eliminates this problem tend to generate anisotropic colonies instead (25).

Many of the drawbacks of discrete spreading directions and discrete displacement distances described above can be surmounted by allowing cell movement with a continuous set of directions and distances. A realistic model of this kind was proposed previously (13) for bacterial colony growth, and it has been used recently for simulation of a multispecies biofilm (14). Bacterial cells are represented as hard spheres, with each cell having, besides variable volume and mass, a set of variable growth parameters. Each cell grows by consuming nutrients and divides when a certain volume is reached. In this model, the pressure buildup due to the increase in biomass is relaxed by minimizing the overlap of cells. The value of this IbM approach is strengthened by the relative ease with which a mechanism for production and spreading 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 polymeric substance production was demonstrated.

The IbM approach seems to be very appealing to microbiologists because it allows individual variability and, by treating bacterial cells as the fundamental units, it can incorporate rare species or rare events, such as mutations. However, the very detailed level of biofilm description can also be a disadvantage when systems with large-scale heterogeneity are modeled. It is evident that IbMs require more computer resources than biomass-based models at the same spatial resolution require. However, a reconciliation between modeling scope and required resources is possible. By allowing the presence of larger biomass particles (e.g., particles that are 10 to 20 μm in diameter), one can give up the cell level biomass representation while still keeping the shoving or pushing principle for biomass redistribution. Below, we call this approach a particle-based modeling approach, which was the focus of this theoretical study.

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

## MATERIALS AND METHODS

Biofilm model description.The new biofilm model, which is capable of describing microbial and 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 implementing a very efficient multigrid solver. In order to facilitate the comparison with 1-D multispecies biofilm models (33, 34), a biofilm 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 liquid compartment having a volume *V* (expressed in cubic meters) is continuously fed at a constant flow rate *Q* (in cubic meters per 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 and continuously withdrawn from the reactor at the same flow rate (

*Q*). The substrate concentration in the bulk liquid is . The biofilm, with average thickness

*L*, grows on a planar support with area

_{F}*A*. The biofilm and bulk liquid compartments are in contact and exchange solutes only by diffusion. The bulk liquid volume (

_{F}*V*) is very large compared with the biofilm volume (

*V*), and thus

_{F}*V*can be considered constant. For simplicity, biological activity and other reactions are not considered in the bulk liquid. Biomass is not present in the inflow.

The computational domain is only a small part of the whole system and includes both a small fraction of the bulk liquid and the biofilm (Fig. 1B). The biofilm grows in a rectangular domain (box) with dimensions *L _{X}* by

*L*by

_{Y}*L*on a planar, inert substratum with the surface at

_{Z}*z*= 0. In 2-D, the size of the domain is

*L*×

_{X}*L*. The processes that contribute to biofilm development that are included in the model are as follows.

_{Z}(i) Biomass growth and decay.It is assumed that the biofilm consists of *N _{B}* types of active biomass (with

*b*= 1,2,…,

*N*) and only one type of inert biomass (with

_{B}*b*= 0) that results from the decay of the active biomass. The separation of active biomass into different categories is usually based on differences in metabolism. For example, the case studies presented here include aerobic ammonium oxidizers, aerobic nitrite oxidizers, and anaerobic ammonium oxidizers. The 3-D biofilm structure is represented by a collection of

*N*nonoverlapping hard biomass spheres, also called biomass particles. Each spherical particle (

_{P}*p*) contains only one type of active biomass and a fraction of inert biomass. The total mass [] (in grams per particle) of particle

*p*is therefore the sum of its active biomass [] and the inert biomass []: (1) It is assumed that all particles have a constant density (ρ

*) (in grams of biomass per cubic meter of particle). When the biomass of the particle changes over time, the volume and radius [*

_{X}*R*

^{(p)}] change accordingly. In 3-D this is: (2) In the 2-D model reduction the biomass particles are in fact cylinders spanning the whole domain (i.e., length [Δ

*x*] and radius): The growth of each biomass type (

*b*) (active or inert) with time is described by an ordinary differential equation representing the mass balance for each biomass particle: (3) The net reaction rates for generation of active biomass (

*r*) and inert biomass (

_{X,b}*r*

_{X,0}) are typically functions of 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*relevant substrates was chosen, then

_{S}*n*= 1,…,

*N*.

_{S}An initial distribution of biomass particles at *t* = 0 must be defined. It is assumed that all *N*_{P,0} initial biomass particles have mass *m*_{0} and, correspondingly, radius *R*_{0}. They are randomly distributed on the planar substratum surface at *z* = 0, and their centers are placed at *z* = *R*_{0}. Therefore, the initial coordinates of biomass particles are (*x*,*y*,*R*_{0}) with *x* ∈ [0,*L _{X}*] and

*y*∈ [0,

*L*]. There are two main differences from the model described by Kreft et al. (14). First, the size of the biomass spheres can be chosen to represent a cluster of similar cells (for example, a microcolony with a diameter of ∼10 to 20 μm) rather than a single bacterial cell (typical size, ∼1 μm). Second, for the sake of simplicity, the metabolic (kinetic and stoichiometric) parameters are the same for all biomass particles of the same type. Consequently, the biomass representation as hard spherical particles used in this study can be considered a simplification and scale-up of the IbM approach (13, 14). This description of biomass as particles can also be seen as an unconventional biomass discretization scheme suitable for differential spreading of the various species in a multispecies biofilm. Biomass spreading in this model does not result in undue dispersal of clonal clusters; the extent of such mixing may affect competition and cooperation within and between species in a biofilm. One can find analogies in other multiparticle models in theoretical physics and biology (2, 10).

_{Y}(ii) Biomass division.If there are sufficient nutrients in the environment, the bacterial mass 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 to a maximum value (*M _{X}*) (in grams of biomass per cubic meter of particle) (

*m*<

_{X}*M*) that is independent of the growth rate for the sake of simplicity (13, 14).

_{X}*M*is conveniently chosen to obtain the desired total biomass density in the biofilm (

_{X}*C*

_{X,max}) (in grams of biomass per cubic meter of biofilm). When the

*M*in a sphere is reached, a new daughter sphere is created, which touches the mother sphere in a randomly chosen direction. Part of the biomass contained in the mother sphere is redistributed to the daughter sphere. There are two significant problems that should be mentioned here, which are related to mass redistribution and to the direction of division. First, if particles are split into precisely equal parts, then synchronized divisions occur until, due to spatial gradients of the substrate, biomass particles grow at different rates and the synchrony is gradually broken (13). Synchrony can be broken from the beginning by prescribing an unequal splitting ratio (e.g., by uniformly randomly redistributing between 40 and 60% of the mother particle). Second, if isotropic spreading of the bacterial colony is desired (producing a round cluster of biomass particles growing at the same rate in all directions), then the new particle must be placed at random around the dividing particle. Although here we consider only this coccoid type of growth, there are cases in which a preferential growth direction prevails (for example, in the development of the filamentous bacteria that form activated sludge flocs).

_{X}(iii) Biomass spreading.Biofilm spreading as a result of biomass growth occurs by shoving of the spheres when they get too close to each other. In this model, the pressure that builds up due to biomass growth is relaxed by minimizing the overlap of spheres. For each particle, the vector sum of all overlap radii with neighboring particles is calculated, and then the position of the particle is shifted in the direction opposite the direction of this vector (similar to the algorithm described by Kreft et al. [13]): *x*_{j,new} = *x _{j}* − δ

*x*, where

_{j}*x*is the coordinate of the current particle

_{j}*p*(i.e.,

*x*=

_{1}*x*,

*x*=

_{2}*y*, and

*x*=

_{3}*z*). The three components of the displacement vector are δ

*x*(or explicitly, δ

_{j}*x*= δ

_{1}*x*, δ

*x*= δ

_{2}*y*, and δ

*x*= δ

_{3}*z*): (4) where

*R*

^{(p)}is the radius of the current particle,

*R*

^{(n)}is the radius of a neighboring cell, is the coordinate of an overlapping neighboring particle, and

*d*

^{(n)}is the distance between particles: (5) The shove radius (

*k*) is a multiplier on the particle radius that allows adjustment of the minimal spacing between particles.

Treatments of domain boundaries constitute special cases of biomass spreading. If the particle's surface extends into the substratum (the particle's center is at *z* < *R*), then the particle is shifted to *z* = *R*. Furthermore, undesired edge effects can simply be avoided by not having edges, that is, by using periodic or wrapped-around boundaries in the directions parallel to the substratum (directions *x* and *y*). This means, for instance, that in the *x* direction, particles close to the domain boundary at *x* = *L _{X}* can also have particles as neighbors at the boundary opposite,

*x*= 0. Consequently, particles shoved out at

*x*=

*L*+ δ

_{X}*x*reenter the computational domain at

*x*= δ

*x*, and those pushed out to

*x*= −δ

*x*are resited at

*x*=

*L*− δ

_{X}*x*.

(iv) Biomass detachment.Detachment is implemented in this model only in a simplistic way, with the sole purpose of keeping the biofilm thickness constant in some of the example cases. The detachment model simply consists of removing every particle which is shifted above an imposed biofilm thickness limit (center *z* > *L _{F}*) due to a shoving step. The removed particle is not relocated but merely considered lost. Examples of biofilm systems growing in similar detachment conditions are the constant-depth film fermentor described by Peters and Wimpenny (20) and, to some extent, the biofilm particles with a surface smoothed due to high abrasion rates in the biofilm airlift reactor during steady state (17, 32).

(v) Mass balances of soluble substrates in the biofilm.The spatial distribution (the field) of the concentrations of a given set of substrates, products, etc. influences the growth rate of a particular biomass type. Conversely, the spatial distribution of bacterial activity affects the substrate concentration fields. Due to growth, division, spreading, and detachment, the spatial distribution of biomass varies over time. This causes temporal variation of the substrate fields as well. For this reason, equations of dynamic state diffusion-reaction mass balances must be written for each of the *n* substrates in the biofilm volume:
(6)
where *r*_{S,n} is the net reaction rate and *D _{n}* is the effective diffusivity of substrate

*n*in the biofilm. The temporal derivative on the left is the rate of substrate accumulation. The right side contains the isotropic 3-D diffusion rate, with uniform diffusivity (a spatial variation of the diffusion coefficient, depending on, for instance, the local bacterial density, can easily be taken into account). The net reaction rate is the algebraic sum of the rates of all individual processes that produce or consume the specific substrate. For example, if we consider a nitrifying biofilm with three microbial populations (see Results and Discussion), the net reaction rate for nitrite is the rate of nitrate production by ammonia oxidizers minus the rates of consumption by nitrite oxidizers and by anammox bacteria. Reaction rates (

*r*) frequently depend on multiple solute concentrations (

_{S}*C*

_{S,n}), and they regularly include nonlinear saturation and inhibition factors.

The boundary conditions for equation 6 were as follows: (i) constant substrate concentrations at the top boundary (the biofilm-bulk liquid interface) (in order to simplify the analysis of model results, the external resistance to mass transfer was neglected; that is, the concentrations of solutes were assumed to be identical to those in the bulk liquid at any point above the plane *z* = *z _{F}*, where

*z*is the maximum biofilm height [Fig. 1B], which increases as the biofilm thickness increases over time) (7) (ii) zero flux of substrate through the impermeable support material (the biofilm-support interface) (8) and (iii) connected lateral boundaries (also called periodic or cyclic boundary type) (9) For the initial state at

_{F}*t*= 0, a uniform distribution of concentrations throughout 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 for dynamic conditions:
(11)
What equation 11 simply accounts for is that the accumulation of dissolved matter in the bulk volume (left side) is the result of (i) flow into and out of the bioreactor (first term on the right side) and (ii) the global conversion rates in the biofilm volume (the triple integral on the right side). The initial condition for equation 11 at *t* = 0 is:
(12)
The solution 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 liquid compartment was kept constant, simulating controlled aeration of the bioreactor.

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

(ii) 2-D and 3-D biofilm model.The dynamics of soluble substrate concentrations due to transport and reaction are very fast compared with the characteristic times of biomass growth, division, and spreading. This observation allowed us to decouple the solution of biomass evolution from that of substrate diffusion-reaction mass balances, as explained in 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 mass

*m*

_{0}and radius

*R*

_{0}]

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

**do**

Biomass dynamics

Biomass division

Biomass spreading (equation 4)

Biomass detachment

Solute dynamics

Mass balances of solutes in the bulk liquid (equation 11)

Mass balances of solutes in the biofilm (equations 6 to 9 solved for steady state)

Time step, *t* ← *t* + Δ*t*

**while** *t* < *t*_{end}

First, a number of inoculum biomass particles are attached to the planar support surface at randomly chosen positions. Values for inlet concentrations of solutes are set in the whole computational domain, discretized in *N _{X}* ×

*N*×

_{Y}*N*cubic elements (or

_{Z}*N*×

_{X}*N*squares in 2-D) of size Δ

_{Z}*x*=

*L*/(

_{X}*N*− 1) (Fig. 1C). Then, for biomass growth and decay, the biomass balances for each particle are integrated directly for one time step (Δ

_{X}*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 according to the shoving mechanism described above is necessary. After some of the particles are removed in the detachment step, the biomass dynamics in the time interval Δ

*t*have been completed; that is, all new positions, masses, and particle radii at

*t*+ Δ

*t*are known. With this new biomass distribution, we now need to compute the new fields of solute concentrations. Because the biomass particles are assumed to be circular and the grid elements of substrate fields are rectangular, the average of the bacterial masses (

*m*) of all cells occupying a certain grid element must be calculated prior to solution of the substrate mass 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 couple a steady-state solute balance for the bulk liquid with the time iterations for biomass and solute balances in the biofilm produced a numerically unstable scheme. The overall solute transformation rates in the biofilm (the triple integral in equation 11) are approximated by summation. The same time step was used for both biomass and bulk solute balances. Second, a steady-state (∂

_{Xi}*C*

_{S,n}/∂

*t*= 0) solution is found for the elliptic form of the partial differential equations (PDEs) for mass balances of soluble substrates in the biofilm (equation 6) with boundary conditions (equations 7 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 PDEs that usually occur in spatial diffusion-reaction problems. Since multigrid algorithms scale linearly with system size (

*N*), they allow a dramatic increase in computational speed for 3-D applications compared with the alternating direction implicit algorithm used in the earlier versions of the model, which scales with

*N*

^{2}(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 biomass detachment (per substratum area) at any moment

*t*was calculated by averaging over the previous 100 time steps (0.5 day).

Case studies. (i) Reaction model.To assess the new biofilm modeling approach against a classic 1-D biofilm model, a multispecies nitrifying biofilm community was chosen as a test case (9). This biofilm contains three active microbial types: aerobic ammonium oxidizers (NH), aerobic nitrite oxidizers (NO), and anaerobic ammonium oxidizers (AN; anammox bacteria). In addition, the inert biomass is a fourth type of particulate material resulting from the decay of all active biomass types. The reaction model includes the substrate dependencies of the various species, which lead to competition among aerobic ammonium oxidizers, aerobic nitrite oxidizers, and anammox bacteria for oxygen, ammonium, and nitrite in the biofilm and production of nitrate and N_{2}. Each of the three microbial types is capable of growth, as well as maintenance (endogenous respiration under aerobic and anaerobic conditions). All model parameters are shown in Table 1. The reaction stoichiometry and kinetics of model processes are summarized in Table 2. The biofilm system analyzed in this study was chosen only to illustrate the model's capability, and it does not represent any real system. Other more realistic but more complex systems can be easily constructed, for example, by including heterotrophic bacteria, as shown by Noguera and Picioreanu (18).

(ii) Test conditions.Three growth regimens were used to evaluate the microbial distribution in biofilms as predicted by multidimensional models. Case I is an oxygen-limited regimen achieved by a low oxygen concentration in the reactor influent (2 g of O_{2} m^{−3}) and a sufficient electron donor concentration (30 g of N-NH_{4} m^{−3} and no other N compounds). In addition, case I* starts with a microbial density eight times lower than that in case I. Case II simulates the biofilm formation at high oxygen levels (10 g of O_{2} m^{−3}) while maintaining the same nonlimiting ammonia concentration. In case III, the influent concentrations are approximately those calculated for the effluent in case II, simulating the behavior of a second biofilm reactor connected in series with the first reactor. In this case, the biofilm forms under nitrogen-limiting conditions [= 4 g of N-NH_{4} m^{−3} and = 3 g of N-NO_{2} m^{−3}]. In all cases, the oxygen concentration in the bulk liquid compartment was assumed to be controlled by aeration, maintaining the constant concentration of the reactor influent.

## RESULTS AND DISCUSSION

The multidimensional biofilm model was evaluated with respect to 1-D models by means of several indicators, such as the spatial distribution of substrate concentrations and fluxes, biofilm structure, total biofilm biomass, dynamics of biomass distribution, dynamics of biomass detachment, dynamics of solute concentrations in the bulk liquid, and steady-state concentrations of biomass and solutes in the biofilm. An example of the spatial distribution of substrate concentrations in 3-D is presented for case III in Fig. 2 and 3. The biofilm geometric structure and total biomass accumulated 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 9 shows the biomass spatial distributions for cases I to III in 2-D simulations. The steady-state concentrations of biomass and solutes in the biofilm are shown in Fig. 10 (case I) and Fig. 11 (case II). Animations of biofilm development can be found at the web page http://www.bt.tudelft.nl/content/ebt/researchgroup_ebt.html.

Spatial distribution of substrate concentrations and fluxes.An essential question is, why do we need spatial (e.g., 2-D or 3-D) biofilm models? What is the relative importance of concentration gradients in directions other than normal to the substratum surface? This computational study and other studies of this kind (7, 23) show that the answer is also a relative one. Depending on the environmental conditions and biofilm structure, for some chemical species there are 3-D gradients, whereas for others the 1-D model description may be sufficient. A 3-D simulation was carried out under the nitrogen-limiting conditions of case III. The spatial distribution of nitrifying bacteria at three moments in time and the corresponding spatial distribution of the intermediate nitrite are shown in Fig. 2. There are obviously regions where NO_{2}^{−} accumulates due to production by ammonium oxidizers and places where NO_{2}^{−} is depleted due to consumption by nitrite oxidizers or anammox bacteria. Magnified representations of the 2-D oxygen and nitrite fluxes and concentration fields in 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 likely to 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 multidirectional gradients are likely to be found. This is clearly shown in Fig. 3b, in which the nitrite flux is oriented along the biofilm surface from the ammonium-oxidizing colonies to the nitrite-oxidizing colonies.

Biofilm geometric structure.Results of multidimensional models can be understood by looking at the differences in geometric biofilm structure predicted for 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 structures is favored under nutrient-limiting conditions. This can be explained as an unavoidable consequence of the existence of strong substrate gradients, which lead to a much greater microbial growth rate at 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 mathematical perspective by Dockery and Klapper (5). Intrinsically, planar 1-D biofilm models assume a planar biofilm surface and concentration boundary layers, as well as 1-D substrate gradients and fluxes. Therefore, an unstable biofilm surface is impossible, and a compact layer with greater total biomass per area develops due to the absence of substrate-depleted valleys.

A sparse initial biomass distribution, as shown previously (24), leads to rugged biofilms with deep vertical channels. This is why 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, corresponding to 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, and 3-D models is an important indicator of biofilm structure. The most direct comparison between the growth curves for the 1-D AQUASIM models and the 2-D and 3-D models can be made under oxygen-limiting conditions because the oxygen gradients are essentially 1-D (case I) (Fig. 5a). If the biofilm thickness is not leveled off by detachment, the multidimensional models predict less biomass than the 1-D models predict. Reduced biomass formation is shown in Fig. 5, which shows that the total biomass was present only until the biofilm thickness reached 500 μm, the critical thickness for the onset of detachment. When the experiment was started with less biomass on the carrier (case I*), not only was the amount of biomass that accumulated lower than the amount that accumulated in case I (and lower than 1-D predictions), but also there was a significant difference between the 2-D and 3-D model results (Fig. 5b). The total biomasses predicted by 1-D, 2-D, and 3-D models do match when compact structures with a flat, horizontal surface form, as they do in the absence of nutrient limitation (case II) (Fig. 7b, inset).

For narrow 2-D systems, in which the length (*L _{X}*) is so small that surface instability cannot develop, a pseudo-1-D system is obtained as a degenerated case of the multidimensional system. The total biomass and bulk concentrations obtained with such a pseudo-1-D model (Fig. 5) are approximately the same as those calculated with 1-D AQUASIM (Fig. 5).

Dynamics of biomass distribution.The development or ecological succession of the total amount of each biomass type attached to the substratum was recorded for all three test cases. The areal biomass concentrations (in grams of chemical oxygen demand [COD] biomass per square meter) were calculated by integration of biomass concentration profiles for the 1-D biofilm:
(13)
and by summation of all particulate biomasses for the 2-D and 3-D models for each *b* = 0, 1, …, *N _{B}*
(14)
Initially, the areal biomass concentrations computed with 1-D and 2-D models match very well (Fig. 6b, 7b, and 8b). However, due to porous biofilm formation in nutrient-limited case I (Fig. 6b), less biomass accumulates in the 2-D model until the detachment begins. During the detachment, the total biomass concentration stays close to 5 g of COD m

^{−2}, but the ratio of the populations changes, progressing to a steady state. The fast-growing population (in this case the ammonium consumers) dominates initially. The thick layer of aerobic bacteria then creates conditions favorable for development of the anaerobic population of anammox bacteria (namely, very low O

_{2}concentrations and very high NH

_{4}

^{+}and NO

_{2}

^{−}concentrations). Gradually, the anammox bacteria displace the ammonium oxidizers, as shown in Fig. 9a. After 57 days, the anaerobic ammonium oxidizers dominate. Note that this is a typical example of ecological succession, where the fast-growing pioneer or opportunist population is replaced by a more slowly growing gleaner population that can sustain a higher standing stock due to, for example, more efficient use of the same resources or changes in resource use (1). Here, the ammonia oxidizers, which use the increasingly limiting oxygen, are replaced by the anammox bacteria, which use nitrite as the alternative resource. The aerobic NO

_{2}

^{−}oxidizers are completely displaced from the biofilm (Fig. 9) because they have to compete on two sides (for oxygen with the aerobic ammonium oxidizers and for nitrite with the anammox bacteria). Good agreement between the biomass ratios of 1-D and 2-D models is shown in Fig. 6b. However, more inert biomass accumulated in the 2-D particle-based model simulations than in the 1-D continuum-based model simulations.

At a high O_{2} concentration (case II; 10 g m^{−3}), the anaerobic layer is much thinner and deeper in the biofilm. Consequently, the anaerobic ammonium oxidizers develop slowly, and the aerobic nitrite oxidizers occupy an important fraction of the biofilm. The initial fast growth of NH_{4}^{+} oxidizers and the catching up by the NO_{2}^{−} oxidizers agree well in 1-D and 2-D models until about 100 days (Fig. 7b). Again, at later stages, formation of inert biomass is much more pronounced in the 2-D model. The 1-D AQUASIM model predicts a steady-state biofilm with very similar proportions of the three populations (Table 3, total biomass in case II), and there is a trend towards a diffuse three-layer distribution with NH_{4}^{+} oxidizers on the top, NO_{2}^{−} oxidizers in the middle, and anammox bacteria on the bottom. However, in the 2-D particle-based model, much of the anammox bacterial layer is replaced by inert biomass (Fig. 9b), because the inert biomass is less mobile in the particle-based model (i.e., it has a lower advective velocity) than in the 1-D continuum model, which reduces the flow of inert biomass to the sheared-off biofilm top. Because of the particle spreading mechanism, the observed trend in all simulations was that the slowest and therefore deepest microbial species were not washed out of the biofilm in the 2-D simulations as quickly as they were in the 1-D simulations.

In N-limited conditions (case III), the whole biofilm development is slowed down. As a result, the maximum thickness (500 μm) and maximum biomass (5 g of COD m^{−2}) are reached much later, only after 200 days. More importantly, the largest biofilm fraction is composed of inert material (Fig. 8b). Active nitrifiers survive only in the top biofilm layers, within the penetration depth of NH_{4}^{+} and NO_{2}^{−} (Fig. 9c), and the bottom biofilm layer is inactive and frozen. This lack of growth and upward pushing of biomass eliminates the reason for the differences between 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 consequence of the fact that the are determined by fluxes to or from the biofilm, which in turn depend mainly on the biomass distribution in the top biofilm layers. Because of the higher substrate concentrations near the biofilm surface, the biomass growth rates are also higher; thus, a stable biomass profile near the surface is established sooner.

Dynamics of biomass detachment.By imposing a maximum biofilm thickness, detachment of biomass particles from the biofilm surface occurs continuously after *z _{F}* = 500 μm has been reached. This results in a gradual filling of biofilm pores with time and, eventually, development of a compact biofilm. The biomass detachment rate during a time step (Δ

*t*) is, in the 2-D model, related to the carrier area (

*L*Δ

_{X}*x*) and results from the removal of

*N*

_{det,Δt}biomass particles: (15) The detachment rate increases during the porous biofilm stage, as the surface area in contact with the

*z*threshold also increases and more biomass particles are pushed out of the upper biofilm boundary (Fig. 6c). The rate of total biomass detachment indicates the rate of overall biomass growth in the biofilm. As expected, under fully aerobic conditions (case II), the flux of detached biomass is higher (0.8 g of COD m

_{F}^{−2}day

^{−1}) (Fig. 7c) than it is in the low-oxygen setting (case I; 0.6 g of COD m

^{−2}day

^{−1}) (Fig. 6c). The fast-growing population of aerobic NH

_{4}

^{+}oxidizers grows preferentially in the top biofilm layer and, consequently, has a high detachment rate. The slowly growing bacteria accumulating in deeper layers must push through the layer of aerobic NH

_{4}

^{+}oxidizers first and, therefore, begin to detach later (anaerobic NH

_{4}

^{+}oxidizers begin to detach 60 days after aerobic NH

_{4}

^{+}oxidizers [Fig. 6c], and aerobic nitrite oxidizers begin to detach 5 days after aerobic NH

_{4}

^{+}oxidizers [Fig. 7c]). Due to the small size of the model biofilm system, noisy oscillations in detachment rates occur due to the synchronized removal of discrete particles.

Steady-state concentrations and fluxes.A perfect steady state with respect to the spatial distribution of biomass and solutes is not reached in the particle-based multidimensional biofilm model. This is due to the stochastic characteristics of the model, as the algorithm uses random numbers on several occasions, including (i) for the initial distribution of biomass particles on the surface, (ii) for the random choice of direction for the placement of daughter particles, and (iii) for uneven cell division, which is a desynchronizing division event in the biofilm.

Four variables that have practical importance, substrate concentration , substrate fluxes (*J*_{S,i}), total biomass in the biofilm (*X _{b}*), and detachment rates , were selected for comparison of the 1-D multispecies model with the 2-D particle-based model at steady state (Tables 3 and 4). Solute concentrations in the bulk liquid were in general very similar in the two models, while a greater difference was observed for nitrite and nitrate in oxygen-limited case I. Profiles for substrate concentrations in the biofilm [

*C*

_{S,i}(

*z*)] are shown in Fig. 10c (case I) and Fig. 11c (case II). The biomass distribution in the biofilm,

*C*

_{X,b}(

*z*), displays some differences (Fig. 8a, 10, and 11). In general, the agreement among biomass profiles is better for the upper biofilm layers, where, due to the presence of higher substrate concentrations, the biomass growth rates and consequently the advective velocities are higher. The substrate profiles computed for the two models agree very well, indicating that they are not very sensitive to the biomass distribution in the biofilm depth. As a consequence, the substrate fluxes to the biofilm (Table 4) are very similar. Finally, the good match in biomass profiles near the biofilm surface is reflected in good agreement for biomass detachment rates at a steady state (Table 3).

These observations suggest that complete knowledge of the 2-D or 3-D distribution of microbial species is probably not needed to calculate overall conversion rates (where the interest is usually from an engineering perspective), but species localization might have a great effect on inherently local microbial competition, cooperation, and selection (where the interest of microbiologists often is).

Conclusions.In this paper we describe a spatially multidimensional (2-D and 3-D) particle-based approach for dynamic modeling of multispecies biofilms growing on multiple substrates. Predictions of the new 2-D and 3-D particle-based models were quantitatively and critically compared with results for equivalent 1-D problems obtained with the classic, 1-D multispecies model (33, 34) implemented in the AQUASIM software (29). For comparison, biofilms containing three 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_{2}) were studied under three growth regimens.

The multidimensional models predicted less biomass per carrier area (i.e., more porous biofilms), especially under nutrient-limited conditions and when biofilm development started from a spatially sparse inoculum.

Chemical species that are either consumed or produced in geometrically homogeneous biofilms mainly form only one-directional gradients of concentration. For intermediates, which are produced and consumed in different locations in the biofilm, strong multidirectional gradients are likely to be found. These gradients can be predicted accurately only by 2-D and 3-D models.

Comparison with the 1-D predictions showed that full knowledge of the 2-D or 3-D distribution of microbial species may not be required for calculation of solute fluxes and overall conversion rates. However, from an ecological and evolutionary point of view, the distribution of species might have great effects on local microbial competition and cooperation, as well as selection processes.

Especially for the slowest growing species and for the inert biomaterial the 1-D and 2-D models predict different dynamics and depth profiles. This is mainly caused by differences in biomass spreading mechanisms (overall advective velocity through the entire 1-D slice versus local pushing of hard spheres in a 2-D or 3-D space), combined with different coupling of inert biomass to active biomass. Experimental observations of the spreading of color-tagged bacteria through growing biofilms should help improve this aspect of biofilm models.

It is believed that IbM of 3-D biofilm structure and the particle-based approach as its scaled-up version, although computationally more demanding than 1-D models, go one step further towards modeling biofilm systems from first principles. Not only physical and chemical processes, such as substrate transport and reactions, but also biomass behavior can be described in a mechanistically meaningful 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 models could explain a large variety of phenomena from microbiology and microbial ecology, retaining at the same time the quantitative features required in engineering applications.

## FOOTNOTES

- Received 20 November 2003.
- Accepted 16 January 2004.
- ↵*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

- American Society for Microbiology