Spatiotemporal modeling of microbial metabolism

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 genome-scale metabolic reconstructions for steady-state 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 genome-scale 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 genome-scale metabolic models with both time and spatial variations and has wide applicability to natural and engineered microbial systems. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0259-2) contains supplementary material, which is available to authorized users.


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 large-scale metabolic networks in model organisms, which paved the way for the extension of FBA to genome-scale metabolic networks [2]. Curated genome-scale metabolic reconstructions are now available for a wide variety of microbial species, with new reconstructions announced on a weekly basis. Because genome-scale 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 steady-state 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 steady-state metabolic models have been addressed through dynamic extensions of stoichiometric models and classical FBA [4][5][6][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 genome-scale 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, genome-scale metabolic reconstructions allow the effects of gene deletions and insertions in mutant strains to be directly investigated. Genome-scale spatiotemporal models have been solved using table lookups of precomputed FBA solutions [11][12][13], lattice based descriptions of nutrient diffusion [14,15] and agent-based 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 genome-scale description of intracellular metabolism and fundamental transport equations for the extracellular environment.

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, ε 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 j-th liquid phase metabolite where v ij is the net flux of metabolite j into the liquid phase from species i, D jL is the liquid-phase 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 gas-phase concentration using Henry's law. The net flux v ij is calculated as the difference between the synthesis rate obtained from the genome-scale metabolic model and the uptake rate calculated from Michaelis-Menten type kinetic expressions [17,18]. The third equation represents a mass balance on the j-th gas phase component where u G is the gas phase velocity, ε G is the gas volume fraction and D jG is the gas-phase 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 liquid-phase metabolites and N P gas-phase 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 right-hand 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 differential-algebraic equation system instead of ODEs with LPs embedded to increase speed. Hierarchical fixed-priority preemptive (lexicographic) optimization is used to determine uniquely all fluxes which appear in the right-hand 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 right-hand side of the ODEs. With lexicographic optimization, the right-hand 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 right-hand side of the ODEs as a unique vector-valued solution with the following order of objectives: appearing in the right-hand 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 vector-valued objective will require at most n + 2 elements to obtain a unique right-hand 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 i-th 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 i-th species was infeasible. This situation occurred when local nutrient uptake rates were insufficient to meet the non-growth ATP maintenance requirements. While this approach had the potential to make the right-hand 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 large-scale production of renewable fuels and chemicals is direct fermentation of Fig. 1 Discretization of the biomass concentration PDE for a single species in Equation 2. The bubble column reactor is divided into sections along the length dimension. Each section is represented by an ODE that has an accumulation term, a source/sink term due to bacterial growth and death, and two convection terms (in/out) 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 model-based techniques for simulating and optimizing these complex multiphase reactors is important to advance syngas fermentation technology.

Bubble column model solution
The bubble column model was formulated by combining a genome-scale 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 liquid-phase concentrations of C. ljungdahlii biomass, ethanol, acetate, CO, H 2 and CO 2 and the gas-phase 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 third-order 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 second-order 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.

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 steady-state solution required only about 8 min running DFBAlab and  Fig. 2. Steady-state 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 frac-tion~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 steady-state 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  (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 ethanol-acetate 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   [26][27][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 ethanol-acetate 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, non-healing 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 tissue-biofilm interface at the bottom of the biofilm and oxygen is primarily available through the biofilm-air 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.

Biofilm model solution
The bacterial biofilm model was formulated by combining a genome-scale 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 genome-scale resolution. Our preliminary FBA calculations showed that the primary byproducts were Fig. 5 Dynamic simulation of the bubble column reactor model for a CO mass transfer coefficient k m,C = 300 h −1 . The first two columns show time resolved predictions at node points in the middle and at the exit of the column, while the third column show spatially resolved predictions for the exit node point at the final time. The H 2 and CO 2 mass transfer coefficients were set to be 2.5k m,C and k m,C , respectively acetate and D-alanine. To obtain better agreement with experimental data [39] showing anaerobic succinate secretion by P. aeruginosa, we constrained the D-alanine 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 liquid-phase 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. Steady-state 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 second-order 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.

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 steady-state 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 large-scale 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 steady-state 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 tissue-biofilm 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 steady-state 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 tissue-biofilm 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 tissue-biofilm 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 steady-state spatial profiles obtained after 1000 h of dynamic simulation for three cases that differ with respect to whether O 2 was supplied at the tissue-biofilm 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 non-homogeneous environments that require metabolic models that account for both temporal and spatial variations. Our spatiotemporal metabolic modeling framework involves combining genome-scale metabolic Oxygen uptake rate Maximize nutrient consumption Figure 6 Dynamic simulation of the bacterial biofilm model at the nominal operating conditions (Table 3) and a width L = 50 μm. Time resolved predictions are shown for nodes points located at the bottom, middle and top of the biofilm 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 genome-scale 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 genome-scale reconstructions to maintain the same uptake-secretion 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.  Spatial profiles obtained at the end of 1000 h dynamic simulations for three scenarios: no O 2 supplied at the tissue-biofilm interface and succinate consumption not allowed ("Nominal"); O 2 supplied at the tissue-biofilm interface by imposing a mass transfer limited boundary condition with a plasma O 2 concentration of 0.05 mmol/L but succinate consumption not allowed ("O2 Tissue"); and O 2 supplied at the tissue-biofilm interface and succinate consumption introduced assuming the same uptake kinetic parameters as used for glucose