Integration of enzymatic data in Bacillus subtilis genome-scale metabolic model improves phenotype predictions and enables in silico design of poly-γ-glutamic acid production strains

Abstract

Background: Genome-scale metabolic models (GEMs) allow predicting metabolic phenotypes from limited data on uptake and secretion fluxes by defining the space of all the feasible solutions and excluding physio-chemically and biologically unfeasible behaviors. The integration of additional biological information in genome-scale models, e.g., transcriptomic or proteomic profiles, has the potential to improve phenotype prediction accuracy. This is particularly important for metabolic engineering applications where more accurate model predictions can translate to more reliable model-based strain design.

Results: Here we present a GEM with Enzymatic Constraints using Kinetic and Omics data (GECKO) model of Bacillus subtilis, which uses publicly available proteomic data and enzyme kinetic parameters for central carbon (CC) metabolic reactions to constrain the flux solution space. This model allows more accurate prediction of the flux distribution and growth rate of wild-type and single-gene/operon deletion strains compared to a standard genome-scale metabolic model. The flux prediction error decreased by 43% and 36% for wild-type and mutants respectively. The model additionally increased the number of correctly predicted essential genes in CC pathways by 2.5-fold and significantly decreased flux variability in more than 80% of the reactions with variable flux. Finally, the model was used to find new gene deletion targets to optimize the flux toward the biosynthesis of poly-γ-glutamic acid (γ-PGA) polymer in engineered B. subtilis. We implemented the single-reaction deletion targets identified by the model experimentally and showed that the new strains have a twofold higher γ-PGA concentration and production rate compared to the ancestral strain.

Conclusions: This work confirms that integration of enzyme constraints is a powerful tool to improve existing genome-scale models, and demonstrates the successful use of enzyme-constrained models in B. subtilis metabolic engineering. We expect that the new model can be used to guide future metabolic engineering efforts in the important industrial production host B. subtilis.

Keywords: Bacillus subtilis; Constraint-based methods; Enzymatic data; Genome-scale metabolic model; Poly-γ-glutamic acid.

Fig. 1

Experimental and predicted fluxes of central carbon reactions for wild-type B. subtilis. The predictions of growth rate, acetate secretion rate, fluxes through the reactions of glycolysis, TCA cycle and pentose phosphate pathway for wild-type strain, using iYO844 and ec_iYO844 (red and green bars, respectively), are compared to the experimental measurements (blue bars). Reaction names are reported according to the annotation used in the original model

Fig. 2

Distribution of prediction errors for B. subtilis mutant strains. Circles represent the prediction error distribution, computed with iYO844 and ec_iYO844, for each mutant strain considered in this work and the red lines indicate the median values

Fig. 3

Experimental and predicted fluxes for B. subtilis mutant strains. Five different mutants strains (∆pgi, ∆mdh, ∆zwf, ∆sdhABC and ∆serA, for which the PGI, MDH, G6PDH, SUCD1 and PGCDr reactions are blocked, respectively) were simulated using iYO844 and ec_iYO844 model (red and green bars, respectively). The predictions of growth rate, acetate secretion rate and fluxes through PGI, G6PDH, PYK, CS reactions are compared with the experimental measures (blue bars). Reaction names are reported according to the annotation used in the original model

Fig. 4

Distribution of predicted errors obtained from enzyme-constrained model with randomized k_cat values on wild-type strain. Bars represent the histograms obtained on 10,000 simulations, dashed line indicates the median value of the prediction errors obtained using randomized k_cat values and solid line represents the prediction error of the ec_iYO844 model. Panels show the results for the k_cat values randomization of all the 17 reactions (a), all the reactions except GAPD_NAD and CS (b) and only GAPD_NAD and CS (c)

Fig. 5

Flux variability analysis on the two models. a Cumulative distribution of all the non-zero FV values. b Identification of reactions and their respective metabolic subsystem with major variability reduction after the enzymatic data integration (red circle)

Fig. 6

Central carbon and γ-PGA production pathways of B. subtilis. The central carbon pathway that includes the pathway of glycolysis (blue), the pentose phosphate pathway (yellow) and the TCA cycle (red) is represented. The pathway of acetate (dark gray) and γ-PGA (green) biosynthesis are also represented. Dashed thin arrows indicate the bypass pathway from akg to succ, with sucsal as intermediate metabolite. The green arrow with dashed edge indicates a set of reactions (not shown) converting akg into glu-L. Metabolite and reaction names are reported in lower and upper case, respectively, according to the annotation used in the original model

Fig. 7

Reaction deletions predicted by MoMA to optimize the production of γ-PGA. Predicted γ-PGA flux and growth rate by iYO844 or ec_iYO844 for the single-reaction deletions and the best double-reaction deletions are reported. The slash indicates alternative double deletion solutions with the same predicted γ-PGA flux and growth rate. See text and Fig. 6 for a description of the reactions reported

Fig. 8

Fermentation experiments. a Growth curves in EM− and EM+ media. b γ-PGA concentration profiles in EM−. c Electrophoretic separation of γ-PGA collected in the experiments with EM−. d γ-PGA concentration profiles in EM−. The reference strain (PB5383) and the two mutants (PB5716 and PB5716) are indicated as A, B and C, respectively. In a, b, d circles represent average data points, error bars represent the 95% confidence intervals of the mean. Solid and dashed lines indicate experiments in EM− and EM+, respectively. In c the size distribution of the polymer can be qualitatively visualized. The intensity of the signal is comparable only among samples in which the polymer shows a similar distribution