 Research article
 Open Access
 Published:
Spatiotemporal modeling of microbial metabolism
BMC Systems Biology volume 10, Article number: 21 (2016)
Abstract
Background
Microbial systems in which the extracellular environment varies both spatially and temporally are very common in nature and in engineering applications. While the use of genomescale metabolic reconstructions for steadystate flux balance analysis (FBA) and extensions for dynamic FBA are common, the development of spatiotemporal metabolic models has received little attention.
Results
We present a general methodology for spatiotemporal metabolic modeling based on combining genomescale reconstructions with fundamental transport equations that govern the relevant convective and/or diffusional processes in time and spatially varying environments. Our solution procedure involves spatial discretization of the partial differential equation model followed by numerical integration of the resulting system of ordinary differential equations with embedded linear programs using DFBAlab, a MATLAB code that performs reliable and efficient dynamic FBA simulations. We demonstrate our methodology by solving spatiotemporal metabolic models for two systems of considerable practical interest: (1) a bubble column reactor with the syngas fermenting bacterium Clostridium ljungdahlii; and (2) a chronic wound biofilm with the human pathogen Pseudomonas aeruginosa. Despite the complexity of the discretized models which consist of 900 ODEs/600 LPs and 250 ODEs/250 LPs, respectively, we show that the proposed computational framework allows efficient and robust model solution.
Conclusions
Our study establishes a new paradigm for formulating and solving genomescale metabolic models with both time and spatial variations and has wide applicability to natural and engineered microbial systems.
Background
Mathematical models of cellular metabolism are a complementary tool to experimentation for analyzing and engineering metabolic function. Over the past several decades, flux balance analysis (FBA) based on stoichiometric descriptions of cellular metabolism has emerged as the dominant approach for microbial metabolic modeling. FBA involves the formulation of stoichiometric equations describing the metabolic network followed by linear program solution of the underdetermined linear equation system subject to an assumed cellular objective such as growth rate maximization [1]. The advent of genome sequencing and bioinformatic technologies has allowed the reconstruction of largescale metabolic networks in model organisms, which paved the way for the extension of FBA to genomescale metabolic networks [2]. Curated genomescale metabolic reconstructions are now available for a wide variety of microbial species, with new reconstructions announced on a weekly basis. Because genomescale modeling is now an established tool, research has increasingly focused on novel ways to use these reconstructions for metabolic systems analysis and design.
Classical FBA methods assume time invariant and spatially homogeneous extracellular conditions and generate steadystate predictions consistent with wellmixed, continuous cultures [3]. Most microbial systems involve time and/or spatially dependent environments that should be incorporated within the metabolic description. The limitations of steadystate metabolic models have been addressed through dynamic extensions of stoichiometric models and classical FBA [4–7]. Dynamic flux balance models are obtained by combining stoichiometric equations for intracellular metabolism with dynamic mass balances on extracellular substrates and products under the assumption that intracellular metabolite concentrations equilibrate rapidly in response to extracellular perturbations [8]. The intracellular and extracellular descriptions are coupled through the cellular growth rate, secretion fluxes and substrate uptake kinetics, which can be formulated to include complex regulatory effects such as growth inhibition by metabolic byproducts. Dynamic flux balance modeling is now an established extension of FBA.
In contrast to the dynamic case, the development of metabolic models that account for spatially varying environments has received little attention. Such problems are very common in natural and engineered microbial systems. For example, naturally occurring microbial biofilms typically exhibit strong spatial gradients due to differential nutrient availability at the biofilm boundaries [9]. Spatial gradients are also present in synthesis gas bubble column reactors because dissolved CO and H_{2} concentrations decrease as the gas flows up the column due to cellular consumption [10]. The incorporation of genomescale metabolic reconstructions within spatiotemporal models that account for both spatial and temporal variations in the environment is desirable to connect genes to metabolic phenotype and system function. For example, genomescale metabolic reconstructions allow the effects of gene deletions and insertions in mutant strains to be directly investigated. Genomescale spatiotemporal models have been solved using table lookups of precomputed FBA solutions [11–13], lattice based descriptions of nutrient diffusion [14, 15] and agentbased simulations [16]. These methods utilize a fixed time step over which the FBA linear program (LP) solution is assumed to remain unchanged and the ordinary differential equations (ODEs) representing the extracellular environment are integrated. By contrast, our approach allows the LP to be directly embedded within the ODEs and to be solved with variable time steps chosen by a stiff integrator. Therefore our computational framework represents an important step towards solving spatiotemporal models that combine a genomescale description of intracellular metabolism and fundamental transport equations for the extracellular environment.
Methods
Model structure
The class of spatiotemporal metabolic models considered below is sufficiently general to encompass a wide variety of potential applications including microbial communities with interacting species and multiphase systems in which the liquid and gas phases move relative to each other. The framework is based on the standard dynamic flux balance modeling assumption that the intracellular metabolism is much faster than the extracellular dynamics, which we do not believe is any more restrictive when the environment exhibits spatial variations. Furthermore, we assume that spatial variations occur only in a single direction z for simplicity. Additional modeling assumptions include that constant gas and liquid phase volume fractions and velocities, constant gas–liquid mass transfer coefficients, constant cell and metabolite diffusion coefficients, and cell incompressibility. The last assumption allows a simple convection term to be used in the species mass balance equations. While cell compressibility could be included in the model if necessary, we expect that this effect would negligible under the low velocity liquid flows encountered in the examples consider here as well as in most practical applications.
Under these assumptions, a general set of model equations can be written as,
The first equation represents a mass balances on the ith microbial species where X _{ i } is the biomass concentration, μ _{ i } is the growth rate obtained from the genomescale metabolic model, μ _{ di } is the death rate, u _{ L } is the liquid phase velocity, \( {\varepsilon}_L \) is the liquid volume fraction and D _{ iX } is the cellular diffusion coefficient that accounts for cell motility. The second equation represents a mass balance on the jth liquid phase metabolite where v _{ ij } is the net flux of metabolite j into the liquid phase from species i, D _{ jL } is the liquidphase metabolite diffusion coefficient, k _{ j } is the gas–liquid mass transfer coefficient, and M _{ j } ^{*} is saturation concentration in the liquid phase calculated from the associated gasphase concentration using Henry’s law. The net flux v _{ ij } is calculated as the difference between the synthesis rate obtained from the genomescale metabolic model and the uptake rate calculated from MichaelisMenten type kinetic expressions [17, 18]. The third equation represents a mass balance on the jth gas phase component where u _{ G } is the gas phase velocity, \( {\varepsilon}_G \) is the gas volume fraction and D _{ jG } is the gasphase diffusion coefficient.
Boundary conditions for these equations are problem specific and can account for the supply/removal of liquid and/or gas phase components at the domain boundaries. Although not discussed here, the general model formulation can be extended to include a moving boundary as would be required for biofilm expansion. In the Appendix, two examples of formulated spatiotemporal metabolic models are presented: (1) a bubble column reactor for bacterial conversion of synthesis gas to ethanol; and (2) a bacterial biofilm associated with chronic wound infections. For the most part, these two models adhere to the general set of equations presented above. However, the species mass balance equation for the biofilm model (Equation 12 in Additional file 1: Appendix) has a slightly different formulation than the general equation to compensate for the lack of biofilm expansion in the model (see Additional file 1: Appendix).
Model solution
Simulation of spatiotemporal metabolic models involves numerically solving a set of nonlinear partial differential equations (PDEs) with embedded linear programs. The efficient and stable solution of such models is a challenging problem at the forefront of microbial metabolic modeling. Our solution approach is based on spatial discretization such that the PDEs are converted into a large set of ordinary differential equations (ODEs) in time with embedded LPs. The spatial domain is discretized with N node points using an appropriate discretization method such as finite difference, finite volume or orthogonal collocation. If the original PDE model contains N _{ X } microbial species, N _{ M } liquidphase metabolites and N _{ P } gasphase components, then the discretization procedure will yield a dynamic FBA model with N _{ X } + N _{ M } + N _{ P } ODEs and N _{ X } LPs at each node point.
Our approach for solving such large discretized models involves the use of DFBAlab [19], a MATLAB code that performs reliable and efficient dynamic FBA simulations. Widespread implementation of dynamic FBA has been hindered by numerical complications resulting from LPs becoming infeasible and having nonunique solution vectors. Infeasible LPs cause simulation failure as the righthand side of the ODEs becomes undefined, and nonunique exchange fluxes cause this same righthand side to become nonunique, producing an ODE system that integrators are unable to solve. These complications are addressed in our previous publication [20].
DFBAlab is a modified MATLAB implementation of our previously developed simulator [20]. DFBAlab reformulates the LP locally as an algebraic system, and it integrates a differentialalgebraic equation system instead of ODEs with LPs embedded to increase speed. Hierarchical fixedpriority preemptive (lexicographic) optimization is used to determine uniquely all fluxes which appear in the righthand side of the ODEs (i.e. exchange fluxes). All other fluxes not optimized lexicographically (i.e. internal fluxes) may still be nonunique, but their values do not affect the righthand side of the ODEs. With lexicographic optimization, the righthand side of the ODEs is guaranteed to be unique, allowing efficient and reliable integration. Finally, DFBAlab uses the Phase I LP of the simplex algorithm combined with lexicographic optimization to avoid infeasibilities.
More specifically, DFBAlab reformulates the FBA LP as a Phase I lexicographic LP to obtain all information required by the righthand side of the ODEs as a unique vectorvalued solution with the following order of objectives:

