SIMULATION OF SINGLE-SPECIES BACTERIAL-BIOFILM GROWTH USING THE GLAZIER-GRANER-HOGEWEG MODEL AND THE COMPUCELL3D MODELING ENVIRONMENT

The CompuCell3D modeling environment provides a convenient platform for biofilm simulations using the Glazier-Graner-Hogeweg (GGH) model, a cell-oriented framework designed to simulate growth and pattern formation due to biological cells’ behaviors. We show how to develop such a simulation, based on the hybrid (continuum-discrete) model of Picioreanu, van Loosdrecht, and Heijnen (PLH), simulate the growth of a single-species bacterial biofilm, and study the roles of cell-cell and cell-field interactions in determining biofilm morphology. In our simulations, which generalize the PLH model by treating cells as spatially extended, deformable bodies, differential adhesion between cells, and their competition for a substrate (nutrient), suffice to produce a fingering instability that generates the finger shapes of biofilms. Our results agree with most features of the PLH model, although our inclusion of cell adhesion, which is difficult to implement using other modeling approaches, results in slightly different patterns. Our simulations thus provide the groundwork for simulations of medically and industrially important multispecies biofilms. 2000 Mathematics Subject Classification. Primary: 92-08, 92C15; Secondary: 65C05, 68U20.

1. Introduction.This paper presents the computational methodology to implement a preliminary model of single-species biofilm growth using the Glazier-Graner-Hogeweg (GGH) model [33,34,37], 1 as implemented in the CompuCell3D [17,18,44] modeling environment. 2Our fairly simple model includes an explicit representation of a nutrient field, but does not include explicit extracellular matrix (ECM).We show that even our simple model reproduces existing results, providing a platform for exploring new phenomenology due to explicit cell shape, ECM properties and multispecies interactions.
A biofilm is an aggregation of microorganisms embedded in a protective and adhesive ECM (also called an extracellular polymeric substance or exopolysaccharide (EPS) matrix) [2].A biofilm is not a multicellular organism but a community of independent individuals.Although biofilms spontaneously form patterns and structure, no evidence suggests that this formation is the result of a genetically controlled developmental program [51,78].Gene expression in biofilms closely correlates with gene expression in stationary-phase bacteria in liquid culture with some degree of anaerobic and iron-limited growth.Biofilms often contain many different species and types of microorganisms, e.g., bacteria, archaea, protozoa, algae, and fungi [20].They may be structurally and genetically heterogeneous, with complex community interactions and with members of each species performing specialized metabolic functions [14].
Biofilms are very common, existing on almost any solid substrate submerged in, or exposed to, water, from river bottoms to the lining of the gut [14].They can also form floating mats on liquid surfaces.Even single-species biofilms exhibit a wide variety of morphologies, making them interesting examples of spontaneous patterning (pattern formation) [13,19,23].Biofilms have no "typical" structure, since, like trees in forests, biofilms can produce all sorts of shapes.These shapes depend on multiple factors including cell motility, the existence of multiple cell phenotypes, the nature of the secreted ECM, flows in the surrounding liquid medium, the properties of the surface on which the biofilm grows, and mass transport in the liquid medium and ECM.Thus biofilms serve as interesting models of the origin of multicellularity, cooperativity, and ecological specialization, making simulations of multispecies biofilms potentially useful in the study of these basic biological questions.
Biofilms contribute to a host of industrial and medical problems, ranging from ship fouling to antibiotic resistance and tooth decay.Biofilms also present a unique opportunity in environmental technologies for cleanup and containment of hazardous waste.Since biofilms are involved in about 80% of microbial infections in humans [1], a realistic understanding of how biofilms grow 3 and how they contribute to antibiotic resistance would be of practical value.
Because of the difficulty of observing and controlling real biofilms, they have attracted much attention from modelers wishing to clarify the roles of different biological and physical mechanisms in determining biofilm properties.Biofilm models generally divide into cell-centered models with discrete cells (cellular automata [8,39,69] or individual-based models [25,52,85]), and continuum models [3,26,28,84].Hybrid models combine discrete cells with a continuum representation of the ECM and/or extracellular fluid [4] or with PDEs for reaction-diffusion [66,67,68].The first continuum, one-dimensional (1D) mathematical models of biofilms considered only gradients perpendicular to the biofilm-liquid interface [83].Later, two-dimensional (2D) models based on physical principles explored the ways lateral gradients create biofilm structures [8,25,28,39,52,66,67,68,69].These models addressed the spontaneous formation of heterogeneous structures, the influence of flow-induced breakage and substrate gradients, and the role of pores in biofilms [82].
However, we still lack both a fundamental theory and accurate experimental observations of biofilm growth [82].In addition, most existing models have certain limitations that may impair their viability in describing complex real biofilms.Many cell-centered and hybrid models describe individual cells as points [66] or spheres [4].However, cell shape can be essential in biological development, 4 and changes in cell shape reflect the elastic properties of cells.
Thus, this paper takes a cell-oriented, hybrid approach to biofilm modeling [56].We treat cells as macroscopic objects with mechanical properties. 5Since the size of substrate molecules is small compared to the size of the biofilm cells, we treat substrates as continuous diffusive fields. 6The main goal of this paper is to validate our GGH-based approach by reproducing known simulation results for single-species biofilms.This validation is essential to future studies, which will apply the GGH model to multispecies biofilms, where differences in cell adhesion and elasticity between various types of cells may greatly affect biofilm development.Because the methods against which we are benchmarking do not include explicit cell shape, we model cells as deformable bodies rather than cells with a fixed shape, as in real bacteria.However, we can easily add explicit cell shapes in future simulations.Treating cells as extended, deformable objects allows us to investigate the role in biofilm formation of cell-cell and cell-ECM adhesion, which are crucial in many types of biological pattern formation [76].Although analyzing the effects of cell adhesion does not require spatially extended objects [12], our approach allows us to introduce 4 Bacteria come in a great range of shapes and sizes, both spherical and nonspherical, and have rigid cell walls.In this paper we study roughly spherical clusters of bacteria, although our approach allows us to describe nonspherical bacteria and control the average shape and rigidity of cells in our biofilm simulations by changing the corresponding parameters in the effective energy (see section 3).We can also introduce different surface energies between various cell types, which drives cell sorting and results in spatially varying species densities within the biofilm (not presented in this paper).
5 Each generalized cell (see section 3) in the model we present in this paper represents a cluster of many real cells and their associated ECM, so it can be deformable while the real cells are rigid.The only mechanical properties of cells we model in this paper are their surface tension and compressibility, although we could easily add a shear modulus and viscosity [6,33]. 6Although the size of substrate molecules is small, the substrate density may, if it is low enough, fluctuate over distances that are large compared to the bacterial length scale.These fluctuations are the reason that bacteria cannot, in general, sense spatial gradients directly [9].We omit these fluctuations since, as we show later, biofilm morphology is sensitive to only one parameter, which relates substrate uptake to biofilm growth rate, and is averaged over a time scale much longer than that of the typical fluctuations in substrate concentration.
elasticity and surface tension directly, in terms of free energies [33].Furthermore, the GGH model allows us to include in a simple and physically correct fashion many other biologically relevant phenomena, e.g., cell polarity and compartmental cell structure [6].
2. The PLH biofilm model.Overall, our choice of biological mechanisms follows the very successful model of Picioreanu, van Loosdrecht, and Heijnen (PLH) [66,67], with the important difference that our cells are extended, elastic, adhesive objects.We do not model ECM explicitly.However, we include it implicitly in the form of a cell-medium adhesion coefficient in the GGH effective energy (see section 3).Thus, our GGH biofilm cells correspond to clusters of real biofilm cells surrounded by a thin layer of ECM.
The hybrid differential-discrete PLH model [66,67] focuses on the development of single-species biofilms containing Nitrosomonas europaea or Nitrobacter agilis, using kinetics derived from [7,24].A 2D rectangular mesh discretizes space.Two variables represent the state of the simulated biofilm: the soluble substrate7 concentration S, and the biomass density C. In addition, the solid component matrix c stores information about the occupation of space.The quantities S and C vary between 0 and 1, while the occupation state c takes either the value 0 (an unoccupied site) or 1 (a site with biomass).The variables S and C are defined as where c X is the biofilm density, c Xm is the maximum biomass density in a colony, c S is the local substrate concentration and c S0 is the substrate concentration in the bulk liquid.The substrate concentration c S in 2D obeys the diffusion equation where D S is the substrate diffusion coefficient, and the substrate consumption rate, r S , depends on both c S and the biomass density in the biofilm, c X .By construction, the PLH model does not include the effective surface energy between cells or between cells and the bulk liquid.
In the absence of detailed experimental data, most biological models, including the PLH model, assume a Monod-like saturation function to describe the substrate consumption rate [67]: where µ m denotes the maximum specific growth rate, Y XS is the growth yield (amount of biomass produced per unit of substrate), m S is the maintenance coefficient and K S is the Monod saturation constant [67].The equation for biomass balance is ∂c where the net rate of biomass formation, r X , is Combining equations ( 4)- (6) gives which is the same as Ċ At low substrate concentrations, this model allows negative biomass growth (shrinkage), while at high substrate concentrations the biofilm growth rate reaches an observable maximal value.To simulate cell division and mechanical elasticity, when the biomass concentration C reaches its maximum value 1 at a given location, it divides into two equal parts.The first stays where it started, while the other moves to a randomly chosen, adjacent grid element, updating the occupation matrix c through a cellular-automaton mechanism [67].Specific boundary conditions govern substrate uptake from the environment.For a biofilm forming on a solid support, the substrate concentration is constant and maximal (S = 1) outside the aqueous medium [66].
In the PLH model, biofilm morphology depends on a single parameter, the nondimensional ratio of the maximum biomass growth rate to the maximum substrate transport rate.In terms of the quantities introduced in the preceding paragraphs, this ratio (called the Growth Group Parameter) is [66] where l z is the maximum biofilm thickness.If G is small, the substrate diffuses fast enough that an appreciable quantity reaches all regions of the growing biofilm, and its isoconcentration lines are flat and horizontal.Therefore, the growing biofilm remains essentially compact and flat.For larger values of G, faster substrate consumption prevents the substrate from reaching all regions of the growing biofilm, inducing competition for the substrate.In this case, both the substrate distribution and the biofilm density become more irregular, and deep vertical channels (valleys) form in the biofilm density in a fingering instability [21,38,60,64].When G is very large, growth is very slow and filamentous structures form.This general behavior holds in both 2D and 3D simulations [66].
3. The GGH biofilm model.A general model of a growing biofilm needs three main elements [26].First, a transport mechanism like diffusion or advectiondiffusion must bring growth substrates (e.g., oxygen, nutrients) to the cells in the biofilm.Second, cells must grow and proliferate while consuming growth substrates.The third and least understood element is biofilm loss due to breakage or cell death.
In addition to these fundamental mechanisms, we could explore the effects of other factors such as cell-cell or cell-ECM adhesion, 8 or the elastic reactions of biofilms to fluid-flow-generated shear stress.
GGH models include objects, interaction descriptions, dynamics and initial conditions.Objects in a GGH model can be generalized cells or fields.Generalized cells are spatially extended domains, which reside on a single cell lattice and may correspond to clusters of biological cells, biological cells, subcompartments of biological cells, or to portions of noncellular materials, e.g., ECM, fluids, solids, etc. [6,33].Generalized cells can carry within them chemical concentrations, state-descriptors (e.g., time since last mitosis), descriptions of biochemical interaction networks, etc.
Fields are continuously variable concentrations, each of which resides on its own lattice.Fields typically represent chemical diffusants, nondiffusing ECM, etc.Multiple fields can be combined to represent materials with textures, e.g., fibers.
Interaction descriptions and dynamics define how the various objects behave both biologically and physically.For generalized cells, these behaviors and interactions are embodied primarily in the effective energy (also known as the CPM Hamiltonian), which determines a generalized cell's shape, motility, adhesion and response to extracellular signals [33].Local minimization of the effective energy [76] drives most pattern dynamics [33,34,37].The effective energy mixes true energies, such as cell-cell adhesion, and terms that mimic energies, e.g., the response of a cell to a chemotactic gradient of a field.Mathematically, N spatially extended cells indexed by σ lie on a 2D or 3D lattice; the value at a lattice site (pixel) i is σ if this site lies in cell σ.A collection of connected lattice sites with the same index represents a cell, as shown in Figure 2. To each cell we associate a cell type τ .Different cells can have the same value of τ .
A typical effective energy in our simulations includes two main terms [33,34,37]: The first term describes the surface energy between cells, and between cells and their environment, in terms of symmetric surface energies J(τ, τ ) = J(τ , τ ) [33,34,37].
To calculate the energy resulting from cell-cell adhesion we consider first-, second-, and third-nearest neighbor lattice sites, which reduces lattice anisotropy effects compared to first-neighbor calculations [43]. 9The second term corresponds to the compressibility of the cells and prevents cell disappearance; V (σ) is the volume in lattice sites of cell σ, V t its target volume, and λ V (τ ) the strength of the volume constraint. 10 of the cell membrane or wall: 2 , where S is the surface area of a cell and St its target surface area [17,44].This term fixes the average shape of the cells.We could also control the shape of the cells explicitly using elongation constraints [57].Another possibility is a term of form (SV −κ − StV −κ t ) 2 (κ = 1/2 in 2D and 2/3 in 3D) to be consistent with the classical theory of elasticity [70].Instead of the surface-adhesion coefficients J(τ, τ ), we can use surface tensions γ(τ, τ ), defined as Negative surface tensions correspond to mixing or dispersing cells, while positive surface tensions correspond to cell sorting.We simulate growth of cells by increasing their V t by a function f which typically depends on the local concentrations of chemicals C that induce growth and the internal state of the cell [71]: If a cell reaches a given doubling volume, it divides and splits along a random axis (with no splitting orientation) into two cells with equal target volumes of one half the parent-cell volume (we can easily add oriented cell division, but it is not included in the PLH model [66,67]).Adding constraints to the effective energy can describe many other cell properties, including osmotic pressure and membrane area [17,44], the average shape of the cells [57,86], chemotaxis [40,54], viscosity and advective diffusion [22], rigid-body motion [6], and cell compartments for the simulation of polarized cells, e.g., in epithelia [11].
The cell lattice evolves through attempts by generalized cells to extend their boundaries into neighboring lattice sites, slightly displacing the generalized cells which currently occupy those sites.These extensions change the effective energy, and we accept each one with a probability that depends on the resulting change of the effective energy according to an acceptance function [33].Thus, the cell lattice evolves to minimize the total effective energy [32,36] (we can translate the effective energy rigorously into a force-based formalism if we choose, but the effective-energy formalism is more convenient).The rearrangement dynamics uses a relaxational Monte-Carlo-Boltzmann-Metropolis dynamics [33,34,37,59], modified so that averaged local movements of cell boundaries obey the overdamped forcevelocity relationship (modified Metropolis dynamics): where µ local is the local mobility of the object experiencing the force and v the local boundary velocity.The cells rearrange their positions to minimize the total effective energy [76].At each step we randomly select a lattice site i and attempt to change its index from σ( i) to the index σ = σ( i ) of an arbitrary lattice site i in its third-order neighborhood (a lattice site which connects to i by a side, edge, or corner).We accept the change with a probability P :11 where ∆H is the difference in effective energy produced by the change, θ is the Heaviside step function, and T is a parameter corresponding to the amplitude of cell fluctuations [61]. 12One Monte Carlo Step (MCS) by default corresponds to n index-change attempts, where n is the total number of lattice sites (however, in our biofilm simulations, we will define a Monte Carlo Step as n 100 index-change attempts; see section 4).
Auxiliary equations describe cells' absorption and secretion of diffusants and extracellular materials (i.e.their interactions with fields), state changes within the cell, mitosis, and cell death [33].Usually, state changes affect generalizedcell behaviors by changing the parameters present in the effective energy (e.g., the surface density of a particular cell adhesion molecule).Fields evolve due to cells' secretion and absorption at their boundaries, and diffusion, reaction and decay under the control of PDEs.The initial condition specifies the initial configurations of all generalized cells and fields, all internal states related to auxiliary equations and any other information required to completely describe a model.
If, as in the PLH papers cited, simulated bacteria are aerobes and we assume that they have plenty of nutrients so their growth is limited only by their oxygen supply, our substrate field (O) represents oxygen.We could also use it to represent a key nutrient to model anaerobic bacteria growing in a nutrient-limited environment.Our single-species biofilm model includes four types of generalized cell, representing: biofilm cells (b), a substrate source (o), a solid bottom (s), and an aqueous medium (m) submerging the growing biofilm. 13The substrate-source and solid-bottom cells are frozen, i.e. they do not exchange pixels with other cells.Fields do not have an excluded volume, i.e. they can co-occupy space with cells.We treat the medium as a large generalized cell without a volume constraint.We can regard the medium as all space in the cell lattice that is not occupied by other generalized cells.Treating the medium as a generalized cell allows us to represent the surface tension of real cells.Figure 3 shows the initial 2D configuration for our simulations. 14The dimension of the domain in pixels is 200 × 100 × 1 (the z-axis is perpendicular to the page).We are free to assign the length scale and time scale in the simulation, which fixes the other parameters.We set 1 pixel to 4µm, so the vertical height of the simulation domain is 400µm, as in [66], and 1 MCS to 3000s, so 1000 MCS corresponds to approximately 35 days.Each biofilm cell initially occupies a 2 × 2 × 1 cube, which is approximately 10 times larger than the real bacterial cells. 15The initial value of V t for biofilm cells is equal to their initial volume.We simulate growth of biofilm cells by applying equation (8) in the equivalent form: In the hybrid PLH model [66,67], cells divide when the local value of C reaches 1.In the GGH model, biofilm cells divide when V reaches the doubling volume V d , which is twice their initial V t .Therefore, both C and V t represent the closeness of cells to division, which justifies equation (15).We can rewrite equation (15) in a form which explicitly depends on G: where Equations ( 3) and (4) give Since c X and V t represent the same object, they are proportional to each other: We rewrite equation ( 18): which shows k's dependence on G.In our simulations, we simplify equation ( 21) by setting V t = 3V d /4, which corresponds to the average cell volume.We also ignore cell loss by setting m S = 0, 16 and do not model cell death.Equation ( 21) has two limiting cases: αS G 1 and αS G 1. In the first case, k does not depend on S, so the second term in the right-hand side of equation ( 20) is proportional to S, which represents substrate decay.In the second case, the decay constant k is proportional to G S , so the second term in the right-hand side of equation ( 20) is independent of S, which represents substrate consumption.We effectively combine both limiting cases by rewriting equation (20) as where and The substrate-source cells, representing the contact of the liquid medium with oxygen in the air, secrete substrate, which diffuses into the medium and biofilm cells.Secretion occurs at a constant rate, and we keep the substrate concentration S everywhere below a threshold value 1 (corresponding to the maximum solubility of substrate in water) by checking whether S > 1 anywhere.If so, we set S = 1 at those points.A constant secretion rate for substrate would implement a constantflux boundary.However, the requirement S ≤ 1 yields a constant-concentration boundary, in agreement with a constant oxygen partial pressure in the surroundings.This algorithm means that the net supply of substrate increases as its consumption by biofilm cells increases and also increases as the separation between the biofilm and the liquid surface decreases, which is realistic for biofilms growing in thin liquid films, where the substrate concentration at the surface is always saturated.It does not describe well a biofilm growing deep under water, where the maximum substrate concentration is well below saturation.The substrate distribution is initially homogeneous and saturated, so S = 1 everywhere.Substrate-secreting generalized cells increase the concentration of substrate at a fixed rate in the pixel corresponding to their center of mass.Biofilm cells consume substrate by secreting negative amounts of substrate at their centers of mass. 17Substrate diffuses on the field lattice using a forward-Euler algorithm run N times per MCS, where N depends on the time and distance scales we have chosen, and the substrate diffusion constant, D S .We use no-flux boundary conditions at the bottom of the simulation domain, and periodic boundary conditions at the sides to simulate a large biofilm.
We must specify all parameters in our simulations in GGH units: pixels, MCS,....In the absence of experimental measurements of these parameters, we choose them in a range which exhibits the general classes of behaviors observed experimentally.The fluctuation amplitude T rescales the surface-adhesion coefficients J(τ, τ ) and the volume-constraint coefficient λ V .Here, we choose these values to prevent cells from disappearing or freezing, and use the hierarchy of surface-adhesion coefficients to keep the cells from desorting (i.e. to keep the bacterial cells from floating off into the liquid medium).Thus, for each pair of generalized-cell types τ and τ we set J(τ, τ ) ≤ J(τ, τ ) and J(τ , τ ) ≤ J(τ, τ ).For the biofilm cells we also have 12 is on the order of the approximate critical surface tension below which cell dissociation occurs (the factor of 12 is the number of neighbors in the 2D third-order neighborhood), which keeps them from dispersing into the medium [34].In the GGH model, stronger adhesion corresponds to lower values of the surface energy J. Since bacterial cells are viscoelastic, we could, in principle, obtain their surface adhesion coefficients from parallel-plate compression experiments [10,29,30,31].
Unless we specify otherwise, we use the following values for the cell-cell and cell-medium adhesion coefficients, which satisfy the above inequalities [35]: The other J(τ, τ ) are irrelevant, so we set them to zero.We also set T = 10 following [34,37].Our choice of parameters is somewhat arbitrary, since only experiments can provide their exact relative values.We choose λ V = 40, which prevents the cells from disappearing and makes them fairly rigid.We take the parameters in equations (1)-( 8) from [66], except for m S , which we set to 0. 18 We control the value of the non-dimensional parameter G in equation ( 9) by changing c S0 .

