A modeling and simulation study of siderophore mediated antagonism in dualspecies biofilms
 Hermann J Eberl^{1}Email author and
 Shannon Collinson^{2}
DOI: 10.1186/17424682630
© Eberl and Collinson; licensee BioMed Central Ltd. 2009
Received: 14 October 2009
Accepted: 22 December 2009
Published: 22 December 2009
Abstract
Background
Several bacterial species possess chelation mechanisms that allow them to scavenge iron from the environment under conditions of limitation. To this end they produce siderophores that bind the iron and make it available to the cells later on, while rendering it unavailable to other organisms. The phenomenon of siderophore mediated antagonism has been studied to some extent for suspended populations where it was found that the chelation ability provides a growth advantage over species that do not have this possibility. However, most bacteria live in biofilm communities. In particular Pseudomonas fluorescens and Pseudomonas putida, the species that have been used in most experimental studies of the phenomenon, are known to be prolific biofilm formers, but only very few experimental studies of iron chelation have been published to date for the biofilm setting. We address this question in the present study.
Methods
Based on a previously introduced model of iron chelation and an existing model of biofilm growth we formulate a model for iron chelation and competition in dual species biofilms. This leads to a highly nonlinear system of partial differential equations which is studied in computer simulation experiments.
Conclusions
(i) Siderophore production can give a growth advantage also in the biofilm setting, (ii) diffusion facilitates and emphasizes this growth advantage, (iii) the magnitude of the growth advantage can also depend on the initial inoculation of the substratum, (iv) a new mass transfer boundary condition was derived that allows to a priori control the expect the expected average thickness of the biofilm in terms of the model parameters.
Background
With but few exceptions, iron is absolutely required for life of all forms, including bacteria. It plays an important role in many biological processes, such as methanogenesis, respiration, oxygen transport, gene regulation and DNA biosynthesis [1]. Iron is abundant in the Earth. However, while in the early ages of life the predominant form of iron was rather soluble, it is now extremely insoluble and, therefore, the bioavailability of this minor nutrient is often low. To overcome iron limitations, some bacteria secrete ironchelation compounds (socalled siderophores) when the environmental iron concentration becomes small. These bind with iron to form a siderophoreiron complex, which is then taken up by the cells and the iron is later liberated internally. This enables the microorganisms to scavenge iron from the environment which, thus, becomes unavailable to other organisms, including hosts.
Under iron limitations, species that produce siderophores and, thus, chelate iron can have a competitive advantage over species that lack this ability [2]. Such siderophore mediated antagonism has been observed in agricultural microbiology [3–5] and in food microbiology for some food spoilage bacteria, e.g. in meat, fish, poultry and dairy [6–10]. In these environments nutrients are often available in abundance, while iron can become growth limiting. The siderophore mediated antagonism is inversely related to the availability of iron [4] in the medium (soil or food); it is not observed if and when iron is not limited [2]. The bacteria that most experimental studies of this phenomenon focus on are pseudomonads, primarily of the Pseudomonas fluorescens  P. putida group, which produce a yellowgreen (under UV light) pigment with high iron binding constant. This is the siderophore pyoverdine.
In the present study we focus on the antagonistic effect against other bacteria, as studied experimentally in [2], but the principle of growth suppression of other microorganisms by iron scavenging from the environment applies also to the control of yeasts; in a medical context this phenomenon has also been suggested as a mechanism to control cancer and other diseases. Because of their antagonistic effect, it is now generally recognized that plant pathogens with this property, in fact, can even have plant growth promoting effect by controlling wilt disease or other root crop diseases. Therefore, such PGPR (plant growth promoting rhizobacteria [5]) have been used for soil inoculation to increase yields.
The majority of experimental studies of iron chelation, as well as the population level modeling studies of pyoverdine production and iron chelation so far have been carried out for suspended cultures. Most bacteria, however, live in biofilm communities and not in suspended cultures. In particular the pseudomonads, which have been most commonly used in iron chelation studies are known to be natural and prolific biofilm formers. While there is increasing evidence that iron chelation can play an important role in biofilms [2, 11–13], no conclusive quantitative studies of siderophore mediated antagonism in biofilms have been conducted so far. Previous laboratory studies of this question in [2] remained inconclusive, because of the affinity of one of the strains involved in the study towards the reactor material. Since the interaction of population and resource dynamics in biofilm communities can be very different from suspended cultures [14], it cannot be answered by straightforward inference from the planktonic case whether or not siderophore production provides a growth advantage. We approach this question by developing a mathematical model, which is then studied in computer simulations. Using a theoretical approach, it becomes possible to focus on the effect at the center of the investigation, without adverse perturbations to which laboratory studies are susceptible, like the ones reported in [2].
Bacterial biofilms are microbial depositions on surfaces and interfaces in aqueous systems. Biofilms form after individual cells attach to the surface, called substratum in the biofilm literature, and begin to produce extracellular polymeric substances (EPS), which form a gellike layer in which the bacteria themselves are embedded. This polymeric layer offers protection against mechanical removal, but also against antimicrobials, that suspended bacteria do not have. One of the most striking differences between life in biofilms and in suspended cultures is that biofilm bacteria live in concentration gradients [14], due to decreased diffusion of dissolved substrates, the spatial organization of the cells, consumption and production of substrates, and biochemical reactions in the EPS matrix. This can lead to spatially structured populations with niches for specialists that cannot be found in suspension. For example, aerobic bacteria close to the biofilm/water interfaces can consume the oxygen in the environment and thus establish anaerobic zones in the deeper regions of the biofilm, closer to the substratum. Similarly, many antimicrobial agents only inactivate the bacteria closest to the biofilm/water interface but do not reach the cells in the deeper layer, which can survive an antibiotic attack virtually unharmed. In environmental systems biofilms are typically considered good, because their sorption and degradation properties contribute to soil and water remediation. Therefore, many environmental engineering technologies are based on biofilm processes, in particular in wastewater treatment, soil remediation, and groundwater protection. In industrial systems, biofilms are responsible for accelerated corrosion (microbially induced corrosion, biocorrosion) and biofouling. Biofilm contamination in food processing plants and hospitals are associated with public health risks [15–17]. In a medical context, biofilms can cause bacterial infections, which are diffiicult to treat with antibiotics, for the reasons indicated above. The list of biofilm originated diseases and infections is long and includes cystic fibrosis pneumonia, periodontitis and dental caries, and native valve endocarditis. A more detailed overview is given in [18]. In order to overcome the limitations of antibiotics in treating biofilm infections other strategies have been suggested recently, such as quorum sensing based methods [18, 19], or iron chelation based methods [12].
Mathematical models for bacterial biofilms have been used for several decades and they have greatly contributed to our understanding of biofilm processes so far. The first generation of biofilm models were continuum models with a focus on population and resource dynamics, formulated under the assumption that a biofilm can be described as a homogeneous layer, cf [20]. In reality, however, biofilms can develop in rather irregular structures, such as clusterandchannel architectures. Homogeneous biofilm layers are primarily obtained under conditions of abundance. Since we are interested in the iron chelation process, we are interested in situations of iron limitations. Therefore, a multidimensional biofilm model is required that supports the formation of clusterandchannel biofilm architectures. In the past decade several such models have been developed [20, 21]. The first group of these models, although utilizing a variety of different mathematical concepts, from individual stochastic based models to stochastic cellular automata, to deterministic continuum models, focused on biofilm growth, population and resource dynamics, i.e. on biofilm processes with typical time scales of days and weeks. This is what we need for our study. The second group of multidimensional models focuses on mechanical aspects of biofilms, such as biofilm deformation and eventual detachment, i.e. on processes on a much shorter time scale. Currently no biofilm model is known that connects both aspects reliably. Therefore, the latter processes are neglected in our model in the same manner as they are neglected in other biofilm growth models.
Mathematical Model
We develop a mathematical model of siderophore production and iron chelation in biofilms by combining the iron chelation model [22, 23], which was originally developed for batch cultures, with the densitydependent diffusion reaction model for biofilm formation that was originally introduced and studied for singlespecies biofilms, both for mathematical and biological interest, in [24–29] and extended to mixedculture systems in [30–32]. Our focus here is on the growth advantage of siderophore producing bacteria over bacteria that lack this ability. Therefore, we formulate the biofilm model for a mixed culture biofilm. A related modeling and simulation study for suspended populations in batch and chemostat like environments was recently conducted in [33], where it was found with a blend of analytical and computational techniques, that iron chelation abilities can greatly affect persistence results in chemostats. Mathematical models of biofilms render the complexity of biofgilm populations. They are essentially more complicated than mathematical models of suspended microbial populations and most mathematical techniques than can be used to study suspended populations cannot be used to study biofilms. In particular, biofilm models do not lend themselves easily to analytical studies but must be investigated in time intensive computer simulations.
Governing equations
Our biofilm formation model is formulated in terms of the dependent variables volume fraction occupied by the siderophore producing species, N, and volume fraction occupied by species that does not produce siderophores, R. We follow the usual approach of biofilm modeling and subsume the EPS that is produced by the bacteria in the biofilm volume fractions. The total volume fraction occupied by the biofilm is then M = N + R.
In our modeling study we focus on siderophore mediated antagonism. Therefore, we assume that iron availability is the only growth limiting factor for the development of the biofilm; all other required nutrients are assumed to be available in abundance. Moreover, we assume that the growth conditions in the medium are not altered by the iron dynamics. Under iron limitations, the chelator produces the siderophore pyoverdine, denoted by P, which binds dissolved iron S and makes it unavailable to the nonchelator. This transformation from dissolved iron S, to chelated iron, Q, is assumed to be 1:1. The dissolved iron diffuses in the aqueous phase and, at a lower rate, in the biofilm. The species R, which does not produce the siderophore, requires dissolved iron, S, for growth, while the siderophore producer's growth is controlled by the total of available iron, dissolved and chelated, S + Q. We assume that pyoverdine and chelated iron do not diffuse into the aqueous environment but are entrapped in the biofilm. The biofilm expands spatially, if the local cell density approaches the maximum cell density, i.e. if it fills up the available volume. It does not expand notably if locally space is available to accommodate new cells. This is described by a densitydependent diffusion mechanism, that shows two nonlinear diffusion effects [25, 34]: (i) it degenerates like the porous medium equation for vanishing biomass densities, and (ii) the diffusion coefficient blows up if the local cell density approaches its maximum value. Effect (i) causes the biofilm/water interface to spread at finite time, i.e., it guarantees a sharp interface between the biofilm and the surrounding aqueous phase. The super diffusion effect (ii) enforces volume filling, i.e. that the maximum cell density is never exceeded. Note that the interplay of both effects is necessary to describe biofilm growth.
Since pyoverdine and chelated iron are associated with the biofilm matrix we assume them to move at the same diffusive rate as the biofilm.
Model parameters used in this study
parameter  symbol  value  unit 