1.
Minimize infeasibilities: If the first objective is equal to zero, the LP is feasible and all other objectives are consistent with the solution of the original FBA LP; otherwise, the objective is positive. If the original FBA LP is infeasible, the reformulated Phase I lexicographic LP still returns values for growth rate and exchange fluxes allowing the integration process to continue. This objective can be integrated to obtain a penalty function. This penalty function can provide useful insights on why and under what conditions the FBA model becomes infeasible.

2.
Maximize growth rate: this is the traditional FBA objective.

3.
Maximize/minimize all of the exchange fluxes appearing in the righthand side of the ODEs. Each one of these objectives involves a linear combination of fluxes that can be minimized or maximized as appropriate. If there are n fluxes appearing in the right hand side of the ODEs, the vectorvalued objective will require at most n + 2 elements to obtain a unique righthand side.
DFBAlab is designed to solve ODE systems; however, it provides a flexible framework that enables the solution of PDE models if the equations can be transformed into ODEs. Consider the following equation that describes the biomass concentration of the ith species in a bubble column reactor:
This PDE can be easily converted into an ODE by discretizing the spatial domain. If a simple backward difference formula is used to approximate the convection term, then the following set of ODEs is obtained for each point j in the spatial domain and each species i:
where L is the length of the spatial domain, ΔL = L/n and n is the number of discretization points. In addition, X _{ i,0 } = X _{feed} = 0 and the outlet biomass concentration of the bubble column reactor is X _{ i,n }. A more detailed explanation for the single species case can be found in Fig. 1. A similar procedure is followed to obtain ODEs corresponding to the discretized PDEs of the liquid and gas phase components in (1). The flexibility of DFBAlab allows for the easy implementation and fast simulation of such discretized PDE systems. To ensure physically meaningful predictions for the two case studies presented below, the species growth rate μ_{i} and the secretion exchange fluxes v _{ ij } were set equal to zero whenever DFBAlab detected that the LP for the ith species was infeasible. This situation occurred when local nutrient uptake rates were insufficient to meet the nongrowth ATP maintenance requirements. While this approach had the potential to make the righthand side of the ODEs discontinuous, we found that DFBAlab had little problem integrating through such points because the growth rates and byproduct secretion rates tended to be very small immediately prior to an infeasibility occurring.
From a biological perspective, the additional objectives involving the exchange fluxes represent lower level cellular strategies than the main objective of growth rate maximization. The choice of these objectives is problem dependent and requires assumptions about the cellular metabolism. We typically assume that the cell regulation machinery is configured to maximize substrate uptake fluxes and minimize byproduct secretion fluxes, which is consistent with the main objective by maximizing the input of carbon containing and electron accepting metabolites and minimizing the output of carbon wasting byproducts. While DFBAlab requires specification of these lower level objectives, they impact the lexicographic optimization only when alternative optima occur. Our experiences with the two examples discussed in the following sections and other problems solved with DFBAlab is that the ordering of these objectives has a negligible impact on spatiotemporal model solutions because alternative optima typically occur only for short periods of simulation time. In other words, DFBAlab allows the integrator to reliably transition across short periods where alternative optima exist.
Simulation codes
All simulations were performed with MATLAB 8.5 (R2015a) using DFAlab for dynamic flux balance model solution and Gurobi 6.0 for linear program solution. DFBAlab is freely available for both education and nonprofit research purposes from https://yoric.mit.edu/dfbalab. Any entity desiring permission to incorporate this software or a work based on the software into commercial products or otherwise use it for commercial purposes should contact Dr. Paul Barton (pib@mit.edu). Simulation codes for the synthesis gas bubble column reactor and bacterial biofilm models can be obtained from www.ecs.umass.edu/che/henson_group/downloads.html.
Results and discussion
Spatiotemporal simulation of a synthesis Gas bubble column reactor model
An emerging route for the largescale production of renewable fuels and chemicals is direct fermentation of waste gas streams and synthesis gas (syngas; mainly comprised of H_{2}/CO/CO_{2}) by specialized CO fermenting microbes. Because syngas can be produced relatively cheaply from a wide variety of biomass feedstocks [21, 22], the bottleneck in this route is the syngas fermentation step. Commercial development efforts are currently focused on bubble column reactors due to their superior gas–liquid mass transfer characteristics and enhanced operational flexibility [10]. Because CO and H_{2} concentrations decrease as the gas flows up the column due to cellular consumption, the column can have strong spatial gradients that affect cellular growth and product synthesis. The development of modelbased techniques for simulating and optimizing these complex multiphase reactors is important to advance syngas fermentation technology.
1. Bubble column model solution
The bubble column model was formulated by combining a genomescale metabolic reconstruction of the syngas fermenting bacterium Clostridium ljungdahlii [23] with uptake kinetics for dissolved gases and reactionconvection–dispersion type equations for gaseous and dissolved substrates and synthesized metabolic byproducts. Our preliminary FBA calculations with the typical maximum growth objective showed that the only metabolic byproducts for growth on CO/H_{2} mixtures were ethanol, acetate and CO_{2}. While other byproducts could be secreted under bubble column operating conditions, we did not attempt to determine or model other byproducts due to our focus on ethanol production. Therefore, the spatiotemporal metabolic model was comprised of 9 PDEs for the liquidphase concentrations of C. ljungdahlii biomass, ethanol, acetate, CO, H_{2} and CO_{2} and the gasphase concentrations of CO, H_{2} and CO_{2} (see Additional file 1: Appendix). Model parameters were obtained from the literature to the extent possible with the remaining parameters specified within reasonable ranges (Table 1). The interested reader is directed to our other paper [24] for additional details about the bubble column model formulation and model sensitivity to various column operating and substrate uptake parameters.
The convection terms were discretized using an upwind finite difference approximation with thirdorder accuracy due to its well established numerical accuracy and stability properties for convection dominated problems [25]. We found that the addition of axial dispersion terms to the liquid phase mass balances greatly improved numerical stability of the model (see Additional file 1: Appendix), as has been well documented in other applications [25]. These dispersion terms were discretized using a central difference approximation with secondorder accuracy. Because the upwind formula was not implementable at the reactor boundaries, a firstorder backward difference approximation was used at these locations. The discretization procedure yielded a set of 9 ODEs at each node point.
The lexicographic optimization objectives required by DFBAlab were specified to reflect the known or expected physiology of C. ljungdahlii (Table 2). We found that the ordering of these objectives had no noticeable effect on predicted metabolic phenotypes and bubble column behavior. Each node point was represented by 9 ODEs for the biomass and biochemical species concentrations, 3 algebraic equations for the local dissolved gas uptake rates and 6 LPs for lexicographic optimization. We typically employed 100 node points to obtain a nearly converged solution using DFBAlab combined with the LP solver Gurobi and the stiff ODE solver ode15s.
2. Prediction of bubble column performance
Our first goal was to investigate the efficiency of DFBAlab for simulating startup of the bubble column reactor with N = 100 node points, which yielded a total of 900 ODEs (9 per point) and 600 LPs (6 per point). Despite the substantial computational complexity of this discretized model, we found that a typical 1000 h dynamic simulation for determining a steadystate solution required only about 8 min running DFBAlab and MATLAB version 7.11 on a Dell XPS laptop with Intel Core i72760QM processor and 8 GB RAM. Time and spatially resolved predictions obtained for reactor startup with a simulation time of 250 h are shown in Fig. 2. Steadystate conditions were achieved approximately 225 h after startup once the rate of biomass production equaled the rate of biomass removal from the top of the column. The gas and liquid phase CO and H_{2} concentrations decreased along the length of the reactor due to gas consumption, while the biomass, acetate and ethanol concentrations increased along the reactor due to liquid flow. The synthesis of CO_{2} was negligible under these nominal operating conditions. Because the feed gas was relatively rich in CO, the H_{2} conversion was 62 % while the CO conversion was only 29 %. As a result of H_{2} being depleted in the first half of the reactor, considerable acetate was produced in the second half of the reactor and the liquid product stream exiting the top of the column contained more acetate than ethanol (ethanol fraction ~40 %). While we are not aware of any published experimental studies that describe the startup dynamics of syngas bubble columns, our model could be a powerful tool for predicting and optimization reactor startup.
To demonstrate that N = 100 node points were sufficient to obtain nearly converged solutions, we performed dynamic simulations for reactor startup with different N values and compared the resulting steadystate solutions obtained at t = 1000 h (Fig. 3). While completely converged solutions appeared to be obtained for 300 node points, this simulation required almost 50 min to complete. For the purposes of this study, we decided that 100 node points provided a suitable compromise between solution accuracy (less than 0.2 % error compared to N = 300) and computational time (~8 min per simulation). All remaining simulations were performed with N = 100.
To demonstrate the power of our computational framework and to gain insights into bubble column reactor dynamics, we performed additional startup simulations with different parameter values. First we changed the feed composition from the nominal 60/40 CO/H_{2} mixture to a 50/50 CO/H_{2} mixture. The column exhibited similar dynamics for this H_{2} rich feed, as the biomass concentration still required approximately 200 h to reach steady state (Fig. 4). However, the increased H_{2} feed concentration produced a more favorable dissolved H_{2} profile along the column, resulting in an enhanced ethanol titer of 102 g/L and a substantially improved ethanolacetate ratio of approximately 3 at the reactor outlet once steady state was reached. The amount of biomass produced was not noticeably changed. Due to the increased H_{2} content of the feed, the H_{2} conversion decreased to 60 % and the CO conversion increased to 35 %. Our model predictions were in qualitative agreement with published experimental studies [26–28] showing that hydrogen rich feeds increase both the ethanol titer and the ethanol/acetate split.
Next we performed a dynamic simulation with the CO gas–liquid mass transfer coefficient changed from the nominal value k _{ m,C } = 80 h^{−1} to a substantially larger value k _{ m,C } = 300 h^{−1}, which could result from syngas microsparging and column internal packing [29]. Consistent with our nominal values (Table 2), we set the H_{2} mass transfer coefficient to be 250 % larger than the CO value and the CO_{2} mass transfer coefficient to equal the CO value. The large increases in gas–liquid mass transfer rates produced faster column dynamics with the biomass concentration requiring only about 150 h to reach steady state (Fig. 5). Once the column reached steady state, the increased mass transfer rates also offered the benefit of increased ethanol titer (120 g/L), a higher ethanolacetate ratio (3.5) and improved CO (33 %) and H_{2} (86 %) conversions compared to the nominal case. Our predictions were in qualitative agreement with published experimental studies [27, 29, 30] showing that enhanced gas–liquid mass transfer increases gas consumption, the ethanol titer and the ethanol/acetate split.
Spatiotemporal simulation of a Bacterial Biofilm
Chronic, nonhealing wounds are a growing medical challenge associated with diabetes and obesity [31]. These wounds are typically colonized by bacterial species such as Pseudomonas growing as biofilms on a complex mixture of wound exudate [32, 33]. Bacteria in biofilms can tolerate antimicrobial agent concentrations 10,000 times higher than the same microbes grown planktonically, making treatment of chronic wound biofilms a major challenge [34]. Carbon sources such as glucose are available only from the exudate through the tissuebiofilm interface at the bottom of the biofilm and oxygen is primarily available through the biofilmair interface at the top of the biofilm. Due to limited diffusion, bacterial biofilms often exhibit strong spatial gradients that affect metabolism, physiology and virulence [35, 36]. The development of predictive metabolic models for simulating these complex spatially structured systems is important to advance understanding and treatment of chronic wound infections.
1. Biofilm model solution
The bacterial biofilm model was formulated by combining a genomescale metabolic reconstruction of the opportunistic human pathogen P. aeruginosa [37] with substrate uptake kinetics and reaction–diffusion equations for substrates and metabolic byproducts. As compared to alternative biofilm modeling approaches based on unstructured intracellular descriptions [38], this model formulation allowed the effects of substrate and byproduct diffusion within the biofilm to be captured with genomescale resolution. Our preliminary FBA calculations showed that the primary byproducts were acetate and Dalanine. To obtain better agreement with experimental data [39] showing anaerobic succinate secretion by P. aeruginosa, we constrained the Dalanine secretion flux to zero such that the only byproducts were acetate and succinate. While other byproducts could be secreted in different biofilm microenvironments, we did not attempt to determine or model other byproducts due to our focus on cellular growth. The spatiotemporal metabolic model was comprised of 5 PDEs for the liquidphase concentrations of P. aeruginosa biomass, glucose, oxygen, acetate and succinate (see Additional file 1: Appendix). Model parameters were obtained from the literature to the extent possible with the remaining parameters specified within reasonable ranges (Table 3). To avoid the complications associated with solving a moving boundary problem, the biofilm was assumed to have a fixed thickness. Therefore, the formulated model was appropriate for predicting the metabolism of P. aeruginosa biofilms of a specified thickness rather than predicting the thickness itself. Model simulations show the spatiotemporal dynamics of cellular metabolism within a fixed spatial domain consistent with growth between two stationary surfaces. Steadystate solutions show the spatial distribution of cell and metabolite concentrations within a biofilm of the prescribed thickness. As expected, we found that growth dynamics were strongly affected and steadystate spatial profiles were less affected by the initial cell concentration.
The diffusion terms were discretized using a central difference approximation with secondorder accuracy, which produced a set of 5 ODEs in time at each node point. The lexicographic optimization objectives were specified to reflect the known or expected physiology of P. aeruginosa (Table 4). We found that the ordering of these objectives had no noticeable effect on predicted biofilm dynamics. Each node point was represented by 5 ODEs for the biomass, glucose, oxygen, acetate and succinate concentrations, 4 algebraic equations for calculating diffusion coefficients as a function of the local biomass concentration [40], and 5 LPs for lexicographic optimization. We used 50 node points for DFBAlab solution with the LP solver Gurobi and the stiff ODE solver ode15s.
2. Prediction of biofilm physiology
We performed a dynamic simulation for a biofilm thickness of 50 μm with N = 50 node points, which produced a discretized model with 250 ODEs (5 per point) and 250 LPs (5 per point). A 1000 h dynamic simulation for determining a steadystate solution required about 15 min on our Dell XPS laptop. This computation time was substantially greater than the 8 min required to simulate the bubble column reactor model over the same time period despite the larger size of the discretized column model (900 ODEs, 600 LPs). While we hypothesize that the increased computation times obtained with the biofilm model were attributable to the diffusion dominated behavior, these results demonstrate the need to better understand the computational complexity of these largescale ODE/LP systems.
Figure 6 shows dynamic simulation results for the 50 μm thick biofilm, where time profiles are presented at the bottom (tissue interface), middle and top (air interface) of the biofilm. The bottom layer was characterized by a high glucose concentration, a very low oxygen concentration and a relatively small biomass concentration with slow dynamics. By contrast, the top layer had a very low glucose concentration, a high oxygen concentration and a relatively large biomass concentration with fast dynamics. Experimental studies [9, 41] also have shown the presence of strong spatial gradients in nutrient (e.g. oxygen) levels within biofilms. Despite having a much lower oxygen concentration, the middle layer exhibited similar dynamic and steadystate behavior as the top layer. Spatially uniform acetate and succinate concentrations were obtained throughout the biofilm due to limited removal of the two byproducts at the tissuebiofilm boundary.
To explore the impact of biofilm thickness on physiology and to further evaluate our modeling framework, we performed a dynamic simulation for a 100 μm thick biofilm (Fig. 7). This thicker biofilm had slower dynamics, with approximately 200 h required to reach a steadystate solution. Major differences between the 50 and 100 μm thicknesses were observed at the top of the biofilm. In particular, the 100 μm thick biofilm exhibited much slower growth dynamics and less biomass accumulation due to the limited glucose diffusion, a prediction in qualitative agreement with experimental data [36, 42] indicating nutrient limited growth in mature biofilms. The thicker biofilm also produced higher levels of the metabolic byproducts, especially succinate, which could potentially serve as a carbon source for aerobic growth in glucose depleted regimes at the top of the biofilm [43].
While the previous simulations were performed assuming the only source of O_{2} was from air at the top of the biofilm, blood plasma has low O_{2} levels [44] that could support limited aerobic growth near the tissuebiofilm interface. To investigate this effect, we modified the boundary condition at z = 0 for the O_{2} mass balance (Equation 13 in Additional file 1: Appendix) from a no flux boundary condition to a mass transfer limited boundary condition with a plasma O_{2} concentration of 0.05 mmol/L. The inclusion of an O_{2} source at this interface resulted in a higher O_{2} level, much faster growth dynamics and more biomass accumulation at the bottom of the biofilm (Fig. 8). The establishment of partially aerobic conditions near the tissuebiofilm interface also reduced the overall level of succinate in the biofilm while the acetate level was unaffected.
Finally we performed a dynamic simulation to investigate the effects of putative succinate reassimilation on biofilm physiology. The thickness was specified as 100 μm and O_{2} was supplied at the tissuebiofilm interface as before. Succinate consumption was included in the model by allowing succinate uptake with the same kinetic parameters as used for glucose (see Equation 11 in Additional file 1: Appendix and Table 3). Figure 9 shows a comparison of steadystate spatial profiles obtained after 1000 h of dynamic simulation for three cases that differ with respect to whether O_{2} was supplied at the tissuebiofilm interface and whether succinate consumption was allowed. If only O_{2} supply at the tissuebiofilm interface was introduced (“O2 Tissue”), the main differences from the nominal case were that more biomass was produced near the interface and lower succinate levels were generated throughout the biofilm. When succinate consumption also was introduced (“Succinate Consume”), then biomass was preferentially accumulated at the top of the biofilm due to succinate oxidation, resulting in a less dense region located in the middle. This prediction was consistent with the known physiology of nutrient limited biofilms [45]. Succinate consumption also resulted in increased acetate levels compared to the other two cases.
Conclusions
Many natural and engineered microbial systems exist in nonhomogeneous environments that require metabolic models that account for both temporal and spatial variations. Our spatiotemporal metabolic modeling framework involves combining genomescale metabolic reconstructions with fundamental transport equations that govern the relevant convective and/or diffusional processes in the extracellular environment. The PDE model is spatially discretized and the resulting system of ODEs with embedded LPs is integrated using DFBAlab [19], a MATLAB code that performs reliable and efficient dynamic FBA simulations. We demonstrated the capabilities of the method by solving large discretized models for a convection dominated syngas bubble column reactor (900 ODEs, 600 LPs) and a diffusion driven bacterial biofilm model (250 ODEs, 250 LPs). The proposed methodology represents an important step towards rigorously solving spatiotemporal models that combine a genomescale description of intracellular metabolism and fundamental transport equations for the extracellular environment. A possible limitation of our modeling framework is computational cost, which depends on the number of microbial species, the number of metabolite uptake and secretion fluxes for each species, and the number of node points used for spatial discretization. Consequently, future research will focus on improving computational efficiency including the reduction of genomescale reconstructions to maintain the same uptakesecretion rate behavior [46] and strategic combination of extracellular byproducts into lumped variables that reduce model size. While our bubble column and biofilm models produce predictions in qualitative agreement with available data, we are currently conducting detailed experimental studies to generate spatially and time resolved data for model validation.
Availability of supporting data
The datasets supporting the results of this article are included within the article and its additional file.
Abbreviations
 DFBAlab:

