Skip to main content
  • Research article
  • Open access
  • Published:

Spatiotemporal modeling of microbial metabolism

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 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.

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 well-mixed, 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 [47]. 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 H2 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 [1113], 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.

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,

$$ \begin{array}{c}\hfill \frac{\partial {X}_i}{\partial t}=\left({\mu}_i-{\mu}_{di}\right){X}_i-\frac{u_L}{\varepsilon_L}\frac{\partial {X}_i}{\partial z}+{D}_{iX}\frac{\partial^2{X}_i}{\partial {z}^2}\hfill \\ {}\hfill \frac{\partial {M}_j}{\partial t}={\displaystyle \sum_{i=1}^N{v}_{ij}{X}_i-\frac{u_L}{\varepsilon_L}\frac{\partial {M}_j}{\partial z}+{D}_{jL}\frac{\partial^2{M}_j}{\partial {z}^2}+\frac{k_j}{\varepsilon_L}\left({M}_j^{*}-{M}_j\right)}\hfill \\ {}\hfill \frac{\partial {P}_j}{\partial t}=-\frac{u_G}{\varepsilon_G}\frac{\partial {P}_j}{\partial z}+{D}_{jG}\frac{\partial^2{P}_j}{\partial {z}^2}-\frac{k_j}{\varepsilon_G}\left({M}_j^{*}-{M}_j\right)\hfill \end{array} $$
(1)

The first equation represents a mass balances on the i-th microbial species where X i is the biomass concentration, μ i is the growth rate obtained from the genome-scale 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 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, \( {\varepsilon}_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 right-hand 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:

  1. 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. 2.

    Maximize growth rate: this is the traditional FBA objective.

  3. 3.

    Maximize/minimize all of the exchange fluxes 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:

$$ \frac{\partial {X}_i}{\partial t}=\left({\mu}_i-{\mu}_{di}\right){X}_i-{u}_L\frac{\partial {X}_i}{\partial z} $$
(2)

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:

$$ \frac{\mathrm{d}{X}_{i,j}}{\mathrm{d}t}=\left({\mu}_i-{\mu}_{di}\right){X}_{i,j}-\frac{u_L}{\varDelta L}\left({X}_{i,j}-{X}_{i,j-1}\right), $$
(3)

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.

Fig. 1
figure 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)

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 non-profit 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 waste gas streams and synthesis gas (syngas; mainly comprised of H2/CO/CO2) 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 H2 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.

