Integrated network reconstruction, visualization and analysis using YANAsquare

Steady state analysis

To increase platform independence and avoid loss of computation time (e.g. by parsing of the Metatool output) we included a Java implementation of the well known Schuster algorithm [2] which computes the EMs through a step-wise satisfaction of the steady state condition for each metabolite. The original version of the algorithm has been improved according to Klamt and Gagneur 2004 [19] by representing EMs during calculation by bit patterns rather than by their fluxes. This is possible due to the existence of a direct mapping of the set of reactions of an EM to the fluxes of these reactions [19]. The most frequently called function during the computation, the test for elementarity of an intermediate EM, is then reduced to a mere bit operation which drastically improves the algorithm's runtime behavior. This implementation is also used by a software package destinated at the computation of chemical organizations in chemical reaction networks [20].

The current version allows computation of the complete set of EMs or only the convex basis, both using either the external Metatool or the newly implemented internal EMA routine. We tested the algorithm thoroughly for consistency with the Metatool results.

KEGG Browser (KGB)

Successful analysis of a metabolic network which is expected to give biologically meaningful results heavily depends upon its accurate reconstruction. Complete sets of all modeled metabolites and enzymes have to be set up, and system boundaries have to be defined carefully. Every enzyme has to be checked for its absence or presence in the actual organism by investigation of genome annotations, homology searches or literature data mining. This task is hindered additionally by the missing or only partially established standardization of compound and enzyme names. Public databases such as Brenda [21], Enzyme [22] or the KEGG encyclopedia [18] have eased this process considerably by providing algorithms and interfaces to find the correct subset of enzymes to work on, their associated pathway structures and check for their occurrences in the organism under study. They further aim to unify the different available naming conventions, but so far the actual network setup in the modeling software had to be done manually. To further accelerate and automate the initial setup of a metabolic network we implemented the KEGG Browser (KGB) module in YANAsquare. It is capable of connecting directly to the KEGG database to browse pathways, reaction lists, metabolites and pathway maps online. Organisms or reference pathways can be selected to only retrieve enzymes which have been annotated as present in the selected organism or to obtain the complete set of available reactions respectively. All associated metabolites are collected automatically, highly abundant metabolites such as H⁺, CO², H₂O, ADP or ATP considered to be well buffered can be filtered out to lower the network complexity in the EMA analysis.

One problem when dealing directly with the KEGG database entries is that metabolites and reaction names are often long and may contain special characters. For a later visualization of the network, for backwards-compatibility to the Metatool input files and most importantly to satisfy the regular expression facets of the SBML2 standard which serves as our main file format, identifiers for enzymes and metabolites have to be standardized. We implemented an automatic abbreviation routine capable of shortening chemical names in a standardized way using a dictionary of keyword abbreviations. Generated abbreviations are validated so that characters illegal in the context of the SBML2 or Metatool format definitions will be eliminated. As enzyme names from ExPASy are generally even shorter, an additional step for the routine is to adapt shorter alternative names from ExPASy, which can again be further shortened by the abbreviation approach if necessary.

When accessing KEGG through a low-bandwidth Internet connection the retrieval of large chunks of reactions can be time consuming. We therefore additionally provide an accelerated query method using a local reaction cache (instead of retrieving all records directly via SOAP). This semi-online mode works significantly faster but the reaction cache has to be updated on a regular basis to assure that the KGB makes use of the latest KEGG definitions and annotations, an operation which can easily be performed through the KEGG Browser's update manager.

Once the desired topology has been assembled, the chosen reactions and metabolites can be imported directly into YANAsquare for further network analysis.

We used version 1.3 of the Apache Axis library, a Java implementation of the Simple Object Access Protocol (SOAP) and Web Services Description Language (WSDL) on top of the HTTP protocol to connect to the database server. We developed a user-friendly graphical interface using the SWING framework to retrieve the information from KEGG (see Fig. 1). Pathways and reactions of interests are specified by the user, all reaction equations are retrieved and listed in a table. A smart editor allows users to edit the results, change metabolite and enzyme names and abbreviations and/or append novel reactions if these are known to the user but not yet reported in KEGG. The whole software can also be easily adopted to extract data from another database, if format requirements are adapted or met, e.g. a private in-house metabolite and enzyme database in the KEGG format.

Visualization

With increasing network complexity the results from a steady-state analysis suffer from the combinatorial explosion of elementary modes. A realistic system comprising of about 100 metabolites and enzymes can easily result in over 100, 000 elementary modes. Interpretation of these data is close to infeasible when the results from the analysis are retrieved in form of a tabular listing of all elementary modes. Typical questions asked here are how the flux distribution for a specific set of active pathways is distributed across the network, which parts of the network are affected by certain sets of elementary modes or how deactivation of a set of pathways changes the overall network flux.