4.
Implementing the GGH mathematical model in CompuCell3D.Now that we have defined our mathematical model of biofilm growth, we must translate it into a simulation.Our simulations use CompuCell3D [15,16,17,18,44], an open-source, multimodel software framework for simulation of the development of multicellular organisms. 19The great advantage of the CompuCell3D implementation of the GGH is that it uses simple CompuCell3D markup language (CC3DML) to specify models, allowing easy model sharing and validation [6].CC3DML is simple to read and write and fairly compact.
To illustrate the translation of our model into CC3DML code, we review the derivation and meaning of the most relevant parts of the model.We provide the complete CC3DML code corresponding to our simulation with G = 10 in Appendix A; all the CC3DML code, the plugins and steppables, and the PIF and TXT files used in our simulations are available for download from http://www.compucell3d.org/Models/biofilm. 18 If we take m S in equation ( 4) equal to its original value in [66], the simulated biofilm detaches from the bottom support (simulations not shown) which is unrealistic in the absence of fluid flows.For m S = 0, our simulated biofilms remain attached but are denser at the bottom than the corresponding simulations in [66]. 19Downloadable from http://www.compucell3d.org/.
The first section of a CC3DML configuration file defines the global parameters of the lattice and the simulation:  14), T = 10.The line Flip2DimRatio 0.01 /Flip2DimRatio tells Com-puCell3D that it should make 0.01 times the number of lattice sites (199 × 99 × 1) index-change attempts per MCS.The purpose of this definition is to model processes with fast time scales, in our case allowing for flexible definition of the diffusion constant of the substrate.In our case, the diffusion constant of the substrate is 100 times larger than the nominal value in the code (see below).The next line FlipNeighborMaxDistance 2.1 /FlipNeighborMaxDistance specifies the range of interactions (in pixels) to be third neighbor.A value of 1 corresponds to nearestneighbor interactions.The longer the interaction range, the more isotropic the simulation and the slower it runs.If the interaction range is comparable to the cell size, apparently nonadjacent cells may contact each other.The CC3DML tags Boundary − x , Boundary − y and Boundary − z impose boundary conditions on the lattice.In our case, the x axis is periodic; so, for example, pixel (0, 1, 1) neighbors pixel (199, 1, 1), while the y and z axis boundaries have no-flux boundary conditions.
The second section of the CC3DML file contains plugins, which CompuCell3D refers to every time it calculates the effective energy in equation (10).The plugins vary from simulation to simulation and define many of the types, behaviors and interactions of objects used in a given simulation.The ability to specify dynamically plugins (and steppables, see below) gives CompuCell3D its flexibility.
The CellType plugin informs CompuCell3D what generalized-cell types we are using in the simulation: Plugin Name="CellType" CellType TypeName="Medium" TypeId="0"/ CellType TypeName="Oxygen" TypeId="1" Freeze=""/ CellType TypeName="Biofilm" TypeId="2"/ CellType TypeName="Support" TypeId="3" Freeze=""/ /Plugin .Each line contains the name of a type that the simulation uses and assigns it to an integer-valued T ypeId.Medium is traditionally assigned a T ypeId = 0 (this type index does not have a volume constraint by construction).In our simulations we create four cell types, Medium, Oxygen, Biofilm, and Support.Oxygen cells represent air above the liquid Medium in which Biofilm cells grow.Like Support cells, they are frozen; that is, they do not move by exchanging pixels with other generalized cells.
The Contact plugin defines the boundary energies J(τ, τ ) between cells of different types and the interaction range (Depth) of the neighborhood used in the contact-energy summation in equation ( 10): Plugin Name="Contact" Energy Type1="Biofilm" Type2="Medium" 16 /Energy Energy Type1="Biofilm" Type2="Biofilm" 8 /Energy Energy Type1="Biofilm" Type2="Support" 4 /Energy ... Depth 2.1 /Depth /Plugin .The VolumeLocal plugin defines the volume-constraint term in equation ( 10): Plugin Name="VolumeLocal" LambdaVolume 40 /LambdaVolume /Plugin .The CenterOfMass plugin enables CompuCell3D to track the center of mass of each cell, for example, to control cell growth (see below): Plugin Name="CenterOfMass"/ .The Mitosis plugin defines the doubling volume at which a Biofilm cell divides into two cells with equal volumes: Plugin Name="Mitosis" DoublingVolume 8 /DoublingVolume /Plugin .The NeighborTracker plugin enables CompuCell3D to track all neighbors of each cell: Plugin Name="NeighborTracker"/ .The third section of the CC3DML file contains steppables.The PIFInitializer steppable, executed only once at the beginning of a simulation, defines the initial conditions of the cell lattice: Steppable Type="PIFInitializer" PIFName Biofilm.PIF /PIFName /Steppable .A PIF file sets the initial cell lattice configuration for the simulations, and for n generalized cells contains at least n lines, each of which consists of: the index of a generalized cell, the cell type, and the range (in pixels) this cell occupies on the lattice (x min , x max , y min , y max , z min , z max ).Using multiple lines/cell allows specification of arbitrary geometries.In our simulations we begin with 50 bacterial clumps distributed uniformly along the bottom (see Figure 3).
CompuCell3D executes the remaining steppables once at the conclusion of every MCS.The FlexibleDiffusionSolverFE steppable introduces fields of chemicals and defines their secretion, consumption, diffusion and decay, in this case, implementing equation ( 18 22) and the SecretionData constant for Biofilm cells to k 2 .The Oxygen.TXT file sets the initial distribution C(x, y, z) of the substrate, with the format: x, y, z, C(x, y, z).In our simulations the initial field is 1 everywhere, representing oxygen saturation in the medium.

