Capturing the response of Clostridium acetobutylicumto chemical stressors using a regulated genome-scale metabolic model

Model comparisons

Model testing

The model was extensively tested to ascertain that it is capable of replicating flux ranges and phenotypes that have been documented for the wild-type (WT) organism and its mutants. The model predicted flux ranges were compared with experimental flux values from ¹³C-metabolic flux analysis (¹³C-MFA) [25]. The ¹³C-MFA data revealed that four reactions (pyruvate carboxylase (PC), fumarate hydratase (FH), pyridoxal phosphate synthase (PLPS), and alanine-glyoxylate (AGT)), which were originally removed to eliminate thermodynamically infeasible cycles, carried flux in the organism, and therefore, they were reinserted in the model. The cycles were instead eliminated by removing three alternate reactions (malate synthase (MS), succinate dehydrogenase (SDH), and malate dehydrogenase (MDH)) and by modifying the directionality of two others: succinyl-coenzyme A synthetase (SCS) was made reversible, and aspartate ammonia-lyase (ASPA) was restricted to the production of fumarate from aspartate. Figure 1C shows one of these cycles. ASPA was initially removed to fix the cycle due to a lack of literature evidence (Figure 1D), however subsequently MFA results indicated that this reaction carried flux whereas the MDH did not. Figure 1E shows how the addition of ASPA (directionally restricted) and removal of MDH avoids the formation of thermodynamically infeasible cycles while agreeing with experimental data.

After these changes in the model, flux variability analysis (FVA) was performed on the core carbon metabolism reactions, and the flux ranges were compared to the values obtained by the ¹³C-MFA analysis [25]. These experiments were carried out with a chemically defined medium that results in slower growth and lower biomass yields. First, all fluxes were normalized for a glucose uptake of 10 mmol gDW⁻¹ h⁻¹. FVA was performed while constraining the growth rate to the WT value of 0.32 h⁻¹ for C. acetobutylicum grown in complex media [37]. The comparison showed that the flux ranges of only four reactions (catalyzed by enolase (ENO), hexokinase (HK), pyruvate kinase (PYK), phosphotransacetylase (PTA), and phosphofructokinase (PFK)) encompassed the reported experimental values, as shown by Figure 2A. The reason for this is that C. acetobutylicum was grown in defined media by Au et al. [25], exhibiting significantly slower growth than in complex media [37]. In addition, the ¹³C-MFA data [25] was collected during the late growth phase with small amounts of solvents being produced, resulting in a reduced growth rate. Matching the FVA results with MFA data, we identified a growth rate value of 0.07 h⁻¹. Upon reapplying FVA with the biomass yield constrained to 0.07 h⁻¹ (see Figure 2B), all reactions except for HK and PC had flux ranges that encompassed experimental values. The two reaction experimental flux values differed from the model predicted range by only a value of 0.02 mmol gDW⁻¹ h⁻¹. It can be observed that, under these slow growing conditions, the TCA cycle reactions carry less flux and lie near the lowest end of the predicted flux range in Figure 2. The remaining flux is directed towards production of acids and solvents through pyruvate. This causes the flux of glycolytic reactions to lie near the high end of the predicted flux ranges (as shown by FVA predictions in Figure 2).

Following the model updates and comparisons with ¹³C-MFA data, the model’s responses to gene knockouts and varying environmental conditions were also tested. The model was used to analyze the effect of increasing the size of the NADH pool on the production of various acids and solvents. It has been shown experimentally that an increase in the level of NADH leads to a concomitant increase in butyrate, solvents, and hydrogen production (Figure 3) [38]. Allowing for the free conversion between NAD and NADH resulted in an increase in their production with the exception of acetate, whose production was, as expected, found to be independent of reducing equivalent availability.

Further model testing was performed by comparing experimental data of solvent yields for a number of C. acetobutylicum mutants [42] with in silico results. Mutants involving gene deletions affecting acid and solvent production in the ABE pathway were used to test iCac802. Biomass was constrained to the reported growth rate values for the respective experiments. Reaction fluxes associated with a deleted gene were set to zero. FVA was performed to determine the possible range of solvent production. FVA was first performed with biomass constrained to the reported growth rate to evaluate the flux ranges for the produced acids and solvents. The identified flux ranges of solventogenic nutrients (glucose, acetate, butyrate, and carbon dioxide) were subsequently calculated by fixing both the growth rate and restricting the acids/solvents to the FVA calculated maximum and minimum values. Yield ranges were determined by evaluating the ratio of acids/solvents to the corresponding minimum nutrient flux. Mutants Δack and Δptb reduce (but do not eliminate) acetate and butyrate production by removing acetate kinase (ACK) and phosphotransbutyrylase (PTB) activities, respectively (Figure 3) [37],[43]. For the two mutant strains, as well as for the WT strain, the model predicts a broad range of yields for the three solvents (butanol, acetone, and ethanol), as shown in the three-dimensional phenotypic solution space (Figure 4). This increased solvent production is also observed in experimental work by Jang et al. along with a reduction in acetate and butyrate production [44]. The study by Jang et al. also demonstrates that the butanol molar yield per glucose mole fed increases by 55% [44]. iCac802 predicts that incorporation of these two knockouts results in an increase in butanol production by 63.6%. An earlier GSM model by Lee et al.[20] predicts an increase in butanol production but by a larger value of 86% due to a lack of internal protons in the model. In the case of mutant strains Δadc and Δhbd, acetone and butanol production is impaired [36],[45] by knocking out acetoacetate decarboxylase (ADC) and hydroxybutyryl-CoA dehydrogenase (HBD), respectively (see Figure 5). In all cases the experimental yield is within the model-based calculated allowable yield for mutant phenotypes.