To circumvent these difficulties we implemented a visualization routine in YANAsquare capable of displaying the metabolic network as a bipartite graph, where the two vertex classes correspond to reactions and metabolites, respectively. We added several field-tested layout algorithms (e.g. spring-embedded layout, radial tree, sugiyama) to automatically arrange the nodes and edges of the network in a user-friendly way. While this visualization greatly facilitates the setup of metabolic networks, it is especially useful after performing a steady-state analysis of the network: YANAsquare is capable of projecting the flux distribution resulting from the set of EM activities onto the network graph, displaying the net flux through a given enzyme next to its vertex and scaling the width of the graph edges accordingly (smaller or bigger blue arrows in Fig. 2). This gives the user the possibility to directly inspect the network graphically for flux bottlenecks, hub or choke points and to see the overall flux distribution across the network at a glance. Individual fluxes through one enzyme for specific modes can also be analyzed and visualized.

Robustness analysis

Experiments in pharmacology, infection biology or genetics aim to modify or eliminate enzymes to study network stability. To further add to YANA's network analysis capabilities for such questions, a robustness test algorithm was implemented, elaborating on earlier ideas [9]. For a stoichiometric m × n matrix S with m metabolites and n enzymes, a matrix of elementary modes E (n × e) of e elementary flux modes and a vector of flux mode activities v the vector of metabolite productions

p = SEv

and is here defined component-wise.

expressed in percent is finally referred to as the overall robustness score of the network against single enzyme deletions.

The average mode loss m is then given via

m = e - t

Rapid network setup

A first example of the murine phospholipid network involved in phagosomal signaling illustrates the setup of a complex metabolic network of interest using YANAsquare. By querying the KEGG pathways "Glycerolipid metabolism (map 00561)", "Glycerophospholipid metabolism (map 00564)" and "Inositol phosphate metabolism (00562)" the complete set of 56 reactions and 64 metabolites available in the mouse metabolism is set up. After sensible setting of the system boundaries this leads directly to 57 elementary modes [see Tab. 4–6, Additional File 1]. Next, typical network extensions and modifications can be analyzed and illustrated. For instance, a reduced network with just 35 metabolites and 24 reactions is more accurate to identify key phospholipid conversions in the phagosome [see Tab. 1–3, Additional File 1]. In particular this reduced network includes as external components a number of phospholipid components which in the meantime have been shown to stimulate or inhibit the phagosome [23].

In this well connected network different conditions can be tested and interesting modes be rapidly visualized applying YANAsquare (as in Fig. 2), e.g. those modes yielding a specific phospholipid or those yielding ATP. Thus, setting ATP as an internal metabolite results in 14 ATP producing actin nucleation modes (from a total of 128 in this network) involving different phospholipids [for network definition and results, see Additional File 1]. Certainly this is a simplistic picture of the much more complex processes in phagosomal signaling. However the prediction of ATP generating actin nucleating modes identified and motivated experiments for different phospholipids, which are in fact in line with the model predictions. In particular, the modes suggest that the phospholipids phosphatidyl-inositol-phosphate and -bisphosphate (PIP and PIP2) as well as sphingosine-1-phosphate (S1P) stimulate the network whereas diacylglycerol (DAG) and phosphatidyl-choline (PC) are not in any of the modes and predicted to be inhibitory. This is in line with experimental data (Kühnel et al., unpublished, [23]) regarding exactly these phospholipids which have been obtained from phagosome actin nucleation assays. In these assays [23], murine phagosomes nucleate GFP labeled actin after phospholipid stimulation and the polymerization is monitored by confocal light microscopy.

Applying the KEGG browser additional phospholipid conversions involving e.g. phosphatidylethanolamine and phosphatidylsphingosine can easily be added. Furthermore, the network can be easily extended to prostaglandines and their conversions [KEGG map No. 00600, 00590, for network definition and results, see Additional File 1].

Genome-scale pathway analysis in Staphylococci

The capabilities of YANAsquare further allow rapid and comparative genome-scale metabolic analyses. The tutorial and examples included in the supplementary material contain helpful hints for the user to achieve this, including network simplification rules to avoid combinatorial explosion of elementary modes. In our illustration example, different Staphylococci strains are compared for their main metabolic capabilities: S. aureus COL, S. aureus N315, S. saprophyticus and S. epidermidis. Note that the comparisons serve illustration purposes and more detailed studies are required for an accurate modeling of these strains. However the examples show how the tool YANAsquare allows a rapid first overview and points out metabolic differences.