5.
Results.In this section, we explore how the morphology of the growing simulated biofilm depends on model parameters, in particular, the dimensionless parameter G in equation ( 9), which Picioreanu, van Loosdrecht, and Heijnen identified as controlling the selection of flat-interface vs. fingered growth in their hybrid model [66,67].Figure 4 shows the terminology we use in this paper to describe the morphology of simulated biofilms and the color code for substrate concentration in Figures 5-16.The leading edge is the vertical position of the highest point in the biofilm.The trailing edge is the average position of the narrowest part of the fingers, or, if fingers are convex, the lowest point occupied by medium.The small white square shows the generalized-cell size.In the following figures, we represent areas containing biofilm cells with white hatch lines, and the biofilm-medium interface as a white boundary line.We do not show the boundaries of individual biofilm cells because they are too small on the scale of the figure.We also show the spatial distribution of the substrate: red-high concentrations (S ∼ 1), blue-low concentrations (S ∼ 0.1), black-almost no substrate (S ∼ 0).For G < 10, the biomass growth rate limits the growth dynamics (growth-limited regime).The substrate penetrates most of the biofilm and reaches most cells.In this case, the simulated biofilms are dense, compact and fast-growing.Growth continues even at the bottom of the valleys between protrusions, resulting in ripples rather than fingers.For G = 1, diffusion is so strong that no valleys form and the biofilmmedium interface remains flat.The isoconcentration lines of the substrate remain essentially horizontal.For G = 2 and G = 5, the biofilm-medium interface becomes rippled.The isoconcentration lines of the substrate curve slightly, following the   biofilm-medium boundary.Figures 5, 6, and 7 show simulations with G = 1, 2, and 5, respectively.Figure 8 shows the time evolution of a simulated biofilm for G = 10.For this value of G, the substrate penetrates several cell diameters into the surface layer of the biofilm.The case G = 10, as we will see below, is a border between two different morphological regimes.As in Figure 9, adapted from [66], the biofilm's structural complexity increases continuously in time.Initially, the substrate penetrates   28)).See equation ( 25) and text for other parameter values.The ripples turn into wide convex fingers (growth-limited regime).Colors represent substrate concentration; see Figure 4(b).
throughout the biofilm colonies, which grow in all directions.As they grow, the cells consume substrate and the substrate concentration develops a gradient, increasing in the vertical direction.As the biofilm thickens, cells in regions which thicken faster due to statistical fluctuations (protrusions) experience higher concentrations of substrate than others.These cells produce new biomass more quickly, while the cells in the valleys (horizontal locations where the interface between biofilm and medium lags significantly behind the highest local vertical position of the biofilm, as in viscous fingering [21,38,60,64]) experience an environment depleted in substrate, begin to starve, and grow slowly or cease to grow, so the space between  growing biofilm clusters fills slowly or not at all with new cells.The valleys deepen and the biofilm becomes more irregular.The isoconcentration lines of the substrate are more curved than for G = 5.This fingering instability generates the fingers observed in experimental biofilms.
For G > 10, we find a substrate-transport-limited regime: the substrate transfer rate, not the biomass growth rate, limits the growth dynamics.The substrate penetrates thin layers of biofilm, but less substrate reaches deeper layers.Therefore, the valleys between biofilm protrusions are wider and deeper and the structure of the biofilm is more finger-like than for G = 10.Cells in the valleys cease to grow and the fingers develop concave necks.The biofilm fingers grow more slowly than for smaller values of G (reaching the top of the simulation domain after a larger number of MCS).The bigger G, the narrower the fingers and the greater the space between them.The neck of the finger also becomes narrower with respect to the width of the finger.Figures 10 and 11 show simulations for G = 20 and G = 50, respectively.At G = 50, the substrate consumption rate relative to biofilm growth is so large that more cells starve (they do not die in our model, but they do stop dividing) and the simulated biofilm produces fewer fingers.As a result, the isoconcentration lines of the substrate curve significantly, because high concentrations of substrate penetrate into the larger valleys.
Figure 12 shows the structures of biofilms and the spatial distributions of substrate for various values of G at the moment when their leading-edge height h reaches 7/8 of the vertical height of the simulation domain (350µm = 86pixels).Our simulated biofilm morphology resembles that in the PLH model, as shown in Figure 13.The substrate penetration length l p describes how far the substrate diffuses into the biofilm before decaying or being consumed: where k is the effective decay constant in equation (21).Note that not all the space covered by white hatch lines (biofilm) is black (near-zero substrate concentration) and the width of the blue layer within the biofilm shows the substrate penetration length, clearly showing that this length is a source of the fingering instability.Our simulations, except for G = 50, agree with the results of the PLH model [66] in both morphology and growth time scale.We do not report detailed numerical comparisons of the morphologies in this paper, which focuses on validation of our computational methods.Instead, we show the general compatibility of GGH and  PLH results and differences which we can attribute to the change in simulation method.Table 1 compares the times (in days) at which our simulated biofilms reach a thickness of 350µm with those of [66].For G = 50, our simulations produce biofilms with fewer fingers than those of [66].This difference results from our including adhesion between biofilm cells.This adhesion generates an effective surface tension, Table 1.The age of the simulated biofilm (in days) when its leading edge reaches a thickness of 350µm, for different values of G, in our simulations (with J(b, m) = 16 and J(b, m) = 8) and in the PLH model [66].which, as in viscous fingering [64], suppresses high-frequency modes in the biofilmmedium boundary, and prevents separate thin structures from forming very close to each other, generating fewer, broader fingers.
Figure 14 shows the effect on biofilm morphology of decreasing J(b, m) from 16 to 8 with G = 10 (or γ(b, m) from 12 to 4), which is equivalent to decreasing the strength of the biofilm's cell-cell adhesion, so that the cells adhere to the medium and to themselves with the same strength (if we decrease J(b, m) to 4, the surface tension of the cell-medium interface, γ(b, m), reaches zero).More fingers form, showing that the effective cell-medium surface tension increases the wavelength of the finger pattern, as occurs in viscous-fluid fingering instabilities [64].Valleys pinch off, leaving pockets of medium surrounded by biofilm as in directional solidification.Figures 15 and 16 show the same effect for biofilms with G = 20 and 50, respectively.The smaller G, the less sensitive the simulated morphology to the strength of cellcell adhesion.Weaker cell-cell adhesion produces fingers with more irregular shapes and rougher surfaces, as we would expect when we eliminate the smoothing effects of surface tension.However, the parameter G still governs the pattern wavelength and overall morphology.
Fractional occupancy (fraction of sites at height y occupied by biofilm cells) c(y) is a convenient parameter to describe biofilm morphology, since it shows the effective average finger shape.Figure 17 shows c(y) at the moment when the biofilm's leading edge reaches 7/8 of the vertical height of the simulation domain (350µm = 86pixels), and compares it to the corresponding results in the PLH model, adapted from [66].Both models agree in curve shapes for smaller values of G, while for larger values, G = 20 and 50, our curves cease to be monotonic, showing that the fingers have concave necks.The differences in positions of curves between the two models that appear at low values of G are not apparent from the biofilm-structure figures.This difference in fractional occupancy occurs because the GGH biofilms always have a base layer of bacterial cells that coat the bottom in a continuous layer.Figure 14 shows also how the sensitivity of the simulated morphology to cellmedium surface tension decreases with decreasing G.The capillary length where γ is the surface tension and ν the kinematic viscosity, defines the critical length below which small structures are suppressed and gives a smoothing rate.If the dimensionless ratio is large, the biofilm structures are thick and smooth.If R cp is small, they are narrow and rough.The exact definition of kinematic viscosity in the GGH model is still debated.Since we are concerned only with relative values of R cp , we take ν = 1 in all of our calculations.To calculate the effective decay constant k in equation (21) we use α = 304 from the PLH model [66], and set S = 0.01 which is on the order of the substrate concentration near the biofilm-medium interface.The profiles for J(b, m) = 8 (R cp ∼ 0.4) are more consistent with the corresponding PLH profiles in [66] than are those for J(b, m) = 16 (R cp ∼ 1), which is the expected result, since the PLH model [66] does not contain surface tension.Further decreasing J(b, m) or increasing J(b, b) produces less stable fingers and, if γ(b, m) < γ c = 5/6, leads to dispersion of biofilm cells into the medium.Overall, however, changes in cell-medium surface tension and hence R cp affect simulated biofilm morphology less than changes in the growth-to-transport ratio G.We leave a detailed exploration of the effects of R cp to a future paper.
Figure 18 shows the leading-edge height h (in µm) of the simulated biofilm as a function of time for different values of G.We measure until the leading edge of the simulated biofilm reaches a height corresponding to 350µm.Initially, the curves in Figure 18 are convex, which corresponds to the initial exponential biofilm growth in a medium with an unlimited supply of substrate.Growth then slows as cells consume substrate, reducing its local concentration and thus local cell-growth rates.Eventually, parts of the biofilm reach the higher concentrations of substrate near the top of the simulated volume and leading-edge growth accelerates again.We expect this acceleration because the net flux of substrate, and hence the maximum possible growth rate, increases as the biofilm nears the surface, which increases the steepness of the substrate gradient.6. Discussion. Figure 19 shows the morphological regimes of biofilms in the GGH simulation.In general, both qualitative and quantitative results of the GGH and PLH models are consistent (see Figures 12 and 13), except for G = 50 where we obtain fewer, smoother fingers than the PLH model, whose fingers look like cacti [66].In our simulations, we restricted parameter variations to G (the only parameter varied in [66]) and J(b, m) and measured only the fractional occupancy c for our simulated biofilms, leaving other quantities studied in [66] for future papers.The difference in fractional occupancy between the two models is significant for high values of G, but also appears for low values of G.This difference results from the base layer of biofilm cells coating the bottom in the GGH biofilms.
At least part of the difference for high G results from our including adhesion between biofilm cells, which creates a cell-medium surface tension, which prevents separate thin structures from forming very close to each other.Decreasing the surface tension increases the number of fingers and their roughness, as expected, making our simulations for high G more similar to those in [66], although we do not reproduce a cactus pattern (characteristic of biofilms with high finger-like structures and without a flat basal layer of cells covering the solid support, like Acinetobacter biofilms).The GGH simulation also shows a new morphology, not found in the PLH model [66], for G ∼ 10 and small surface tension, in which valleys tend to pinch off (see Figure 14). 20We need to repeat our simulations for generalized cells of the size of real biofilm cells to check if this morphology is a feature specific to cell-cell adhesion.
The lack of cacti in the GGH simulation, which is the main morphological difference between the GGH and PLH biofilm models, partially results from the inclusion of surface tension in the GGH model.However, when we decreased the strength of the biofilm's cell-cell adhesion to J(b, m) = 4 (when the cell-medium surface tension is zero), cacti still did not form and the simulated morphology did not change much, except that some cells dispersed into the medium (results not shown).In this case, the suppression of the very-short-wavelength instabilities seen in the PLH model may result from the cell length scale itself.By construction, the GGH model cannot form structures with length scales smaller than a cell diameter.In this case, this length cutoff eliminates morphologies intrinsic to continuum and hybrid models, as discussed in [56].Since our simulated cells are ten times the size of real cells, we need to check how this smoothing changes as we change the simulated cell size.We will report studies of the effects of the cell size in future papers.We also need to explore the role of cell adhesion in determining GGH biofilm morphologies by running simulations for more values of J(b, m).
Another difference between the biofilm structures in the GGH and PLH models is that all GGH biofilms have a flat base layer of biofilm, on top of which the fingers form.This base layer is absent in the PLH model at higher values of G and the  cacti grow up from an initially empty plane.All our simulations have a base layer of biofilm cells covering the supporting surface.Since cell loss is irreversible in the GGH model, we did not set m S in equation ( 4) equal to its original value in [66] (for which GGH-simulated biofilms detach from the bottom support, see footnote 18), but set m S = 0.As a result, our simulated biofilms remained attached but were denser at the bottom (having a base layer) than the corresponding simulations in [66].
Since mass transfer, as we have confirmed in this paper, is a key determinant of biofilm morphology, we must consider the effects of boundary conditions when we compare GGH and PLH simulation results.The PLH simulations used a constantconcentration boundary at a fixed distance from the leading edge of the growing biofilm, while our GGH simulations use a constant-concentration boundary at a fixed position (which greatly simplifies the computation).Thus the distance between the source of the substrate and the leading edge of the biofilm decreases in time in our GGH simulations, while in the PLH model the distance between the substrate source and the trailing edge of the biofilm increases in time.Since the biofilm patterns in both models are nonstationary and the definition of the length scale l z is somewhat arbitrary, experimental conditions should determine which of these (or another choice entirely) is most appropriate.Since l z , the thickness of the fluid layer, enters directly into G (equation ( 9)), changing the z-dimension of the simulation, while keeping all other parameters constant, changes both the growth rate and the biofilm morphology, as Figure 20 shows.However, if we vary l z and keep G constant by changing c S0 , the morphology does not change greatly, as Figure 21 shows.Thus the net effect of the difference in boundary conditions is small and should not significantly change the phase diagram in Figure 19.The greater distance between the leading edge of the growing biofilm and the constantconcentration boundary at early times in the GGH model may, however, explain more subtle differences from PLH results, like the greater spacing between fingers which we observe in our GGH simulations for a given value of G.
This paper is a computational methods paper, where our goal is to show how to translate the PLH biofilm model into the GGH model using CompuCell3D.For this reason, we chose a generalized cell size of ten times the size of biological cells, which significantly sped up our simulations.For the same reason we did not discuss  While several research groups study the fundamental mechanisms controlling the structure and growth dynamics of bacterial biofilms, including the Montana State University Center for Biofilm Engineering and the biofilm group at the Delft University of Technology, as far as we know, none has measured the typical values of the growth group parameter G in equation ( 9) experimentally.However, experiments have observed biofilm morphologies that correspond to all those we obtain in our simulations with 1 < G < 50.We lack experimental measurements of the kinetics of the thickness of the growing biofilm (our Figure 18).Experimental tests of whether biofilm morphologies and dynamics depend on G would be valuable.7. Summary.We have shown that the GGH model can successfully simulate the growth of a single-species biofilm on a solid surface in a liquid medium, reproducing qualitatively the results of previous continuum and hybrid models.Our simple 2D, single-species biofilm model includes cell growth, proliferation, and division, diffusion of substrate, and cell-cell and cell-medium adhesion.We neglected the kinetics of initial cell attachment, recruiting of cells from the medium, breaking off of chunks of the biofilm and cell death.We coarse-grained to a cell size approximately ten times the real value.Instead of modeling ECM explicitly, we treated cells and their ECM surface layer as a generalized cell and used the generalized-cell-generalized-cell adhesion coefficients to represent the aggregate interactions among cells and ECM.Our GGH simulations using the kinetics of the PLH model [67] reproduced the structures obtained in [66], and showed that biofilm size and morphology depend mainly on the non-dimensional parameter G, the ratio of the maximum biomass growth rate to the maximum substrate transport rate.
Cell-cell adhesion significantly affects growth and morphology for larger values of G, in the substrate-transport-limited regime, since adhesion produces additional surface-tension forces which smooth the biofilm.Finite cell size has similar effects, suppressing high-frequency instabilities.Surface adhesion will be even more important in multispecies biofilms, where adhesion will vary greatly within and between species.
The great appeal of the GGH model is its extensibility.Adding new biological mechanisms that depend on cell-cell and cell-field interactions is as easy as adding new potential energies or constraints to the GGH effective energy.For example, including variable adhesion among different bacterial species in CompuCell3D requires only a few lines of CC3DML code, making CompuCell3D a natural framework for simulating multispecies biofilms.Each CompuCell3D simulation presented in this paper took at most a few hours, using an earlier version of CompuCell3D.Recent improvements to CompuCell3D have increased its speed by a factor of ten.Steppable Type="TargetVolume" InitialTargetVolume Type="Biofilm" 4 /InitialTargetVolume Growth Type="Biofilm" Field="Oxygen" GrowthRate="1" Saturation="0.033"/ /Steppable /CompuCell3D