Dynamic flux balance analysis laboratory
 FBA:

Flux balance analysis
 LP:

Linear program
 ODE:

Ordinary differential equation
 PDE:

Partial differential equation
References
 1.
Papin JA, Price ND, Wiback SJ, Fell DA, Palsson BO. Metabolic pathways in the postgenome era. Trends Biochem Sci. 2003;28(5):250–8.
 2.
Price ND, Papin JA, Schilling CH, Palsson BO. Genomescale microbial in silico models: the constraintsbased approach. Trends Biotechnol. 2003;21(4):162–9.
 3.
Palsson B. Systems biology: properties of reconstructed networks. Cambridge: Cambridge University Press; 2006.
 4.
Hanly TJ, Henson MA. Dynamic flux balance modeling of microbial cocultures for efficient batch fermentation of glucose and xylose mixtures. Biotechnol Bioeng. 2011;108(2):376–85.
 5.
Hjersted JL, Henson MA, Mahadevan R. Genomescale analysis of Saccharomyces cerevisiae metabolism and ethanol production in fedbatch culture. Biotechnol Bioeng. 2007;97(5):1190–204.
 6.
Varma A, Palsson BO. Stoichiometric flux balance models quantitatively predict growth and metabolic byproduct secretion in wildtype Escherichiacoli W3110. Appl Environ Microbiol. 1994;60(10):3724–31.
 7.