1. 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 reaction-convection–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/H2 mixtures were ethanol, acetate and CO2. 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, H2 and CO2 and the gas-phase concentrations of CO, H2 and CO2 (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.

Table 1 Parameter values for the synthesis gas bubble column reactor model

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 first-order 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.

Table 2 Cellular objective functions used for C. ljungdahlii metabolism

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 steady-state solution required only about 8 min running DFBAlab and MATLAB version 7.11 on a Dell XPS laptop with Intel Core i7-2760QM 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. 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 H2 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 CO2 was negligible under these nominal operating conditions. Because the feed gas was relatively rich in CO, the H2 conversion was 62 % while the CO conversion was only 29 %. As a result of H2 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.

Fig. 2
figure 2

Dynamic simulation of the bubble column reactor model at the nominal operating conditions (Table 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

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.

Fig. 3
figure 3

Effect of the number of discretization node points (N) on biomass and ethanol concentration spatial profiles (top) and on biomass and ethanol concentrations exiting the reactor (bottom). The chosen value of N = 100 is indicated by the dashed lines

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/H2 mixture to a 50/50 CO/H2 mixture. The column exhibited similar dynamics for this H2 rich feed, as the biomass concentration still required approximately 200 h to reach steady state (Fig. 4). However, the increased H2 feed concentration produced a more favorable dissolved H2 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 H2 content of the feed, the H2 conversion decreased to 60 % and the CO conversion increased to 35 %. Our model predictions were in qualitative agreement with published experimental studies [2628] showing that hydrogen rich feeds increase both the ethanol titer and the ethanol/acetate split.

Fig. 4
figure 4

Dynamic simulation of the bubble column reactor model for a CO/H2 feed composition of 50/50. 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

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 H2 mass transfer coefficient to be 250 % larger than the CO value and the CO2 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 H2 (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.

Fig. 5
figure 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 H2 and CO2 mass transfer coefficients were set to be 2.5k m,C and k m,C , respectively

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.

1. 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 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 steady-state spatial profiles were less affected by the initial cell concentration.

Table 3 Parameter values for the P. aeruginosa biofilm model

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.

Table 4 Cellular objective functions used for P. aeruginosa metabolism

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 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.

Figure 6
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

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].

Fig 7
figure 7

Dynamic simulation of the bacterial biofilm model at the nominal operating conditions (Table 3) and a width L = 100 μm. Time resolved predictions are shown for nodes points located at the bottom, middle and top of the biofilm

Fig. 8
figure 8

Dynamic simulation of the bacterial biofilm model at the nominal operating conditions (Table 3), width L = 50 μm and a mass transfer limited boundary condition with a plasma O2 concentration of 0.05 mmol/L imposed at the tissue-biofilm interface. Time resolved predictions are shown for nodes points located at the bottom, middle and top of the biofilm

While the previous simulations were performed assuming the only source of O2 was from air at the top of the biofilm, blood plasma has low O2 levels [44] that could support limited aerobic growth near the tissue-biofilm interface. To investigate this effect, we modified the boundary condition at z = 0 for the O2 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 O2 concentration of 0.05 mmol/L. The inclusion of an O2 source at this interface resulted in a higher O2 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 O2 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 O2 was supplied at the tissue-biofilm interface and whether succinate consumption was allowed. If only O2 supply at the tissue-biofilm 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.

Fig. 9
figure 9

Spatial profiles obtained at the end of 1000 h dynamic simulations for three scenarios: no O2 supplied at the tissue-biofilm interface and succinate consumption not allowed (“Nominal”); O2 supplied at the tissue-biofilm interface by imposing a mass transfer limited boundary condition with a plasma O2 concentration of 0.05 mmol/L but succinate consumption not allowed (“O2 Tissue”); and O2 supplied at the tissue-biofilm interface and succinate consumption introduced assuming the same uptake kinetic parameters as used for glucose

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 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.

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 post-genome era. Trends Biochem Sci. 2003;28(5):250–8.

    Article  CAS  PubMed  Google Scholar 

  2. Price ND, Papin JA, Schilling CH, Palsson BO. Genome-scale microbial in silico models: the constraints-based approach. Trends Biotechnol. 2003;21(4):162–9.

    Article  CAS  PubMed  Google Scholar 

  3. Palsson B. Systems biology: properties of reconstructed networks. Cambridge: Cambridge University Press; 2006.

    Book  Google Scholar 

  4. Hanly TJ, Henson MA. Dynamic flux balance modeling of microbial co-cultures for efficient batch fermentation of glucose and xylose mixtures. Biotechnol Bioeng. 2011;108(2):376–85.

    Article  CAS  PubMed  Google Scholar 

  5. Hjersted JL, Henson MA, Mahadevan R. Genome-scale analysis of Saccharomyces cerevisiae metabolism and ethanol production in fed-batch culture. Biotechnol Bioeng. 2007;97(5):1190–204.

    Article  CAS  PubMed  Google Scholar 

  6. Varma A, Palsson BO. Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia-coli W3110. Appl Environ Microbiol. 1994;60(10):3724–31.

    CAS  PubMed  PubMed Central  Google Scholar 

  7. Mahadevan R, Edwards JS, Doyle FJ. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys J. 2002;83(3):1331–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Hjersted JL, Henson MA. Steady-state and dynamic flux balance analysis of ethanol production by Saccharomyces cerevisiae. IET Syst Biol. 2009;3(3):167–79.

    Article  CAS  PubMed  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  11. Cole JA, Kohler L, Hedhli J, Luthey-Schulten Z. Spatially-resolved 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.

    Article  CAS  PubMed  Google Scholar 

  13. Fang Y, Scheibe TD, Mahadevan R, Garg S, Long PE, Lovley DR. Direct coupling of a genome-scale microbial in silico model and a groundwater reactive transport model. J Contam Hydrol. 2011;122(1–4):96–103.

    Article  CAS  PubMed  Google Scholar 

  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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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 co-culture: predicting and optimizing ethanol production from glucose/xylose mixtures. Biotechnol Biofuels. 2013;6:44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Hanly TJ, Urello M, Henson MA. Dynamic flux balance modeling of S. cerevisiae and E. coli co-cultures for efficient consumption of glucose/xylose mixtures. Appl Microbiol Biotechnol. 2012;93(6):2529–41.

    Article  CAS  PubMed  Google Scholar 

  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.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Höffner K, Harwood SM, Barton PI. A reliable simulator for dynamic flux balance analysis. Biotechnol Bioeng. 2013;110(3):792–802.

    Article  PubMed  Google Scholar 

  21. Kirkels AF, Verbong GPJ. Biomass gasification: Still promising? A 30-year global overview. Renew Sust Energ Rev. 2011;15(1):471–81.

    Article  CAS  Google Scholar 

  22. McKendry P. Energy production from biomass (part 3): gasification technologies. Bioresour Technol. 2002;83(1):55–63.

    Article  CAS  PubMed  Google Scholar 

  23. Nagarajan H, Sahin M, Nogales J, Latif H, Lovley DR, Ebrahim A, et al. Characterizing acetogenic metabolism using a genome-scale 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.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Finlayson BA. Numerical Methods for Problems with Moving Fronts. Incorporated: Ravenna Park Publishing; 1992.

    Google Scholar 

  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.

    Article  CAS  Google Scholar 

  29. Munasinghe PC, Khanal SK. Biomass-derived syngas fermentation into biofuels: Opportunities and challenges. Bioresour Technol. 2010;101(13):5013–22.

    Article  CAS  PubMed  Google Scholar 

  30. Bredwell MD, Srivastava P, Worden RM. Reactor design issues for synthesis-gas fermentations. Biotechnol Prog. 1999;15(5):834–44.

    Article  CAS  PubMed  Google Scholar 

  31. Kirketerp-Moller 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:S84-S.

  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.

    Article  PubMed  Google Scholar 

  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.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Stewart PS. Diffusion in biofilms. J Bacteriol. 2003;185(5):1485–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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.

    Article  CAS  PubMed  Google Scholar 

  37. Oberhardt MA, Puchalka J, Fryer KE, Martins dos Santos VA, Papin JA. Genome-scale metabolic network analysis of the opportunistic pathogen Pseudomonas aeruginosa PAO1. J Bacteriol. 2008;190(8):2790–803.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Horn H, Lackner S. Modeling of Biofilm Systems: A Review. Productive Biofilms. 2014;146:53–76.

    Google Scholar 

  39. Eschbach M, Schreiber K, Trunk K, Buer J, Jahn D, Schobert M. Long-term anaerobic survival of the opportunistic pathogen Pseudomonas aeruginosa via pyruvate fermentation. J Bacteriol. 2004;186(14):4596–604.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Beyenal H, Tanyolac A, Lewandowski Z. Measurement of local effective diffusivity in heterogeneous biofilms. Water Sci Technol. 1998;38(8–9):171–8.

    Article  CAS  Google Scholar 

  41. Hall-Stoodley L, Costerton JW, Stoodley P. Bacterial biofilms: From the natural environment to infectious diseases. Nat Rev Microbiol. 2004;2(2):95–108.

    Article  CAS  PubMed  Google Scholar 

  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.

    Article  CAS  PubMed  Google Scholar 

  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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Pittman RN. Oxygen gradients in the microcirculation. Acta Physiol. 2011;202(3):311–22.

    Article  CAS  Google Scholar 

  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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Erdrich P, Steuer R, Klamt S. An algorithm for the reduction of genome-scale metabolic network models to meaningful core models. BMC Syst Biol. 2015;9.

Download references

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

Authors and Affiliations

Authors

Corresponding author

Correspondence to Michael A. Henson.

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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Chen, J., Gomez, J.A., Höffner, K. et al. Spatiotemporal modeling of microbial metabolism. BMC Syst Biol 10, 21 (2016). https://doi.org/10.1186/s12918-016-0259-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12918-016-0259-2

Keywords