The KEGG browser allows easy establishment of the complete network of interest. Regarding the growth equation of S. aureus [24], it is evident that we require the metabolism of central sugars, amino acids, and fatty acids as well as the nucleotide metabolism. This is easily compiled using the KEGG browser. Two different strains are readily accessed using KEGG: S. aureus COL and S. aureus N315. The central metabolism of S. aureus COL including carbohydrates, nucleotides, amino acids and lipid metabolism is rapidly established applying YANAsquare. Using the KEGG maps 00010, 00020, 00030, 00251, 00252, 00260, 00271, 00272, 00290, 00300, 00330, 00340, 00400, 00620, 00720, a system with 134 metabolites and 3957 modes is obtained [see Additional File 1, Tab. 10–15; for rules for abbreviation and network simplification see tutorial]. YANAsquare shows also quickly that regarding central metabolism, strain N315 is rather similar to S. aureus COL (same number of modes and metabolites). Clear differences in central metabolism become readily apparent changing to S. epidermidis and S. saprophyticus instead. Reading in the KEGG maps (the serial numbers are the same as in S. aureus) we now obtain a system with 134 metabolites and 11910 modes (S. epidermidis) as well as 132 metabolites and 11844 modes (S. saprophyticus), respectively [for example modes see Additional File 1; the complete list of modes is available from the authors on request].

A first comparison of the network models obtained predicts that S. aureus has lower ATP production than S. epidermidis, since its tricarboxylic acid cycle (TCA) is not complete: one involved enzyme – malate dehydrogenase (MDH; EC: 1.1.1.37) – is absent in S. aureus strains. It is a well-known NAD-dependent enzyme, which synthesizes oxaloacetate and NADH, the latter metabolite will contribute 3 moles of ATP through oxidative phosphorylation. A second enzyme, membrane-associated malate: quinone-oxidoreductase (EC: 1.1.99.16) can be expressed at certain stages. As this is then highly expressed in S. aureus, this partially compensates for the lack of MDH. However, the corresponding reaction requires external supply of Vitamin K (the S. aureus strain lacks the capability of producing it). The model for this can be readily established, e.g. by switching the menaquinol and menaquinone metabolites from internal to external in the software to simulate the synthesis of metabolites if Vitamin K is present in the medium. Comparing the resulting pathways between the different strains further, we find that the S. aureus strain can bypass the above TCA gap via malic enzyme (EC: 1.1.1.38), phosphoenolpyruvate-protein phosphotransferase (EC: 2.7.3.9) and phosphoenolpyruvate carboxylase (EC: 4.1.1.49) to produce oxaloacetate which can take part in the next round of the TCA cycle or be used to synthesize amino acids. However, this alternative way consumes ATP and requires the energy-consuming PTS system (EC: 2.7.1.69, EC: 2.7.3.9).

S. aureus has more functional enzymes and modes regarding amino acid production, whereas S. saprophyticus has fewer modes for certain amino acids. For instance, the absence of threonine aldolase (EC 4.1.2.5) disrupts the connection between threonine and glycine pathways. Even its urea cycle is not complete (e.g., S. saprophyticus modes no: 4 - 8 etc.), as for instance carbamoyl phosphokinase (EC: 2.7.2.2) is not present. S. saprophyticus lacks the capability of producing carbamoyl phosphate directly. Instead, this is supplied from amino acids (e.g., S. saprophyticus modes no: 4 - 8 etc.). All missing metabolic capabilities (e.g., compared to the modes no. 3, 19, 21, 81, 91, 100 in S. aureus) should potentially reduce the robustness, even though the system can still synthesize all the essential amino acids with an appropriate medium. S. epidermidis lacks arginase (EC: 3.5.3.1) and its urea cycle is incomplete. Arginine has to be converted into citrulline by arginine deiminase (EC: 3.5.3.6), forming a shunt (e.g. modes no. 22). An advantage shared by S. saprophyticus and S. epidermidis is that their native malate dehydrogenase (EC: 1.1.1.37) is present (e.g. modes no: 21, 22, 24 in S. saprophyticus and modes no. 16, 17, 18 in S. epidermidis), so that their citric acid cycle is more efficient and produces more ATP. These are only some strain-specific differences for illustration and exploiting the new software to show how the full scale metabolic models may be easily and quickly compared. In practice, after this initial analysis in most cases a more detailed study of individual sub-networks has to follow, often in combination with further genome re-annotation as well as experimental tests (e.g. [17]).

Robustness analysis

Furthermore, our software provides a systematic test routine to get an estimate of the robustness of a metabolic network. The principle is to validate the flux preservation when a certain involved enzyme is knocked out from the system. Our program simulates this evolutionary or genetic event by removing enzymes one at a time and re-inspecting all generated flux modes. These survival modes after a single enzyme gene knock-out are used to predict the number of metabolites which can be still produced in the restricted system. Consequently all the enzymes are examined and an average metabolite production is calculated. By dividing the number through the maximum number of external metabolites which can be produced by the complete system, we achieve a percentage, referred to as the robustness score, which is an important feature of the system (see implementation). Table 1 illustrates robustness differences among these three systems. From the results (the third row), S. saprophyticus achieves the highest score (92.84), it is definitely the most robust organism among the three Staphylococci compared. This is potentially caused by its flexible living environments, not only in vivo as a pathogen but also outside any host. In contrast, S. aureus has strictly adapted to the host environment so that it requires more complex medium in order to acquire a recognizable growth rate.