Mahadevan R, Edwards JS, Doyle FJ. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys J. 2002;83(3):1331–40.
 8.
Hjersted JL, Henson MA. Steadystate and dynamic flux balance analysis of ethanol production by Saccharomyces cerevisiae. IET Syst Biol. 2009;3(3):167–79.
 9.
Burmolle M, Ren DW, Bjarnsholt T, Sorensen SJ. Interactions in multispecies biofilms: do they actually matter? Trends Microbiol. 2014;22(2):84–91. doi:DOI 10.1016/j.tim.2013.12.004.
 10.
Daniell J, Kopke M, Simpson SD. Commercial Biomass Syngas Fermentation. Energies. 2012;5(12):5372–417.
 11.
Cole JA, Kohler L, Hedhli J, LutheySchulten Z. Spatiallyresolved metabolic cooperativity within dense bacterial colonies. BMC Syst Biol. 2015;9.
 12.
Jayasinghe N, Franks A, Nevin KP, Mahadevan R. Metabolic modeling of spatial heterogeneity of biofilms in microbial fuel cells reveals substrate limitations in electrical current generation. Biotechnol J. 2014;9(10):1350–61.
 13.
Fang Y, Scheibe TD, Mahadevan R, Garg S, Long PE, Lovley DR. Direct coupling of a genomescale microbial in silico model and a groundwater reactive transport model. J Contam Hydrol. 2011;122(1–4):96–103.
 14.
