Original Research ARTICLE
Regulation-Structured Dynamic Metabolic Model Provides a Potential Mechanism for Delayed Enzyme Response in Denitrification Process
- Pacific Northwest National Laboratory, Richland, WA, United States
In a recent study of denitrification dynamics in hyporheic zone sediments, we observed a significant time lag (up to several days) in enzymatic response to the changes in substrate concentration. To explore an underlying mechanism and understand the interactive dynamics between enzymes and nutrients, we developed a trait-based model that associates a community's traits with functional enzymes, instead of typically used species guilds (or functional guilds). This enzyme-based formulation allows to collectively describe biogeochemical functions of microbial communities without directly parameterizing the dynamics of species guilds, therefore being scalable to complex communities. As a key component of modeling, we accounted for microbial regulation occurring through transcriptional and translational processes, the dynamics of which was parameterized based on the temporal profiles of enzyme concentrations measured using a new signature peptide-based method. The simulation results using the resulting model showed several days of a time lag in enzymatic responses as observed in experiments. Further, the model showed that the delayed enzymatic reactions could be primarily controlled by transcriptional responses and that the dynamics of transcripts and enzymes are closely correlated. The developed model can serve as a useful tool for predicting biogeochemical processes in natural environments, either independently or through integration with hydrologic flow simulators.
Microbes in terrestrial and aquatic ecosystems are primary drivers of biogeochemical processes, including carbon and nitrogen cycling (Gougoulias et al., 2014). Predicting ecosystem functions and dynamics in a varying environment therefore requires a mechanistic understanding of how microorganisms interact with each other to catalyze biogeochemical functions. Microbial communities are complex adaptive systems; their dynamics and emergent properties are difficult to understand and predict (Konopka et al., 2015). A long time lag in the microbial responses to changes in local environment is an example of complex community dynamics commonly observed in the field and experimental studies, but the underlying reasons are poorly understood (Wood et al., 1995; Shade et al., 2012). Metabolic time lags in microbes can affect the fate and transport of microbial substrates, particularly when the timescales of both metabolic lag and transport are comparable (Nilsen et al., 2012). Currently, there are no fully mechanistic approaches for incorporating metabolic lags into reactive transport models. Existing models account for metabolic lags by using temporal convolution integrals (Nilsen et al., 2012) or by introducing exposure time as an additional dimension (or coordinate; Wood et al., 1995). These models, however, do not address the fundamental mechanisms that lead to such delayed metabolic responses; thus, limiting their use for prediction.
In a recent experimental study (Li et al., 2017), we reported the complex dynamics of microbial communities that drive denitrification in hyporheic zone sediments from the Hanford Reach of the Columbia River in Washington state. The hyporheic zone is a biologically active domain where (nitrogen-limited, but carbon-abundant) surface water and (carbon-poor, but nitrogen-enriched) groundwater mix (Stegen et al., 2016). Denitrification in hyporheic zone sediments is mediated by microbial communities through multiple stages of reduction reactions from to N2, with each step catalyzed by distinct enzymes. Our experimental data showed a significant time lag in in the synthesis of enzymes that catalyze reduction to . Microbial communities continued to express catalyzing enzymes even after was depleted. Enzyme concentrations did start to decrease only several days after reduction was completed and were maintained at base levels thereafter.
Understanding the processes governing delayed enzyme synthesis and catalysis and its translation into quantitative modeling are essential for predicting biogeochemical dynamics, particularly in dynamic environments such as the hyporheic zone where chemical conditions can change rapidly in response to hydrologic exchange. For these purposes, we developed a new concept of modeling that accounts for the regulated synthesis of enzymes through transcriptional and translational processes. The developed model collectively describes biogeochemical dynamics in terms of functional enzymes synthesized in a community, without directly considering the dynamics of individual species or their guilds that synthesize those enzymes. Using this model, we explored how the dynamic interplay between regulatory molecules such as transcripts, enzymes, and internal resource molecules affects biogeochemical functions. The resulting model not only provided an excellent fit to data (both nutrient and enzyme dynamics), but also suggested a complex regulatory machinery in microbes as a plausible mechanism for the observed time lags in enzyme responses. The model also provided novel insights into the biological roles of transcripts as a control point of dynamic metabolic regulation.
The microbial community model developed in this article provides several new conceptual components for microbial community modeling, which we summarize in the following section as an essential pre-requisite for understanding the idea of the proposed model framework and the implications of the simulation results.
Functional Enzyme-Based Formulation
Modeling environmental microbial communities is challenging, primarily due to immense complexity in their compositional diversity (Song et al., 2014). A common idea to alleviate this structural complexity is trait-based modeling, which maps organism-based information into a functional space (Allison, 2012; Boon et al., 2014). The trait-based modeling often considers grouping organisms that share certain metabolic functional similarities (i.e., traits) into a fewer number of species guilds or functional guilds (Taffs et al., 2009; Jin and Roden, 2011; Bouskill et al., 2012). This guild-based grouping can be less effective, however, in the case where microorganisms can perform multiple functions (functional versatility) that are partially overlapping with one another (functional degeneracy) in a dynamically changing environment (Whitacre and Bender, 2010; Song et al., 2015). As a consequence, the number of functional guilds (and accordingly, the number of parameters) can become significantly large.
As an alternative, the complexity of microbial communities can be reduced by associating a community's metabolic traits with functional enzymes. This enzyme-based approach views a microbial community as a collective assembly of metabolic capabilities as opposed to a complex network of individual organisms or their groups (e.g., functional guilds), therefore requiring no a priori knowledge of the functional roles of individual organisms.
Figure 1 illustrates how the enzyme- and guild-based formulations can be different, particularly with respect to model structure and parameterization. Figure 1A shows two cases where a microbial community is composed of (1) single-function organisms only (top) or (2) a combination of single- and multi-function organisms (bottom). These functions denote the microbes' capability of expressing enzymes that catalyze specific biogeochemical reactions such as those in denitrification. Therefore, single-function organisms (specialists) express one functional enzyme, while multi-function organisms (generalists) express multiple enzymes. The functional guild-based approach groups organisms based on the type and number of metabolic functions each guild catalyzes. For the community composed of specialists only, this approach can successfully reduce nine species into three guilds that express the same enzymes. In this case, the enzyme-based formulation is not much different from the species guild-based approach. If the community contains member organisms expressing multiple functional enzymes, the guild-based reduction results in the expanded number of guilds, implying the increase in the number of model parameters. In contrast, the enzyme-based approach reduces the model complexity by focusing on the functional enzymes independent of associated guilds. As a result, the functional enzyme-based approach only requires the identification of the dynamics of three enzyme classes regardless of the metabolic characteristics of organisms. Due to its focus on enzyme pool, the enzyme-based model carries no information on whether a given function is occurring within one or multiple organisms.
Figure 1. A schematic representation of trait-based microbial community modeling strategies: (A) functional guild- vs. functional enzyme-based formulations, (B) enzyme-based modeling of denitrification. In (A), communities are considered as being composed of specialists (that can express one functional enzyme) only, or being mixed with specialists and functionally versatile organisms (that can express multiple enzymes). In (B), the latter case is illustrated using denitrification, which is composed of four steps of reduction catalyzed by distinct sets of genes and enzymes: NapA and NarG for reduction to , NirS and NirK for to NO, Nor for NO to N2O, and NosZ1 and NosZ2 for N2O to N2.
Figure 1B provides a schematic of the functional enzyme-based approach for denitrification, which is composed of four catalytic steps that reduce to N2. The functionally versatile organisms (generalists) can contribute to more than one catalytic step as implied by the multiple arrows. Since the enzyme functions are directly modeled, the enzyme-based formulation does not have to track which organisms are involved in each step; thereby reducing the number of parameters to describe the community dynamics in comparison to the guild-based formulation.
The Cybernetic Description of Metabolic Regulation
The enzyme-based formulation describes an entire microbial community as a single unit that responds to the environment by regulating its metabolism through expression of genes and enzymes. As a key advantage, this approach facilitates the application of existing dynamic modeling platforms that have been developed to predict metabolic behaviors of single organisms (Song et al., 2014; Vasilakou et al., 2016). In this regard, the cybernetic approach developed by Ramkrishna and co-workers (Ramkrishna, 1983; Ramkrishna and Song, 2012; Young, 2015) provides an ideal platform for the enzyme-based microbial community modeling due to its rational description of metabolic regulation through the control of enzyme syntheses and activities.
The cybernetic model postulates organisms as teleonomic systems that regulate metabolism toward achieving a certain metabolic objective (e.g., maximization of nutrient uptake rate or growth rate). This postulate removes the need to account for detailed mechanistic information on regulation, which is generally unknown for environmental organisms. The term “cybernetics” denotes a goal-seeking aspect of metabolism (Wiener, 1961; Mayr, 1988). The predictive power of cybernetic approach has been successfully demonstrated over the past few decades through a variety of challenging case studies of modeling metabolic switches in varying environments (Dhurjati et al., 1985; Kompala et al., 1986; Turner et al., 1988; Alexander and Ramkrishna, 1991; Baloo and Ramkrishna, 1991; Straight and Ramkrishna, 1991; Varner, 2000; Kim et al., 2008; Young et al., 2008; Song et al., 2009, 2013b; Song and Ramkrishna, 2010, 2011; Franz et al., 2011), dynamic response to genetic perturbations (Young et al., 2008; Song and Ramkrishna, 2012), and nonlinear behavior such as multiple metabolic states (Namjoshi et al., 2003; Kim et al., 2012; Song and Ramkrishna, 2013).
Accounting for the Dynamic Interplay among Regulatory Molecules
In a previous study, we applied the cybernetic approach to model the dynamics of denitrification process driven by single organisms (Song and Liu, 2015). The focus of this previous work was to simulate sequential utilization of multiple electron acceptors associated with denitrification and did not include enzyme measurements for model parameterization and validation. Due to the relatively simple description of the regulatory process, this form of cybernetic model is inappropriate for predicting a long enzymatic lag we observed in a denitrifying microbial community. Therefore, in this work, we significantly extended the previous model by elaborating cellular regulation program by accounting for the dynamic interplay among transcripts, enzymes, and internal resources. Following the central dogma of molecular biology (Crick, 1970), we modeled the regulation process in microbes based on transcription and translation: (1) the synthesis of transcripts from DNA and other internal resource molecules (such as ATP, NADH, DNA polymerases, RNA polymerases, ribosomes, etc.), and (2) the subsequent synthesis of enzymes from transcripts and internal resources. Note that the coupling between transcripts and enzymes is not occurring in a linear fashion due to the involvement of internal resources in both transcription and translation. We have termed this approach the Regulation-Structured Cybernetic Model (RSCM). Figure 2 illustrates the control structure of the RSCM for two competing reactions catalyzed by distinct enzymes.
Figure 2. A schematic representation of the control structure in the RSCM. The variable ui controls the synthesis of resource (ρi) based on the contribution of individual reactions to the postulated metabolic objective (i.e., return-on-investment pi, e.g., pi = ri); the synthesis of transcript (τi) is subsequently determined by ρi; the synthesis of ei is determined by both τi and ρi; the enzyme ei, in turn, catalyzes reaction (ri).
The idea of elaborating the cybernetic regulation model at the transcriptional and translational levels was also considered previously in the literature (Varner, 2000), the regulation structure of which is, however, much more complex than the RSCM due to the direct consideration of the control of inductive synthesis of both transcripts and enzymes. In contrast, the RSCM indirectly realizes those controls only at the resource generation stage (see Figure 2 and model equations in the next section). As another key distinction, the RSCM accounts for the dynamics of internal resources and their dynamic interplay with transcripts and enzymes, which was neglected in the previous work.
We structured and parameterized the denitrification model based on the experimental data collected by the Pacific Northwest National Laboratory (PNNL)'s subsurface biogeochemistry research group. As all experimental details on data collection and enzyme assay are published elsewhere (Li et al., 2017), we here provide only a brief summary. Batch denitrification experiments under anaerobic conditions were performed using the Columbia River hyporheic zone sediments collected from the Hanford Reach. These sediment samples were mixed with synthetic groundwater. Experimental data measured for model parameterization include the concentrations of gas phase nitrous oxide (N2O), dissolved nitrate (), nitrite (), acetate, dissolved inorganic carbon (DIC), and functional enzymes. Concentrations of functional enzymes were quantified using the PNNL-developed signature peptide-based method. The targeted quantification of functional enzymes include dissimilatory reductases (NapA and NarG) and dissimilatory N2O reductase (NosZ), which represent the first and last steps in denitrification. Other enzymes associated with intermediate steps (such as NirS/NirK and Nor) were not measured. Nutrient concentrations were determined by analyzing three replicates. The enzyme (NarG, NapA, NosZ1, and NosZ2) concentrations were obtained as the average of their two signature peptide concentrations: each signature peptide was determined from two replicates. Dissimilatory nitrate reduction to ammonium (DNRA) is another potential reaction pathway for the anaerobic reduction of nitrate, but the experimental data by Li et al. (2017) provide several evidences indicating that denitrification is a major reaction pathway active in their experiments. First of all, the ammonia assay detected no formation of ammonium, which is the end product of DNRA. From the same samples, however, the formation of N2O (an intermediate during denitrification) was detected. While it is known that N2O could also be released as a byproduct of DNRA (Stevens et al., 1998; Kraft et al., 2011), this could also be an evidence with no ammonium detected. Finally, Li et al. measured NosZ1 and NosZ1 enzymes that catalyze the reduction of N2O to N2, the final reduction step in denitrification. The temporal profiles of NosZ1/NosZ2 enzyme concentrations were qualitatively similar to those of NarG/NapA enzymes, indicating their coupling as the first and last steps of denitrification.
Denitrification Reaction Network and Stoichiometry
Based on the experimental data by Li et al. (2017) described above, we developed a biogeochemical model with a focus on denitrification. Denitrification is an anaerobic process in which is reduced to N2 gas. The reduction occurs in four steps, as written below for one electron equivalent (e− eq) (Rittmann and McCarty, 2001):
We normalized the above reactions with respect to one mole of dissolved organic carbon (CH2O), so that the coefficients will represent the stoichiometric relationships between reactants and products for a unit mole of consumed organic carbon, as follows:
Accurate parameterization of some of the intermediate reactions was difficult because (1) the concentration of N2O was very low (i.e., of the order of 10−6–10−3 mM), indicating its consumption is faster than generation, and (2) no NO measurement was available. Therefore, we simplified the denitrification process into two step reactions by assuming relatively fast dynamics of NO and N2O, i.e.,
We also consider the synthesis of biomass (C5H7O2N) as follows:
Microbes obtain energy through anaerobic respiration pathways, reactions (3) or (4) or both, depending on which electron acceptors ( or ) are available in the environment. To account for this coupling, we combined equations (3) and (4) with (5) as follows:
where f1 and f2 denote the fractions of energy production associated with and . Note that reactions and share the same electron donor (CH2O), but involve different electron acceptors (i.e., and ), thus representing alternative pathways for the production of biomass (C5H7O2N). In the above reactions, we excluded H+ and H2O to focus on carbon and nitrogen balances. For simplicity, hereafter we use DOC (dissolved organic carbon), DIC (dissolved inorganic carbon) and BM (biomass) to denote CH2O, CO2 and C5H7O2N.
Mass Balances of Nutrients and Biomass
In a homogeneous batch reactor, mass balances of substrates and biomass can be written based on stoichiometric equations (6) and (7) as follows:
In the above equation, x denotes concentrations of reactants and products, and r1 and r2 are reaction rates associated with and , respectively. We modify the mass balances of BM and DOC to account for biomass degradation and the impact on the formation of organic carbon as follows:
where kdeg denotes the rate of biomass degradation. The coefficient factor “5” on the right hand side of Equation (10) denotes the carbon-based stoichiometric relationship between BM (i.e., C5H7O2N) degradation and the resulting increase of DOC (i.e., CH2O).
Description of Metabolic Regulation
It is often that simple Monod-type kinetic equations do not properly describe biogeochemical reactions (Tang and Riley, 2013), which might be ascribed to the lack of description for the regulated metabolic behaviors of microbial communities. Regulation is a hallmark of microbial metabolism, which therefore should be a key component of biogeochemical modeling. Microbes regulate metabolism through the control of enzyme syntheses (and their activities). To account for those regulatory processes, we modeled the DOC uptake rate through the reaction i (ri) as being regulated by the enzyme concentration (ei), i.e.,
where denotes the kinetic form of (unregulated) reaction rate. The synthesis of enzyme requires the energy and material resources (i.e., internal resources) such as ATP, NADH, DNA polymerases, RNA polymerases, ribosomes, and so on. Due to their limited availability, microbes allocate resources among different enzymes to selectively activate certain reactions. The cybernetic model provides a rule for this resource allocation based on “the return-on-investment” concept, i.e., how much microbes will gain profit (such as nutrient uptake rates or growth rates) by investing internal resources for the synthesis of certain enzymes.
The RSCM applies the same rule for the resource allocation, but in a more structured way. That is, the model accounts for the resource balance in combination with the two steps of enzymes synthesis: transcription (i.e., the information flow from genes into transcripts) and translation (i.e., from transcripts to enzymes). With a focus on the dynamics of transcripts and enzymes, we treated the genes as part of the resource pool to minimize the number of parameters. The resulting model contains mass balance equations for the three major components associated with regulation: the resource pool (denoted by R), transcripts (T), and enzymes (E). The equations are given as follows:
where ρi and τi denote the concentrations of the resource pool and transcripts associated with the reaction i; rR,i, rT,i, and rE,i represent the kinetics of synthesis of resource pools, transcripts and enzymes, respectively; αR,i is the constitutive synthesis rate of resource pool i; βR,i, βT,i, and βE,i denote the parameters of degradation rates of resource pool Ri, transcript Ti, and enzyme Ei, respectively. The three terms on the right hand side of Equation (12), respectively, denote the constitutive synthesis rate, inductive synthesis rate, and the degradation rate of resource pool associated with the reaction i, where the variable ui denotes the control of inductive synthesis of resource. Similarly, the two terms on the right hand side of Equations (13) and (14) denote the synthesis and degradation rates of transcripts and enzymes. Note that the resource allocation is determined by the cybernetic variable ui in Equation (12), which subsequently dictates the syntheses of transcripts and enzymes as in Equations (13) and (14).
The control variable ui in Equation is determined by the cybernetic control law termed the matching law (Young and Ramkrishna, 2007) as follows:
where pi is the return-on-investment (or profit), which reflects the expected metabolic benefit derived from the allocation of internal resources for the syntheses of transcripts and enzymes and for the subsequent catalysis of specific biogeochemical reactions (see Figure 2). In other words, the variable pi denotes the contribution that the reaction i makes to the postulated metabolic objective function of the microbial community. The typical choice of this objective includes maximization of carbon uptake rate and growth rate; we set pi to ri (i.e., carbon uptake rate) in this work. The functional form of ui in Equation (15) was heuristically determined in earlier cybernetic models, but was later theoretically derived from an optimal control theory (Young and Ramkrishna, 2007).
The typical cybernetic modeling formulation adds the post-translational control of enzyme activities as well, which was not accounted for in the current form of RSCM due to its relatively low impact on the prediction as discussed in our previous denitrification modeling analysis (Song and Liu, 2015). However, the enzyme activity control is an important regulatory mechanism in general (Oliveira and Sauer, 2012) and may need to be incorporated into the model in certain circumstances, e.g., in the case that the system exhibits fast nonlinear regulatory dynamics over a short time period, because it enables the prompt control of metabolic shifts through allosteric hindrance of enzyme activities.
We used Monod-type equations to describe the unregulated kinetic rates (i.e., ) of and in Equation (11), i.e.,
A similar form of kinetic equations was considered for describing the synthesis rates of resources, transcripts, and enzymes, i.e.,
where the k'S and K's denote reaction rate and half-saturation constants associated with specific reactions as implied by their subscripts. It is important to note that (1) any alternative forms of kinetic models can be used in the cybernetic modeling formulation, and (2) even in the case where simple Monod-type kinetics were used for (as above), the actual kinetics in Equation (11) show much more complex dynamics due to the regulatory processes represented in Equations (12–15).
We determined key parameters of the denitrification model through the data fit such that the sum of the squared errors between simulations and measured data is minimized. Experimental data used for this purpose included temporal profiles of substrates (, , DOC, and DIC) and enzymes (NarG, NapA, NosZ1, and NosZ2). Specifically, the catalysis of the reduction from to in Equation (6) was parameterized in association with the dynamics of NarG/NapA enzyme concentrations. The reaction from to N2 is originally composed of three reduction steps catalyzed by NirS/NirK, Nor, and NosZ1/NosZ2, respectively (Figure 1B), but as mentioned earlier, we associated it only with the dynamics of NosZ1/NosZ2 enzymes that catalyze the last step of reduction, due to the lack of other enzyme measurements. The fitted parameters include the rate constants of carbon uptake (k1, k2), the rate constants of the syntheses of internal resources (kR,1, kR,2), transcripts (kT,1, kT,2), and enzymes (kE,1, kE,2), the biomass degradation rate (kdeg), and the fraction of energy production associated with reduction to (f1). Due to the limited availability of experimental data, the parameters associated with the dynamics of internal resources and transcripts (kR,1, kR,2, kT,1, and kT,2) were indirectly identified by fitting the profiles of nutrients and enzymes. Other parameters were fixed based on the literature values (i.e., half saturation constants) or determined by manual adjustment (f2, α 's, and β 's). For fitted parameters, we used the bootstrapping technique to quantify the uncertainty in the estimates based on 95% confidence intervals. Bootstrapping was performed by estimating the best-fit parameter values for 500 experimental data sets, where each data set was created by randomly sampling a uniform distribution of the measured values at each time point. The sampled values varied within ±standard deviation from the mean experimental value. The resulting parameter sets containing values below the 2.5th and above the 97.5th quantile for at least one parameter were removed. The lowest and highest values of each parameter in the remaining parameter sets were taken as the lower and upper 95% confidence limits, respectively. The resulting parameter values (with lower and upper bounds) are summarized in Table 1. The initial concentrations of all variables including nutrients, internal resources, transcripts, and enzymes are shown in Table 2.
Table 2. Initial condition used for the simulation of denitrification dynamics: DOC = dissolved organic carbon, DIC = dissolved inorganic carbon, and BM = biomass; ρ, τ and e denote the concentrations of internal resource, transcript and enzyme, respectively; the subscripts 1 and 2 denote the association with and reduction pathways, respectively.
Fitting Results and Prediction of Time Lags in Enzyme Response
The batch denitrification experiment by Li et al. (2017) with only one initial condition (Table 2) did not provide sufficient data for the accurate determination of all parameters in the RSCM. Thus, we fixed some of the parameters that might not be directly determinable from simple batch data, including half-saturation constants and the parameters associated with constitutive synthesis and degradation of regulatory molecules. This strategy led us to reproducibly determine the remaining parameters through optimal data fit (see Table 1 and Figure 3).
Figure 3. Model fits and predictions of the dynamic changes in (A) substrates and (B) enzymes associated with reduction to (top) and the subsequent reduction to N2 (bottom). Symbols in (A,B) represent the average concentrations of nutrients (DOC, DIC, , and ) and enzymes (NarG, NapA, NosZ1, and NosZ2) obtained from three and two biological replicates, respectively. Error bars denote the standard deviations of each measurement.
The resulting model successfully fitted the overall dynamics of nutrients (Figure 3A) and enzymes (Figure 3B). The model also captured the non-intuitive trend of DOC profile, i.e., the gradual increase of DOC after day 6 (when no further denitrification was actively taking place), by attributing it to biomass degradation (see Equation 10). The model also predicted N2 production on the assumption that all is converted into N2 through the stoichiometric relation given in Equation (7). Figure 3B shows that enzyme responses to the change of and were significantly delayed. Data on the top panel of Figure 3B shows that the reduction rate reached its maximal value around day 4, while the maximal concentrations of catalyzing enzymes (i.e., NarG and NapA) occurred about 4 days later. The RSCM simulated the maximal reduction rate and enzyme concentrations to occur slightly earlier, but its prediction of the time delay in enzyme expression was aligned with experimental observation. The delayed enzymatic response was also inferred for reduction from the observation that the maximal concentrations of NosZ1 and NosZ2 enzymes appeared about 1 day after the depletion of . However, it was difficult to accurately estimate the time delay in this case because of (1) the uncertainty about the exact time when the reduction rate reached its maximum (due to the non-monotonic profile shape), and (2) the absence of the measurement of enzymes other than NosZ1 and NosZ2, which catalyze the last step of denitrification. Continued synthesis of enzymes even after substrate depletion is an interesting phenomenon, which we may describe as dynamic momentum in metabolic control. In addition to delayed enzymatic responses, the data showed that enzyme concentrations (for both and reduction reactions) were maintained from day 10 up to day 20. The RSCM predicted this persistence of enzyme levels as the result of dynamic balance between constitutive enzyme synthesis and degradation rates, as in Equations (12–14).
Dynamic Control of Denitrification Reactions
Figure 4 shows the effects of control realized by the cybernetic variables on the and reduction rates and the biomass yield. On the top panel, the value of u1 is initially high, and gradually decreases. This means that the initial allocation of the resource pool is directed primarily toward promoting reduction; its portion decreases as the amount of increases as implied by the increasing value of u2. The middle panel shows the rates of and reduction as their concentrations change with time. The RSCM predicted these two competitive reactions can take place simultaneously until is depleted. This is an interesting finding from a modeling point of view because the conventional cybernetic control that considers the control of the inductive enzyme synthesis directly in the enzyme balance equation usually predicts competitive reactions (such as consecutive and reduction reactions) to occur in a sequential manner as observed in our previous denitrification model. Prediction of the simultaneous occurrence of competitive reactions by the conventional cybernetic model often requires the consideration of additional reactions, which is not necessary for the RSCM that indirectly accounts for the control of inductive enzyme synthesis through the allocation of internal resources. On the bottom panel, we showed the prediction for dynamic changes in the community's biomass “yield” during denitrification. The biomass yield at each time instant was calculated as the ratio between the biomass production rate and the DOC (or carbon) uptake rate, i.e.,. That is, the biomass yield is a dimensionless quantity representing the amount of produced biomass [mM] per the unit amount of consumed DOC [mM]. The biomass synthesis initially is high due to the higher contribution of reduction (because f1 < f2; see section Methods); gradually decreases for the first 3 days; maintains an almost constant level until day 5.7 where the relative ratio of and reductions does not change much (as shown by the flattened profiles of u1 and u2. The reduction continues until day 6.8 while the biomass yield becomes very small after day 5.7, which is because the high value of f2 is close to 1 (i.e., 0.99) in Equation (7).
Figure 4. Model predictions of the dynamic changes in cybernetic variables (top), reaction rates (middle), and biomass yield (bottom). In the bottom panel, Y1 and Y2 denote the stoichiometric biomass yields associated with and reductions, respectively.
Transcripts Playing a Role as a Dynamic Control Point
In Figure 5, we provided (a) the dynamic profiles of resources, transcripts, and enzymes, and (b) time lags quantified based on the differences between the peak times of each variable (as marked by upside-down triangles in Figure 5A). For the reduction of to , the figure shows that a long lag in enzyme response is fundamentally caused by the interplay of enzymes with transcripts and internal resources, which do not quickly decay even after the substrates are depleted. The same interpretation also applies to the conversion of to N2, but the level of time delay is less pronounced, particularly when considering the fact that the reduction was associated with the dynamics of NosZ1/NosZ2 enzymes that catalyze the last step of reduction (see section Methods). This result—i.e., the prediction of a longer time delay in enzymatic response for the reduction (i.e., 3.4 day) in comparison to the reduction (0.9 day)—agrees with experimental observation, as shown in Figure 3.
Figure 5. Model predictions of (A) the dynamic changes in resources, transcripts and enzymes and (B) time delays from resource generation to transcriptional and enzymatic responses. In (A), the upside-down triangles denote the time where each variable reaches its maximal value.
For the reduction of to , the model predicted a long delay between the resource generation and the transcriptional response (i.e., about 3 days), while predicting a relatively short delay between transcriptional and enzymatic responses (0.5 day) (top panel of Figure 5B). This indicates that the transcriptional response might be a key control point in reduction. However, this was not the case for the reduction where time delays between three processes were comparable (i.e., 0.5 and 0.4 days, respectively) (bottom panel of Figure 5B). Interestingly, time delays between transcriptional and enzymatic responses were predicted to be similar for (0.5 day) and reduction (0.4 day).
This result also reveals the complex interplay between resources, transcripts, and enzymes in denitrification. The dynamic profiles of enzymes (bottom panel of Figure 5A) bear a qualitative resemblance to the profiles of reaction rates shown on the middle panel of Figure 4. That is, profiles in both figures show that initially reduction is dominant but is subsequently exceeded by reduction when is depleted. The initial trend was also the same for resource abundance profiles (top panel of Figure 5A). In contrast, the transcript abundance profiles (middle panel of Figure 5A) are fundamentally different; the abundance of transcripts associated with reduction is relatively low throughout the entire simulation, despite being at similar abundance to reduction transcripts initially. Prediction of these distinct features of transcripts is an outcome of data fit in order to match the overall denitrification dynamics where reduction is more dominant than reduction, while the corresponding enzyme concentrations that catalyze them are comparable. This implies that transcripts might play a role of controlling the overall dynamics of regulation at the interface between resources and enzymes.
Highly Correlated Dynamics of Transcripts and Enzymes
We also analyzed how the abundances of resources, transcripts and enzymes were dynamically inter-related. To remove the effect of time delay, we appropriately adjusted the time scales of the profiles of these three variables to draw phase diagrams between all possible pairs of three variables. Figures 6A,B provide the relationships between resources and transcripts, and between resources and enzymes, respectively; we found the reasonably close relationships between those pairs for reduction, while there was no such correlation for reduction, i.e., transcript and enzyme levels did not proportionally change to the abundance of resources. In Figure 6C, we consistently observed, however, the close relationships between transcripts and enzymes for both and reduction; surprisingly, their relationships were almost linear along both increasing and decreasing phases.
Figure 6. Phase diagrams of resources, transcripts and enzymes: the dynamic relationships (A) between the concentrations of resources and transcripts, (B) between the concentrations of resources and enzymes, and (C) between the concentrations of transcripts and enzymes. The time scales of all data are adjusted to remove the effect of time delay based on the quantification in Figure 5B.
Key Features of the RSCM
We provided a new conceptual platform (termed RSCM) for modeling complex microbial communities. There are several unique features of the RSCM as highlighted in the following. First, the RSCM is built upon the functional enzyme-based approach, which enables the consistent model formulation of complex microbial communities. This enzyme-based parameterization is particularly useful in circumstances where metabolic functions of individual organisms are difficult to uniquely characterize (and therefore difficult to group). Second, the RSCM elaborates regulatory control by accounting for transcription and translation, the products of which can be quantified. Considering the difficulty in obtaining reliable estimates of guild-specific biomass for the function guild-based approach, the ability to directly quantify transcripts and enzyme concentrations is likely an advantage of the enzyme-based approach. As another key distinction, the RSCM describes the dynamic microbial regulation based on the cybernetic approach. In the current practices of biogeochemical modeling, metabolic regulation is often accounted for based on empirical inhibition kinetics, which not only increases the number of parameters to identify, but also may lead to less accurate predictions (Tang and Riley, 2013). In contrast, the cybernetic model can provide reliable predictions over a wider range of conditions with a fewer number of parameters as demonstrated in the past case studies of metabolic modeling (Ramkrishna and Song, 2012; Song et al., 2013a). Data requirement for parameter identification is accordingly lesser for the cybernetic model due to the absence of inhibition-associated parameters, which was necessary for empirical formulation. Based on the principle of Occam's razor (Gauch, 2003), the simplest model is preferred over over-parameterized models, if performance is the same. Accounting for metabolic control based on the cybernetic approach holds a significant advantage in modeling environmental biogeochemical processes, not only that experimental data required for parameter identification is usually insufficient, but also that a priori knowledge is unavailable for inhibition control in a complex system.
Insights into Metabolic Lags and the Role of Transcripts
The RSCM showed that enzymatic time lags observed in denitrification experiments can be caused by the complex interplay among internal resource molecules, transcripts and proteins. The delayed enzymatic response could be ascribed to microbes' regulatory machinery that controls energetically-expensive enzyme synthesis to ensure its survival in frequently varying environment (Ramkrishna et al., 1987). This postulate is rational because immediate and frequent shifts in enzyme settings in response to environmental variation can lead to the rapid depletion of energy required for microbial growth and maintenance due to the cost for protein synthesis (Wessely et al., 2011; Noor et al., 2016).
The RSCM also provided insights in regard to the potential role of transcripts in regulatory dynamics of denitrifying microbial communities. First, it suggests that the overall time delay in denitrification may be primarily controlled by the transcriptional response. Model predictions showed that transcriptional response can be slower than or comparable to enzymatic responses depending on the type of reactions. This implies that accurate simulation of denitrification requires an appropriate consideration of transcriptional dynamics in modeling. Second, the abundances of transcripts may not necessarily reflect the level of reaction activities. The simulation showed that the transcripts associated with reduction were less abundant than those associated with reduction even though the former reaction was dominant until the depletion of . This suggests that within a given reaction, changes in transcript abundances (e.g., through time) may indicate changes in the rate of that reaction, but that such a comparison cannot be made across reactions. In this regard, it would be beneficial to consider collecting additional complementary data (e.g., gene abundances, metabolomic profiles, etc.) and/or performing advanced tracer experiments using isotopically labeled substrates. Third, the model suggests that transcripts could be used as a proxy for enzymes due to their highly inter-correlated dynamics. This finding is of practical importance considering the relatively higher cost and technical difficulties in directly measuring enzyme concentrations from environmental samples. While we were able to use quantitative enzyme data thanks to the advanced analytical method recently developed by Li et al. (2017), this method is costly and labor intensive in relative to transcript quantification.
Concluding Remarks and Future Considerations
Using the RSCM framework, we explored how the dynamic interplay between internal resources, transcripts, and enzymes can be considered as a potential mechanism responsible for the non-intuitive biogeochemical dynamics observed in our denitrification experiments. The good agreement we obtained between data and the model indicated that the proposed mechanism is plausible, while we could not provide a direct experimental proof. Further experimental studies would be required for a more conclusive understanding, particularly on the role of transcripts as a control point of regulatory dynamics. One could consider other biological or physical processes that may potentially cause our observed phenomena. For example, the slow process of transcription and translation could be ascribed to nutrient transport across the cell membrane mediated by the membrane-localized proteins (permeases) (Stephanopoulos et al., 1998). As another example, enzyme adsorption to minerals could be a potential process explaining their persistent activities in the sediments (Zimmerman and Ahn, 2010). Despite these caveats, the predictions of RSCM provide fresh insights and novel hypotheses that cannot be derived by simple structured biogeochemical models. Thus, the simulation results provide a reasonable basis for a deeper understanding of the role of microbes in regulating biogeochemical functions. In future applications, the RSCM will be further tested under various other environmental conditions where the importance of accounting for cellular regulation becomes more pronounced, including the transition from oxic to anoxic conditions. We expect the RSCM to serve not only as an independent simulator of denitrification dynamics, but also as a key biogeochemical modeling component of multi-scale reactive-transport models. As will be reported in a near future, the RSCM is currently incorporated into multi-scale simulation software, PFLOTRAN (Lichtner et al., 2015) to investigate hydro-biogeochemical processes and to predict the long-term impact of dam operation on thermo-hydro-biogeochemical dynamics.
All authors contributed to the design of the research. HS developed the modeling concept and performed simulations. DT, XS, and XC contributed to the parameter identification. HS and DT drafted out the article, which was edited by JS and JF. All authors read and approved the final manuscript.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This research was supported by the U.S. Department of Energy (DOE), Office of Biological and Environmental Research (BER), as part of Subsurface Biogeochemistry Research Program's Scientific Focus Area (SFA) at the Pacific Northwest National Laboratory (PNNL). The PNNL is operated for DOE by Battelle Memorial Institute under contract DE-AC05-76RL01830.
Boon, E., Meehan, C. J., Whidden, C., Wong, D. H., Langille, M. G., and Beiko, R. G. (2014). Interactions in the microbiome: communities of organisms and communities of genes. FEMS Microbiol. Rev. 38, 90–118. doi: 10.1111/1574-6976.12035
Bouskill, N. J., Tang, J., Riley, W. J., and Brodie, E. L. (2012). Trait-based representation of biological nitr fication: model development testing, and predicted community composition. Front. Microbiol. 3:364. doi: 10.3389/fmicb.2012.00364
Dhurjati, P., Ramkrishna, D., Flickinger, M. C., and Tsao, G. T. (1985). A cybernetic view of microbial-growth - modeling of cells as optimal strategists. Biotechnol. Bioeng. 27, 1–9. doi: 10.1002/bit.260270102
Franz, A., Song, H. S., Ramkrishna, D., and Kienle, A. (2011). Experimental and theoretical analysis of poly(beta-hydroxybutyrate) formation and consumption in Ralstonia eutropha. Biochem. Eng. J. 55, 49–58. doi: 10.1016/j.bej.2011.03.006
Gougoulias, C., Clark, J. M., and Shaw, L. J. (2014). The role of soil microbes in the global carbon cycle: tracking the below-ground microbial processing of plant-derived carbon for manipulating carbon dynamics in agricultural systems. J. Sci. Food Agric. 94, 2362–2371. doi: 10.1002/jsfa.6577
Kim, J. I., Song, H. S., Sunkara, S. R., Lali, A., and Ramkrishna, D. (2012). Exacting predictions by cybernetic model confirmed experimentally: steady state multiplicity in the chemostat. Biotechnol. Prog. 28, 1160–1166. doi: 10.1002/btpr.1583
Kim, J. I., Varner, J. D., and Ramkrishna, D. (2008). A hybrid model of anaerobic E. coli GJT001: combination of elementary flux modes and cybernetic variables. Biotechnol. Prog. 24, 993–1006. doi: 10.1002/btpr.73
Kompala, D. S., Ramkrishna, D., Jansen, N. B., and Tsao, G. T. (1986). Investigation of bacterial-growth on mixed substrates - experimental evaluation of cybernetic models. Biotechnol. Bioeng. 28, 1044–1055. doi: 10.1002/bit.260280715
Li, M., Gao, Y., Qian, W.-J., Shi, L., Liu, Y., Nelson, W. C., et al. (2017). Targeted quantification of functional enzyme dynamics in environmental samples for microbially mediated biogeochemical processes. Environ. Microbiol. Rep. doi: 10.1111/1758-2229.12558. [Epub ahead of print]
Lichtner, P., Karra, S., Hammond, G., Lu, C., Bisht, G., Kumar, J., et al. (2015). PFLOTRAN User Manual: A Massively Parallel Reactive Flow and Transport Model for Describing Surface and Subsurface Processes. Los Alamos National Laboratory (LANL).
Noor, E., Flamholz, A., Bar-Even, A., Davidi, D., Milo, R., and Liebermeister, W. (2016). The protein cost of metabolic fluxes: prediction from enzymatic rate laws and cost minimization. PLoS Comput. Biol. 12:e1005167. doi: 10.1371/journal.pcbi.1005167
Oliveira, A. P., and Sauer, U. (2012). The importance of post-translational modifications in regulating Saccharomyces cerevisiae metabolism. FEMS Yeast Res. 12, 104–117. doi: 10.1111/j.1567-1364.2011.00765.x
Shade, A., Peter, H., Allison, S. D., Baho, D. L., Berga, M., Burgmann, H., et al. (2012). Fundamentals of microbial community resistance and resilience. Front. Microbiol. 3:417. doi: 10.3389/fmicb.2012.00417
Song, H.-S., Morgan, J. A., and Ramkrishna, D. (2009). Systematic development of hybrid cybernetic models: application to recombinant yeast co-consuming glucose and xylose. Biotechnol. Bioeng. 103, 984–1002. doi: 10.1002/bit.22332
Song, H.-S., and Ramkrishna, D. (2011). Cybernetic models based on lumped elementary modes accurately predict strain-specific metabolic function. Biotechnol. Bioeng. 108, 127–140. doi: 10.1002/bit.22922
Song, H.-S., and Ramkrishna, D. (2013). Complex nonlinear behavior in metabolic processes: global bifurcation analysis of Escherichia coli growth on multiple substrates. Processes 1, 263–278. doi: 10.3390/pr1030263
Song, H.-S., Ramkrishna, D., Pinchuk, G. E., Beliaev, A. S., Konopka, A. E., and Fredrickson, J. K. (2013b). Dynamic modeling of aerobic growth of Shewanella oneidensis. Predicting triauxic growth, flux distributions, and energy requirement for growth. Metab. Eng. 15, 25–33. doi: 10.1016/j.ymben.2012.08.004
Song, H.-S., Renslow, R. S., Fredrickson, J. K., and Lindemann, S. R. (2015). Integrating ecological and engineering concepts of resilience in microbial communities. Front. Microbiol. 6:1298. doi: 10.3389/fmicb.2015.01298
Stegen, J. C., Fredrickson, J. K., Wilkins, M. J., Konopka, A. E., Nelson, W. C., Arntzen, E. V., et al. (2016). Groundwater-surface water mixing shifts ecological assembly processes and stimulates organic carbon turnover. Nat. Commun. 7:11237. doi: 10.1038/ncomms11237
Stevens, R. J., Laughlin, R. J., and Malone, J. P. (1998). Soil pH affects the processes reducing nitrate to nitrous oxide and di-nitrogen. Soil Biol. Biochem. 30, 1119–1126. doi: 10.1016/S0038-0717(97)00227-7
Taffs, R., Aston, J. E., Brileya, K., Jay, Z., Klatt, C. G., McGlynn, S., et al. (2009). In silico approaches to study mass and energy flows in microbial consortia: a syntrophic case study. BMC Syst. Biol. 3:114. doi: 10.1186/1752-0509-3-114
Tang, J. Y., and Riley, W. J. (2013). A total quasi-steady-state formulation of substrate uptake kinetics in complex networks and an example application to microbial litter decomposition. Biogeosciences 10, 8329–8351. doi: 10.5194/bg-10-8329-2013
Turner, B. G., Ramkrishna, D., and Jansen, N. B. (1988). Cybernetic modeling of bacterial cultures at low growth-rates - mixed-substrate systems. Biotechnol. Bioeng. 32, 46–54. doi: 10.1002/bit.260320108
Vasilakou, E., Machado, D., Theorell, A., Rocha, I., Noh, K., Oldiges, M., et al. (2016). Current state and challenges for dynamic metabolic modeling. Curr. Opin. Microbiol. 33, 97–104. doi: 10.1016/j.mib.2016.07.008
Wessely, F., Bartl, M., Guthke, R., Li, P., Schuster, S., and Kaleta, C. (2011). Optimal regulatory strategies for metabolic pathways in Escherichia coli depending on protein costs. Mol. Syst. Biol. 7:515. doi: 10.1038/msb.2011.46
Yan, S., Liu, Y., Liu, C., Shi, L., Shang, J., Shan, H., et al. (2016). Nitrate bioreduction in redox-variable low permeability sediments. Sci. Total Environ. 539, 185–195. doi: 10.1016/j.scitotenv.2015.08.122
Young, J. D., Henne, K. L., Morgan, J. A., Konopka, A. E., and Ramkrishna, D. (2008). Integrating cybernetic modeling with pathway analysis provides a dynamic, systems-level description of metabolic control. Biotechnol. Bioeng. 100, 542–559. doi: 10.1002/bit.21780
Keywords: nitrogen cycle, denitrification, hyporheic zone, environmental microbial community, trait-based modeling, the cybernetic approach, functional enzyme-based modeling
Citation: Song H-S, Thomas DG, Stegen JC, Li M, Liu C, Song X, Chen X, Fredrickson JK, Zachara JM and Scheibe TD (2017) Regulation-Structured Dynamic Metabolic Model Provides a Potential Mechanism for Delayed Enzyme Response in Denitrification Process. Front. Microbiol. 8:1866. doi: 10.3389/fmicb.2017.01866
Received: 26 April 2017; Accepted: 12 September 2017;
Published: 29 September 2017.
Edited by:Kelly Wrighton, The Ohio State University Columbus, United States
Reviewed by:Nick Bouskill, Lawrence Berkeley National Laboratory, United States
Karrie A. Weber, University of Nebraska Lincoln, United States
Copyright © 2017 Song, Thomas, Stegen, Li, Liu, Song, Chen, Fredrickson, Zachara and Scheibe. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hyun-Seob Song, email@example.com
†Present Address: Minjing Li, School of Environmental Studies, China University of Geosciences, Wuhan, China
Chongxuan Liu, School of Environmental Science and Engineering, Southern University of Science and Technology, Shengzhen, China