Individual-based modelling of biofilms

Jan-Ulrich Krefta,1,2, Cristian Picioreanu2, Julian W. T. Wimpenny1 and Mark C. M. van Loosdrecht2

Cardiff School of Biosciences, Cardiff University, PO Box 915, Cardiff CF10 3TL, UK1
Kluyver Institute of Biotechnology, Delft University of Technology, Julianalaan 67, 2628 BC Delft, The Netherlands2

Author for correspondence: Jan-Ulrich Kreft. Tel: +44 29 2087 5278. Fax: +44 29 2087 4305. e-mail: kreft{at}cardiff.ac.uk


   ABSTRACT
TOP
ABSTRACT
INTRODUCTION
THE INDIVIDUAL-BASED MODEL...
ANALYSIS OF MODEL OUTPUT
TESTING (VALIDATION) OF THE...
RESULTS
DISCUSSION
REFERENCES
 
Understanding the emergence of the complex organization of biofilms from the interactions of its parts, individual cells and their environment, is the aim of the individual-based modelling (IbM) approach. This IbM is version 2 of BacSim, a model of Escherichia coli colony growth, which was developed into a two-dimensional multi-substrate, multi-species model of nitrifying biofilms. It was compared with the established biomass-based model (BbM) of Picioreanu and others. Both models assume that biofilm growth is due to the processes of diffusion, reaction and growth (including biomass growth, division and spreading). In the IbM, each bacterium was a spherical cell in continuous space and had variable growth parameters. Spreading of biomass occurred by shoving of cells to minimize overlap between cells. In the BbM, biomass was distributed in a discrete grid and each species had uniform growth parameters. Spreading of biomass occurred by cellular automata rules. In the IbM, the effect of random variation of growth parameters of individual bacteria was negligible in contrast to the E. coli colony model, because the heterogeneity of substrate concentrations in the biofilm was more important. The growth of a single cell into a clone, and therefore also the growth of the less abundant species, depended on the randomly chosen site of attachment, owing to the heterogeneity of substrate concentrations in the biofilm. The IbM agreed with the BbM regarding the overall growth of the biofilm, due to the same diffusion-reaction processes. However, the biofilm shape was different due to the different biomass spreading mechanisms. The IbM biofilm was more confluent and rounded due to the steady, deterministic and directionally unconstrained spreading of the bacteria. Since the biofilm shape is influenced by the spreading mechanism, it is partially independent of growth, which is driven by diffusion-reaction. Chance in initial attachment events modifies the biofilm shape and the growth of single cells because of the high heterogeneity of substrate concentrations in the biofilm, which again results from the interaction of diffusion-reaction with spreading. This stresses the primary importance of spreading and chance in addition to diffusion-reaction in the emergence of the complexity of the biofilm community.

Keywords: biofilm structure, nitrification, spatial heterogeneity, chance, complexity

Abbreviations: 2D, 3D, two-, three-dimensional; IbM, individual-based model/modelling; BbM, biomass-based model/modelling; CV, coefficient of variation

a Present address: Abteilung Theoretische Biologie, Botanisches Institut, Universität Bonn, Kirschallee 1, D-53115 Bonn, Germany.


   INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
THE INDIVIDUAL-BASED MODEL...
ANALYSIS OF MODEL OUTPUT
TESTING (VALIDATION) OF THE...
RESULTS
DISCUSSION
REFERENCES
 
Biofilms are microbial assemblages that occupy an interface, creating a chemically distinct microenvironment. They develop at any interface with conditions at least temporally favourable for microbial growth. Important examples where biofilms occur are teeth (dental plaque), catheters, stones in riverbeds, trickling filters, plant leaves, soil, marine and freshwater sediments, ship hulls, water and sewage pipelines, and bioreactors (for reviews see Bryers & Characklis, 1982 ; Costerton et al., 1995 ; Stoodley et al., 1999 ; Tolker-Nielsen & Molin, 2000 ).

Biofilm models have undergone an evolution of increasing complexity (Noguera et al., 1999a ), prompted by a growing body of experimental evidence of biofilm heterogeneity. The first models, in the 1970s, described biofilms as uniform steady-state films of a single species with one-dimensional mass transport and biochemical reactions. In the 1980s, stratified dynamic models of multisubstrate/multispecies biofilms were developed. However, they could not generate the characteristic biofilm morphology but had to use a given biofilm structure as input into the model. With the aim to explain new experimental findings, and facilitated by the advances in both computational power and numerical methods, 2D and three-dimensional (3D) models were developed in the 1990s. They incorporate the whole range of transport processes as well as biofilm growth and detachment to various extents (Wimpenny & Colasanti, 1997 ; Picioreanu et al., 1998a , b ; Hermanowicz, 1999 ; Noguera et al., 1999b ; Eberl et al., 2000 , 2001 ). Such models will be referred to as biomass-based models (BbMs) in this paper. In these models, biofilm structure is an emergent property rather than the model input because they follow a bottom-up approach where the complex community emerges as a result of the actions and interactions of the biomass units with each other and the environment. They can also be classified as spatially structured population models.

Individual-based modelling (IbM) is also a bottom-up approach as it attempts to model a population or community by describing the actions and properties of the individuals comprising the population or community (Huston et al., 1988 ; DeAngelis & Gross, 1992 ; Grimm, 1999 ). In contrast to BbM, IbM allows individual variability and treats organisms, in our case bacterial cells, as the fundamental entities. This has important consequences regarding biomass spreading. When such an entity changes position, its biomass, fixed and variable properties (such as its genome, state of differentiation, etc.) and its cell number (one) are displaced together.

All bottom-up models have the potential to address questions about the relationship of microscopic and macroscopic properties. How do simple microscopic units (cells) give rise to complex macroscopic structures (biofilms)? How does chance affect these structures?

IbM is particularly suited to address questions about the effects of individual variability. Before we can take full advantage of this potential, however, we must first emphasize the prediction of ‘normal’ biofilm structure.

We have applied IbM to a nitrifying biofilm with one ammonia- and one nitrite-oxidizing species as our model system, because it provides a simple example of a food chain and, moreover, nitrification is an important process in natural environments, sewage treatment plants and agriculture.

The aims of this paper were the development and validation of the first IbM of biofilms, based on an IbM of Escherichia coli colony growth (Kreft et al., 1998 ), and the comparison of results with a BbM adapted from Picioreanu et al. (1998a , b ) by using the same kinetics, parameters and initial conditions.


   THE INDIVIDUAL-BASED MODEL BacSim
TOP
ABSTRACT
INTRODUCTION
THE INDIVIDUAL-BASED MODEL...
ANALYSIS OF MODEL OUTPUT
TESTING (VALIDATION) OF THE...
RESULTS
DISCUSSION
REFERENCES
 
