Most natural microbial systems have evolved to function in environments with temporal and spatial variations. A major limitation to understanding such complex systems is the lack of mathematical modelling frameworks that connect the genomes of individual species and temporal and spatial variations in the environment to system behaviour. The goal of this review is to introduce the emerging field of spatiotemporal metabolic modelling based on genome-scale reconstructions of microbial metabolism. The extension of flux balance analysis (FBA) to account for both temporal and spatial variations in the environment is termed spatiotemporal FBA (SFBA). Following a brief overview of FBA and its established dynamic extension, the SFBA problem is introduced and recent progress is described. Three case studies are reviewed to illustrate the current state-of-the-art and possible future research directions are outlined. The author posits that SFBA is the next frontier for microbial metabolic modelling and a rapid increase in methods development and system applications is anticipated.
Metabolic flux balance analysis
Whole cell metabolic modelling is highly challenging due to the large-scale and interconnected nature of microbial metabolic pathways. Rigorous modelling of metabolism requires descriptions of enzyme synthesis, kinetics and regulation at the individual reaction level. With the notable exception of primary metabolic pathways in model organisms [1–3], the development of such kinetic models is currently impractical due to lack of in vivo data on enzyme kinetics and regulation. As a result, a less-detailed modelling approach commonly termed flux balance analysis (FBA), based purely on metabolite mass conservation and reaction stoichiometry, has emerged as the dominant methodology for describing whole cell metabolism [4,5]. A stoichiometric model can be constructed by incorporating reactions derived from an annotated genome and augmented with additional ‘gap filling’ reactions required to complete otherwise incomplete pathways [6–9]. Due to the key role of the annotated genome in reaction identification, the resulting model is usually called a genome-scale metabolic (GSM) reconstruction. GSMs are now available for dozens of microbial organisms [10–12] and new GSMs are published on a monthly basis.
Given the stoichiometric matrix and a specified set of available nutrient uptake rates, the goal is to solve the GSM for the unknown intracellular reaction rates (i.e. intracellular fluxes) and the metabolic byproduct secretion rates (i.e. secretion fluxes). Because the stoichiometric equations invariably contain more unknown fluxes than mass balance equations, the GSM has an infinite number of solutions and the flux distribution cannot be uniquely determined. To overcome this limitation, cell metabolism is assumed to be regulated to achieve some type of cellular objective [13,14]. The most common objective is maximal growth, which requires the incorporation of a biomass stoichiometric equation in terms of biomass precursors (e.g. RNA, DNA, proteins, carbohydrates, etc.) . The combination of the stoichiometric equations and the growth rate maximization objective produces a linear programming (LP) problem, which can be efficiently solved to generate predictions of the intracellular flux distribution, uptake and secretion (exchange) fluxes and the growth rate (Figure 1). FBA and other so-called constraint-based computational tools [16,17] have been widely used to analyse the metabolism of both wild-type [18–22] and engineered [23–27] microbial strains.
Dynamic flux balance analysis
FBA is based on assumptions that both intracellular metabolism and the extracellular environment are time invariant. Whereas removal of the intracellular steady-state assumption requires the incorporation of enzyme kinetics, the steady-state assumption on the extracellular environment is more easily relaxed through an extension of FBA commonly termed dynamic FBA (DFBA) [28,29]. A DFBA model is formulated by combining a GSM of intracellular metabolism with kinetic expressions for the uptake rates of growth limiting nutrients and dynamic mass balance equations for cellular biomass, limiting nutrients and secreted metabolic by-products under the assumption that cells rapidly equilibrate to environmental changes (Figure 2). Although difficult to verify, this assumption has been accepted for many systems including batch biochemical reactors where nutrient concentrations vary due to cellular consumption and possibly nutrient feeding . A major advantage of DFBA is that the model outputs include biomass and extracellular metabolite concentrations rather than just the growth rate and exchange fluxes as in FBA . Furthermore, the extracellular concentrations as well as intracellular fluxes are predicted with temporal resolution.
A DFBA model comprises of ordinary differential equations (ODEs) describing the extracellular environment and a LP describing intracellular metabolism. As compared with FBA, a notable disadvantage of DFBA is the computational complexity of this hybrid ODE–LP system. Initial attempts to solve DFBA models were based on sequential strategies in which the LP was repeatedly solved according to a pre-specified time step and the ODEs were integrated between time steps using the most current LP solution [28,31]. Due to the unpredictable accuracy and stability of these sequential methods, more current solution techniques are typically based on simultaneous strategies in which the LP is embedded with the ODE solution [32,33]. DFBA has emerged as a standard tool for the analysis [34–39] and optimization [32,40–42] of individual species as well as microbial communities [43–45] in dynamic environments. Although some DFBA models can be effectively solved by a straightforward combination of ODE and LP solvers, more advanced techniques may be required due to complications associated with the LP unpredictably becoming infeasible or producing alternative optima characterized by non-unique exchange fluxes. Sophisticated simulation methods that explicitly address and effectively overcome these complications are now available [46,47].
Spatiotemporal flux balance analysis
Although DFBA accounts for the effects of extracellular dynamics on intracellular metabolism, the method is based on the assumption that the extracellular environment is well mixed and spatially homogeneous. Many engineered microbial systems, such as biochemical reactors, are carefully designed to achieve spatial homogeneity through liquid mixing . However, most naturally occurring microbial systems exist in spatially heterogeneous environments that also exhibit time variations. The presence of spatial heterogeneity plays an essential role in the evolution and function of natural microbial species [49–52]. Perhaps the most common and important example of spatially heterogeneous environments is that established due to biofilm formation and development (Figure 3a) [53–56]. Concentration gradients in key nutrients due to limited diffusion establish unique metabolic niches within the biofilm that produce spatial variations in biomass density in the case of single-species biofilms  and, additionally, spatial partitioning of species in the case of biofilm consortia . The development of metabolic models that captures such spatial and temporal variations is important to analyse and manipulate complex microbial systems.
Spatial heterogeneities in microbial biofilms
Mathematical models that account for both spatial and temporal variations are commonly termed spatiotemporal models [58,59]. Therefore, we refer to the extension of DFBA to include spatial heterogeneity in the environment as spatiotemporal FBA (SFBA). The author believes that SFBA is the next frontier for microbial metabolic modelling and that the nascent activities to date will soon be followed by a rapid increase in methods development and system applications.
A SFBA model is formulated by replacing the time-varying ODEs in DFBA with partial differential equations (PDEs) expressed in terms of time and some spatial co-ordinate(s) as independent variables . The PDEs represent extracellular mass balance equations for biomass, metabolite and possibly other chemical species concentrations and account for the transport mechanisms that induce the spatial variations, which commonly include metabolite diffusion and liquid/gas phase convection. Boundary conditions must be imposed at the boundaries of the spatial domain to ensure that the PDEs are well posed. As in DFBA, the GSM and the extracellular mass balance equations are linked through nutrient uptake kinetics.
Some systems, such as microbial biofilms, involve mixed boundary conditions when nutrients enter and/or by-products exit the biofilm at different locations [60,61]. Furthermore, the spatial domain may shrink and/or grow with time due to cellular growth and/or death . For the sake of illustration, consider a single species biofilm (Figure 3b) with a fixed thickness L in which a single growth limiting nutrient is available at the bottom of the biofilm (z=0) and a single synthesized byproduct exits the top of the biofilm (z=L). The PDEs describing biofilm diffusional processes can be written as follows assuming that spatial variations occur only in the axial direction z of the biofilm,
Here X(z,t), S(z,t) and P(z,t) represent the biomass, substrate and byproduct concentrations at location z and time t. The growth rate μ, substrate uptake rate vS and byproduct synthesis rate vP are obtained by solution of the GSM. The substrate and byproduct diffuse through the biofilm with efficient diffusion coefficients DS and DP respectively, whereas the biomass is assumed to be non-motile. The second column lists the boundary conditions at the bottom of the biofilm, which are based on the assumptions that substrate is available at a concentration S0 and the biomass and byproduct do not flux across this boundary. The boundary conditions at the top of the biofilm listed in the third column are based on the assumptions that the biomass and substrate do not flux across this boundary whereas the byproduct removal rate is sufficiently high that the boundary concentration is zero. Finally, the initial conditions in the fourth column are based on the assumption that the biofilm is spatially uniform at t=0.
A SFBA model comprises PDEs describing the extracellular environment and LPs describing intracellular metabolism . Pure species systems have a single LP , whereas multispecies systems require solution of a LP for each species . Because no computational methods exist to directly solve such hybrid PDE–LP models, the PDE transport behaviour must be approximated in some manner to generate a solvable model. To date, three alternative methods have been proposed for implementing this approximation: (M1) table lookups of computed FBA solutions combined with integration of the PDEs on a coarse spatial grid [61,63,65]; (M2) real-time FBA solution combined with lattice-based descriptions of metabolite transport [62,64]; and (M3) spatial discretization of the PDEs followed by time integration of the resulting ODE–LP system (Figure 4) [60,66]. These methods differ according to whether the LP solutions are pre-computed (M1) or generated in real-time (M2, M3) and whether the PDEs are discretized directly (M1, M3) or approximated indirectly (M2). Although exhaustive studies of these methods have not yet been performed, the most appropriate method is very likely to be problem dependent.
Genome-scale SFBA 
To illustrate the current state-of-the-art, three representative studies that represent the breadth of the spatiotemporal modelling methods, described above, are reviewed. In the first study , method M1 was used to study electricity generation in a microbial fuel cell where the anaerobic bacterium Geobacter sulfurreducens formed a biofilm on the anode. A metabolic model was developed to analyse the biofilm that oxidized acetate (the electron donor) to carbon dioxide with electron transfer to Fe(III) (the electron acceptor) mimicking the anode. Cellular growth was assumed to create an axial velocity that drove biofilm expansion and along with acetate diffusion created spatial heterogeneities across the biofilm. The steady-state assumption was invoked such that time was eliminated as an independent variable and the transport equations only involved spatial derivatives in the axial direction, thereby simplifying the spatiotemporal model to a time-invariant spatial model. The spatial model consisted of a published G. sulfurreducens GSM , Michaelis–Menten uptake kinetics for acetate and Nerst kinetics for Fe(III) reduction and reaction–convection–diffusion type equations for the fractions of active, respiring and inert biomass, acetate concentration, current density, local potential and axial velocity.
The G. sulfurreducens GSM was solved off-line at different values of the acetate uptake and Fe(III) reduction rates to generate a 20×20 lookup table of LP solutions. The boundary value ODEs representing the reaction–diffusion–convection equations were solved using a continuation method  with the LP lookup table queried as needed to resolve intracellular metabolism at different points in the biofilm. Unfortunately, details concerning the continuation solution method and the spatial grid used were not reported. The spatial metabolic model was solved for varying extracellular conditions (i.e. the presence/absence of NH4) and maintenance energy demands to investigate their effects on biofilm thickness and electrical current production. Simulations combined with validation experiments provided key predictions into biofilm and fuel cell behaviour including that: (1) limited acetate diffusion through the biofilm induced low local acetate concentrations that restricted biomass formation and current generation; and (2) respiring cells that did not grow but produce current located in acetate limited regions had a substantial impact on fuel cell performance.
In the second study  based on method M2, a spatiotemporal modelling framework called Computation of Microbial Ecosystems in Time and Space (COMETS) was developed to investigate the emergent behaviour of two- and three-member synthetic communities. A 2D lattice was defined, with each box within the lattice representing a distinct spatial location. GSMs of the participating species were used to solve an independent DFBA problem for each box assuming that each species had the potential to consume available carbon sources according to the same uptake kinetics. A sequential solution strategy was used where the species LPs were solved to generate local growth rates and fluxes and then the extracellular ODEs were integrated over a fixed time step with these constant LP solutions to generate local biomass and metabolite concentrations. Before resolving the LPs to start the next iteration, 2D diffusion equations were solved over the same time step to allow species cell mass and metabolites to diffuse between boxes.
COMETS was used to investigate engineered metabolite cross-feeding in a two-species synthetic community consisting of Salmonella enterica and Escherichia coli. The spatiotemporal model was able to reproduce experiments showing the two species robustly co-existed and converged to a population fraction of 79% E. coli regardless of the initial fractions. Similar results were obtained for a considerably more complex three-species synthetic community consisting of S. enterica, E. coli and Methylobacterium extorquens. The two-species system was used to investigate the impact of spatial structure. The model correctly predicted that community growth would diminish as the two species were inoculated further apart from each other due to diffusional limitations of the cross-fed metabolites. Moreover, the model reproduced experiments showing community growth would be enhanced by placement of a second colony of engineered S. enterica between the S. enterica and E. coli populations but that placement of a wild-type S. enterica colony which did not require the cross-fed metabolite from E. coli would reduce community growth.
In the third study , method M3 was used to predict synthesis gas (CO, H2) conversion into ethanol by the anaerobic bacterium Clostridium ljungdahlii in a vertical bubble column reactor. Synthesis gas and liquid media feed streams were introduced into the bottom of the column and flowed up the column with different velocities, producing large spatial gradients due to cellular growth and gas depletion. The metabolic by-products ethanol and acetate were recovered in the liquid phase stream exiting the top of the column. The spatiotemporal metabolic model consisted of a published C. ljungdahlii GSM , Michaelis–Menten uptake kinetics for CO and H2 and reaction–convection type PDEs for C. ljungdahlii biomass, liquid phase CO, H2, ethanol and acetate, and gas phase CO and H2. The PDEs were discretized with 100 spatial node points and lexicographic optimization  with six LPs at each node point was used to ensure unique exchange fluxes. The large discretized model consisting of 900 ODEs in time and 600 LPs was efficiently solved within MATLAB (MathWorks) using the DFBAlab (DFBA laboratory) tool (Figure 4) .
Dynamic simulations performed for a wide range of column operating conditions and nutrient uptake parameters generated predictions about column behaviour and bottlenecks to ethanol production consistent with available experiments including: (1) typical CO rich syngas will produce substantial acetate due to H2 depletion in the upper part of the column, suggesting that H2 augmentation of the syngas feed may be beneficial; (2) efficient gas–liquid mass transfer is critical to achieve high ethanol production and high conversions, demonstrating the need for continued development of advanced bubble column designs that achieve very high gas–liquid mass transfer rates; and (3) enhanced H2 uptake rates substantially increase the ethanol titre and the ethanol–acetate ratio, suggesting that C. ljungdahlii engineering efforts should focus on increasing H2 uptake rates.
Spatially heterogeneous environments are the rule rather than the exception for naturally occurring microbial systems. Our continued ability to fundamentally understand and rationally manipulate microbial systems will depend on the development of predictive modelling frameworks that connect the genetic capabilities of individual species and the spatiotemporal characteristics of growth environments to system function. Over the past five years, spatiotemporal metabolic modelling techniques that leverage the increasing availability of genome-scale reconstructions have been developed and evaluated. Despite their limitations with regard to intracellular organization, dynamics and regulation, these SFBA methods hold great promise for analysis of natural microbial systems as well as the design of engineered systems that exploit the evolved capabilities of microbes to optimally function in highly dynamic and spatially heterogeneous environments [70,71]. Possible research directions for the nascent SFBA field are myriad and include: (1) continued identification and study of microbial systems such as the human microbiome [72,73] and lignocellulosic degrading communities [74,75] that would benefit from spatiotemporal analysis; (2) incorporation of higher fidelity intracellular models that extend beyond just reaction stoichiometry [76,77]; (3) development of alternative methods for formulating SFBA models that are more amenable to efficient numerical solution; (4) development and testing of general purpose software such as DFBAlab  for the solution of the large-scale differential equation and LP systems that result from SFBA models; and (5) experimental testing of SFBA model predictions through the collection of omics data with both temporal and spatial resolution [78–80]. Based on the overarching importance of the problem and the initial successes reported in this review, SFBA can be expected to become the next major frontier for microbial metabolic modelling.
The author wishes to acknowledge assistance from Jin Chen in compiling the references.
This work was supported by the National Institutes of Health [grant numbers 1U01EB019416-01]; and the National Science Foundation [grant number CBET 1511346].
computation of microbial ecosystems in time and space
dynamic flux balance analysis
dynamic flux balance analysis laboratory
flux balance analysis
genome-scale metabolic reconstruction
ordinary differential equation
partial differential equation
spatiotemporal flux balance analysis
Metabolic Pathways Analysis 2015: Held at Bom Jesus, Braga, Portugal., 8–12 June 2015.