Chubiz LM, Granger BR, Segre D, Harcombe WR. Species interactions differ in their genetic robustness. Front Microbiol. 2015;6.
 15.
Harcombe WR, Riehl WJ, Dukovski I, Granger BR, Betts A, Lang AH, et al. Metabolic Resource Allocation in Individual Microbes Determines Ecosystem Interactions and Spatial Dynamics. Cell Rep. 2014;7(4):1104–15.
 16.
Biggs MB, Papin JA. Novel Multiscale Modeling Tool Applied to Pseudomonas aeruginosa Biofilm Formation. PLoS One. 2013;8(10).
 17.
Hanly TJ, Henson MA. Dynamic metabolic modeling of a microaerobic yeast coculture: predicting and optimizing ethanol production from glucose/xylose mixtures. Biotechnol Biofuels. 2013;6:44.
 18.
Hanly TJ, Urello M, Henson MA. Dynamic flux balance modeling of S. cerevisiae and E. coli cocultures for efficient consumption of glucose/xylose mixtures. Appl Microbiol Biotechnol. 2012;93(6):2529–41.
 19.
Gomez JA, Höffner K, Barton PI. DFBAlab: a fast and reliable MATLAB code for dynamic flux balance analysis. BMC Bioinformatics. 2014;15:409.
 20.