growth rate, chelator  μ _{1}  12.3  d ^{1} 
growth rate, nonchelator  μ _{2}  12.3  d ^{1} 
halfsaturation concentration, chelator  k _{1}  3.7  μM 
halfsaturation concentration, nonchelator  k _{2}  3.7  μM 
decay rate, chelator  d _{1}  0.49  d ^{1} 
decay rate, nonchelator  d _{2}  0.49  d ^{1} 
yield coefficient, chelator  Y _{1}  0.6003   
yield coefficient, nonchelator  Y _{2}  0.6003   
maximum biomass density, chelator  N _{∞}  10^{4}  gm ^{3} 
maximum biomass density, nonchelator  R _{∞}  10^{4}  gm ^{3} 
chelation rate  β  1.92 

pyoverdine production rate  δ  2.56  OD _{ P } d ^{1} 
pyoverdine inhibition concentration  S∞  0.3762  μM 
pyoverdine inhibition exponent  n  3   
biomass motility parameter  ϵ  10^{12}  m ^{2} d ^{1} 
biomass interface exponent  a  4   
biomass threshold exponent  b  4   
diffusion coefficent of iron in water  d(0)  8.64·10^{4}  m ^{2} d ^{1} 
diffusion coefficent of iron in biofilm  d(1)  7.776·10^{4}  m ^{2} d ^{1} 
The iron chelation reaction terms in the biofilm model (1) are a slightly generalized from those that have been proposed and identified for the suspended batch culture population model in [22]. In the latter the saturation that is described by Monod kinetics is not relevant for practical purposes since always Q ≪ k_{1} and S ≪ k_{1} after a very short initial transient phase. Therefore, first order reactions could be assumed.
This is not necessarily the case in the biofilm setting, depending on the amount of iron supplied to the system, where the iron concentration can be very different between locations close to the substratum and at the biofilm/water interface. Therefore, an extension of the model to Monod kinetics became necessary. Analytical results for densitydependent diffusionreaction models with degeneracy and super diffusion effects as implied by (3) can be found in [24, 27–29, 32, 35]. These include existence results and for singlespecies models also uniqueness results, as well as studies on long term behaviour and stability. The study [34] gives a derivation of this deterministic, fully continuous model from a discretespace model that is based on local behavioural rules similar to cellular automata models for biofilm growth, e.g [36–38]. Moreover, the underlying prototype biofilm model [25] can also be derived with hydrodynamic arguments similar to those used in [39] but under weaker assumptions, (cf Frederick et al, "A mathematical model of quorum sensing in patchy biofilm communities with slow background flow", submitted). Biological systems that have been previously described using this modeling approach include disinfection of biofilms with antibiotics [26, 35], competition between species for shared substrates [30, 40, 41], and amensalistic control [31, 32, 42].
Initial and boundary conditions
In order to close model (1) above, suitable initial and boundary conditions must be specified.
Initial conditions
In laboratory experiments the inoculation of the substratum, i.e. the sites at which the cells initially attach to the substratum, is difficult to control and appears random. In most of our simulations (except where noted) below, we will mimic this by choosing the actual sites of attachment at the substratum randomly. However, in order to ensure comparability across simulations we specify the initial number of colonies of both bacterial species as input data. Moreover, the volume fraction occupied by biomass in these inoculation sites is chosen randomly (uniform) between a given minimum and maximum value.
The initial biomass densities N and R are thus positive in the attachment sites on the substratum and 0 everywhere else. Initially, we choose a constant dissolved iron concentration S(0, ·) = S_{0} = 2 μM below the half saturation concentrations k_{1} and k_{2} but higher than the pyoverdine inhibition concentration S_{∞} that triggers the chelation process. Both, the concentration of chelated iron Q and the pyoverdine concentration P, are assumed to be 0 initially.
While (1) represents a completely deterministic model, this choice of inoculation adds a stochastic element.
It is naturally expected that different inoculation sites lead to different local biofilm morphologies and, hence, to different substrate distributions, but it is not clear a priori whether this also affects global, lumped results such as bacterial population sizes, mass conversion rates etc. For example, in [31] an amensalistic biofilm control strategy was investigated where the actual initial distribution of the control agent relative to the pathogen determines success or failure of the control strategy. Other studies, such as [26, 40] showed no or only little quantitative and no qualitative effect of inoculation sites on global measurements. The modeling studies of the impact of inoculation sites on biofilm processes conducted so far allow the conclusion that it depends on (i) the type of interaction between species (e.g., competition, amensalism), (ii) the response to limiting substrates (e.g., growth, disinfection), and (iii) density (vs. sparsity) of the inoculation. Since the effect of inoculation on the overall biofilm performance is a priori not clear, it is advisable to run simulation experiments in the form of trials with several replicas, cf also [40]. This is the approach that we take in this study.
Boundary conditions
A so far only unsatisfactorily solved, open problem in mesoscopic biofilm modeling is the specification of boundary conditions for the dependent variables. While it is relatively straightforward to prescribe boundary conditions for biomass and biomass associated components of the biofilm, formulating appropriate boundary conditions for dissolved, growth limiting substrates requires more thought.
where ∂_{ n }denotes the outer normal derivative at the boundary of the domain. This ensures that no biomass or biomass associated components leave or enter the domain across the boundary. For the part of the boundary of the domain that consists of the substratum this is the natural boundary condition. For the lateral boundaries these are symmetry conditions, which enable us to view the small simulation section as a part of a much larger system.
In other words, in order to obtain enough biomass for a noteworthy biofilm community, the domain Ω must be huge relative to the desired biofilm size. Otherwise, all iron will be immediately consumed before a biofilm can develop. Hence, since for computational reasons the domain size Ω must be restricted, the boundary conditions must include a mechanism that describes replenishment of the consumed substrate, even if it is expected to become limited eventually.
Usually this problem is dealt with by prescribing the concentration of the dissolved growth promoting substrate on some part of the boundary (Dirichlet condition), often the boundary opposing the substratum, while noflux conditions are specified everywhere else, which can be interpreted in the same manner as above for the biomass associated components. When the biofilm grows, the substrate concentration inside the domain decreases due to consumption. However, since under these boundary conditions the concentration is fixed along the Dirichlet boundary, this leads to an increased substrate gradient into the domain there, and, thus, to an increasing diffusive flux into the domain as the biofilm grows. Hence, if Dirichlet conditions are specified to model substrate replenishment, biofilm growth implies an increased supply of growth limiting substrate. Since we are here interested in studying biofilms under substrate limitations, which trigger the chelation mechanism, this is not appropriate for our application. In order to alleviate the effect of increasing substrate supply in response to biofilm growth we propose here two alternative boundary conditions to describe substrate replenishment.
Iron boundary condition I
at the substratum, y = 0, and the biofilm/water interface, y = L_{ f }. These boundary conditions have two parameters, the bulk concentration S_{0} and the thickness of the concentration boundary layer L_{ BL }, i.e., one parameter more than the traditional Dirichlet condition. Note that this concentration boundary layer is an abstract, not experimentally observable concept. It is qualitatively related to the bulk flow velocity, in the sense that a small bulk flow velocity implies a thick concentration boundary layer, while a thin concentration boundary layer represents fast bulk flow. However, a quantitative corelation between these two quantities is yet unknown [20].
Thus the boundary condition for iron is a mixed boundary condition consisting of a homogeneous Neumann boundary and a Robin boundary. Compared to the traditional Dirichlet boundary condition discussed above it has the effect that the growing biofilm not only lowers the substrate concentration inside the domain, but also on the boundary. While the diffusive flux into the system still increases with increasing biofilm size, it is bounded by d(0)S_{0}/L_{ BL }. In the case of the Dirichlet condition, on the other hand, it grows unbounded. Thus iron replenishment will be slower under (6) than under the usually used Dirichlet conditions.
Iron boundary condition II
where ∂Ω_{ S }denoted the part of the boundary of Ω that forms the substratum. The parameter λ is the average thickness that a completely compressed biofilm would have, i.e. a biofilm for which R ≡ 1 in Ω_{2}.
If, as in our simulations and in the vast majority of biofilm modeling studies in general, Ω is rectangular, and if the substrate flux is applied on the opposite side of the substratum, then the integral terms in (8) cancel out.
When using this boundary condition we will specify it in terms of λ, rather than the actual substrate flux. Unlike the previous boundary condition (6) and the more traditional Dirichlet boundary condition discussed above, the nonhomogeneous Neumann boundary condition allows us not only to estimate but to control the size that the biofilm will eventually have.
Note that (5) together with a boundary condition for S, such as (6) or (7) suffices. Since the solutions of the diffusionreaction system (1) are understood in the weak sense, no internal boundary conditions must be specified across the biofilm/water interface to close the model, which, however, are necessary for other biofilm models, such as [44].
Parameters
The model parameters used in this study are summarized in Table 1. In [22, 23] a set of model parameters of the chelation process was determined from laboratory experiments with batch cultures of Pseudomonas fluorescens. In the absence of measurements for the biofilm setting, this is also what we use here. The remaining parameters for the biofilm growth model were chosen in the usual parameter range, cf. [20] and [25]. In order to ensure that competition effects are entirely due to differences in the strains' ability to utilize chelated iron, we choose that both species have the same specific growth rate μ_{1} = μ_{2}, half saturation constant k_{1} = k_{2}, yield coefficient Y_{1} = Y_{2} and decay rate d_{1} = d_{2}. Thus, we assumed that X_{2} is a genetic modification of X_{1}, which switches off iron chelation but leaves the growth kinetics unaffected.
Computational realisation
The mathematical model (1) is discretized on a regular grid using an nonstandard finite difference scheme for time integration and a second order finite difference based finite volume discretization. This is a straightforward adaptation of the method that has been introduced in [43] for single species biofilms and extended to mixedculture systems in [31]. The main difference between (1) and other mixedculture applications of the nonlinear diffusionreaction biofilm model is that P and Q are controlled by the degeneratesingular diffusion operator, which, however, does not depend on P and Q directly. Thus, in the discretization these two equations behave essentially like semilinear equations which to incorporate into the simulation algorithm does not pose any new problems. In every timestep, five sparse linear systems need to be solved, one for each dependent variable. This is the computationally most expensive part of the simulation code and was prepared for parallel execution on multiprocessor/multicore computers using OpenMP; cf [41] for a more detailed discussion of this aspects, where this approach was applied to a dualspecies biofilm system that plays a role in groundwater protection. For the simulations presented here usually four threads were used on a SGI Altix 330 system. The visualisation of simulation results shown here were created using the Kitware ParaView visualisation package.
Numerical experiments
Simulations illustrate siderophore mediated growth advantage
where Z = 0 (only nonchelators, no siderophore producers) is depicted in yellow and Z = 1 (only siderophore producers, no nonchelator) in dark green.
The biofilm grows throughout the simulation experiment, despite the maximum concentration of dissolved iron being clearly smaller than the half saturation constant, i.e. despite growth limitations. The simulation starts out from twelve small initial colonies. As these colonies grow bigger they grow closer together and eventually neighboring colonies merge into bigger colonies. At t = 2d, we observe three mixedculture colonies, three clearly siderophore producer dominated colonies and one nonchelating colony. At t = 4d the nonchelating colony remains separated from the other colonies which now merge into two large mixedculture colonies, which at t = 6d merge into one large clearly siderophore producer dominated mixedculture colony. Also the nonchelating colony continues growing and the interfaces of the nonchelator and the mixedculture colony collide at the substratum. For t = 8d and t = 10d we notice siderophore producers slowly invading the nonchelating colony. While the larger chelator dominated biofilm colony keeps growing toward the iron source, i.e. the top boundary, the nonchelator colony cannot grow further due to a severe limitation of dissolved iron S.
Initially, S took the bulk concentration value S_{0} = 2.0 μM but continuously decreases due to biofilm growth. By the end of the simulation the maximum concentration of dissolved iron in the system (attained at the upper boundary, where the replenished iron enters the system) drops to S ≈ 0.23 μM. The iron concentration S is smaller in the chelator dominated colonies than in the nonchelator colony. In the larger chelator dominated biofilm colonies, the iron concentration S drops below S_{∞} and chelation starts. Thus, in addition to dissolved iron being directly consumed by chelators and nonchelators alike it is scavenged from the environment and transformed into chelated iron Q by the chelator. This leads to a diffusive flux of dissolved iron from the nonchelator colony into the chelator dominated colony. Hence iron does not only enter the mixedculture colony from the top boundary but also laterally.
Initially, up to t ≈ 1d, as long as iron is not limited, both species grow at about the same rate. After that, the growth of the species that does not produce siderophores lags behind the siderophore producer's growth, indicating the expected growth advantage. Eventually, at about t ≈ 4, the population that is not able to chelate declines, while the chelating population continues growing throughout the simulation, albeit at a decreased rate. The simulation stops at t = 10d, where, as indicated already before, all the nonchelator is accumulated in a single colony that is not yet notably invaded by the siderophore producer.
Simulations with controlled inoculation show that the effect of siderophore mediated antagonism is sensitive to initial attachment sites
In order to investigate the effect of the competition between siderophore producing and nonproducing species further we conduct a small simulation experiment, in which the initial biomass distribution is controlled in the following manner. Initially, the substratum is only inoculated by two colonies of identical, semispherical shape. One is situated at the left end of the simulation domain and one at the right end of the simulation domain.
The simulations are carried out on a grid of 300 × 200 cells covering a computational domain of size L × H = 0.75 mm × 0.5 mm.
We differentiate between the following four cases
(a) Two simulations are conducted. In one of them, both colonies are siderophore producers with an initial biomass density N_{0} = 0.3 (R_{0} = 0.0). In the second simulation both colonies are formed by the nonchelating species, R_{0} = 0.3 (N_{0} = 0). The concentration boundary layer thickness is set at L_{ BL }= 500 μm.
(b) The same as (a) but with a thicker concentration boundary layer L_{ BL }= 1000 μm, implying reduced rate of iron replenishment.
(c) A simulation in which one of the colonies is a singlespecies siderophore producer colony with initial biomass volume fraction N_{0} = 0.3, R_{0} = 0, the other colony is a singlespecies colony that is not able to produce siderophores, with R_{0} = 0.3 and N_{0} = 0. The concentration boundary layer is as in (b), L_{ BL }= 1000 μm.
(d) A simulation in which both colonies are identical, occupied by equal parts of each species, N_{0} = R_{0} = 0.15. The concentration boundary layer is as in (b), L_{ BL }= 1000 μm.
Thus, in all four scenarios the total amount of biomass, and, hence, iron consumers is initially the same, which enables comparison across simulations.
In the competition case (c), where two separate single species colonies are considered, we observe again that initially, as long as iron is not severely limited both colonies grow at the same rate. Eventually the growth of R lags behind the growth of N. At t ≈ 7d the amount of R starts decreasing as a consequence of iron limitation, while the siderophore producer colony keeps growing, albeit at a decreased rate. The competition between both colonies in this case is nonlocal in the sense that both species are spatially separated. Fickian diffusion of dissolved iron is the facilitator of this competition. As discussed in the context of cases (a) and (b), the chelation mechanism is a sink for dissolved iron. Therefore, the iron concentration S in the chelating colony is lower than in the nonchelating colony (while the concentration of chelated iron is higher). This leads to a diffusive flux from the nonchelator to the chelator. The iron gained by the chelator in this way is converted into chelated iron and utilized for growth. Eventually the iron available to the nonchelator drops below the levels that are required to sustain growth.
In the competition case (d), where two mixed culture colonies are considered, the overall picture is qualitatively similar as in case (c), however with big quantitative differences. The maximum amount of biomass that is reached by the nonsiderophore producing species R remains below the one of case (c) and is attained earlier. On the other hand the chelator N grows faster and to larger population levels. The competition for iron between both species is direct and local. Both have access to the same amount of dissolved iron, however, the chelating species scavenges some or most of it and transforms it into chelated iron that cannot be utilized by its competitor, which, therefore, eventually stays behind in its development and is outcompeted.
In summary, the simulations of cases (a)  (d) not only clearly show the absolute growth advantage of the chelator compared to the nonchelator but also indicate that the competition can be due to two different effects, namely direct competition for dissolved iron and the advantage that the chelator gains by scavenging, as well as indirect competition that is facilitated by mass transfer of the growth controlling substrate, dissolved iron S, from regions, in which no chelation takes place, to regions, in which chelation takes place. However, the simulations imply, by comparing (c) with (d), that local competition is stronger than the nonlocal diffusion facilitated competition.
Moreover, this simulation experiment indicates that the overall competition for dissolved iron depends quantitatively on the initial distribution of chelator and nonchelators relative to each other. In mixedspecies colonies competition is direct and fiercer than between separated single species colonies. Initial attachment sites are difficult to control in biofilm experiments and appear to some extent stochastic. In order to investigate the effect that this random inoculation of both species has on the iron chelation process, we conduct the following simulation experiment.
Simulations with uncontrolled inoculation show that sensitivity to initial attachment sites is due to substrate diffusion
We conduct two simulation experiments, one for each of the iron replenishment mechanisms described above. One simulation trial consists of ten biofilm growth simulation with different randomly chosen biomass inoculations, where across the trial the number of inoculations sites is kept constant and also the minimum and maximum values between the volume fractions in the inoculation sites are chosen constant across the simulations in one trial.
Boundary layer concentration prescribed (Robin conditions)
The computations were carried out on a grid with 600 × 200 cells and size L × H = 1.5 mm × 0.5 mm. Initially the substratum is inoculated in 6 randomly chosen sites each for the chelator N and the nonchelator R, with initial volume fraction in these sites randomly chosen between 0.2 and 0.4. Iron replenishment is in these simulation described by Robin boundary conditions with concentration boundary layer thickness L_{ BL }= 1 mm.
The specific biofilm morphology that develops is of course clearly affected by the initial inoculation sites. The population sizes themselves show the same qualitative behavior for all simulations in the trial. However, there are quantitative differences. For example, in addition to the difference in the time at which R_{ avg }starts declining there are also differences observed in the maximum population size, ranging between R_{ avg }≈ 0.05 and R_{ avg }(t) ≈ 0.09, i.e. by a factor of 80%. A trend seems to be that the population at maximum is the higher, the later it is achieved, although there are some exceptions. The simulations were stopped at t = 10d, at which time the chelator population size varied between N_{ avg }= 0.2 and N_{ avg }= 0.3, i.e. by a factor of 50%. This suggests that the initial inoculation sites can affect the outcome of the simulation experiment quantitatively.
Iron flux prescribed (nonhomegenous Neumann conditions)
The boundary conditions (7) are designed to grow a biofilm with an a priori specified average biofilm thickness in the sense defined above. The previous simulation experiments have shown rather high variations in biofilm size as a consequence of variations in substratum inoculation. Furthermore, as discussed above, the substrate flux into the domain increases when the substrate concentration in the domain decreases. The chelation process is an iron sink and thus increases substrate supply to the system under this boundary condition. The nonhomogenous Neumann boundary condition (7) instead of (6) avoids this diffusion effect on the competition for iron, because the substrate flux is prescribed as a constant and independent of the biofilm and the substrate utilization itself. The anticipated biofilm size is prescribed as an input parameters of the boundary condition.
Eventually, after the biofilm has grown to a certain size and after the excess iron is depleted, the growth slows down and is only controlled by iron replenishment. The duration of this initial transient phase does not depend on the inoculation sites. The overall qualitative behavior with respect to population growth is in all cases as in the previous simulations: the species that cannot produce siderophores first grows but eventually declines and dies out, while the siderophore producing species survives. Since the maximum sustainable biofilm size is specified as an input parameter, in all simulations in one trial the population size converges to the same value, namely λ/H, where it reaches a plateau. In the transient phase between the initial growth period and the plateau phase, variations of population counts are observed across trials. The maximum populations sizes of the nonchelator vary between approximately 0.033 and 0.047 [λ = 100 μm; 42% variation], between approximately 0.055 and 0.068 [λ = 150 μ; 24% variation] and between approximately 0.070 and 0.088 [λ = 200 μm; 25% variation]. While inoculation induced variations across trials were observed in all cases, they seem smaller than in the case of the Robin boundary condition above. Similarly, the variations across population size of the chelator in one trial appear much smaller than in the simulation experiments in which the concentration boundary layer was prescribed instead of the diffusive flux. These observations indicate that in addition to a priori and explicitly controlling the size of the biofilm that is expected, using the nonhomogeneous Neumann boundary conditions leads to population counts that are less sensitive to variations in the substratum inoculation.
Conclusions
We presented a mathematical model for siderophore mediated competitive advantage in dualspecies biofilm systems under iron limitations. The model is based on previously developed building blocks, namely a densitydependent diffusionreaction model for biofilm growth and a mathematical model for pyoverdine production and iron chelation. The model is studied in simulation experiments. So far most laboratory studies of iron chelation and siderophore mediated antagonism have been conducted with suspended cultures, while the biofilm setting has only received little attention. A notable exception is the study [2], in which, however, the siderophore mediated competition between Pseudomonas fluorescens and Bacillus cereus was also found to be influenced by the material of the substratum that was used in the experiments. Therefore, this study could not answer the question whether siderophore mediated anatagonism in biofilm communities is possible. While, for this lack of quantitative experimental data, our model predictions cannot be quantitatively validated against existing experiments, it allows us to draw the following conclusions.
1. Under iron limitations, siderophore producing bacteria have a growth advantage in biofilm systems, compared to bacteria that lack this ability. In contrast to the much better understood situation in suspended populations, this growth advantage is manifested in two aspects: (i) Direct competition between between bacteria that locally share the same space: The siderophores bind iron, which the siderophore producers can utilize later, while it becomes unavailable to their competitors. (ii) The chelation process is a local iron sink that induces Fickian diffusion toward siderphore producing colonies from colonies that do not produce siderophores. This increased supply of iron constitutes a further growth advantage for the siderophore producing species. Also in mixed colonies, siderophore producers benefit from Fickian diffusion of iron to the biofilm more than nonsiderophore producers, because they convert the increased supply of dissolved iron to an increased supply of chelated iron. Comparing local, direct siderophore mediated antagonism (i) with the nonlocal mass transfer facilitated effects (ii), our results indicate that iron chelation provides a greater growth advantage in mixedspecies colonies than between single species cultures of of siderophore producers and species that lack this ability. In suspended cultures, (i) is the only competition effect, while (ii) is specific to the biofilm setup. Indeed, in [14] it is argued that diffusion induced substrate concentration gradients is in many respects one of the major differences between life in biofilm communities and living in a suspended culture. Our results show that this applies to siderophore mediated antagonism under iron limitation as well.
2. The location where bacteria initially attach on the substratum are difficult to control or predict in experimental and natural biofilm systems. We mimicked this by stochastic inoculation. Our simulations indicate, however, that siderophore mediated antagonism and competition for iron can be greatly affected by the initial spatial distribution of biomass. Generally one expects that this sensitivity to attachment sites is smaller in densely inoculated systems (where one can expect that all colonies are soon mixed) than in sparsely inoculated systems (where singlespecies colonies can persist). The former is the case in soils and root systems where bacteria are found in abundance, while the latter is the case in food safety applications, where pathogens are hopefully scarce.
3. Boundary conditions for the dependent variables are an important part of every biofilm model. While formulating boundary conditions along physical boundaries is often unproblematic, the crux in biofilm modeling is that, due to computational limitations, all simulation studies are constrained to open subdomain of physical systems. In this case, the boundary conditions connect the computational domain with the physical processes outside. In biofilm modeling the major issue is to describe the replenishment of consumed required resources through the boundaries of the computational domain. The boundary condition that is probably most frequently used to describe replenishment of consumed substrates is to prescribe the concentration value along some part of the boundary. Such Dirichlet conditions imply, as a consequence of diffusion of consumed substrates from the boundary toward the biofilm, that the amount of material that becomes available to the biofilm grows unbounded as the biofilm grows. This is not an appropriate description for modeling studies that focus on substrate limitations. We presented here two alternate boundary conditions that dampen or alleviate this effect. (i) The first one is a Robin boundary condition which introduces, in accordance with the boundary conditions used in traditional onedimensional biofilm models, the abstract concentration boundary layer thickness L_{ BL }as an additional parameter. The traditional Dirichlet condition is a special case of this condition for L_{ BL }= 0. Introducing this concentration boundary layer puts an upper bound on the substrate flux in to the system but still implies that increasing biomass increases the flux into the system. (ii) The second boundary condition tested is a nonhomogeneous Neumann condition. The substrate flux into the system can be simply correlated with the target biofilm size. Thus, this boundary conditions allows an easy control over the expected biomass accumulation in simulation experiments. Since the target biofilm size is a priori specified, the eventually prevailing siderophore producing species is less less sensitive to stochastic uncertainty effects in the longterm than in case of boundary conditions (i). Also the uncertainty effects in the growth curves of the nonsiderophore producer, albeit still notable, are smaller than in case (i). The boundary conditions chosen not only affect simulation results quantitatively, but also qualitatively. Which boundary condition is more appropriate for a particular application depends on the reactor and the environmental conditions modeled. Therefore, sensitivity of the chelation process with respect to initial attachment can depend on the reactor type.
4. The model presented here is a barebone model that focuses solely on the competition effect. It can serve as a qualitative tool and guide experimental studies of the phenomenon. As it is the case with all multidisciplinary research, this requires the collaboration of researchers with complimentary skill sets and research infrastructure. Such a cooperation between theoreticians and experimentalists will allow the experiments to feed back into and improve the theory. The mathematical model framework used here has been used to study several other biofilm systems before. It has been shown that this modeling concept as well as the computational techniques employed to study it are fairly flexible to incorporate additional biofilm processes from population and resource dynamics, such as competition for substrate, amensalism, biofilm response to biocides, convective transport of dissolved substrates, or quorum sensing. Therefore, the model presented here can serve as a starting point for model refinements and future improvements that may be revealed by experimental results to be necessary to account for a qualitatively and quantitatively accurate description of siderophore mediated microbial antagonism.
Declarations
Acknowledgements
This study was conducted as part of the project "Bacteria, Biofilms, and Foods", funded the Advanced Foods and Material network (AFMNet), a Network of Centers of Excellence (NCE). HJE wishes to thank Hedia Fgaier (University of Guelph) and Robin McKellar (now retired from Agriculture and Agrifood Canada) for many discussions on the topic of siderophore producing and iron chelating bacteria.
Authors’ Affiliations
References
 Andrews SC, Robinson AK, RodriguezQuinones F: Bacterial iron homeostasis. FEMS Microbiology Reviews. 2003, 27: 215237. 10.1016/S01686445(03)00055X.View ArticlePubMedGoogle Scholar
 Simoes M, Simoes LC, Pereira MO, Viera MJ: Antagonism between Bacillus cereus and Pseudomonas fluorescens in planktonic systems and in biofilms. Biofouling. 2008, 24 (5): 339349. 10.1080/08927010802239154.View ArticlePubMedGoogle Scholar
 Kloepper JW, Leong J, Teintz M, Schroth MN: Enhanced plant growth by siderophores produced by plant growthpromoting bacteria. Nature. 1980, 286: 885886. 10.1038/286885a0.View ArticleGoogle Scholar
 Loper JE, Buyer JS: Siderophores in microbial interactions on plant surfaces. Molec PlantMicrobe Interactions. 1991, 4 (1): 513.View ArticleGoogle Scholar
 O'Sullivan DJ, O'Gara F: Traits of fluorescent Pseudomonas spp. involved in suppression of plant root pathogens. Microbiol Reviews. 1992, 56 (4): 662676.Google Scholar
 Chen CM, Doyle MP, Luchansky JB: Identification of Pseudomonas fluorescens strains isolated from raw pork and chicken that produce siderophores antagonistic towards foodborne pathogens. J Food Protection. 1995, 58 (12): 13401345.Google Scholar
 Freedman DJ, Kondo JK, Willrett DL: Antagonism of foodborne bacteria by Pseudomonas spp.: a possible role for iron. J. Food Protection. 1989, 52 (7): 484489.Google Scholar
 Gram L: Inhibitory effect against pathogenic and spoilage bacteria of Pseudomonas strains isolated from spoiled and fresh fish. Appl Env Microbiol. 1993, 59 (7): 21972203.Google Scholar
 Gram L, Melchiorsen J: Interaction betwee fish spoilage bacteria Pseudomonsas sp. and Shewanella putrefaciens in fish extracts and on fish tissue. J Appl Bacteriol. 1996, 80: 589595.View ArticlePubMedGoogle Scholar
 Gram L, Ravn L, Rasch M, Bruhn JB, Christensen AB, Giskov M: Food spoilage  interactions between food spoilage bacteria. Int J Food Microbiol. 2002, 78: 7997. 10.1016/S01681605(02)002337.View ArticlePubMedGoogle Scholar
 Banin E, Vasil ML, Greenberg EP: Iron and Pseudomonas aeruginosa biofilm formation. PNAS. 2005, 102 (31): 1107611081. 10.1073/pnas.0504266102.PubMed CentralView ArticlePubMedGoogle Scholar
 Greenberg EP, Banin E: Ironing out the biofilm problem. Control of Biofilm Infections by Signal Manipulation. Edited by: Balaban N. 2008, Springer, BerlinHeidelbergGoogle Scholar
 Visca P, Imperi F, Lamont L: Pyoverdine siderphores: from biogenesis to biosignificance. Trends in Microbiology. 2006, 15 (1): 2230. 10.1016/j.tim.2006.11.004.View ArticlePubMedGoogle Scholar
 Stewart PS: Diffusion in biofilms. J Bacteriol. 2003, 185 (5): 14851491. 10.1128/JB.185.5.14851491.2003.PubMed CentralView ArticlePubMedGoogle Scholar
 Chmielewski RAN, Frank JF: Biofilm formation and control in food processing facilities. Compr Revs Food Sc and Food Saf. 2003, 2: 2232. 10.1111/j.15414337.2003.tb00012.x.View ArticleGoogle Scholar
 Hota S, Hirji Z, Stockton K, Lemieux C, Dedier H, Wolfaardt GM, Gardam M: An outbreak of multidrugresistant Pseudomonas aeruginosa colonization and infection secondary to imperfect intensive care unit room design. Infect Control Hosp Epidemiol. 2009, 30: 2533. 10.1086/592700.View ArticlePubMedGoogle Scholar
 Verran J: Biofouling in food processing: Biofilm or biotransfer potential?. Trans. IChemE. 2002, 80 (C): 292298.Google Scholar
 Costerton JW, Stewart PS, Greenberg EP: Bacterial Biofilms: A Common Cause of Persistent Infections. Science. 1999, 241: 13181322. 10.1126/science.284.5418.1318.View ArticleGoogle Scholar
 Anguige K, King JR, Ward JP: Modelling antibiotic and quorum sensing treatment of a spatiallystructured Pseudomonas aeruginosa population. J Math Biol. 2005, 51: 557594. 10.1007/s0028500503168.View ArticlePubMedGoogle Scholar
 Wanner O, Eberl H, Morgenroth E, Noguera D, Picioreanu C, Rittmann B, van Loosdrecht M, (IWA Task Group on Biofilm Modeling): Mathematical modeling of biofilms. 2006, IWAPublishing, LondonGoogle Scholar
 Klapper I, Dockery J: Mathematical description of microbial biofilms. SIAM Review, to appear, electronic preprint.http://www.math.montana.edu/klapper/home/PSman/review4.pdfhttp://www.math.montana.edu/klapper/home/PSman/review4.pdf
 Fgaier H, Feher B, McKellar RC, Eberl HJ: Predictive modeling of siderphore production by Pseudomonas fluorescens under iron limitation. J Theoretical Biology. 2008, 251 (2): 348362. 10.1016/j.jtbi.2007.11.026.View ArticleGoogle Scholar
 Fgaier H, Eberl HJ: Parameter identification and quantitative comparison of differential equations that describe physiological adaptation of a bacterial population under iron limitation. DCDS Supplements. 2009, 2009: 230239.Google Scholar
 Duvnjak A, Eberl HJ: Timediscretisation of a degenerate reactiondiffusion equation arising in biofilm modeling. El Trans Num Analysis. 2006, 23: 1538.Google Scholar
 Eberl HJ, Parker DF, van Loosdrecht MCM: A new deterministic spatiotemporal continuum model for biofilm development. J Theor Medicine. 2001, 3 (3): 161175. 10.1080/10273660108833072.View ArticleGoogle Scholar
 Eberl HJ, Sudarsan R: Exposure of biofilms to slow flow fields: the convective contribution to growth and disinfections. J Theor Biol. 2008, 253 (4): 788807. 10.1016/j.jtbi.2008.04.013.View ArticlePubMedGoogle Scholar
 Efendiev MA, L Demaret L: On the structure of attractors for a class of degenerate reactiondiffusion systems. Adv Math Sci & Appls. 2008, 18: 105118.Google Scholar
 Efendiev MA, Eberl HJ, Zelik SV: Existence and longtime behavior of solutions of a nonlinear reactiondiffusion system arising in the modeling of biofilms. RIMS Kokyuroko. 2002, 1258: 4971.Google Scholar
 Efendiev MA, Zelik SV, Eberl HJ: Existence and longtime behavior of a biofilm model. Comm Pure Appl Analysis. 2009, 8 (2): 509531. 10.3934/cpaa.2009.8.509.View ArticleGoogle Scholar
 Eberl HJ, Schraft H: A diffusionreaction model of a mixed culture biofilm arising in food safety studies. Mathematical modeling of biological system, Birkhäuser. Edited by: Deutsch A, et al. 2007, II: 109120.Google Scholar
 Eberl HJ, Khassehkhan H, Demaret L: A mixedculture model of a probiotic biofilm control system. Comp. Math. Meth. Med., accepted, electronic preprinthttp://www.informaworld.com/smpp/contentcontent=a915661754db=alljumptype=rsshttp://www.informaworld.com/smpp/contentcontent=a915661754db=alljumptype=rss
 Khassehkhan H, Efendiev MA, Eberl HJ: A degenerate diffusionreaction model of an amensalistic probiotic biofilm control system: Existence and simulation of solutions. Discr Cont Dyn Sys B. 2009, 12 (2): 371388. 10.3934/dcdsb.2009.12.371.View ArticleGoogle Scholar
 Fgaier H, Eberl HJ: A Competition Model Between Pseudomonas fluorescens and Pathogens Via Iron Chelation. J Theor Biol. 2008,Google Scholar
 Khassehkhan H, Hillen T, Eberl HJ: A nonlinear master equation for a degenerate diffusion model of biofilm growth. LNCS. 2009, 5544: 735744.Google Scholar
 Demaret L, Eberl HJ, Efendiev MA, Lasser R: Analysis and simulation of a mesoscale model of diffusive resistance of bacterial biofilms to penetration of antibiotics. Adv Math Sci Appls. 2008, 18: 269304.Google Scholar
 Picioreanu C, van Loosdrecht MCM, Heijnen JJ: A new combined differentialdiscrete cellular automaton approach for biofilm modeling: Application for growth in gel beads. Biotech Bioeng. 1998, 57 (6): 718731. 10.1002/(SICI)10970290(19980320)57:6<718::AIDBIT9>3.0.CO;2O.View ArticleGoogle Scholar
 Pizarro G, Griffeath D, Noguera D: Quantitative cellular automaton model for biofilms. J Env Eng. 2001, 127 (9): 782789. 10.1061/(ASCE)07339372(2001)127:9(782).View ArticleGoogle Scholar
 Pizarro GE, Texiera J, Sepulveda M, Noguera DR: Bitwise implementation of a 2D cellular automata biofilm model. J Comp Civ Eng. 2005, 19 (3): 258268. 10.1061/(ASCE)08873801(2005)19:3(258).View ArticleGoogle Scholar
 Dockery J, Klapper I: Finger formation in biofilm layers. SIAM J Appl Math. 2002, 62: 853869. 10.1137/S0036139900371709.View ArticleGoogle Scholar
 Eberl HJ, Muhammad N, Sudarsan R: Computing Intensive Simulations in Biofilm Modeling. Proc. 22nd Int. High Performance Computing Symposium (HPCS2008), Quebec City, IEEE Conference Series. 2008, 132138.Google Scholar
 Muhammad N, Eberl HJ: OpenMP Parallelization of a Mickens TimeIntegration Scheme For a MixedCulture Biofilm Model and its Performance on Multicore and Multiprocessor Computers. LNCS. 2010,Google Scholar
 Khassehkhan H, Eberl HJ: Modeling and simulation of a bacterial biofilm that is controlled by pH and protonated lactic acids. Comp Math Meth Med. 2008, 9 (1): 4767. 10.1080/17486700701797922.View ArticleGoogle Scholar
 Eberl HJ, Demaret L: A finite difference scheme for a degenerated diffusion equation arising in microbial ecology. El J Diff Equs CS. 2007, 15: 7795.Google Scholar
 Cogan NG: Twofluid model of biofilm disinfection. Bull Math Biol. 2008, 70 (3): 800819. 10.1007/s1153800792803.View ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.