This fully quantitative IbM, called BacSim, is the much improved and extended version 2 of the original BacSim model, which was the first spatially explicit IbM of bacterial growth (Kreft et al., 1998 ). BacSim is an extension of the Objective-C program Gecko, which in turn is based on the Swarm toolkit (Kreft et al., 1998 ). In addition, BacSim and Gecko have recently been ported to Java by Ginger Booth (Yale University, CT, USA) and the Java versions of BacSim and Gecko can now be run as applets within a web browser (see http://www.cf.ac.uk/biosi/staff/kreft/bacsim.html).

BacSim essentially consists of two parts: one deals with the simulation of the growth and behaviour of individual bacteria as autonomous agents; the other deals with the simulation of substrate and product diffusion and reaction (Fig. 1). Since biofilm growth is usually a much slower process than diffusion of substrate into the biofilm, the diffusion process can be simulated assuming the growth process to be frozen, and the growth process can be simulated assuming the diffusion process to be in a pseudo-steady state (Picioreanu et al., 1999 ). This decoupling of processes based on differences in their time scale is standard for BbM and the main improvement from the previous version of BacSim (Kreft et al., 1998 ). Although the numerous extensions of the features of BacSim make it impossible to compare execution times of the two versions of BacSim exactly, as a rough guideline, runs are now about 500 times faster.



View larger version (36K):
[in this window]
[in a new window]
 
Fig. 1. IbM (BacSim) model structure. The two components of the model, the ‘World’ and the ‘Bacteria’, and the program flow, alternating between diffusion and growth steps, is shown. The mutatable parameters are shown in bold type.

 
The biofilm is simulated in a Cartesian space with width (lx), depth (ly) and height (lz) (Fig. 2). For the substrate diffusion-reaction, a two-dimensional (2D) space is used, but for the bacteria, a continuous 3D space is used to allow the same degrees of freedom for cellular movement as in reality. However, it is essentially a 2D model because the extra dimension used for the bacterial movements was kept very small at about two cell diameters.



View larger version (63K):
[in this window]
[in a new window]
 
Fig. 2. Schematic of the system geometry. Periodic boundaries are in place for cells and substrates in the x (width) and y (depth) dimensions, whereas the z (height) dimension is non-periodic. Located at the bottom is the inert substratum and at the top the bulk liquid. The boundary layer at z=zb separates the bulk liquid with constant substrate concentrations (infinite reservoir) from the biofilm and its surrounding liquid where substrate transport takes place by diffusion of substrate from the bulk liquid into the biofilm. The boundary layer moves up with biofilm growth, since the distance to the top of the biofilm at z=zt is kept constant. The y dimension is also referred to as the extra dimension in this paper. Note that this system is essentially 2D, but sketched wider (10 µm), shorter (20 µm) and lower (20 µm) than usual for visualization purposes.

 
The assumptions of the model are as follows. Only diffusion, reaction, growth (including cell growth, cell division and cell spreading) and chance are considered. The bacteria are discrete spherical and immotile cells with kinetic and yield parameters, which may vary randomly around a given mean. Initially, the bacteria are randomly placed close to the inert, impermeable substratum. The top of the biofilm is at a constant distance from the bulk liquid. This distance constitutes the mass transfer boundary layer through which the chemical species are transported by diffusion only. The spreading of cells only occurs when cells get too close to each other. It is a deterministic, continuous process in time and space.

Simulation of bacteria. Each bacterium is individually simulated as a sphere of variable volume in a continuous, 3D space.

The single cell uptake rate for the nitrogen-containing substrate vN is described by a Monod equation with double substrate limitation for both the electron donor (denoted N for nitrogen-containing substrate, i.e. ammonia or nitrous acid, depending on the species) and acceptor (denoted O for oxygen, for both species), as well as a substrate inhibition term for the respective electron donor substrate:

(1)

where MX is the amount of biomass (dry mass), Vmax,N is the maximum specific uptake rate for the nitrogen-containing substrate, CN and CO are the concentrations of the respective substrates, KS is the Monod saturation constant and Ki is the substrate inhibition constant. The uptake rate for oxygen, vO , is related to vN by the yield ratio:

(2)

where YN is the growth yield for the electron donor and YO is the growth yield for the electron acceptor.

Although the nitrogen-containing compounds are ionized at neutral pH, the saturation and inhibition constants given in Table 1 are expressed in terms of the neutral species because NH3 and HNO2 are the real substrates. This makes the constants independent of temperature and pH. Given temperature and pH, these kinetic constants, K, can be converted to an expression for the ionized species:


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

 

(3)


(4)

where T is the absolute temperature in Kelvin (Wiesmann, 1994 ).

The observed growth rate of the cell is given by:

(5)

where µ is the specific growth rate and m is the maintenance rate, expressed on the basis of rate of biomass decrease (Herbert model; see Kreft et al., 1998 ).

The size of a cell is important, since it affects the cell’s metabolic rates and starvation survival. To model the well known dependence of cell size on growth rate, we have used the Donachie model for cell division (Donachie & Robinson, 1987 ) in a descriptive form (for a detailed discussion see Kreft et al. (1998) . Donachie and coworkers found that cells divide a constant time after initiation of DNA replication. Faster growing cells will gain more volume during this time interval than slower growing cells, therefore they will be larger at division. Donachie and coworkers also established an empirical relationship for the dependence of cell volume on growth rate, which we have used here. The volume at cell division Vd, is a function of the number of generations per hour, g:

(6)

The minimal volume at division is a function of the maximal volume at division:


(7)

The spreading of the growing biomass is simulated by maintenance of a minimum distance between neighbouring cells. For each cell, the vector sum of all positive overlap radii with neighbouring cells is calculated and then the position of the cell is shifted in the direction opposite to this vector. The overlap radius, Ro, is given by:



(8)

where Rc is the radius of the cell, Rn is the radius of the neighbouring cell, d is the euclidean distance between the centres of the cell and its neighbour and k is a multiplier on the cell radius that allows adjustment of the minimal spacing between cells (shove radius). The shove radius multiplier prevents a closer packing of cells than is physically possible (Kreft et al., 1998 ). If the cell’s surface extended into the substratum, the overlap with the substratum was added to the vector sum of overlap radii. This spreading mechanism is illustrated by a movie available at our website (http://www.cf.ac.uk/biosi/staff/kreft/biofilms.html) and also at Microbiology Online (http://mic.sgmjournals.org).

Boundary conditions. The x and y dimensions had periodic boundaries, that is, cells that become shifted beyond such a boundary plane re-enter the domain through the opposite boundary plane. In other words, the sides are wrapped around. This avoids edge effects, since there are, effectively, no edges. The z dimension was non-periodic, with bacteria allowed to spread towards the bulk liquid but not into the solid substratum (Fig. 2).

Initial conditions. If not mentioned otherwise, the system was inoculated with 10 bacteria of each species, placing them at randomly chosen locations on the surface of the substratum at time zero.

Sorting bacteria by location. A tree data structure was used to sort bacteria by location in the vertical xz plane to find neighbours efficiently. Each leaf of the tree stores a reference to one bacterium. Traversing this tree will visit each bacterium in turn and therefore all bacteria in a particular sequence that is determined by the tree structure. The tree is constructed from a list of bacteria by calculating the bounding box around all area patches occupied by the list members. This box is then subdivided into quadrants and the list members are sorted into these quadrants according to their location. This procedure is repeated recursively for each sublist of quadrant members until each list contains fewer than 11 members. Resorting the tree is done every 10 steps by reorganizing all entries in the tree into a single list in a way that inverts the order of this new list compared to the original list from which the tree was constructed. (This inversion is exact only if no list members were added due to cell divisions or removed due to deaths.) This new inverted list is then used to reconstruct the tree in the same way as above. If list members have changed their location during the last 10 steps, they may become sorted into a different quadrant, and the quadrants may also change in area and location.

Random variation of cell parameters. If not mentioned otherwise, both the maximal uptake rate and the volume-at-division of each bacterium were varied by random draws from a Gaussian distribution with a coefficient of variation (CV) of 10%. This CV is a value typical for the few cases were this has been studied (for a discussion see Kreft et al., 1998 ). Draws resulting in parameters outside their default ±2{sigma} were skipped, as were changes of sign. This is necessary as the Gaussian distribution has an interval of [-{infty},+{infty}], which would lead to draws of parameter values that are physically impossible or beyond the range allowed by biological constraints.

Simulation of substrate fields
The concentrations of the dissolved substrates and products (here oxygen, ammonia, nitrite and nitrate), result from a solution of their mass balances, including transport and reaction terms. Here, we only consider diffusive transport. The continuous diffusion-reaction equation for each compound i is given by:

(9)

where D denotes the diffusion coefficient and r the reaction rate of the given compound.

Note that the uptake rates are calculated by each bacterium, but this information has to be transferred to a grid for computing the substrate fields. Since the specific uptake rates of different bacteria of the same species were allowed to vary randomly, a rectangular lattice for storing the momentary reaction rates must be constructed by querying all cells for their uptake rates at given substrate concentrations. Each cell’s reaction rates are apportioned into the lattice elements covered by this cell on an area percentage basis (this procedure of accurately gridding the reaction rates is left out of equations 10–13 below for the sake of simplicity). This process has to be repeated for every iteration of the substrate field relaxation method, because reaction rates depend on substrate concentrations. The reaction rates for the different compounds in grid element p with volume Vp relate to the single-cell uptake rates of all the bacteria j of species 1 and all the bacteria k of species 2 in grid element p as follows (assimilation of nitrogen not considered):


(10)


(11)


(12)


(13)

For the substrate fields, the space is discretized with a rectangular mesh with Nx1xL elements, numbered from 0 to N-1, 0 and L-1 for the dimensions x, y and z, respectively. Note that this mesh was only one layer deep in the y direction but the depth of this single layer was variable. Unless stated otherwise, 100x1x100 volume elements were used with a length of 2 µm for all sides (lengths in x and y were always equal), giving a computational domain of 200x2x200 µm.

Boundary conditions. The zero-flux boundary condition ({partial}Ci/{partial}z=0) is assumed at the substratum, z=0. At the top of the system, z=lz, the bulk liquid with constant substrate concentration (infinite reservoir) was located. This bulk liquid phase extended down to the mass transfer boundary layer, the top of which was approximated by a plane located at a constant distance above the top of the changing biofilm front (Fig. 2). Here, we use a typical value for the boundary layer thickness of 40 µm (Picioreanu et al., 1998b ), thereby modelling the effect of convective flow rather than modelling convective flow directly. For the bulk liquid, a fixed-value boundary condition was used with the substrate concentration being set to the external substrate concentration throughout the bulk liquid phase. For the other dimensions, periodic boundaries were used.

Initial conditions. The substrate concentration was set to the external substrate concentration throughout.

Solution. The system of non-linear second order partial differential equations (9–13) was solved for the steady-state solution using the alternating direction implicit method (Peaceman & Rachford, 1955 ; Ames, 1977 ). The discretization used for each grid point was a five-point centred finite difference scheme also used in Picioreanu et al. (1998b ).

Comparison with the BbM
The chosen BbM was essentially the model developed by Picioreanu et al. (1998a , b ), tailored to the simulation of a two-species nitrifying biofilm in two dimensions. The two models are based on the same algorithms for the solution of substrate diffusion-reaction mass balances. Regarding the biomass, the IbM and BbM differed in three fundamental aspects: (a) the space for the biomass was discrete for the BbM but continuous for the IbM; (b) the biofilm consisted of uniform ‘bricks’ of biomass for the BbM but individual cells with mutable parameters for the IbM; and (c) biomass spreading algorithms. In the BbM, when ‘bricks’ of biomass divide, one of the halves is randomly shifted into the neighbourhood with a preference for an unoccupied site. In the IbM, the spherical cells shove each other every time step to minimize overlap. Both models use the same equations and parameters (Table 1) for the growth of the biomass and the same system geometry and initial conditions. The maximal uptake rate and the size-at-division were randomly varied in the IbM (with the mean equal to the setting in the BbM). However, this hardly affected results (see below).

Biomass spreading in the BbM was as follows. If the maximum biomass density was reached in a given lattice element, the biomass density in that lattice element was halved and an equal amount of biomass was placed into a randomly chosen free grid element in the eight-connected neighbourhood. If none of the grid elements in this neighbourhood was free, a randomly chosen neighbour was displaced. The search for a free grid element was then recursively started over again with the displaced neighbour as the new starting point. Based on these cellular automata rules, three different versions of spreading were derived that treated the mixing of the two species differently. In the first scheme, called ‘apart’, different species occupied different grid elements and stayed apart during spreading. In the second scheme, called ‘coupled’, different species could jointly occupy the same grid element and were redistributed as a unit. In the third scheme, called ‘mixed’, different species could occupy the same grid element as in ‘coupled’, but upon reaching the maximum density allowed per grid element, only the biomass of the species constituting the larger fraction was redistributed.

It was essential that both models only differed in those respects that we intended to compare and be equivalent otherwise. Therefore, in addition to using the same equations and parameters for diffusion, reactions and growth, the following points had to be taken care of as well. The initial random distribution of bacteria along the substratum in the IbM was converted into a matrix of initial biomass distribution that was read into the BbM to ensure equal but random initial distribution of biomass in both models. Also, the mean biomass density of the biofilm resulting from the IbM was used to calculate the maximum biomass density for the BbM, using a previously established ratio of mean to maximum biomass density for the BbM. Since both models differed in the way the maintenance energy requirements were modelled, the maintenance energy was set to zero for the sake of comparison. The effective diffusion coefficient was assumed to be equal to the diffusion coefficient in water for the same reason.


   ANALYSIS OF MODEL OUTPUT
TOP
ABSTRACT
INTRODUCTION
THE INDIVIDUAL-BASED MODEL...
ANALYSIS OF MODEL OUTPUT
TESTING (VALIDATION) OF THE...
RESULTS
DISCUSSION
REFERENCES
 
All analysis procedures were programmed in MATLAB (The MathWorks, Natick, MA, USA) with the Image Processing and Statistics Toolboxes.

Biofilm structure. The aim was to establish biologically relevant quantitative descriptors of biofilm structure that allowed comparison of various biofilms, simulated or real. Many such descriptors have been proposed (Zhang & Bishop, 1994a , b ; Gibbs & Bishop, 1995 ; Murga et al., 1995 ; Hermanowicz et al., 1995 , 1996 ; Kreft et al., 1998 ; Kwok et al., 1998 ; Picioreanu et al., 1998b ; Lewandowski et al., 1999 ; Heydorn et al., 2000 ; Yang et al., 2000 ), but their relative merit has not yet been assessed by a comparative study. See Table 2 for a list of measures used.


View this table:
[in this window]
[in a new window]
 
Table 2. Deviations of growth, shape and texture parameters among repeated runs

 
When applying such measures taken from the literature, the methods used to calculate them have to be modified to take the geometry of the biofilm system into account. (a) Vertical and horizontal cross-sections may have to be treated differently. The vertical cross-section is anisotropic, as it is made along the vertical and one horizontal dimension, whereas the horizontal cross-section is isotropic as it is made along the two horizontal dimensions. (b) The slice through the system to which the 2D image analysis procedure is applied must be thin enough to avoid artefacts of projection of the 3D biofilm surface onto a plane, or the image analysis procedure must be suitable for 3D images (Russ, 1994 ). (c) The biofilm image should be cropped down to the bounding box of the biofilm features so as to avoid including bulk fluid in the image analysis to various extents, which affects most measures drastically. (d) The substratum, if impermeable to the substrate, should be treated as a part of the solid phase rather than as a second biofilm surface when calculating measures relating to mass transfer, such as diffusion distance. (e) For the same reason, interior holes should be ignored for the purpose of calculation of mass-transfer-related measures, since they do not contribute significantly to mass transfer. (f) As a special feature of simulated biofilms, the periodic boundaries have to be considered in all calculations involving padding the image border, e.g. diffusion distance. Instead of padding with zeros, the borders have to be padded with the last column or row from the opposite side of the image. (g) The measures should ideally be independent of image resolution and orientation. (h) Last, but not least, the descriptors should be unaffected by short-term fluctuations in the biofilm structure. These may result from the shoving of cells from one time step to the next, illustrated by a movie available at our website (http://www.cf.ac.uk/biosi/staff/kreft/biofilms.html) and also at Microbiology Online (http://mic.sgmjournals.org). It is unclear whether these fluctuations occur in natural biofilms in a similar manner or whether they are artefacts of the simulation. The variation of the shape measures step by step from 29975 to 30025 min (1 step min-1) was calculated. Since some of the measures, e.g. biomass and solids hold-up, showed a clear trend, the trend was removed by differentiation before the standard deviation was calculated. Taking the CV as the criterion for estimating the degree of short-term fluctuations, all growth and shape measures varied by less than 1%, apart from contrast, which varied by 5%.

Discussion of biofilm structure analysis. All measures can be used to quantify the structures of model biofilms, since they are not too affected by fluctuations from time step to time step. The dependence of the shape measures on the resolution, 1 or 2 µm, of the analysis of the same run (IbM with 1 µm spatial resolution; data not shown) was very strong, particularly for measures of heterogeneity, although image analysis shape measures are supposed to be independent of resolution. Only surface enlargement and solids hold-up were the same at the different resolutions. This renders comparisons of biofilm structures sampled at different resolutions meaningless. Contrast and heterogeneity are both measures of some aspects of heterogeneity, with contrast showing higher sensitivity, differentiating between the mixed and the other two versions of the BbM as well as the IbM (see Fig. 8). Also, contrast was the only measure that varied within a given biofilm structure (data not shown) as would be expected of a good measure of heterogeneity. Heterogeneity, however, is highest close to the biofilm surface, and low and level within the biofilm (data not shown). Therefore, contrast seemed to be the better measure of internal heterogeneity. By this token, the IbM, compared with the BbM, appeared to have a more homogeneous internal structure and a less homogeneous structure close to the surface. We feel that the objective quantification of biofilm structure remains an important challenge, despite recent progress, since the utility of shape measures can only be assessed when they have been applied to the full diversity of biofilm shapes, whether real or simulated. This study is only a small contribution to this task.



View larger version (26K):
[in this window]
[in a new window]
 
Fig. 8. Comparison of biofilm shape parameters (see Table 2) in the IbM (——) with the three versions of the BbM (......, apart; -·-·-·, coupled; -----, mixed). (a) Surface enlargement, (b) roughness, (c) solids hold-up, (d) diffusion distance, (e) heterogeneity and (f) contrast.

 
Statistical analysis. With computer simulations, statistical significance of the differences between results from different runs, as marked by a 95% confidence interval, is almost always obtained due to the high precision of the raw data. This does not invalidate these procedures in any way, but renders them less useful. Therefore, a criterion for the decision whether statistically significant differences between runs are biologically relevant or significant must be introduced. We define, for our purposes, biological significance as a divergence in the trends of the time series of different runs of more than 5%. Divergence is not given if trends are interwoven or largely parallel while starting from the same point. If both statistical and biological significance are given, the runs will be called considerably different, and if not, their differences will be called negligible in this paper. Since traces for different runs often deviated only after some time, the analysis should be able to distinguish different phases in the time series and treat them individually.

The following statistical method was found to be appropriate. The difference between a given and a reference dataset, divided by the reference dataset, was computed (called ‘relative residuals’ later on). To test for a significant difference of the mean of these relative residuals from zero, a T-test can be used. Since a divergence of trends may only occur for a part of the trace, say the last half, we have used a moving T-test with a window of 10 values. If a stretch of significant deviations from random scatter around zero was found, and this stretch was long enough to allow an estimate of the trends for this time interval (five data points), we determined the range of the deviation of the trends in this window of significance. The trends were estimated by robust smoothing (3RSR, followed by two Hannings; see Tukey, 1977 ) for less noisy traces, or a stronger averaging procedure (moving median with a window of three data points to remove spikes, followed by a moving mean with a window of five data points to dampen) for noisier traces. (The trend in the data could be removed by one or two differentiations.) The relative difference (divergence) between these trends was plotted and a moving mean curve (five data point window) drawn. If the range of this curve was below 0·05, the divergence of trends was defined as negligible.


   TESTING (VALIDATION) OF THE MODEL
TOP
ABSTRACT
INTRODUCTION
THE INDIVIDUAL-BASED MODEL...
ANALYSIS OF MODEL OUTPUT
TESTING (VALIDATION) OF THE...
RESULTS
DISCUSSION
REFERENCES
 
Biomass spreading. The original colony model (Kreft et al., 1998 ) uses a unilateral shoving mechanism (called ‘classic’) where the shoving algorithm accesses all bacteria in a particular sequence, computing the vector sum of the overlaps of the current bacterium with all its neighbours and then shoving only the current bacterium by this vector. This method is biased, since the sequence of accessing the bacteria affects who is shoved by whom. To minimize the effects of this bias, the sequence of who shoves whom was inverted every 10 steps. Here, a mutual shoving mechanism (called ‘mutual’) was implemented that is free of this bias, since the current bacterium is not only shoved by its neighbours, but in addition also shoves the neighbours. This makes shoving independent of the sequence with which the bacteria are accessed. The results of runs with different shoving schemes differ considerably, in growth and more so in shape parameters, with the notable exception of surface enlargement. All these differences were caused by the difference in biomass density attained with the two shoving schemes (Fig. 3). The biomass density was lower in the ‘mutual’ scheme and fluctuated more. This indicates that the ‘mutual’ scheme is better in avoiding overlaps and can achieve a relocation of cells that got stuck. The lower biomass density achieved implies that the biofilm is more stretched out. Since this will reduce the proportion of cells close to the substratum, the overall growth of the biofilm will be enhanced.



View larger version (19K):
[in this window]
[in a new window]
 
Fig. 3. Biomass density within BbM or IbM biofilms grown at 1 mg oxygen l-1. For the IbM, the ‘classic’ ({circ}) and the ‘mutual’ ({bullet}) versions of the shoving mechanism are compared. All versions of the BbM are similar, therefore only the ‘mixed’ version is shown here ({square}).

 
Sequential execution biases. Since the program carries out activities of the bacteria sequentially, which in reality are simultaneous, the possibility of biases due to the sequence of execution arises. The bacterial activities can be divided into two groups. First, those that change the size and number of cells (growth, division and death) do not directly affect the outcome of the same activities in the other bacteria. Second, those that change the location of bacteria (shoving) do affect the other bacteria. Three types of bias can be expected. If only one pass through the bacteria is performed, and the list of activities of each bacterium is executed as a single group, some bacteria will change in size or position, or divide or die, before others, yet they will influence the shoving of bacteria that are updated only later. This bias can be avoided by making two passes through the bacteria, in the first pass executing the group of activities that change size and number and then carrying out the ‘mutual’, unbiased shoving of cells in the second pass. A second bias comes from where in the calling sequence daughter cells arising from cell division are inserted. This will determine whether both the mother and daughter cells will be called at the time step of the cell division event or not. (Clearly, mother and daughter cells cannot exist at the same time in reality, but they do co-exist in the model during the time step in which the division takes place.) There are three alternatives: (a) inserting them at the start of the sequence; (b) inserting them at the end of the sequence; (c) performing shoving before division and death. Comparing the results of these three different schemes showed that this second bias is negligible. The third bias comes from the sequence of shoving as the repositioning of a cell is affected by the previous repositioning of other cells. This third bias can be avoided only by a ‘mutual’ shoving algorithm (see above). The effect of the first two of these biases was negligible in practice and, for performance reasons, usually only a single pass through the bacteria was performed. Inverting the sequence of bacteria every 10 growth steps probably kept the first two biases very low. Also, a higher time resolution (see below) reduces the effect that any sequential execution biases may have.

The extra dimension for bacteria. The extra dimension (y), the depth of the biofilm (Fig. 2), was kept small to avoid biomass gradients and for reasons of computational efficiency. To check for the absence of a biomass gradient across the total extent of this dimension, we simulated biofilms with depths of 1, 2, 3, 4, 5 and 10 µm (data not shown). No gradient spanned the extent of the extra dimension, but small-scale local gradients reflecting arrangement of cells into layers were found.

Spatial resolution of the lattices for the diffusion-reaction scheme. The size of the lattice cells of all the lattices used in the diffusion-reaction scheme was varied, using 1, 2, 4, 5 and 10 µm for the sides of the squares. Initially, all of the growth curves of the majority species deviated strongly from the reference run (1 µm resolution), but growth curves converged again later (Fig. 4). Specifically, growth curves were back to less than 5% relative difference when the biofilm reached a height of about 6–8 grid cells, e.g. 16 µm at 2 µm resolution and 30 µm at 5 µm resolution. In other words, the thinner the biofilm, the higher the resolution has to be for accurate growth rates. If a 5% relative difference is chosen as the acceptable maximal difference, then the resolution should be at least 1/8 of the biofilm height. Since we have used a constant resolution of usually 2 µm, growth results are not very accurate initially, differing by up to 12%.



View larger version (20K):
[in this window]
[in a new window]
 
Fig. 4. Relative differences of growth curves at 2 ({circ}), 4 ({square}), 5 ({lozenge}) and 10 µm ({triangledown}) spatial resolution, using the highest resolution of 1 µm as the reference. Growth of the majority species, the ammonia oxidizer, is shown in terms of the number of grid elements occupied by the biofilm (relative maximal height). Time is implicit in this plot, since the biofilm height is a function of time; data points (symbols) are shown at 1000 min intervals. The growth curves first diverge and then converge again; the differences fall below 5% when the biofilm has reached a thickness of 6–8 grid elements.

 
No considerable deviations were found at 2 µm resolution for any of the biofilm shape parameters; however, for the lower resolutions of 4, 5 and 10 µm, several of these parameters deviated considerably. Therefore, we chose a resolution of 2 µm as generally suitable. We plan to use an adaptive multigrid method which is free of this problem in the future.

Accuracy of gridding biomasses, reaction rates and substrate concentrations. The procedure of doling out the biomass of each spherical cell, according to the cell’s position, into the square biomass grid was carried out in a simple and also in a more accurate way. With the simple procedure, all of the cell’s biomass was added to the biomass grid in which the cell’s centre was located. With the more accurate procedure, the biomass was doled out into the eight neighbouring biomass grid elements in addition to the central grid, in proportion to the area coverage of an equivalent-volume cube representing the spherical cell. The reaction rates were allocated and the substrate concentrations were averaged in the same manner. The growth of the majority species as well as biofilm shape parameters did not deviate considerably, but the growth of the minority species was clearly affected by the biomass apportioning. Therefore, the more accurate scheme was used for all runs. Apparently, the growth of the minority species is more sensitive to slight local changes in substrate concentration, since it only exists in a few different locations in a small number. It must be assumed that the growth of individual cells can be strongly affected by these slight changes, but for the growth of larger populations spread over many different locations, these changes are averaged out.

Time resolution for growth. The growth time step of the simulation was set to 1, 2, 5 or 10 min and the results compared. At a time step of 10 min, many of the parameters deviated; at a time step of 5 min, growth of both species deviated; at a time step of 2 min, results agreed well. Note that these results depend on the growth rate of the species, hence, time resolution should be higher for faster growing bacteria.

Discussion of model validation. For validation of the IbM, we checked that effects due to (a) code optimization (data not shown), (b) implementing the algorithms in two different programming languages (Objective-C and Java; data not shown) and (c) sequence of execution of the bacterial activities were negligible. Further, the biomass distribution in the third dimension was level. Sufficiently accurate partitioning of the reaction rates of the cells into a reaction rate matrix required distribution of rates into an eight-connected neighbourhood on the basis of area coverage. The spatial and time resolution of the simulation was fine enough and the convergence criterion for diffusion was stringent enough (data not shown).

The new ‘mutual’ shoving scheme was better than the ‘classic’ unilateral shoving scheme (Kreft et al., 1998 ), since the biomass was more spread out at the same setting for the minimum distance between cells (shove radius). Also, the ‘mutual’ scheme is free of bias due to the sequence of shoving.


   RESULTS
TOP
ABSTRACT
INTRODUCTION
THE INDIVIDUAL-BASED MODEL...
ANALYSIS OF MODEL OUTPUT
TESTING (VALIDATION) OF THE...
RESULTS
DISCUSSION
REFERENCES
 
In all our simulations, the growth of the nitrite oxidizer was only marginal, since the ammonia oxidizer has higher uptake rates (Table 1) under conditions of low oxygen tension which were due to the high concentrations of the nitrogen substrates (Table 1). The low abundance of the nitrite oxidizer was useful for comparing different simulations and therefore substrate concentrations were not altered to facilitate the growth of the nitrite oxidizer.

Effect of random variations
Stochastic daughter cell placement upon division. There are three uses of random numbers in the simulations. One is varying parameters of the bacteria (see section on individual variability below), another one is in choosing the initial location of cells within given bounds (see section on stochastic initial cell locations below) and the last one is in choosing the directions of displacement of the daughter cells from the mother’s position prior to division. To evaluate the effect of the latter on biofilm growth and shape, we made 50 repeats of a simulation starting with identical positions of the bacteria and drawing random numbers from one non-repeating series. Each run was compared with the mean of all runs (Table 2, Fig. 5a, b). The growth of the dominant species (and therefore the biofilm) did not deviate considerably in any of the runs (Fig. 5a). But the growth of the less abundant species deviated considerably in two runs (Fig. 5b). Almost all of the biofilm shape parameters agreed amongst the 50 runs. Therefore, the placement of daughter cells did not have a major effect in the system studied here.



View larger version (16K):
[in this window]
[in a new window]
 
Fig. 5. Comparison of growth curves of 50 runs starting with identical (a, b) or different (c, d) random initial cell locations. The median is enveloped by the 5 and 95% percentiles. (a, c) Growth of the ammonia oxidizer; (b, d) growth of the nitrite oxidizer -·-·-·, 5%; ——, 50%; -----, 95%.

 
Stochastic initial cell locations. As above, we compared the results of 50 runs of the simulation, drawing from one series of random numbers, but starting with randomly chosen locations for the bacteria (Table 2, Fig. 5c, d). The variability among these runs was far higher than for the set of 50 runs with identically located initial attachment sites. In many cases, the total growth of the biofilm and of the majority species was divergent and the growth of the minority species was at the mercy of the growth conditions set by the majority species, and therefore almost always different.

Individual variability. We compared a simulation where bacterial parameters were kept constant with one where both Vmax,N and volume-at-division were varied randomly. No differences between the two runs were found apart from temporary deviations for contrast.

Comparison of IbM with BbM
The IbM was compared with the BbM using two different sets of random initial locations for the bacteria. These sets of locations were generated with the IbM, converted into biomass matrices for the two species, saved to file and read into the BbM. Since the original version of the BbM, called ‘apart’, did not allow the coexistence of biomass of both species in the same grid element, the initial biomass matrices for the two species should not have overlapping entries. However, the other two versions, called ‘coupled’ and ‘mixed’, allow such coexistence. Since the two species do not mix in the ‘coupled’ version, they will stay apart if they were apart at the beginning, and results will be the same as with the ‘apart’ version. Therefore, to compare all versions of the BbM with the IbM, one set of non-overlapping and one set of partially overlapping initial locations were used.

The biomass density of the BbM version was set to the biomass density of the IbM, since this density is not known a priori (Fig. 3). It emerges from the interaction of cell sizes and growth rates with the spreading mechanism of the IbM, which uses the ‘shove radius’ parameter (Table 1).

The growth of the majority species, and therefore the growth of the biofilm, agreed reasonably well between the IbM and all the versions of the BbM (Fig. 6a). Although the trend of some of the BbM versions diverged from the IbM considerably, the differences between some of the BbM versions were more pronounced than the difference to the IbM. This is also true for the growth of the minority species (Fig. 6b), but here the discrepancies between some of the runs were dramatic.



View larger version (15K):
[in this window]
[in a new window]
 
Fig. 6. Comparison of growth of (a) the ammonia oxidizer and (b) the nitrite oxidizer in the IbM (——) with the three versions of the BbM (......, apart; -·-·-·, coupled; -----, mixed).

 
The IbM biofilm developed into three fingers that corresponded to similar but more irregular and less-rounded fingers in the BbM biofilms (Fig. 7). The shape of the IbM biofilm differed from the BbM biofilms in all measures (Fig. 8). The IbM biofilm was more compact and less porous, the cell clusters were more confluent and therefore larger (higher diffusion distance), and the surface was smoother (roughness) and less convoluted (surface enlargement). Regarding measures of heterogeneity, results are somewhat inconsistent, as the IbM biofilm had the highest heterogeneity, but the lowest contrast.



View larger version (100K):
[in this window]
[in a new window]
 
Fig. 7. Images of IbM and BbM biofilms with contours of oxygen concentration after 67 d (96000 min) of growth. (a) BbM apart, (b) BbM coupled, (c) BbM mixed and (d) IbM. The oxygen limitation of growth and the height of the mass transfer boundary layer were evident from the contours; the highest concentration was 1 mg oxygen l-1. Ammonia oxidizers (light grey) were more abundant than nitrite oxidizers (darker grey); coexistence of both leads to intermediate shades of grey. Note the wrapped around periodic boundary; there are only three fingers in the biofilms shown. Colour images and movies visualizing the growth of these four biofilms are available at our website (http://www.cf.ac.uk/biosi/staff/kreft/biofilms.html) and also at Microbiology Online (http://mic.sgmjournals.org).

 
The shapes of the biofilms of all three versions of the BbM were similar, apart from solids hold-up, roughness and diffusion distance. The ‘apart’ and ‘coupled’ versions were particularly close (Figs 7 and 8), even though the initial conditions were not completely identical in the runs shown, since there was one case of overlap in the initial biomass lattices. (If the initial conditions are non-overlapping, the two versions will give identical results.) The ‘mixed’ version of the BbM might be expected to be more similar to the IbM, since the spreading of the two species is independent in both, but it is furthest from the IbM regarding the growth of the minority species, solids hold-up, roughness and contrast.

Finer textures have a higher contrast (Weszka et al., 1976 ). This would indicate a coarser texture for the IbM.

Results with the non-overlapping set of initial conditions were similar to the results with the partially overlapping initial conditions (data not shown).


   DISCUSSION
TOP
ABSTRACT
INTRODUCTION
THE INDIVIDUAL-BASED MODEL...
ANALYSIS OF MODEL OUTPUT
TESTING (VALIDATION) OF THE...
RESULTS
DISCUSSION
REFERENCES
 
The nitrite oxidizers were less abundant than the ammonia oxidizers as a result of the high nitrogen substrate concentrations relative to the oxygen concentration in the simulations, which led to a strong oxygen limitation of growth. This choice of substrate concentrations, exemplifying the treatment of sewage with a high nitrogen load, was found suitable for the purpose of comparing the two models, since it resulted in poor growth of the nitrite oxidizers. They are also referred to as the minority species in this paper. For this marginal species, the effect of substrate concentration heterogeneity on growth was more pronounced because the small number of cells makes the variability more obvious.

The effect of random variation of the two key growth parameters, Vmax,N and volume-at-division, was negligible in contrast to the colony model (Kreft et al., 1998 ), due to the influence of the very high heterogeneity of substrate concentration in the biofilm. The growth and shape of biofilms was independent of random variations in the placement of daughter cells upon division.

The interplay of the stochastic process of initial attachment with the heterogeneity of substrate concentration in the biofilm drastically altered the growth of single cells into clones (Fig. 9), and therefore of less abundant species with a small number of individuals in the population. By analogy, stochastic attachment events after formation of the biofilm would have a similar, though smaller effect. Furthermore, in the BbM, the biomass spreading mechanisms were themselves stochastic, explaining why the different versions of biomass spreading of the BbM resulted in strong differences in the growth of the minority species. It should be noted that this implies that modelling the growth of rare species requires careful consideration of the spreading mechanism appropriate for this species.



View larger version (95K):
[in this window]
[in a new window]
 
Fig. 9. Growth of clones from the 20-cell inoculum. Each of the 10 clones of the ammonium oxidizer is shown in a different colour, but the 10 clones of the nitrite oxidizer, which mostly do not grow very much, are all shown in the same colour (magenta; see the foot of the leftmost finger for the largest clone of nitrite oxidizers). The length of the substratum was 200 µm.

 
Comparison of the BbM and IbM shows that both models agreed regarding the overall growth of the biofilm, due to the same diffusion-reaction process. But the biofilm shape was different and, to a degree, independent of diffusion-reaction since it was also affected by the biomass spreading mechanism in combination with chance. The IbM biofilm was more confluent, compact and smoother, due to the steady, deterministic and directionally unconstrained shoving of the spherical bacteria. Smoother, more confluent biofilm shapes, similar to the IbM biofilms, were also obtained with a spreading mechanism based on diffusion of the biomass (diffusion coefficient increases sharply as biomass reaches critical density), as might be expected from a diffusion process (Eberl et al., 2001 ). This diffusional spreading mechanism is also operating steadily, deterministic and directionally unconstrained, but the biomass is continuous rather than discrete as in the IbM. In the IbM, the diffusion-reaction process shaped the biofilm into fingers, as in other bottom-up biofilm and colony models (Ben-Jacob et al., 1994 ; Matsushita, 1997 ; Wimpenny & Colasanti, 1997 ; Picioreanu et al., 1998b ; Hermanowicz, 1999 ; Noguera et al., 1999b ), but the effect is lessened by the more fluid biomass spreading mechanism, placing a higher emphasis on the biological processes. All these models use growth parameters from pure liquid cultures, where bacteria are not exhibiting a biofilm phenotype, but predict realistic biofilm structures. However, different phenotypic characteristics of cells in biofilms may affect the spreading of these cells. The importance of biomass spreading is also highlighted by the strikingly different structures formed by various closely related strains of Pseudomonas growing under identical conditions (Tolker-Nielsen et al., 2000 ; Heydorn et al., 2000 ).

The BbM biofilm shapes were clearly more related to one another than to the IbM (Fig. 8), since all versions use the same principle of biomass spreading. The observed differences in the biofilm shapes of the three versions of the BbM can be explained by the stochastic nature of the spreading mechanism which let different fingers grow better in different runs of the same version (data not shown) or in runs of the different versions (Fig. 7). Since the spreading of the minority species is independent of the spread of the biofilm in the mixed version, results of this version were particularly stochastic. Also, the minority species had higher chances to be shifted towards higher oxygen concentrations, resulting in the best growth of all models (Fig. 6b).

BbM biofilms from runs with a spatial resolution of 1 µm were more confluent and rounded than at 2 µm (data not shown) and therefore more similar to the IbM, which gave biofilms that were even more rounded than the 1 µm BbM biofilms. The spatial resolution of the biomass is fixed in the IbM, since it corresponds to the mean size of the bacteria, which were about 1 µm in diameter. Due to the stochastic nature and the resolution dependence of the spreading mechanism of the BbM, its biofilms often grew into different patterns of fingers at different resolutions and the growth of the minority species varied drastically. At the higher resolution, the species became more mixed in the ‘mixed’ version; therefore, growth of the minority species was much better. But the growth of the majority species was still almost identical.

In conclusion, both the BbM and IbM results agreed in principle, due to modelling the same physical processes, but differed in details of biofilm shape and growth of minority species. Therefore, if one wants to focus on rare species or rare events, the choice of model should depend on which biomass spreading mechanism is supported by experimental results. Until enough is known about this, the choice of model is rather arbitrary and one might want to use both. The BbM is clearly better suited to model larger scale systems, which is particularly important when the system is characterized by larger scale heterogeneities.

Biofilm growth, driven by reaction-diffusion, interacts with biofilm shape, influenced by the biomass spreading mechanism. Chance events modify the biofilm shape and the growth of single cells, because of the high heterogeneity of substrate concentrations in the biofilm, which again results from the interaction of diffusion-reaction with spreading. This stresses the primary importance of spreading and chance in addition to diffusion-reaction in the emergence of the complexity of the biofilm community.


   ACKNOWLEDGEMENTS
 
We would like to thank Ginger Booth (Yale University) and Hermann Eberl (GSF, Munich) for fruitful discussions about biomass spreading and Ginger also for porting BacSim to Java. Financial support from the BBSRC (Biotechnology and Biological Sciences Research Council, UK, grant 72/E09748) is gratefully acknowledged.


   REFERENCES
TOP
ABSTRACT
INTRODUCTION
THE INDIVIDUAL-BASED MODEL...
ANALYSIS OF MODEL OUTPUT
TESTING (VALIDATION) OF THE...
RESULTS
DISCUSSION
REFERENCES
 
Ames, W. F. (1977). Numerical Methods for Partial Differential Equations 2nd edn. London:Thomas Nelson & Sons.

Ben-Jacob, E., Schochet, O., Tenenbaum, A., Cohen, I., Czirok, A. & Vicsek, T. (1994). Generic modeling of cooperative growth-patterns in bacterial colonies. Nature 368, 46-49.[Medline]

Bergey, D. H., Buchanan, R. E., Gibbons, N. E. & Cowan, S. T. (1974). Bergey’s Manual of Determinative Bacteriology 8th edn. Baltimore:Williams & Wilkins.

Bryers, J. D. & Characklis, W. G. (1982). Processes governing primary biofilm formation. Biotechnol Bioeng 24, 2451-2476.

Costerton, J. W., Lewandowski, Z., Caldwell, D. E., Korber, D. R. & Lappin-Scott, H. M. (1995). Microbial biofilms. Annu Rev Microbiol 49, 711-745.[Medline]

DeAngelis, D. L. & Gross, L. J. (1992). Individual-Based Models and Approaches in Ecology: Populations, Communities, and Ecosystems New York:Chapman and Hall.

Donachie, W. D. & Robinson, A. C. (1987). Cell division: parameter values and the process. In Escherichia coli and Salmonella typhimurium: Cellular and Molecular Biology, pp. 1578–1593. Edited by F. C. Neidhardt and others. Washington, DC: American Society for Microbiology.

Eberl, H. J., Picioreanu, C., Heijnen, J. J. & Loosdrecht, M. C. M. (2000). A three-dimensional numerical study on the correlation of spatial structure, hydrodynamic conditions, and mass transfer and conversion in biofilms. Chem Eng Sci 55, 6209-6222.

Eberl, H. J., Parker, D. F. & Loosdrecht, M. C. M. (2001). A new deterministic spatio-temporal continuum model for biofilm development. J Theor Med 3, 161–176.

Gibbs, J. T. & Bishop, P. L. (1995). A method for describing biofilm surface roughness using geostatistical techniques. Water Sci Technol 32(8), 91–98.

Grimm, V. (1999). Ten years of individual-based modelling in ecology: what have we learned and what could we learn in the future? Ecol Model 115, 129-148.

Haralick, R. M., Shanmugam, K. & Dinstein, I. (1973). Textural features for image classification. IEEE Trans Syst Man Cybern SMC-3, 610-621.

Hermanowicz, S. W. (1999). Two-dimensional simulations of biofilm development: Effects of external environmental conditions. Water Sci Technol 39(7), 107–114.

Hermanowicz, S. W., Schindler, W. & Wilderer, U. (1995). Fractal structure of biofilms: new tools for investigation of morphology. Water Sci Technol 32(8), 99–105.

Hermanowicz, S. W., Schindler, W. & Wilderer, U. (1996). Anisotropic morphology and fractal dimensions of biofilms. Water Res 30, 753-755.

Heydorn, A., Nielsen, A. T., Hentzer, M., Sternberg, C., Givskov, M., Ersbøll, B. K. & Molin, S. (2000). Quantification of biofilm structures by the novel computer program COMSTAT. Microbiology 146, 2395-2407.[Abstract/Free Full Text]

Huston, M. A., DeAngelis, D. L. & Post, W. (1988). New computer models unify ecological theory. BioScience 38, 682-691.

Kreft, J.-U., Booth, G. & Wimpenny, J. W. T. (1998). BacSim, a simulator for individual-based modelling of bacterial colony growth. Microbiology 144, 3275-3287.[Abstract]

Kwok, W. K., Picioreanu, C., Ong, S. L., van Loosdrecht, M. C. M., Ng, W. J. & Heijnen, J. J. (1998). Influence of biomass production and detachment forces on biofilm structures in a biofilm airlift suspension reactor. Biotechnol Bioeng 58, 400-407.[Medline]

Lewandowski, Z., Webb, D., Hamilton, M. & Harkin, G. (1999). Quantifying biofilm structure. Water Sci Technol 39(7), 71–76.

Matsushita, M. (1997). Formation of colony patterns by a bacterial population. In Bacteria as Multicellular Organisms , pp. 366-393. Edited by J. A. Shapiro & M. Dworkin. New York:Oxford University Press.

Murga, R., Stewart, P. S. & Daly, D. (1995). Quantitative analysis of biofilm thickness variability. Biotechnol Bioeng 45, 503-510.

Noguera, D. R., Okabe, S. & Picioreanu, C. (1999a). Biofilm modeling: Present status and future directions. Water Sci Technol 39(7), 273–278.

Noguera, D. R., Pizarro, G., Stahl, D. A. & Rittmann, B. E. (1999b). Simulation of multispecies biofilm development in three dimensions. Water Sci Technol 39(7), 123–130.

Peaceman, D. W. & Rachford, H. H.Jr (1955). The numerical solution of parabolic and elliptic differential equations. J Soc Indust Appl Math 3, 28-41.

Picioreanu, C., van Loosdrecht, M. C. M. & Heijnen, J. J. (1997). Modelling the effect of oxygen concentration on nitrite accumulation in a biofilm airlift suspension reactor. Water Sci Technol 36(7), 147–156.

Picioreanu, C., van Loosdrecht, M. C. M. & Heijnen, J. J. (1998a). A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads. Biotech Bioeng 57, 718-731.

Picioreanu, C., van Loosdrecht, M. C. M. & Heijnen, J. J. (1998b). Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach. Biotech Bioeng 58, 101-116.

Picioreanu, C., van Loosdrecht, M. C. M. & Heijnen, J. J. (1999). Discrete-differential modelling of biofilm structure. Water Sci Technol 39(7), 115–122.

Russ, J. C. (1994). Fractal Surfaces. New York: Plenum Publishing Corporation.

Shuler, M. L., Leung, S. K. & Dick, C. C. (1979). A mathematical model for the growth of a single bacterial cell. Ann NY Acad Sci 326, 35-55.

Stoodley, P., Boyle, J. D., De Beer, D. & Lappin-Scott, H. M. (1999). Evolving perspectives of biofilm structure. Biofouling 14, 75-90.

Tolker-Nielsen, T. & Molin, S. (2000). Spatial organization of microbial biofilm communities. Microbiol Ecol 40, 75-84.[Medline]

Tolker-Nielsen, T., Brinch, U. C., Ragas, P. C., Andersen, J. B., Jacobsen, C. S. & Molin, S. (2000). Development and dynamics of Pseudomonas sp. biofilms. J Bacteriol 182, 6482-6489.[Abstract/Free Full Text]

Tukey, J. W. (1977). Exploratory Data Analysis. Reading, MA: Addison-Wesley.

Weszka, J. S., Dyer, C. R. & Rosenfeld, A. (1976). A comparative study of texture measures for terrain classification. IEEE Trans Syst Man Cybern SMC-6, 269-285.

Wiesmann, U. (1994). Biological nitrogen removal from wastewater. Adv Biochem Eng Biotechnol 51, 113-154.[Medline]

Wimpenny, J. W. T. & Colasanti, R. (1997). A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models. FEMS Microbiol Ecol 22, 1-16.

Yang, X. M., Beyenal, H., Harkin, G. & Lewandowski, Z. (2000). Quantifying biofilm structure using image analysis. J Microbiol Methods 39, 109-119.[Medline]

Zhang, T. C. & Bishop, P. L. (1994a). Density, porosity and pore structure of biofilms. Water Res 28, 2267-2277.

Zhang, T. C. & Bishop, P. L. (1994b). Evaluation of tortuosity factors and effective diffusivities in biofilms. Water Res 28, 2279-2287.

Received 10 January 2001; revised 9 July 2001; accepted 20 July 2001.