Figure 2 .
Figure 2. A typical GGH configuration in 2D showing portions of nine generalized cells.The numerals indicate index values.The colors indicate cell types.A generalized cell is a collection of connected lattice sites (squares) with the same index value.The number of lattice sites in a cell is its volume and the number of lattice links on its boundary (interfaces with other indices) is its surface area.

Figure 3 .
Figure 3.The initial lattice configuration, with biofilm cells (colonists) at the bottom, sitting on a flat solid, and generalized cells producing substrate at the top.The space between the cells is aqueous medium with an initially homogeneous concentration of substrate.The vertical height corresponds to 400µm and the horizontal length to 800µm.The inset inside the simulation domain shows an enlarged fragment of the solid bottom and a few biofilm cells.
CompuCell3D Potts Dimensions x="200" y="100" z="1"/ Steps 20000 /Steps Temperature 10 /Temperature Flip2DimRatio 0.01 /Flip2DimRatio Boundary − x Periodic /Boundary − x Boundary − y NoFlux /Boundary − y Boundary − z NoFlux /Boundary − z FlipNeighborMaxDistance 2.1 /FlipNeighborMaxDistance /Potts .The line Dimensions x="200" y="100" z="1"/ declares the dimensions of the lattice to be 200 x 100 x 1, i.e. the lattice is 2D and extends in the xy plane.The basis of the lattice is 0 in each direction, so the lattice sites in the x and y directions have indices ranging from 0 to 199 and from 0 to 99, respectively.The lines Steps 20000 /Steps and Temperature 10 /Temperature tell Compu-Cell3D that the simulation should last 20000 MCS, with the fluctuation amplitude in equation (

Figure 4 .
Figure 4. (a) Terminology describing the morphology of simulated biofilms.(b) Color code for substrate concentrations in Figures 5-16.

(a) 80 Figure 8 .
Figure 8. Convex-finger biofilm (white hatch lines) growth in a 2D biofilm simulation with G = 10, J(b, m) = 16 and R cp = 1.06 (see equation (28)).See equation (25) and text for other parameter values.The fingers are usually convex and the valleys do not pinch off.Colors represent substrate concentration; see Figure 4(b).

Figure 12 .
Figure 12.Biofilm (white hatch lines) morphologies and spatial distributions of substrate as a function of G in the GGH model, all with J(b, m) = 16 (see equation (25) and text for other parameter values), observed when the leading edge of the simulated biofilms reach 7/8 of the vertical height of the simulation domain.Colors represent substrate concentration; see Figure 4(b).

Figure 13 .
Figure 13.Biofilm morphologies as a function of G in the 2D PLH model.Adapted from [66].

(a) 80 Figure 14 .
Figure 14.Rough-finger biofilm (white hatch lines) growth in a 2D biofilm simulation with G = 10, J(b, m) = 8 and R cp = 0.35 (see equation (28)).See equation (25) and text for other parameter values.The weaker cell-medium surface tension results in more, narrower, fingers, with relatively narrower valleys, than in the corresponding high-surface-tension G = 10 simulation in Figure 8.The valleys pinch off.Note the greater irregularity of the finger shapes and their rougher surfaces.Colors represent substrate concentration; see Figure 4(b).

Figure 15 .
Figure 15.Rough-finger biofilm (white hatch lines) growth in a 2D biofilm simulation with G = 20, J(b, m) = 8 and R cp = 0.37 (see equation (28)).See equation (25) and text for other parameter values.The weaker cell-medium surface tension results in more, narrower, fingers, with relatively narrower valleys, than in the corresponding high-surface-tension G = 20 simulation in Figure 10.Note the greater irregularity of the finger shapes and their rougher surfaces.Colors represent substrate concentration; see Figure 4(b).

Figure 16 .
Figure 16.Concave rough-finger biofilm (white hatch lines) growth in a 2D biofilm simulation with G = 50, J(b, m) = 8 and R cp = 0.39 (see equation (28)).See equation (25) and text for other parameter values.The weaker cell-medium surface tension results in more, narrower, fingers, with relatively narrower valleys, than in the corresponding high-surface-tension G = 50 simulation in Figure 11.Note the greater irregularity of the finger shapes and their rougher surfaces.Colors represent substrate concentration; see Figure 4(b).

50 Figure 17 .
Figure 17.Averaged finger shapes, shown by the fractional occupancy c as a function of the vertical coordinate y for different G, observed when the leading edge of the simulated biofilm reaches 7/8 of the vertical height of the simulation domain.We have rescaled the data so that the relative position y of the top of the finger (the leading edge) is 1 and the trailing edge has y = 0. Squares: J(b, m) = 16.Diamonds: J(b, m) = 8.Rectangles: the corresponding data in the PLH model, from [66].

Figure 18 .
Figure 18.The leading-edge height h of simulated biofilms as a function of time for different values of G. Initial exponential biofilm growth in medium with an unlimited supply of substrate slows as cells consume substrate and reduce its local concentration.When parts of the biofilm reach the higher concentrations of substrate near the top of the simulation domain, leading-edge growth accelerates again.

Figure 19 .
Figure 19.Schematic representation of the morphological regimes of GGH-simulated biofilms.The borderline between ripples and fingers is at about G = 10.Valleys between convex fingers pinch off at low surface tensions.Valleys between concave fingers do not pinch off in our simulations.Concave fingers become rough at low surface tensions.

Figure 20 .
Figure 20.Biofilm growth for different vertical heights of the simulation domain and thus different G.All other parameters as in Figure 8.

Figure 21 .
Figure 21.Biofilm growth for different vertical heights of the simulation domain with fixed G = 10, achieved by setting c S0 = c S0