Modeling metabolic stressors using the CoreReg

iCac802 is a metabolic model and does not include any regulatory information. This section describes model regulation under various conditions by using transcriptomic data. Regulation was implemented in order to better describe the metabolism of C. acetobutylicum under butyrate and butanol stress and to pinpoint the reactions where flux changes are needed to explain the observed impact on biomass formation (that is, growth inhibition). Regulatory constraints on the iCac802 model were imposed using the transcriptomic data from Wang et al. [3] in the form of modified reaction flux bounds for each of the stress conditions using the CoreReg method (see Methods section for full description).

When regulation was imposed on the model, the biomass yield decreased by approximately 70% for all stress conditions except for the low-level butyrate stress, where the biomass yield was unaffected. For each one of the stress conditions we identify the reactions (core set) for which the application of regulatory constraints is sufficient to reproduce the flux ranges seen in the model when all regulatory constraints are imposed. Core sets of reactions were identified for each of the stress conditions by comparing flux bounds of the regulated model with the imposed regulatory constraints (Step 4 in Figure 6). The core sets represent likely nexus points of regulation under stress conditions, as they can broadly propagate the regulatory effect to the stress affected pathways through model stoichiometry. When regulatory bounds were imposed exclusively on the core set of reactions, the flux ranges of all reactions matched those of the model with all regulatory constraints. Subsequent core sets were found for the various stress levels by excluding the regulatory constraints on previously identified core sets (primary, secondary, tertiary sets, and so on). These subsequent core sets consisted of reactions whose regulatory constraints restrict the fluxes from wild-type distribution, and represent additional reactions that may be focal points of regulation. All these core sets are listed in Tables 5 and 6. In most cases, the same core sets of reactions were shared among the different levels of butanol stress. Three of the four reactions that made up these core sets (ornithine carbamoyltransferase (OCBT), arginosuccinate lyase (ARGSL), and arginosuccinate synthase (ARGSS)) belonged to arginine metabolism. OCBT was present in the primary core set of all levels of butanol stress. The final reaction, N-acetyl-gamma-glutamyl-phosphate reductase (AGPR), which is associated with amino acid metabolism, was present in the primary high-level butanol stress core set. The arginine metabolism genes identified in the core set for butanol stress are regulated by ArgR, the arginine repressor [3]. Expression of the genes corresponding to these identified arginine metabolic reactions was strongly down-regulated under butanol stress [3],[5] with a corresponding effect on biomass formation (growth inhibition). Identification of reactions in the arginine metabolism using the regulated model and its corroborative evidence from transcriptomic studies confirms the key role of arginine metabolism in response to stress and its subsequent effect of growth and metabolism. Furthermore, apart from the arginine metabolism, these genes are also involved in the biosynthesis of proline and lysine, which further emphasizes their role in regulating the amino acid metabolism and hence growth and biomass formation.

The primary medium level butyrate stress core set contained a reaction from pyrimidine metabolism, sulfate adenylyltransferase (SAT). The presence of this reaction can be related to the regulation of the DNA replication and repair mechanism which is initiated to protect the DNA from any damage owing to the oxidative stress component of the butyrate stress [3]. However, the subsequent core sets contained reactions involved in arginine metabolism, such as ARGSS in the high-level butyrate stress core set. Under butyrate stress, the effect that the regulatory constraints had on biomass yield was small. In comparison to butanol stress, under butyrate stress (low and medium), there is a strong up-regulation of genes in the arginine metabolism [3]. The addition of butyrate has a direct effect on the induction of solventogenesis, as the formation of solvents is directly related to reassimilation of butyrate from the medium (Figure 3). Jones et al. [46] have reported induction of the genes in arginine metabolism during the onset of solventogenesis, and this suggests that up-regulation of arginine genes under low and medium butyrate stress is associated with the induction of solventogenesis. The observation of arginine metabolism in the core set of high butyrate stress can be linked to the role of arginine metabolism as the acid response (AR3) system [3]. This ability of the model and regulatory modeling framework CoreReg to explicitly delineate the effect of two different metabolite stresses (at various levels) demonstrates the robustness and discriminatory capabilities of the model.