Höffner K, Harwood SM, Barton PI. A reliable simulator for dynamic flux balance analysis. Biotechnol Bioeng. 2013;110(3):792–802.
 21.
Kirkels AF, Verbong GPJ. Biomass gasification: Still promising? A 30year global overview. Renew Sust Energ Rev. 2011;15(1):471–81.
 22.
McKendry P. Energy production from biomass (part 3): gasification technologies. Bioresour Technol. 2002;83(1):55–63.
 23.
Nagarajan H, Sahin M, Nogales J, Latif H, Lovley DR, Ebrahim A, et al. Characterizing acetogenic metabolism using a genomescale metabolic reconstruction of Clostridium ljungdahlii. Microb Cell Fact. 2013;12.
 24.
Chen J, Gomez JA, Hoffner K, Barton PI, Henson MA. Metabolic modeling of synthesis gas fermentation in bubble column reactors. Biotechnol Biofuels. 2015;8:89.
 25.
Finlayson BA. Numerical Methods for Problems with Moving Fronts. Incorporated: Ravenna Park Publishing; 1992.
 26.
Mohammadi M, Mohamed AR, Najafpour GD, Younesi H, Uzir MH. Kinetic Studies on Fermentative Production of Biofuel from Synthesis Gas Using Clostridium ljungdahlii. Sci World J. 2014:8. doi:10.1155/2014/910590
 27.
