We present the RAVEN (Reconstruction, Analysis and Visualization of Metabolic Networks) Toolbox: a software suite that allows for semi-automated reconstruction of genome-scale models. It makes use of published models and/or the KEGG database, coupled with extensive gap-filling and quality control features. The software suite also contains methods for visualizing simulation results and omics data, as well as a range of methods for performing simulations and analyzing the results. The software is a useful tool for system-wide data analysis in a metabolic context and for streamlined reconstruction of metabolic networks based on protein homology. The RAVEN Toolbox workflow was applied in order to reconstruct a genome-scale metabolic model for the important microbial cell factory Penicillium chrysogenum Wisconsin54-1255. The model was validated in a bibliomic study of in total 440 references, and it comprises 1471 unique biochemical reactions and 1006 ORFs. It was then used to study the roles of ATP and NADPH in the biosynthesis of penicillin, and to identify potential metabolic engineering targets for maximization of penicillin production.
Genome sequencing projects have in recent years contributed enormously to our understanding of the metabolic capabilities of cellular systems. Functional annotation of the gene products allow for reconstruction of genome-scale metabolic models (GEMs) that summarize these metabolic capabilities in a consistent and compact way [1], [2]. A number of mathematical tools, including sampling of available metabolic states [3], [4] and methods borrowed from computational geometry [5], have been developed to analyze the resulting networks and to gain insight into the complex interactions that give rise to the metabolic capabilities. GEMs have also been used extensively for simulation of metabolism, particularly for metabolic engineering purposes [6], [7]. Since these models connect metabolites, proteins, and genes they are particularly well suited for the integration of metabolomics, proteomics, and genomics which is, in a sense, the goal of systems biology [8], [9]. The foundation of a GEM is the functional annotation of the genes. The first GEMs were primarily for model organisms for which direct evidence exists in the literature for a large proportion of the genomically encoded functions [10], . However, as the number of genome sequencing projects increases there is a growing demand of GEMs for less well known organisms. These models must by necessity be built largely relying on protein homology to more well-characterized organisms [12]. This, together with the large amount of manual work that is involved in a strict bottom-up reconstruction, has sparked interest in more automated approaches to model reconstruction. There are now a number of tools available for automated annotation of genes [13], [14]. However, the annotated genes must be linked to metabolic reactions in a way so as to generate a functional metabolic model. This includes the addition of spontaneous reactions and non-carrier mediated transport across membranes as well as sub-cellular localization of enzymes. Most importantly, the model must also be constructed in a way so that all reactions are balanced and well-connected [15]. This tends to become a problem if the gene-reaction relationship is automatically inferred from databases, partly due to differences in metabolite naming, but mainly because of how complex carbohydrates and complex lipids are represented. Because of the aforementioned issues it makes sense to use previously reconstructed models as templates for new GEMs. Here we present the RAVEN Toolbox, which allows the user to input GEM (s) for one or more template organisms, their corresponding protein sequences, and the protein sequences of the target organism. A GEM for the target organism is then constructed based on orthology between the protein sequences of the target organism and the organisms of the template models. Metabolic functions not present in the template models can obviously not appear in the new model, and to account for these missing reactions the RAVEN Toolbox also includes a functionality that matches proteins to KEGG Orthology (KO) categories [16] by using Hidden Markov models to capture the representative amino acid pattern in each KO. The resulting metabolic network can be used for automatic or manual gap filling, or it can be used on its own as a draft network. Several approaches which also aim at generating GEMs from either a template model or from a general database have been published [17]–[20]. Table 1 summarizes the capabilities of the RAVEN Toolbox compared to some other published approaches when it comes to automatic reconstruction. However, the largest difference is maybe not in the approaches taken, but in that the RAVEN Toolbox is a complete software for all tasks involving reconstruction and simulation of GEMs. In this aspect the RAVEN Toolbox is more similar to the COBRA Toolbox, but with extensive reconstruction capabilities [21]. Even though the RAVEN Toolbox can be used for fully automated reconstruction, in a manner similar to Model SEED, the intended purpose is to make use of the extensive quality control and gap identification/gap filling features for increasing the quality of reconstructions, as well as for decreasing the time needed for reconstructing high-quality models. The RAVEN Toolbox was evaluated for its ability to reconstruct a GEM for the well studied yeast Saccharomyces cerevisiae, and the resulting GEM was compared with a manually reconstructed model. Thereafter the RAVEN toolbox was used for the reconstruction of a GEM of the industrially important mold Penicillium chrysogenum. The Penicillium genus encompasses species of great economical, medical, and environmental importance [22]. Members of the Penicillium genus serve important roles in the food industry, both as some of the main spoilers of fresh vegetables and as essential actors in the production of blue cheeses. Most importantly though, they are sources of major antibiotics, particularly penicillin and griseofulvin. The industrial production of β-lactam antibiotics, such as penicillins and cephalosporins, is one of the success stories of biotechnology. Today the β-lactams represent one of the largest biotechnological products in terms of value, with sales of about USD 15 billion [23]. The industrial P. chrysogenum strains have been subjected to 50 years of directed evolution to increase the yields and titers of penicillin, with great cost reduction and productivity gain, but the yields are still far from the theoretical maximum [24]. A GEM of P. chrysogenum could aid in identification of metabolic bottlenecks as well as in elucidating the underlying reason for the significantly better performance of industrial strains compared to low producing strains.
A software suite named the RAVEN Toolbox (Reconstruction, Analysis, and Visualization of Metabolic Networks) was developed. The toolbox is a complete environment for reconstruction, analysis, simulation, and visualization of GEMs and runs within MATLAB. The software imports and exports models in two formats: the widely used Systems Biology Markup Language (SBML) format [25] and a Microsoft Excel model representation. Both these formats allow for extensive annotation of model components, such as International Chemical Identifier strings (InChI) [26] for metabolites or database identifiers for reactions and genes. The native model format for the RAVEN Toolbox follows the format of the yeast consensus metabolic network [27], but models in the COBRA Toolbox format can also be imported [21]. The Microsoft Excel representation enables the user to set simulation parameters such as bounds and objective function coefficients directly in the spread sheets. This simplifies the modeling process for users not comfortable with working within a scripting environment, as well as providing a simpler, but less rigorous, model format compared to SBML. The software, together with a manual, a set of tutorials and a detailed description of the supported file formats is available through the BioMet Toolbox [28] (http: //www. sysbio. se/BioMet). Figure 1 summarizes the capabilities of the RAVEN Toolbox. The software has three main foci: 1) automatic reconstruction of GEMs based on protein homology, 2) network analysis, modeling and interpretation of simulation results, 3) visualization of GEMs using pre-drawn metabolic network maps. In order to assign metabolic functions to the genes present in the P. chrysogenum genome, sequence alignment analysis was performed. Three fungi from the Aspergillus genus (A. oryzae, A. niger and A. nidulans) were selected for sequence comparison based on being closely related fungi outside of the Penicillium genus and on having previously reconstructed GEMs (see Figure S1). Table 2 shows some genome characteristics of the Aspergilli in comparison with P. chrysogenum. Initially pairwise comparison was done by similarity searching of the protein sequences of P. chrysogenum against the protein sequences known to be involved in the metabolism of the three Aspergillus species as described in the Methods section. With a chosen threshold of the E-value, identity, and alignment length, a list of inferred metabolic functions was generated. The results are summarized in Table 2. Pairwise comparison shows that A. oryzae has the highest number of sequence homologues of proteins with metabolic functions with P. chrysogenum (915 sequences). This result suggests that metabolism of A. oryzae is probably closer related to P. chrysogenum than A. nidulans and A. niger which have less sequence homologues of 576 and 563, respectively. Upon completion of the similarity searching, the results suggest that 1143 genes in P. chrysogenum could be assigned as orthologous metabolic genes from the three Aspergillus species used for comparison. The large number of metabolic orthologues indicates that the existing GEMs for closely related species could be a sound foundation upon which to reconstruct the new model. Using the RAVEN Toolbox the metabolic network of P. chrysogenum metabolic network was reconstructed. The metabolic network comprises 1471 unique metabolic reactions in four sub-cellular compartments; extracellular, cytosolic, mitochondrial, and peroxisomal (Table 3). 1006 ORFs are associated to the reactions, 89 of which participate in one of 35 protein complexes. In parallel to the automatic reconstruction, an extensive literature study was performed. In total 440 cited articles provide experimental evidence for the majority of the reactions. All model components were extensively annotated to adhere to the MIRIAM standard for biological models [41]. The model was validated with respect to 76 important metabolic functions (the supplementary file simulations. xls is an input file to checkTasks, which was used to performed the validation). There are 30 reactions in the model which cannot carry flux if all uptakes are allowed, i. e. dead-end reactions (see Table S7). According to the naming conventions for metabolic networks the presented model is denoted as iAL1006 [42]. Table 3 shows the division of model elements between the four compartments. Figure 5 summarizes the literature support for the reactions in the model and shows a classification of the ORFs in the model based on the KEGG pathways. The full list of reactions, metabolites, and genes are supplied in Microsoft Excel format and SBML format in Dataset S1. Both these formats are compatible with the RAVEN Toolbox. To illustrate the metabolic network, and to aid in interpretation of gene expression data and simulation results, a map of the full model was drawn in CellDesigner and annotated so as to be compatible with the visualization functions in the RAVEN Toolbox. The CellDesigner file is available in Dataset S1. Even though the map is drawn for Penicillium metabolism it can be used as a template for generation of maps for other organisms as well. Lastly, the P. chrysogenum GEM has also been added to the model repository in the BioMet Toolbox, which allows for a variety of analyses and simulations to be carried out. To evaluate the similarity between the reconstructed network and the template networks, the networks were compared with respect to identical reactions and involved metabolites. Only A. oryzae iWV1314 and A. niger iMA871 were used in the comparison since only a small number of reaction were inferred from A. nidulans iHD666. Figure 6 illustrates the results. 534 reactions were unique to iAL1006. The large discrepancy between the models is primarily because of differences in how lipid metabolism is formulated and due to differences in localization. There are also differences in how reactions catalyzed by protein complexes are described, where one reaction for each subunit is formulated rather than lumping the reactions. The difference in metabolic capabilities between the models is therefore smaller than what is indicated by the Venn diagrams (Figure 6). The unique capabilities of the P. chrysogenum model are mainly in penicillin metabolism and transport (data not shown). In general, the reactions that were inferred from A. oryzae but not from A. niger are predominantly involved in co-factor synthesis and in sugar polymer metabolism. The reactions inferred from A. niger iMA871 but not from A. oryzae iWV1314 are mainly involved in lipid metabolism. The key statistics of the reconstructed P. chrysogenum network compared to those of other fungal networks is available in Table S8. Growth is described as production of biomass, which in turn is regarded as drain of the macromolecules and building blocks that constitute the cellular components. The demand of each component is estimated based on published data on the biomass composition. The main components and their content within the biomass are listed in Table 4 (see also Table S9 for a detailed description). The cost of biomass production does not only include synthesis of precursors and polymerization of macromolecules, but also factors such as maintaining turgor pressure, transport costs, protein turnover, and membrane leakage. These costs are summarized as an ATP requirement for non-growth associated maintenance, mATP, and for growth associated maintenance, KxATP, i. e. ATP costs not directly associated with biomass synthesis but associated with cell growth (can be maintenance of membrane potentials across an expanding cell). These parameters were determined by linear regression to glucose-limited chemostat experiments in the presence of phenoxyacetic acid (POA) [43]. The values were hereby estimated to be 4. 14 mmol ATP/g DW/h for mATP and 104 mmol ATP/g DW for KxATP. The growth-associated ATP cost is significantly higher than for the template organisms (64 mmol ATP/g DW/h in A. niger iMA871). This could possibly be an effect of the presence of phenoxyacetic acid, which is added to the fermentation medium under industrial penicillin producing conditions. It is believed that phenoxyacetic acid, being a lipophilic weak acid, acts as a proton un-coupler which would manifest itself as a high ATP maintenance cost [44]. The P/O ratio is fitted by assigning the number of cytosolic protons needed to synthesize one ATP by the F0F1-ATPase. This is a small simplification since the number of protons pumped across the mitochondrial membrane might also differ between organisms. This parameter was estimated to 3. 75 (3. 88 in A. niger iMA871). Figure S2 shows the agreement of model simulations with experimental fermentation data after parameter fitting. A genome-scale metabolic model is a powerful tool that can be used for exploring the metabolic capabilities of the cell, as well as being used as a scaffold for integrative data analysis. Here we present two case studies to illustrate the use of the reconstructed P. chrysogenum model. The first case is a study of penicillin yields and in particular the relative importance of ATP and NADPH provision during penicillin production. In the second study we show how the model can be used to integrate fermentation data with transcriptome data using a recently published sampling algorithm to aid in the interpretation of high-throughput data [4].
The RAVEN Toolbox, a software suite for semi-automated reconstruction and simulation of genome-scale metabolic models was developed. The RAVEN Toolbox is the first software that contains methods both for model reconstruction and for a wide variety of simulation approaches. A visualization feature for simulation results and a feature that allows the user to manipulate metabolic models and set simulation parameters via Microsoft Excel are provided in order to make the software easy to use. The RAVEN Toolbox was evaluated for its ability to reconstruct GEMs by generating a model for S. cerevisiae. The reconstructed model compares well with a manually reconstructed model. This demonstrates that the RAVEN Toolbox is useful for reconstruction of novel models, in particular eukaryotic models, due to its feature for automatic assignment of sub-cellular localization. We used the RAVEN Toolbox to reconstruct a GEM for P. chrysogenum by using three models for closely related fungal species. Extensive manual validation of the model was performed; both to validate the reconstruction method and to ensure a high-quality model. The resulting P. chrysogenum model consists of 1471 reactions, 1235 metabolites and 1006 genes. 440 cited articles provide experimental evidence for the majority of the reactions. Considerable efforts were spent on standardizing and annotating the template models in order to adhere to MIRIAM standards. The standardized template models and the reconstructed P. chrysogenum model are available through the BioMet Toolbox. This collection of fungal models, together with the complementary method of generating metabolic networks based on the KEGG database, constitutes an excellent platform for the reconstruction of metabolic networks for other eukaryotic organisms.
Three GEMs for other filamentous fungi, A. nidulans iHD666 [52], A. niger iMA871 [53], and A. oryzae iWV1314 [54], were used as template models for the reconstruction of a P. chrysogenum model. Efforts were taken in order to standardize the template models to facilitate the automatic reconstruction. This standardization primarily involved metabolite naming, but also how to represent more complex aspects of metabolism such as polymers and lipids. As part of the standardization effort a large majority of the metabolites were assigned database identifiers and chemical structure information. This annotation step allowed for verification that all reactions were elementally balanced, which in turn led to a number of inconsistencies in the template models being corrected. The revised models for the three Aspergillus species are available as up-dates in the BioMet Toolbox (http: //www. sysbio. se/BioMet) [28]. A draft GEM for P. chrysogenum was then constructed based on bidirectional best hits of BLASTp between the template model proteins and their orthologues in P. chrysogenum using the RAVEN Toolbox. Proteins were regarded as orthologues if E-value <1e-30, identity >40%, sequence coverage (>50%) and alignment length (>200 amino acids). The protein sequences of P. chrysogenum Wisconsin 54-1255 (annotation, version 1) were obtained from the EMBL database (http: //www. ebi. ac. uk/embl/). The protein sequences of A. nidulans FGSC A4 (annotation, version 4) were taken from the Broad Institute database (http: //www. broadinstitute. org/annotation/genome/aspergillus_group). The protein sequences of A. oryzae RIB40 (annotation, version 1) were taken from the DOGAN database (http: //www. bio. nite. go. jp/dogan/project/view/AO). The protein sequences of A. niger ATCC1015 (annotation, version 3) were taken from the JGI database (http: //genome. jgi-psf. org/Aspni5/Aspni5. home. html). The first draft model based on homology to template models contained gaps due to incorrect annotation in the template models and lacked reactions in parts of metabolism that were unique to P. chrysogenum (Figure 8). Therefore, another draft model was generated from KEGG using the RAVEN Toolbox. An E-value <1e-50 was used as cut off in the gene assignment. This model was used for filling gaps in the draft network and for suggesting metabolic pathways that were not included based on the template models. No reactions were included based solely on presence in the KEGG model, and gene assignments were only included after careful manual validation against the NCBI RefSeq database. Gaps in the draft metabolic network were identified using the gap finding capabilities of the RAVEN Toolbox (Figure 8). The initial network was rather disconnected and biomass production from glucose was not possible. Firstly, the software was used to identify which metabolites had to be connected in order to produce biomass. This resulted in the addition of some spontaneous reactions from the template models. The second step was to ensure that as many metabolites as possible could be produced. When non-connected metabolites were identified, the KEGG model was queried for candidate reactions which could connect that metabolite. A targeted literature search was then conducted to find evidence for the presence of such a reaction. In only few cases this resulted in the addition of transport reactions based solely on connectivity issues (<2%, see Figure 5). The final step was to use the RAVEN Toolbox to identify reactions that could not carry a flux during growth on any of the available carbon sources. There are, however, sets of reactions that are included in the model even though they cannot currently carry a flux. One such example is the synthesis and loading of tRNA, which is included to allow for possible future extension of the model to cover protein synthesis. The model has four compartments: extracellular space, cytosol, mitochondrion, and peroxisome. Reactions with unknown localization, or where the real localization is not represented by one of the compartments in the model, were assigned to the cytosol. Enzymes reported to be present in the cell wall were assigned to the extracellular space, and those present in the mitochondrial membrane were mainly assigned to the mitochondria. The peroxisome was included primarily because of its role in penicillin metabolism and β-oxidation of fatty acids. Transport reactions between compartments were inferred mainly from the template fungal models and backed up with literature evidence. However, there were situations where a transport reaction had to be included in order to have a functional network, even when no literature evidence could be found. For enzymes reported to have isoenzymes in several compartments, the ORF assignments to each compartment were based on localization predictions from CELLO [36] and pTarget [55]. Unless otherwise specified, simulations were carried out using unlimited uptake of oxygen, phosphate, sulfate, NH3, thiamin and pimelate. The carbon source was glucose. Excretion of all exchange metabolites was allowed. Biosynthesis of L-cysteine was only allowed through the transsulfuration pathway. Since the energy content of NADH and NADPH is similar it is possible that cycles convert one into the other. These cycles normally take the form of reversible reactions that can utilize either NADH or NADPH. Since such futile cycles make it difficult to study NADPH/NADH metabolism using GEMs, and since they are probably not active in the cell, they were identified and blocked in analogy to what has been done in previous models [10]. The way by which they were identified was minimizing/maximizing for the flux through an artificial reaction NADH+NADP (+) < = >NAD (+) +NADPH while not allowing for uptake of any carbon source. Reactions were then identified and deleted until this reaction could no longer carry a flux. The following reactions were deleted: D-mannitol: NAD+ 2-oxidoreductase (r0181); ethanol: NADP+ oxidoreductase (r0019); ethanol: NAD+ oxidoreductase (r0020); (S) -3-hydroxybutanoyl-CoA: NADP+ oxidoreductase (r0081). None of these reactions can be expected to be active during the studied conditions. Deletion of D-mannitol: NAD+ 2-oxidoreductase does not inactivate the NADPH regenerating mannitol cycle, which has a role in NADPH regeneration in many fungal species [56]. The penicillin yields were calculated by setting the glucose uptake rate to 1. 0 mmol/gDW/h and maximizing for penicillin production. Free ATP was simulated by including an artificial reaction in the form ADP+Pi = >ATP+H2O. A random sampling algorithm was applied in order to identify transcriptionally regulated metabolic bottlenecks [4]. Flux data and gene expression levels for aerobic, glucose-limited chemostat fermentation of DS17690 and Wis 54-1255 were used as input to the algorithm [57]. The exchange fluxes were fitted to the reported values using a quadratic fitting. 5000 sampling iterations were performed for each of the two strains. The expression data set of the study were retrieved from GEO database (GSE9825) as CEL format then normalized together with PLIER workflow (http: //media. affymetrix. com/support/technical/technotes/plier_technote. pdf). Two way ANOVA were employed to evaluate the differentially expressed genes with respect to the strain (DS17690 and Wis 54-1255) with multiple correction following [58]. The Reporter algorithm [47] was employed to integrate the transcriptome data with the reconstructed GEM to identify key metabolites in the network.
Genome-scale models (GEMs) are large stoichiometric models of cell metabolism, where the goal is to incorporate every metabolic transformation that an organism can perform. Such models have been extensively used for the study of bacterial metabolism, in particular for metabolic engineering purposes. More recently, the use of GEMs for eukaryotic organisms has become increasingly widespread. Since these models typically involve thousands of metabolic reactions, the reconstruction and validation of them can be a very complex task. We have developed a software suite, RAVEN Toolbox, which aims at automating parts of the reconstruction process in order to allow for faster reconstruction of high-quality GEMs. The software is particularly well suited for reconstruction of models for eukaryotic organisms, due to how it deals with sub-cellular localization of reactions. We used the software for reconstructing a model of the filamentous fungi Penicillium chrysogenum, the organism used in penicillin production and an important microbial cell factory. The resulting model was validated through an extensive literature survey and by comparison with published fermentation data. The model was used for the identification of transcriptionally regulated metabolic bottlenecks in order to increase the yield in penicillin fermentations. In this paper we present the RAVEN Toolbox and the GEM for P. chrysogenum.
systems biology
biology
computational biology
genetics and genomics
The RAVEN Toolbox and Its Use for Generating a Genome-scale Metabolic Model for Penicillium chrysogenum