The addition of butyrate in the media leads to earlier onset of solventogenesis with higher butanol production [47],[48], which is due to the corresponding up-regulation of the genes involved in solvent production and notably those of the sol operon (CAP0162-CAP0164, adhe2-ctfA-ctfB) [3],[5]. The CoreReg method was able to simulate increased flux ranges for these reactions involved in solvent production during butyrate stress (Additional file 4).

Model construction

The general principles of the metabolic model reconstruction process have been previously described [49]-[51]. Construction of iCac802 entailed the following processes: 1) identification of biotransformations using previous models and homology searches; 2) assembly of reaction sets into a genome-scale metabolic model and subsequent conversion into a computations-ready format; 3) identification and removal of thermodynamically infeasible cycles; and 4) evaluation and improvement of model performance when compared to in vivo growth phenotypes.

We began metabolic model reconstruction by parsing the existing genome-scale models for C. acetobutylicum ATCC 824 [20]-[22] and the automated model generated on the Model SEED website [23]. We also made use of MetRxn, a knowledgebase that includes standardized metabolite and reaction descriptions drawn from multiple databases and genome-scale models [52], to parse and compare the contents of these models. We converted the gene associations in the Model SEED to the same open reading frame (ORF) naming scheme as the other models and available experimental data (CACxxxx or CAPxxxx) using the start and stop sites and in the genome annotation at the Model SEED website and those in the TIGR gene annotation [53]. We examined and updated the elemental and charge balancing of all reactions by making use of the chemical formulas and charges provided by the SEED database (molecular charge values calculated at neutral pH) [23]. All reactions were checked to verify that the reactants and products shared the same total numbers of different atoms and the same total charge. Reactions that did not pass these tests were replaced with equivalent balanced reactions from Model SEED (Additional file 3). In the case of no alternate Model SEED reaction, the reactions underwent rebalancing, where protons, hydroxide, and/or water molecules were added to balance atoms and/or charges. Metabolites that were not fully specified due to the use of generic R groups as side chains, or that were oligomeric with a non-specified number of repeat units, were flagged for manual examination and balancing. One-to-one equivalency of metabolites with such generic R groups and metabolites without generic R groups was verified to ensure balanced reactions. The biomass equation was adapted from the McAnulty model [22], the only modifications being the use of charged and uncharged tRNAs to account for amino acid utilization, the doubling of the coefficients for solutes in the equation, and a slight decrease in the coefficient for protein use from 3.1 to 3.

Model curation began with the removal of thermodynamically infeasible cycles present in the model. Flux variability analysis (FVA) was performed with no biomass constraints to identify unbounded reaction fluxes. The number of unbounded reactions was reduced by restricting the directionality of certain reactions using Model SEED’s calculated values of the reaction’s Gibbs free energy. If the entire range, including error, was more than 4 kcal/mol removed from zero, the reaction was restricted to the direction specified by the free energy [54]. The method described by Schellenberger et al. [55] was then used to identify the core set of thermodynamically infeasible cycles which form the basis of all such possible cycles. This method relies on the observation that all possible thermodynamically infeasible cycles form the null space of the stoichiometric matrix [55]. During the identification of the basis of the null space, only the set of reactions whose fluxes hit the bounds during FVA were used to decrease processing time. The null space basis of this set was evaluated by LU decomposition of the stoichiometric matrix [55].

Of the cycles found, those involving only two reactions were examined first. These were composed of equivalent reactions, that is, reactions which both converted a given set of reactants to products or vice versa using the same co-factors, varying only in their directionality. One of the two identified equivalent reactions was subsequently removed based on a preponderance of literature or annotation information concerning its directionality. Cycles containing more than two reactions were formed due to the presence of pathways that recycled metabolites consumed or produced by an alternate route with a zero net metabolic cost. Reactions in these cycles were reviewed, and those with no literature or annotation reference were either removed or had their directionality restricted. Figure 1A is an example of such a cycle. The 2-acetolactate methylmutase (ACLM) reaction is more accurately represented by two separate reactions, (R)-2, 3-dihydroxy-3-methylbutanoate: NADP + oxidoreductase (R-DMBO) and 2, 3-dihydroxy-3-methylbutanoate: NADP + oxidoreductase (DMBO). The removal of this aggregate reaction eliminates the thermodynamically infeasible cycle, as shown in Figure 1B. Figure 1 (C-E) represents another cycle and its subsequent removal using experimental ¹³C-metabolic flux analysis (¹³C-MFA) data.

