iBsu1103: a new genome-scale metabolic model of Bacillus subtilis based on SEED annotations

A new and validated genome-scale metabolic model of Bacillus subtilis 168, iBsu1103, is presented that has significantly improved completeness and accuracy.


Background
Bacillus subtilis is a naturally competent, Gram-positive, sporulating bacterium often used in industry as a producer of high-quality enzymes and proteins [1]. As the most thor-oughly studied of Gram-positive and sporulating bacteria, B. subtilis serves as a model cell for understanding the Grampositive cell wall and the process of sporulation. With its similarity to the pathogens Bacillus anthracis and Staphylococ-cus aureus, B. subtilis is also important as a platform for exploring novel medical treatments for these pathogens. Moreover, the natural competence of B. subtilis opens the way for simple and rapid genetic modification by homologous recombination [2].
For all these reasons, B. subtilis has been the subject of extensive experimental study. Every gene essential for growth on rich media is known [3]; 60 gene intervals covering 49% of the genes in the genome have been knocked out and the resulting phenotypes analyzed [4]; 13 C experiments have been run to explore the cell response to mutations in the central carbon pathways [5]; and Biolog phenotyping experiments [6] have been performed to study the ability of B. subtilis to metabolize 271 different nutrient compounds [7].
As genome-scale experimental datasets begin to emerge for B. subtilis, genome-scale models of B. subtilis are required for the analysis and interpretation of these datasets. Genomescale metabolic models may be used to rapidly and accurately predict the cellular response to gene knockout [8,9], media conditions [10], and environmental changes [11]. Recently, genome-scale models of the metabolism and regulation of B. subtilis have been published by Oh et al. [7] and Goelzer et al. [12], respectively. However, both of these models have drawbacks and limitations. While the Goelzer et al. model provides regulatory constraints for B. subtilis on a large scale, the metabolic portion of this model is limited to the central metabolic pathways of B. subtilis. As a result, this model captures fewer of the metabolic genes in B. subtilis, thereby restricting the ability of the model to predict the outcome of large-scale genetic modifications. While the Oh et al. metabolic model covers a larger portion of the metabolic pathways and genes in B. subtilis, many of the annotations that this model is based upon are out of date. Additionally, both models lack thermodynamic data for the reactions included in the models. Without these data, the directionality and reversibility of the reactions reported in these models is based entirely on databases of biochemistry such as the Kyoto Encyclopedia of Genes and Genomes (KEGG) [13,14]. Hence, directionality is often over-constrained, with a large number of reactions listed as irreversible (59% of the reactions in the Goelzer et al. In this work, we introduce a new genome-scale model of B. subtilis based on the annotations generated by the SEED Project [15][16][17]. The SEED is an attractive source for genome annotations because it provides continuously updated annotations with a high level of accuracy, consistency, and completeness. The exceptional consistency and completeness of the SEED annotations are primarily a result of the subsystems-based strategy employed by the SEED, where each individual cellular subsystem (for example, glycolysis) is annotated and curated across many genomes simultaneously. This approach enables annotators to exploit comparative genomics approaches to rapidly and accurately propagate biological knowledge.
During the reconstruction process for the new model, we applied a group contribution method [18] to estimate the standard Gibbs free energy change of reaction (Δ r G'°) for each reaction included in the model. We then developed new extensions to an existing methodology [19][20][21] that uses these estimated Δ r G'° values along with the reaction stoichiometry to predict the reversibility and directionality of every reaction in the model. The Δ r G'° values reported for the reactions in the model may also be of use in applying numerous forms of thermodynamic analysis now emerging [22][23][24] to study the B. subtilis metabolism on a genome scale.
Once the reconstruction process was complete, we applied a significantly modified version of the GrowMatch algorithm developed by Kumar and Maranas [25] to fit our model to the available experimental data. In the GrowMatch methodology, an optimization problem is solved for each experimental condition that is incorrectly predicted by the original model, in order to identify the minimal number of reactions that must be added or removed from the model to correct the prediction. As a result, many equivalent solutions are generated for correcting each erroneous model prediction. We propose new solution reconciliation steps for the GrowMatch procedure to identify the optimal combination of GrowMatch solutions that results in an optimized model. We also propose significant alterations to the objective function of the GrowMatch optimization to improve the quality of the solutions generated by GrowMatch.

Reconstruction of the Core iBsu1103 model
We started the model reconstruction by obtaining the annotated B. subtilis 168 genome from the SEED. This annotated genome consists of 2,691 distinct functional roles associated with 3,257 (79%) of the 4,114 genes identified in the B. subtilis 168 chromosome. Of the functional roles included in the annotation, 50% are organized into SEED subsystems, each of which represents a single biological pathway such as histidine biosynthesis. The functional roles within subsystems are the focus of the cross-genome curation efforts performed by the SEED annotators, resulting in greater accuracy and consistency in the assignment of these functional roles to genes. Reactions were mapped to the functional roles in the B. subtilis 168 genome based on three criteria: match of the Enzyme Commission numbers associated with the reaction and the functional role; match of the metabolic activities associated with the reaction and the functional role; and match of the substrates and products associated with the reaction and functional role [26]. In total, 1,263 distinct reactions were associated with 1,032 functional roles and 1,104 genes. Of these reactions, 88% were assigned to functional roles included in the highly curated SEED subsystems, giving us a high level of confidence in the annotations that form the basis of the B. subtilis model.
Often genes produce protein products that function cooperatively as a multi-enzyme complex to perform a single reaction. To accurately capture the dependency of such reactions on all the genes encoding components of the multi-enzyme complex, we grouped these genes together before mapping them to the reaction. We identified 111 such gene groups and mapped them to 199 distinct reactions in the B. subtilis model. Reactions were mapped to these gene groups instead of individual genes if: the functional roles assigned to the genes indicated that they formed a complex; multiple consecutive non-homologous genes were assigned to the same functional role; or the reaction represented the lumped functions of multiple functional roles associated with multiple genes.
The metabolism of B. subtilis is known to involve some metabolic functions that are not associated with any genes in the B. subtilis genome. During the reconstruction of the B. subtilis model, 71 such reactions were identified. While 19 of these reactions take place spontaneously, the genes associated with the remaining reactions are unknown. These reactions were added to the model as open problem reactions, indicating that the genes associated with these reactions have yet to be identified (Table S3 in Additional data files 1 and 2).
Data from Biolog phenotyping arrays were also used in reconstructing the B. subtilis model. The ability of B. subtilis to metabolize 153 carbon sources, 53 nitrogen sources, 47 phosphate sources, and 18 sulfate sources was tested by using Biolog phenotyping arrays [7]. Of the tested nutrients, B. subtilis was observed to be capable of metabolizing 95 carbon, 42 nitrogen, 45 phosphate, and 2 sulfate sources. Transport reactions are associated with genes in the B. subtilis 168 genome for only 94 (51%) of these proven nutrients. Therefore, 73 open problem transport reactions were added to the model to allow for transport of the remaining Biolog nutrients that exist in our biochemistry database (Table S3 in Additional data files 1 and 2).
In total, the unoptimized SEED-based B. subtilis model consists of 1,405 reactions and 1,104 genes (Table 1). We call this model the Core iBsu1103, where the i stands for in silico, the Bsu stands for B. subtilis, and the 1,103 stands for the number of genes captured by the model (one gene is lost during the model optimization process described later). In keeping with the modeling practices first proposed by Reed et al. [27], protons are properly balanced in the model by representing all model compounds and reactions in their charge-balanced and mass-balanced form in aqueous solution at neutral pH [28].

Construction of a biomass objective function
In order to use the reconstructed iBsu1103 model to predict cellular response to media conditions and gene knockout, a biomass objective function (BOF) was constructed. This BOF was based primarily on the BOF developed for the Oh et al. genome-scale model of B. subtilis [7]. The 61 small molecules that make up the Oh et al. BOF can be divided into seven categories representing the fundamental building blocks of biomass: DNA, RNA, lipids, lipoteichoic acid, cell wall, protein, and cofactors and ions. In the Oh et al. BOF, all of these components are lumped together as reactants in a single biomass synthesis reaction, which is not associated with any genes involved in macromolecule biosynthesis. In the iBsu1103 model, we decomposed biomass production into seven synthesis reactions: DNA synthesis; RNA synthesis; protein synthesis; lipid content; lipoteichoic acid synthesis; cell wall synthesis; and biomass synthesis. These abstract species produced by these seven synthesis reactions are subsequently consumed as reactants along with 22 cofactors and ionic species in the biomass synthesis reaction. This process reduces the complexity of the biomass synthesis reaction and makes the reason for the inclusion of each species in the reaction more transparent. Additionally, this allows the macromolecule synthesis reactions to be mapped to macromolecule biosynthesis genes in B. subtilis. For example, genes responsible for encoding components of the ribosome and genes responsible for tRNA loading reactions were all assigned together as a complex associated with the protein synthesis reaction. Two new biomass precursor compounds were added to the biomass synthesis reaction of the iBsu1103 model to improve the accuracy of the gene essentiality predictions: coenzyme A (CoA) and acyl-carrier-protein (ACP). Both of these species are used extensively as carrier compounds in the metabolism of B. subtilis, making the continuous production of these compounds essential. The biosynthesis pathways for both compounds already existed in the iBsu1103, and two of the steps in these pathways are associated with essential genes in B. subtilis: ytaG (peg.2909) and acpS (peg.462). If these species are not included in the BOF, these pathways become nonfunctional, and the essential genes associated with these pathways are incorrectly predicted to be nonessential.
The coefficients in the Oh et al. BOF are derived from numerous analyses of the chemical content of B. subtilis biomass [29][30][31][32][33]. We similarly derived the coefficients for the iBsu1103 model from these sources. While no data were available on the percentage of B. subtilis biomass represented by our two additional biomass components CoA and ACP, we assume these components to be 0.5% of the net mass of cofactors and ions represented in the BOF.

Results of automated assignment of reaction reversibility
The group contribution method [18] was used to estimate standard Gibbs free energies of formation (Δ f G'°) for 948 (83.3%) of the metabolites and Δ r G'° for 1,372 (97.4%) of the reactions in the unoptimized iBsu1103 model. Estimated Δ r G'° values were used in combination with a set of heuristic rules (see Materials and methods) to predict the reversibility and directionality of each reaction in the model under physiological conditions ( Figure 1). Based on these reversibility rules, 635 (45%) of the reactions in the model were found to be irreversible. However, when the directionality of the irreversible reactions was set according to our reversibility criteria, the model no longer predicted growth on LB or glucoseminimal media. This result indicates that the direction of flux required for growth under these media conditions contradicted the predicted directionality for some of the irreversible reactions in the model. Six reactions were identified in the model that met these criteria ( Table 2). In every case, these reactions were irreversible in the reverse direction because the minimum Gibbs free energy change ( ) of each reaction was greater than zero. However, all of these reactions involve uncommon molecular substructures for which few experimental thermodynamic data are available [18]. Thus, in combination with the strong experimental evidence for the activity of these reactions in the direction shown in Table 2, we assumed that the Δ r G'° values of these reactions were overestimated by the group contribution method and that these reactions are, in fact, reversible.

Results of the model optimization procedure
The unoptimized model was validated against a dataset consisting of 1,500 distinct experimental conditions, including gene essentiality data [3], Biolog phenotyping data [7], and gene interval knockout data [4] (Table 3). Initially, 85 errors arose in the gene essentiality predictions, including 58 false positives (an essential gene being predicted to be nonessen- tial) and 27 false negatives (a nonessential gene being predicted to be essential). The annotations of all erroneously predicted essential and nonessential genes were manually reviewed to identify cases where the prediction error was a result of an incorrect gene annotation. Of the essential genes that were predicted to be nonessential, 30 were mapped to essential metabolic functions in the model. However, these essential genes all had homologs in the B. subtilis genome that were mapped to the same essential metabolic functions (Table S4 in Additional data files 1 and 2). Three explanations exist for the apparent inactivity of these gene homologs: they are similar to the essential genes but actually perform a different function; they are nonfunctional homologs; or the regulatory network in the cell deactivates these genes, making them incapable of taking over the functions of the essential genes when they are knocked out. In order to correct the essentiality predictions in the model, these 30 homologous genes were disassociated from the essential metabolic functions.
We then applied our modified GrowMatch model optimization procedure (see Materials and methods) in an attempt to fix the 116 remaining false negative predictions and 39 remaining false positive predictions ( Figure 2). First, the gap filling algorithm was applied to identify existing irreversible reactions that could be made reversible or new reactions that could be added to correct each false negative prediction. This step produced 686 solutions correcting 78 of the false negative predictions. The gap filling reconciliation algorithm was used to combine the gap filling solutions into a single solution that corrected 45 false negative predictions and introduced five new false positive predictions. Next, the gap generation algorithm was applied to identify reactions that could be removed or made irreversible to correct each false positive prediction. The gap generation algorithm produced 144 solutions correcting 32 of the false positive predictions. The gap generation reconciliation algorithm combined these solutions into a single solution that corrected 11 false positive predictions without introducing any new false negative predictions. Overall, two irreversible reactions were made reversible, 35 new reactions were added to the model, 21 reversible reactions were made irreversible, and 3 reactions were removed entirely from the model (Table S5 in Additional data files 1 Distribution of reactions conforming to reversibility rules Figure 1 Distribution of reactions conforming to reversibility rules. (a) The distribution of reactions in the iBsu1103 model conforming to every possible state in the proposed set of rules for assigning reaction directionality and reversibility is shown. This distribution indicates that most of the irreversible reactions in the model were determined to be irreversible because the Δ r G ' max value calculated for the reaction was negative. (b) The distribution of reactions in the iBsu1103 model involving the compounds used in the reversibility score calculation is also shown. These compounds are prevalent in the reactions of the iBsu1103 model, with 64% of the reactions in the model involving at least one of these compounds.  and 2). As a result of these changes, the model accuracy increased from 89.7% to 93.1%.

Model overview
The final optimized version of the iBsu1103 model consists of 1,437 reactions, 1,138 metabolites, and 1,103 genes (Table 1). Based on the reversibility rules and the estimated thermodynamic data, 653 (45.0%) of the model reactions were determined to be irreversible. All data relevant to the model are provided in the Additional data files, including metabolite structures (Additional data file 3), metabolite data (Table S1 in Additional data files 1 and 2), reaction data (Table S2 in Additional data files 1 and 2), estimated thermodynamic data (Table S2 in Additional data files 1 and 2), model stoichiometry in SBML format (Additional data file 4), and mappings of model compound and reaction IDs to IDs in the KEGG and other genome-scale models (Tables S1 and S2 in Additional data files 1 and 2).
The reactions included in the optimized model were categorized into ten regions of B. subtilis metabolism ( Figure 3a; Table S2 in Additional data files 1 and 2). The largest category of model reactions is 'fatty acid and lipid biosynthesis'. This is due to the explicit representation of the biosynthesis of every significant lipid species observed in B. subtilis biomass as opposed to the lumped reactions used in other models. The explicit representation of these pathways has numerous advantages: Δ f G'° and Δ r G'° may be estimated for every species and reaction; every species has a distinct structure, mass, and formula; and the stoichiometric coefficients in the reactions better reflect the actually biochemistry taking place. The other most significantly represented categories of model reactions are carbohydrate metabolism, amino acid biosynthesis and metabolism, and membrane transport. These categories are expected to be well represented because they represent pathways in the cell that deal with a highly diverse set of substrates: 20 amino acids, more than 95 metabolized carbon sources, and 244 transportable compounds.
Reactions in the model were also categorized according to their behavior during growth on Luria-Bertani (LB) media ( Figure 3b; Table S2 in Additional data files 1 and 2). Of the model reactions, 300 (21%) were essential for minimal growth on LB media. These are the reactions fulfilling essential metabolic functions for B. subtilis where no other pathways exist, and they form an always-active core of the B. subtilis metabolism. Another 697 (49%) of the model reactions were nonessential but capable of carrying flux during growth on LB media. While these reactions are not individually essential, growth is lost if all of these reactions are simultaneously knocked out. The reason is that some of these Model optimization procedure results  reactions represent competing pathways for performing an essential metabolic function. Another 229 (16%) of the reactions cannot carry flux during growth on LB media. These reactions are on the periphery of the B. subtilis metabolism involved in the transport and catabolism of metabolites not included in our in silico representation of LB media. Moreover, 210 (14%) of the model reactions are disconnected from the network, indicating that these reactions either lead up to or are exclusively derived from a dead end in the metabolic network. Presence of these reactions indicates miss-annotation or overly generic annotation of the gene associated with the reaction, or a gap in the metabolic network. Thus, these reactions represent areas of the metabolic chemistry where more experimental study and curation of annotations must occur.

Comparison with previously published models of B. subtilis
We performed a detailed comparison of the Oh et al. and iBsu1103 models to identify differences in content and elucidate the conflicts in the functional annotation of genes (Table  1). Our comparison encompassed the reactions involved in the models, the genes involved in the models, the mappings between genes and reactions in the models, and the gene complexes captured by the models (Figure 4). Our comparison revealed significant overlap in the content of the two models. Of the 1,020 total reactions in the Oh et al.  (Figure 3). This indicates a significantly more complete handling of complexes in the iBsu1103 model.
All of the additional content in the iBsu1103 model translates into a significant improvement in the accuracy of the gene knockout predictions, the Biolog media growth predictions, and the gene interval knockout predictions (Table 3). Even before optimization, the iBsu1103 model is 0.7% more accurate than the Oh et al. model. After optimization, the iBsu1103 model is 4.1% more accurate. In addition to the improvement in accuracy, the improved coverage of the genome by the iBsu1103 model also allows for the simulation of 337 additional experimental conditions by the model. In this work we also demonstrate new extended reversibility criteria for consistently and automatically assigning directionality to the biochemical reactions in genome-scale metabolic models. The extended criteria enabled us to identify 306 additional irreversible reactions that are missed when using existing methodologies alone [19][20][21]. However, we also found that even with the extended criteria, the predicted reversibility was not correct for every reaction in the model. In order for model predictions to fit available experimental observations, the predicted reversibility had to be adjusted for 29 (2%) of the model reactions. Some possible explanations for these exceptions to the reversibility criteria include: the estimated Δ r G'° may be too high or too low; the reactant or product con-centrations may be tightly regulated to levels that prohibit reactions from functioning in certain directions; or the reactions involve additional/alternative cofactors not accounted for in current reversibility calculations. These exceptions to the reversibility rules emphasize the importance of using a model correction method to adjust predicted reversibility based on experimental data. While these rules were very effective with the iBsu1103 model, they still need to be validated with a wider set of organisms and models. The extended version of GrowMatch presented in this work was also demonstrated to be a highly effective means of identifying and correcting potential errors in the metabolic network that cause errors in model predictions. This method is driven entirely by the available experimental data, requiring manual input only in selecting the best of the equivalent solutions generated by the solution reconciliation steps of the method. The reconciliation steps we introduced to the GrowMatch method also proved to be effective for identifying the minimal changes to the model required to produce the optimal fit to the available experimental data. The reconciliation reduced 830 distinct solutions involving hundreds of changes to the model to a single solution that combined 62 model modifications to fix 51 (33%) of the 155 incorrect model predictions.
Overall, we demonstrate the iBsu1103 model to be the most complete and accurate model of B. subtilis published to date. The identification and encoding of gene complexes, the removal of lumped reactions and compounds, and the refinements of the biomass objective function make this model especially applicable to thermodynamic analysis and gene knockout prediction. This model will be a valuable tool in the ongoing efforts to genetically engineer a minimal strain of B. subtilis for numerous engineering applications [2,4]. The thermodynamic data published with this model will be invaluable in the application of the model to numerous emerging forms of thermodynamic analysis [22][23][24]. Additionally, the new extensions that we have proposed for methods of automatically predicting reaction reversibility and automatically correcting model errors are valuable steps towards the goal of automating the genome-scale model reconstruction process [34,35].

Validation of the B. subtilis model using flux balance analysis
Flux balance analysis (FBA) was used to simulate all experimental conditions to validate the iBsu1103 model. FBA defines the limits on the metabolic capabilities of a model organism under steady-state flux conditions by constraining the net production rate of every metabolite in the system to zero [36][37][38][39]. This quasi-steady-state constraint on the metabolic fluxes is described mathematically in Equation 1: In Equation 1, N is the m × r matrix of the stoichiometric coefficients for the r reactions and m metabolites in the model, and v is the r × 1 vector of the steady-state fluxes through the r reactions in the model. Bounds are placed on the reaction fluxes depending on the reversibility of the reactions: -(CDW = cell dry weight). When simulating a gene knockout, the bounds on the flux through all reactions associated exclusively with the gene being knocked out (or associated exclusively with a protein complex partially encoded by the gene being knocked out) were reset to zero. When simulating media conditions, only nutrients present in the media were allowed to have a net uptake by the cell. All other transportable nutrients were allowed only to be excreted by the cell. Details on conditions for all FBA simulations performed are provided in Table S8 in Additional data files 1 and 2.

Prediction of reaction reversibility based on thermodynamics
The reversibility and directionality of the reactions in the iBsu1103 model were determined by using a combination of thermodynamic analysis and a set of heuristic rules based on knowledge of metabolism and biochemistry. In the thermodynamic analysis of the model reactions, Δ r G'° was estimated for each reaction in the model by using the group contribution method [40][41][42]. The estimated Δ r G'° values were then used to determine the minimum and maximum possible values for the absolute Gibbs free energy change of reaction (Δ r G ' ) using Equations 4 and 5, respectively: In these equations, x min is the minimal metabolite activity, assumed to be 0.01 mM; x max is the maximum metabolite activity, assumed to be 20 mM; R is the universal gas constant; T is the temperature; n i is the stoichiometric coefficient for species i in the reaction; U r is the uncertainty in the estimated Δ r G'°; and ΔG Transport is the energy involved in transport of ions across the cell membrane. Any reaction with a negative maximum Gibbs free energy change of reaction ( ) was assumed to be irreversible in the forward direction, and any reaction with a positive was assumed to be irreversible in the reverse direction. These criteria form the basis of many existing methods for predicting reaction reversibility [19][20][21].
However, in our work with the iBsu1103 model we found that and alone are insufficient to exhaustively identify every irreversible reaction in a model. Many reactions that are known to be irreversible have a negative and a positive due primarily to a lack of knowledge of true metabolite concentration ranges. To identify every irreversible reaction in the iBsu1103 model, we developed and applied a set of three heuristic rules based on common categories of biochemical reactions that are known to be irreversible: carboxylation reactions, phosphorylation reactions, CoA and ACP ligases, ABC transporters, and reactions utilizing ATP hydrolysis to drive an otherwise unfavorable action. We applied our new heuristic rules to identify any irreversible reactions that were missed by previous methods based only on and .
The first reversibility rule is that all ABC transporters are irreversible. As a result of the application of this rule, ATP synthase is the only transporter in the iBsu1103 model capable of producing ATP directly. The second reversibility rule is that any reaction with a milli-molar Gibbs free energy change (Δ r G 'm ) that is less than 2 kcal/mol and greater than -2 kcal/ mol is reversible. The Δ r G 'm is calculated by using Equation 6: Δ r G 'm is preferred over Δ r G'° when assessing reaction feasibility under physiological conditions because the 1-mM reference state of Δ r G 'm better reflects the intracellular metabolite concentration levels than does the 1-M reference state of Δ r G'°.
The final reversibility rule uses a reversibility score, S rev , calculated as follows: In this equation, n x is the number of molecules of type x involved in the reaction, Pi represents phosphate, Ppi represents pyrophosphate, and λ i is a binary parameter equal to 1 when i is a low-energy substrate and equal to zero otherwise.  to the final reversibility rule, if the product of S rev and Δ r G 'm is >2 and Δ r G 'm is <0, the reaction is irreversible in the forward direction; if the product of S rev and Δ r G 'm is >2 and Δ r G 'm is >0, the reaction is irreversible in the reverse direction. All remaining reactions that fail to meet any of the reversibility rule criteria are considered to be reversible.

Model optimization procedure overview
We applied an extended version of the GrowMatch procedure developed by Kumar et al. [25] to identify changes in the stoichiometry of the iBsu1103 model that would eliminate erroneous model predictions. The procedure consists of four steps applied consecutively (

Model optimization step one: gap filling
The gap filling step of the model optimization process, originally proposed by Kumar et al. [43], attempts to correct false negative predictions in the original model by either relaxing the reversibility constraints on existing reactions or by adding new reactions to the model. For each simulated experimental condition with a false negative prediction, the following optimization was performed on a superset of reactions consisting of every balanced reaction in the KEGG or in any one of ten published genome-scale models [7,12,20,27,[44][45][46][47][48][49]: Objective: Subject to: The objective of the gap filling procedure (Equation 8) is to minimize the number of reactions that are not in the original model but must be added in order for biomass to be produced under the simulated experimental conditions. Because the gap filling is run only for conditions with a false negative prediction by the original model, at least one reaction will always need to be added.
In the gap filling formulation, all reactions are treated as reversible, and every reversible reaction is decomposed into separate forward and reverse component reactions. This decomposition of reversible reactions allows for the independent addition of each direction of a reaction by the gap filling, which is necessary for gaps to be filled by the relaxation of the reversibility constraints on existing reactions.   If v max, i corresponds to a reaction associated with a knocked-out gene in the simulated experiment, this v max, i is set to zero. If v max, i corresponds to the uptake of a nutrient not found in the media conditions being simulated, this v max, i is also set to zero. Equation 11 constrains the flux through the biomass reaction in the model, v bio , to a nonzero value, which is necessary to identify sets of component reactions that must be added to the model in order for growth to be predicted under the conditions being simulated.
Each solution produced by the gap filling optimization defines a list of irreversible reactions within the original model that should be made reversible and a set of reactions not in the original model that should be added in order to fix a single false negative prediction. Recursive mixed integer linear programming (MILP) [50] was applied to identify the multiple gap filling solutions that may exist to correct each false negative prediction. Each solution identified by recursive MILP was implemented in a test model and validated against the complete set of experimental conditions. All incorrect predictions by a test model associated with each gap filling solution were tabulated into an error matrix for use in the next step of the model optimization process: gap filling reconciliation.

Model optimization step two: gap filling reconciliation
The gap filling step in the model optimization algorithm produces multiple equally optimal solutions to correct each false negative prediction in the unoptimized model. While all of these solutions repair at least one false negative prediction, they often do so at the cost of introducing new false positive predictions. To identify the cross-section of gap filling solutions that results in an optimal fit to the available experimental data with minimal modifications to the original model, we apply the gap filling reconciliation step of the model optimization procedure. In this step, we perform the following integer optimization that maximizes the correction of false negative errors, minimizes the introduction of new false positive errors, and minimizes the net changes made to the model: Objective: Subject to: Once again, recursive MILP was applied to identify multiple equivalently optimal solutions to the gap filling reconciliation problem, and each solution was validated against the complete set of experimental data to ensure that the combination of multiple gap filling solutions did not give rise to additional false positive predictions. The solutions that resulted in the most accurate prediction of growth in all experimental conditions were manually curated to identify the most physiologically relevant solution. This solution was then implemented in the original model to produce the gap-filled model.

Model optimization step three: gap generation
The gap-filled model produced by the gap filling reconciliation step not only will retain all of the false positive predictions generated by the original model but also will generate a small number of new false positive predictions that arise as a result of additions and modifications made during the gap filling process. In the gap generation step of the model optimization procedure we attempt to correct these false positive predictions either by removing irreversible reactions or by converting reversible reactions into irreversible reactions. For each simulated experimental condition with a false positive prediction by the gap-filled model, the following optimization was performed: Objective: Subject to: The objective of the gap generation procedure (Equation 17) is to minimize the number of component reactions that must be removed from the model in order to eliminate biomass production under conditions where the organism is known not to produce biomass. Equations 20 and 21 define the dual constraints associated with each flux in the primal FBA formulation. In these con- v bio no growth N v gapfilled growth Equation 23 is the constraint that sets the original primal FBA objective (maximization of biomass production) equal to the dual FBA objective (minimization of flux slack). This constraint ensures that every set of v no-growth fluxes that satisfies the constraints in Equations 20 to 23 represents an optimal solution to the original FBA problem that maximizes biomass production. Therefore, if the biomass flux is set to zero, as is done in Equation 24, this is equivalent to stating that the binary use variables z i must be set in such a way that the maximum production of biomass when simulating the false positive experimental conditions must be zero.
With no additional constraints, the gap generation optimization would produce solutions recommending the knockout of component reactions that cause the loss of biomass production under every experimental condition instead of just the false positive conditions. Constraints are required to ensure that only solutions that eliminate biomass production under the false positive conditions while preserving biomass production in all other conditions will be feasible. These constraints are defined by Equations 25, 26, and 27, which represent the FBA constraints simulating an experimental condition where the organism being modeled is known to grow. When the false positive condition being simulated by the v max, no-growth, i values is the knockout of an essential gene or interval, the v max, growth, i values in Equation 26 simulate the same media conditions with no reactions knocked out. When the false positive condition being simulated is an unviable media, the v max, growth, i values simulate a viable media. Because the binary z i variables are shared by the 'no growth' and 'growth' FBA constraints, z i will be set to zero only for those reactions that are not essential or coessential under the 'growth' conditions but are essential or coessential under the 'no growth conditions'. To further reduce the probability that a gap generation solution will cause new false negative predictions, we identified the component reactions in the gapfilled model that were essential for the correct prediction of growth in at least three of the experimental conditions prior to running the gap generation optimization. The z i variables associated with these essential component reactions were fixed at one to prevent their removal in the gap generation optimization.
As done in previous steps, recursive MILP was used to identify up to ten equally optimal solutions that correct each false positive prediction error in the gap-filled model. Each solution was implemented and validated against the complete set of experimental data, and the accuracy of each solution was tabulated into a matrix for use in the final step of the model optimization procedure: gap generation reconciliation.

Model optimization step four: gap generation reconciliation
Like the gap filling step, the gap generation step of the model optimization process produces multiple equally optimal solutions to correct each false positive prediction in the gap-filled model, and many of these solutions introduce new false negative prediction errors. To identify the cross-section of gap generation solutions that results in the maximum correction of false positive predictions with the minimum addition of false negative predictions, we perform one final optimization step: gap generation reconciliation. The optimization problem solved in the gap generation reconciliation step is identical to the gap filling reconciliation optimization except that the constraints defined by Equations 14 and 15 are replaced by the constraints defined by Equations 29 and 30: Equation 29 is written for any experimental condition with a false positive prediction by the gap-filled model. This constraint states that at least one gap generation solution that corrects the false positive prediction must be implemented for the condition to be correctly predicted by the optimized model. Equation 30 is written for any experimental condition where the original model correctly predicts that growth will occur. This constraint states that implementation of any gap generation solution that causes a new false positive prediction will result in a new incorrect prediction by the optimized model. All of the variables and constants used in Equations 29 and 30 have the same meaning as in Equations 14 and 15.
Although the objective, remaining constraints, and remaining variables in the gap generation reconciliation are mathematically identical to the gap filling reconciliation, some variables take on a different physiological meaning. Because gap generation solutions involve the removal (not the addition) of reactions from the gap-filled model, the reaction use variable z i is now equal to 1 if a reaction is to be removed from the gapfilled model and equal to zero otherwise.
The gap generation reconciliation was solved repeatedly by using recursive MILP to identify multiple solutions to the gap generation reconciliation optimization, and each solution was implemented in a test model and validated against the complete set of experimental data. The solutions associated with the most accurate test models were manually examined to identify the most physiologically relevant solution. The selected solution was then implemented in the gap-filled model to produce the optimized iBsu1103 model.

Additional data files
The following additional data are available with the online version of this paper: an Excel file containing all supplementary data associated with the iBsu1103 model as Tables S1 to S11 (Additional data file 1); a zip archive containing tabdelimited text files for all of the supplementary tables included in Additional data file 1 (Additional data file 2); a zip archive containing data on the structure of every molecule in the model in molfile format (Additional data file 3); an SBML version of the model, which may be used with the published COBRA toolbox [18] to run FBA on the model (Additional data file 4).
Additional data file 1 Tables S1 to S11 Tables S1 and S2 contain all compound and reaction data associ-ated with the model, respectively; Table S3 lists all of the open problem reactions in the model; Table S4 lists all of the essential genes that have nonessential homologs in the B. subtilis genome; Table S5 lists all of the changes made to the model during the model optimization process; Table S6 lists the reactions in the Oh et al. model that are not in the iBsu1103 model; Table S7 shows simula-tion results for all 1,500 experimental conditions; Table S8 pro-vides the details on the media formulations used for each FBA simulation; and Tables S9, S10, and S11 show all data on the genes, functional roles, and subsystems in the B. subtilis SEED annota-tion. Click here for file Additional data file 2 Tab-delimited text files for all of the supplementary tables included in Additional data file 1 The text files may be copied into any spreadsheet program of choice to visualize the data for the iBsu1103 model. Click here for file Additional data file 3 Data on the structure of every molecule in the model in molfile for-mat These molfiles reflect the structure of the predominant ionic form of the compounds at neutral pH as predicted using the Marvin-Beans software [28]. These structures were used with the group contribution method [40][41][42] to estimate the Δ f G'° and Δ r G'° for the compounds and reactions in the model. Click here for file Additional data file 4 SBML version of the model This may be used with the published COBRA toolbox [18] to run FBA on the model. Click here for file