Liew FM, Köpke M, Simpson SD. Gas Fermentation for Commercial Biofuels Production. INTECH Open Access Publisher. 2013. doi:10.5772/52164
 28.
Younesi H, Najafpour G, Mohamed AR. Ethanol and acetate production from synthesis gas via fermentation processes using anaerobic bacterium. Clostridium ljungdahlii Biochem Eng J. 2005;27(2):110–9.
 29.
Munasinghe PC, Khanal SK. Biomassderived syngas fermentation into biofuels: Opportunities and challenges. Bioresour Technol. 2010;101(13):5013–22.
 30.
Bredwell MD, Srivastava P, Worden RM. Reactor design issues for synthesisgas fermentations. Biotechnol Prog. 1999;15(5):834–44.
 31.
KirketerpMoller K, Zulkowski K, James G. Chronic Wound Colonization, Infection, and Biofilms. Biofilm Infections. 2011. p. 11–24.
 32.
Rani SA, Hoon R, Najafi R, Wang L, Debabov D. What Is the Antimicrobial Activity of Wound and Skin Cleansers at Nontoxic Concentrations? J Wound Ostomy Cont. 2013;40:S84S.
 33.
James GA, Swogger E, Wolcott R, Pulcini E, Secor P, Sestrich J, et al. Biofilms in chronic wounds. Wound Repair Regen. 2008;16(1):37–44.
 34.