After reviewing all of the cycles in the null basis with up to five participating reactions, the null basis was reevaluated. After fixing the smaller cycles, the same review process was extended to larger cycles. This reevaluation process was iterated until no more cycles were detected using FVA. A total of 101 reactions were removed from the model, and 58 had their directionality restricted to stop the involvement of 404 reactions in thermodynamically infeasible cycles.

Model simulations and analysis

where S _ij is the stoichiometric coefficient for metabolite i in reaction j. v _j represents the flux of reaction j, while $v_j^{\mathit{min}}$ and $v_j^{\mathit{max}}$ denote the minimum and maximum flux bounds on reaction j. N and M represent the overall sets of metabolites and reactions, respectively.

No constraints were placed on the biomass in order to identify all thermodynamically infeasible cycles. This analysis was performed iteratively for all metabolites.

Model testing

The model predictions were tested against experimental ¹³C flux data [25], experimental fermentation data [3],[35],[38],[40],[41], and in vivo gene knockout data from the literature [37],[43],[45],[57].

While the acidogenic phase was simulated in the model solely through the inclusion of required nutrients, additional constraints were required for the solventogenic phase. The solventogenic phase of Clostridia is characterized by the uptake of acetate and butyrate along with the reduction of carbon flux towards amino acids [58]. Thus, the solventogenic phase was simulated with acetate and butyrate as additional nutrients, and constraints on the export of acids and excess amino acids.

Incorporation of regulation using CoreReg

Regulation was incorporated into the model by a stepwise procedure that modified the minimum and maximum flux bounds of reactions based on fold change values of corresponding gene expression values by a new method called CoreReg. Gene expression data for each stress condition (GSE48031 and GSE48039), obtained from Wang et al. [3], were used to calculate the fold change from unstressed conditions using significance analysis of microarrays (SAM) [3],[59]. Transcriptomic data were collected for three concentrations each of butyrate and butanol. We refer to these as low (30 mM butyrate, 30 mM butanol), medium (40 mM butyrate, 60 mM butanol), and high (50 mM butyrate, 90 mM butanol) stress conditions. The fold change for each reaction under each stress condition was calculated from gene expression fold changes under stressed conditions by using GPRs. In the case of multiple enzyme subunits, the minimum expression value for the genes associated with the subunits was considered for calculating the reaction fold change. In the case of isozymes the total transcript level, obtained by summation of all isozyme transcripts, was considered. The unregulated model reaction bounds without any biomass constraint represent the minimum and maximum possible bounds of each reaction. Thus, a further increase in these bounds does not affect any maximum yield calculations, as the bounds are not active. Thus, only down-regulated genes were considered for evaluating the regulated model. The procedure for implementing regulation can be divided into five steps as explained below (see Figure 6).

Step 1: FVA is performed on the unregulated (UR) model to obtain lower $v_j^{\mathrm{L},\mathrm{U}\mathrm{R}}$ and upper $v_j^{\mathrm{U},\mathrm{U}\mathrm{R}}$ reaction flux bounds. The value of k is set equal to one to indicate primary core set.

Note that the updated lower bound is non-zero only for reversible reactions, effectively lowering the maximum possible flux value in the reverse direction.

Step 3: FBA is performed on the unregulated and the regulated model to obtain maximum biomass yields $v_{\mathit{biom}}^{\max ,\mathrm{U}\mathrm{R}}$ and $v_{\mathit{biom}}^{\max ,\mathrm{R}}$, respectively. If $v_{\mathit{biom}}^{\max ,\mathrm{R}}$ varies from $v_{\mathit{biom}}^{\max ,\mathrm{U}\mathrm{R}}$ by less than 2%, then the process is terminated, because the effect of the remaining regulation in the model is too small to cause any significant changes in metabolism as exemplified by the max biomass yield. Therefore, no additional regulatory core sets are extracted.

Step 4: FVA is performed at max biomass $v_{\mathit{biom}}^{\max ,\mathrm{R}}$ on the regulated (R) model to obtain lower $v_j^{\mathrm{l},\mathrm{R}}$ and upper $v_j^{\mathrm{u},\mathrm{R}}$ reaction flux bounds. These bounds are next compared with the imposed regulatory constraints from Step 2.

Step 5: Reactions j whose flux bounds are equal to the regulatory constraints (that is, $v_j^{\mathrm{u},\mathrm{R}}\equiv$$v_j^{\mathrm{U},\mathrm{R}}$ or $v_j^{\mathrm{l},\mathrm{R}}\equiv$$v_j^{\mathrm{L},\mathrm{R}}$) are assembled into the core set of reactions (of order k). To identify secondary, tertiary, and higher order core sets, the fold change values (c _j ) for the previously determined k core sets are set to one, thus removing their regulatory role in the model. The process is repeated from Step 2 with the value of k increased by one.