Folsom JP, Richards L, Pitts B, Roe F, Ehrlich GD, Parker A, et al. Physiology of Pseudomonas aeruginosa in biofilms as revealed by transcriptome analysis. BMC Microbiol. 2010;10:294.
 35.
Stewart PS. Diffusion in biofilms. J Bacteriol. 2003;185(5):1485–91.
 36.
Stewart PS. A review of experimental measurements of effective diffusive permeabilities and effective diffusion coefficients in biofilms. Biotechnol Bioeng. 1998;59(3):261–72.
 37.
Oberhardt MA, Puchalka J, Fryer KE, Martins dos Santos VA, Papin JA. Genomescale metabolic network analysis of the opportunistic pathogen Pseudomonas aeruginosa PAO1. J Bacteriol. 2008;190(8):2790–803.
 38.
Horn H, Lackner S. Modeling of Biofilm Systems: A Review. Productive Biofilms. 2014;146:53–76.
 39.
Eschbach M, Schreiber K, Trunk K, Buer J, Jahn D, Schobert M. Longterm anaerobic survival of the opportunistic pathogen Pseudomonas aeruginosa via pyruvate fermentation. J Bacteriol. 2004;186(14):4596–604.
 40.
Beyenal H, Tanyolac A, Lewandowski Z. Measurement of local effective diffusivity in heterogeneous biofilms. Water Sci Technol. 1998;38(8–9):171–8.
 41.
HallStoodley L, Costerton JW, Stoodley P. Bacterial biofilms: From the natural environment to infectious diseases. Nat Rev Microbiol. 2004;2(2):95–108.
 42.
Okabe S, Yasuda T, Watanabe Y. Uptake and release of inert fluorescence particles by mixed population biofilms. Biotechnol Bioeng. 1997;53(5):459–69.
 43.
Tralau T, Vuilleumier S, Thibault C, Campbell BJ, Hart CA, Kertesz MA. Transcriptomic analysis of the sulfate starvation response of Pseudomonas aeruginosa. J Bacteriol. 2007;189(19):6743–50.
 44.
Pittman RN. Oxygen gradients in the microcirculation. Acta Physiol. 2011;202(3):311–22.
 45.
Woods J, Boegli L, Kirker KR, Agostinho AM, Durch AM, Pulcini ED, et al. Development and application of a polymicrobial, in vitro, wound biofilm model. J Appl Microbiol. 2012;112(5):998–1006.
 46.
Erdrich P, Steuer R, Klamt S. An algorithm for the reduction of genomescale metabolic network models to meaningful core models. BMC Syst Biol. 2015;9.
Acknowledgements
JC and MAH wish to acknowledge financial support from the National Science Foundation (CBET 1511346) and National Institutes of Health (Award U01EB019416). The authors wish to thank Harish Nagarajan (UCSD) for providing the C. ljungdahlii metabolic reconstruction and Prof. Jason Papin (U. Virginia) for providing the P. aeruginosa metabolic reconstruction.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
JC, PP and MAH conceived the study and developed the model. JC, JAG, KH, PIB and MAH developed the model solution method. JC and MAH performed the simulations and analyzed the results. JC, JAG, KH, PP, PIB and MAH prepared the manuscript. All authors read and approved the final manuscript.
Additional file
Additional file 1: Appendix
The Appendix contains detailed equations for the two spatiotemporal metabolic models studied in this paper: (1) a bubble column reactor for bacterial conversion of synthesis gas to ethanol; and (2) a bacterial biofilm associated with chronic wound infections. (DOCX 45 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Metabolic modeling
 Bioreactors
 Biofuels
 Gas fermentation
 Biofilms