Volume 17 Supplement 4
From highthroughput structural bioinformatics to integrative systems biology: Selected Works from the 14th International Workshop on Network Tools and Applications in Biology (NETTAB 2014)
Multiplex methods provide effective integration of multiomic data in genomescale models
 Claudio Angione^{1}Email author,
 Max Conway^{2} and
 Pietro Lió^{2}
DOI: 10.1186/s1285901609121
© Angione et al. 2016
Published: 2 March 2016
Abstract
Background
Genomic, transcriptomic, and metabolic variations shape the complex adaptation landscape of bacteria to varying environmental conditions. Elucidating the genotypephenotype relation paves the way for the prediction of such effects, but methods for characterizing the relationship between multiple environmental factors are still lacking. Here, we tackle the problem of extracting networklevel information from collections of environmental conditions, by integrating the multiple omic levels at which the bacterial response is measured.
Results
To this end, we model a large compendium of growth conditions as a multiplex network consisting of transcriptomic and fluxomic layers, and we propose a multiomic network approach to infer similarity of growth conditions by integrating layers of the multiplex network. Each node of the network represents a single condition, while edges are similarities between conditions, as measured by phenotypic and transcriptomic properties on different layers of the network. We then fuse these layers into one network, therefore capturing a global network of conditions and the associated similarities across two omic levels. We apply this multiomic fusion to an updated genomescale reconstruction of Escherichia coli that includes underground metabolism and new geneproteinreaction associations.
Conclusions
Our method can be readily used to evaluate and crosscompare different collections of conditions among different species. Acquiring multiomic information on the topology of the space of experimental conditions makes it possible to infer the position and to build conditionspecific models of untested or incomplete profiles for which experimental data is not available. Our weighted network fusion method for genomescale models is freely available at https://github.com/maxconway/SNFtool.
Keywords
Genomescale models Multiplex networks Multiomics Network fusion Transcriptomics Fluxomics Flux balance analysisBackground
As the cost of collecting omic data is likely to decrease in the coming years, methods to integrate different types of biological data are likely to become increasingly important. Biological data is often expressed as network data. As a result, methods to tackle networks with different layers, each representing a different omic, are likely to play a key role to integrate these data types. A common approach in this case is a layerwide analysis of each layer separately, finally integrating the separate results into a single aggregate measure.
However, integrating results collected on single layers separately may cause loss of information, especially in cases where the interlayer interactions are nonnegligible. Intuitively, by analyzing different types of data in isolation we miss the information that is implicit in the coordinated activity of the different layers. For instance, in an application to networks of patients to investigate cancer, this approach resulted in different clusters of patients depending on which dataset was used, therefore causing incorrect subdivision of patients or samples in different molecular cancer subtypes [1].
Genomescale models of biological organisms are often used to analyze the metabolic potential and to identify the metabolic interventions required to produce metabolites of interest. In this postgenomic era, more than 90 genomescale metabolic reconstructions have been published to date [2], for both prokaryotes and eukaryotes. Computational methods on genomescale models make it possible to explore the reaction network and find solutions that take into account noncomparable objectives, while satisfying all the given constraints.
With the growing availability of multiomic databases, modelbuilding and data integration techniques, a number of methods for integration of omic data have been proposed recently (for a comprehensive review and current challenges, the reader is referred to Saha et al. [3] and Bordbar et al. [4]). The wide variation in response of a single organism across differing conditions is facilitated by complex metabolic networks, which often include multiple levels of redundancy and adaptability. This is evident in the heterogeneity at the transcriptomic, proteomic and fluxomic levels across conditions [5]. In this regard, evaluating and aggregating in a single network the response to growth condition on different omic layers, while accounting for the dependence between the metabolic response at different levels (e.g. transcriptomic and fluxomic), is a highly desired feature, as it increases the reliability of any comparative analysis performed on the conditions, and provides overall topological information.
Our idea is to integrate multiple data types, collected from a comprehensive set of environmental conditions in which Escherichia coli was grown. These conditions are modeled as a multiplex network. Specifically, we focus on the transcriptomic layer (microarray data) and on the fluxomic layer (reaction rates at steady state, which we consider as the minimum proxy for the phenotype). The idea is that analyzing growth conditions by focusing only on a single omic level will likely miss complementary information and therefore lead to incorrect evaluations. For instance, considering only gene expression profiles as a response to conditions will miss the actual metabolic response of the bacterium, and will not provide insights into the resulting cellular behavior.
Two main methods have been recently proposed to address the task of integrating multiple data types in an aggregate layer of a multiplex network. iCluster, developed by Shen et al. [6], proposes a joint latent variable model for integrative clustering of multiple genomic datasets. More specifically, tumor subtypes are modeled as unobserved latent variables, which are then estimated from the multiple data types available. Since iCluster does not scale to the entire set of available measurements and therefore needs gene preselection steps, an alternative approach was developed by Wang et al. [7], which is instead based on networks of samples as a starting point for the creation of an aggregate layer. Both methods have been only proposed on genes and, since a genotypephenotype link is missing, they are not able to target and classify growth conditions.
Unlike the existing methods for network aggregation, we propose a weighted network fusion that takes into account the importance of each layer. This allows us to apply our method to multiomic genomescale models where nodes represent environmental conditions and layers represent transcriptomic and phenotypic information. The weight measures the relative importance of each layer, which is quantified by the reliability of the genomescale model in predicting correct flux rates starting from gene expression data. As a result, our multiplex network approach is able to integrate and elucidate condition similarity in multiomic genomescale models, and to be calibrated to the quality of different data types available for a given organism. To test our method, we build a genomescale model of Escherichia coli, starting from the iJO1366 reconstruction [8] and including recently discovered geneprotein associations and a comprehensive set of underground metabolic reactions.
Methods
An improved model of Escherichia coli with underground metabolism and new geneproteinreaction associations
To obtain the model used in this study, we started from the iJO1366 E. coli reconstruction [8]. We then added the full set of underground reactions reported by Notebaart et al. [9]. These underground reactions occur at low rates compared to the rest of the metabolic network, and new metabolites are usually detected with low abundances. Interestingly, the effect of the underground reactions is promptly detected in the newly built genomescale model, since many underground reactions share metabolites with existing reactions or constitute new pathways leading to precursors for biomass production.
New geneproteinreaction (GPR) associations added in the E. coli model augmented with underground metabolic reactions, according to the experimental studies in [10]
Reaction name  GPR in iJO1366  GPR in this study 

Aspartate transaminase  aspC  aspC OR tyrB 
Tyrosine transaminase  aspC OR tyrB  aspC OR tyrB OR ilvE 
Acetylornithine transaminase  argD OR astC  argD OR astC OR gabT OR puuE 
Succinyldiaminopimelate transaminase  argD  argD OR astC OR gabT OR puuE 
Citrate synthase  gltA  gltA OR prpC 
We remark that adding the underground metabolic reactions and the new GPR associations did not generate detrimental intermediates or reduced flux in those pathways supporting generation of biomass. The inclusion of underground metabolic reactions also increased flexibility and fitness under various growth environments. The resulting model, provided as Additional file 1, contains 1380 genes, 3027 reactions (including exchange reactions), and 2151 metabolites.
Steadystate analysis of genomescale models
Various techniques for modelling and simulating metabolic networks have been developed in recent years to better understand bacterial metabolism, as well as to support rational design for reprogramming microorganisms and overproducing biochemical compounds [4]. Arguably, the most successful technique for predicting flux distributions at steady state is fluxbalance analysis (FBA) [12]. FBA is based on two main assumptions: (i) homeostatic assumption, i.e. the organism has reached a steady state where the metabolite concentrations are constant and a set of nutrients are being constantly converted to generate biomass; (ii) (multilevel) optimality, i.e. in each state the organism tends to maximize one or multiple objectives, usually related to growth, excretion of biotechnologicallyrelevant metabolites and important energycarrying molecules. FBA is suitable for analyzing the flow of metabolites through a metabolic network (e.g. their formation and degradation, transport and cellular utilization).
Let the metabolic network be composed of m metabolites with concentration x _{ i }, i=1,…,m and n reactions with flux rates v _{ j }, j=1,…,n. The derivative of the concentration of each metabolite can be computed through a linear combination of the input and output reaction fluxes, which produce and consume (respectively) that metabolite. The balance that metabolite concentrations x _{ i } must satisfy is \(\frac {dx_{i}}{dt} = \sum _{j=1}^{n} A_{\textit {ij}}v_{j}, \quad i=1,\ldots,m\), where A _{ ij } is the stoichiometric coefficient of the ith metabolite in the jth reaction. Under steady state conditions \(\frac {dx_{i}}{dt} = 0, \ \forall i\), the balance for the ith metabolite is \(\sum _{j=1}^{n} A_{\textit {ij}}v_{j} = 0\) (homeostatic assumption). Therefore, at steady state, the balance equation is A v=0, where A is the stoichiometric matrix (m rows and n columns), and v is the vector of the flux rates (metabolic and transport fluxes).
Quantitative gene expression levels in genomescale models
Recently, FBA has been integrated with regulatory constraints taken from gene expression data. For a comprehensive evaluation of these methods, the reader is referred to the paper by Machado and Herrgård [13]. The most widely used approach for linking gene expression level and FBA models consists of removing those reactions in the model that are linked to genes whose expression is below a specific threshold. This approach allows a fast implementation of the effect of gene expression variations on the model, but does not provide the necessary level of information needed for a comprehensive and predictive steadystate biological model. Furthermore, it forces the gene expression levels to be mapped onto a binary domain {0,1}, depending on whether the corresponding reaction is active or switched off. Finally, the introduction of the on/off threshold restricts the possible optimal configurations to a few points in a discrete variable space [14]. A further limitation of these methods is that they consider only one objective, or a linear combination of objectives (usually encoded in the biomass reaction); this does not allow one to fully explore the metabolic potential of the organism, where competing goals may lead to a maximized production of a given chemical while simultaneously ensuring high growth rate.
Note that the expression level of nested gene sets is computed by applying these rules recursively.
Multilevel conditionspecific metabolic networks
In the present study, each condition and the corresponding gene expression profile are used as indicators for the activity of the associated reactions in the model, and are therefore mapped onto the model and finally associated with a point in a multidimensional phenotypic space. We assume that the genetic level is slower than the metabolic one; the time that chemical reactions within the cell need in order to reach a steady state is usually less than a minute [17], and therefore the steady state is reached faster than the variation of enzyme concentrations caused by changes in the gene expression profile. Furthermore, we account for the fact that the importance of a gene for the metabolism is negatively correlated with the variance of that gene across the experimental conditions [18]. Essential genes are in fact more tightly regulated and evolve slowly, as a high variance of gene expression level of essential genes would affect all the downstream genes and is more likely to be less tolerated in a large interaction network [19, 20]. For instance, essential genes in the metabolic network are those coding for an enzyme controlling a reaction which is upstream of many others (e.g. upstream of the TCA cycle). To this end, we take into account the inverse of the variance of a gene across conditions as a multiplicative factor for the lower and upper bound of the flux regulated by the gene.
(and φ(Θ _{ i })=1 if Θ _{ i }=1) maps the gene set expression value onto the metabolic model. The sign operator is defined as sgn(Θ _{ i }−1)=(Θ _{ i }−1)/Θ _{ i }−1. \({\sigma _{i}^{2}}\) is the variance of the ith gene set, computed from the variance of the genes involved using the rules (1) defined to map the gene expressions to the gene set expressions. γ is the weight for the variance, and constitutes the reliability of the variance as an indicator of the importance of the genes in the model.
The reasons for choosing this mathematical structure are as follows. Logarithmic and multilayer processes are not uncommon in biology [21, 22], where many processes display a multiple layer architecture in order to produce amplification in cascades (e.g. blood clotting [23] and MAP kinase cascades [24]). Although the correlation between gene expression and metabolic phenotype is still a matter of debate, our assumption of a logarithmic map φ is supported by recent evidence that, although protein synthesis rate increases with increasing mRNA abundance, the rate of increase is lower for high mRNA abundance [25]. Furthermore, φ provides full compatibility with the abovementioned Boolean on/off approaches, as it approaches zero when the expression level approaches zero. Being an approximation of a linear function around 1 (which in our method represents the wildtype gene expression level for E. coli), we also capture the property of roughly linear relation between gene expression and enzyme activity for the wildtype E. coli [26, 27].
The definition of φ is also useful when searching for optimal gene expression values in a given condition [28]. Such algorithm would not be encouraged to search for unrealistically high values of gene expression levels, since this would not be converted into weak constraints (as it would happen with, e.g., a linear map). As introduced above, the parameter γ is quantifying the role attributed to the variance as an indicator of the importance of a gene (we assume that low variance across different conditions means high importance). By increasing the parameter γ, we increase the ability of the gene expression values to fine tune the final reaction rates. The method is robust with respect to perturbations of γ, while variations of orders of magnitude lead to an increase of the metabolic sensitivity to the different environmental conditions.
Here we use the logarithmic map (3) to set constraints for the metabolic model, then we solve the linear program (2) to find the flux distribution. While suggesting a logarithmic map, we remark that φ is fully adjustable depending on the type of model and on further data available on protein abundance, which would allow a reactionspecific definition of φ. Given existing evidence for a logarithmlike regulation of reaction rates, when incorporating additional experimental data into the model, we suggest keeping the logarithmic map and calibrating the base of the logarithm and γ as adjustable parameters. With good accuracy, our method enables the prediction of the conditionspecific growth rate measured experimentally. It also outperforms existing methods based on linear associations between gene expression and reaction flux bounds. Overall, it shows high correlation between experimental and predicted rates (see Results ).
Multiplex networks to model omics response to environmental conditions
Multilayer networks are particular types of networks used to model interactions between components in a system where different channels of connectivity are explicitly taken into account. A channel is associated with a layer, and each node may have different interactions with the other nodes, depending on the layer where these interactions take place [29, 30]. Our focus is on investigating similarity between environmental growth conditions by analyzing the E. coli response from a multiomic perspective. To this end, we link the two layers through a geneexpression driven trilevel linear program. This multilevel structure accounts for the fact that the actual bacterial response to any environmental condition is highly dependent on multiple cellular objectives that the bacterium is required to meet [31, 32].
Since we are interested in analyzing conditions and finding condition similarity from the transcriptomic and fluxomic viewpoints, we construct a multiplex with these two layers, where nodes are environmental conditions. The edges between nodes in each layer represent a measure of similarity between conditions in terms of gene expression profile and metabolic flux profile respectively. In order to obtain a global picture of similarity between conditions, we aggregate both layers in a singlelayer network through a weighted network fusion approach. After aggregating the system into a single network, we adapt measures on graphs with the aim of identifying clusters of similar conditions, where the similarity takes into account the two omics (transcriptomic and fluxomic) and the metabolic network.
Similarity network fusion (SNF)
Motivation
In recent years, complex networks have come to increased prominence as a dataset type. Many types of information, particularly in biology, are best represented in terms of pairwise interactions (edges) between entities (nodes), which can be a particularly compact form for large systems where only a small proportion of elements have interesting interactions.
As more datasets are gathered and stored as complex networks, it becomes increasingly common to find multiple network datasets (edge sets) that refer to the same entities (node sets). We might suppose that sometimes these network datasets are different reflections of some common underlying network, and we may wish to integrate the known dataset to discover the structure of the underlying, unknown, dataset. This is the purpose of Similarity Network Fusion (SNF) [7].
Definition of W
W are the base similarity networks that need to be integrated. In each layer, if the entries of W (edge weights) are already representing similarity between nodes of the network, then they only need to be normalized to be positive and in [ 0,1], e.g. by squaring and dividing by the maximum. For instance, social networks are often natively similarity networks, where edge weights are counts of various friendly events between individuals.
Definition of P _{0}
which is more robust to numerical instabilities.
Definition of S
Core operation
where v is the index of each of the m layers, k ranges over all layers except for the one under consideration, n is the iteration number, \(\mathbf {P}^{(v)}_{n}\) is the similarity matrix, and \(\mathbf {S}^{(v)}_{n}\) is the local similarity matrix, both referred to the vth layer at the nth iteration.
For a set of networks (layers of a multivariate network), this iteration repeatedly makes each network more similar to the others until they converge to a single integrated network.
Code listings
For those readers who find source code more illuminating than equations, the file https://github.com/maxconway/SNFtool/blob/master/R/SNF.R is a good place to start.
Biased similarity network fusion
Since our multiomic model and trilevel formulation provide phenotypic flux rates from gene expression profiles associated with growth conditions, one may object that the transcriptomic layer should carry no weight in the fusion process. However, we remark that the fluxomic layer is a result of predictions of the metabolic model, and therefore it should not be considered as the only indicator of the response to a given condition. In fact, in applications to genomescale models (as in this study), we suggest setting the weight as the level of confidence in the model itself.
where \(d_{\textit {ij}}^{(v)}\) denotes the Euclidean distance between the two arrays representing conditions i and j in the vth layer.
Let P be the multiplex network of similarity between environmental conditions. P is computed, using (10) on m=2 omic levels, i.e., from the distance between gene expression arrays in the transcriptomic layer, and from the distance between flux rate arrays in the phenotype layer. The central equation of standard SNF is Eq. 9, which is iterated to actually conduct the fusion process, namely a finite number of message passing steps in which the m layers coevolve. Eq. 9 describes the update step for each of the m status matrices P ^{(v)}, representing similarities of conditions in each layer. These matrices are initialized as a normalized form of the similarity matrices of the m networks. S ^{(v)} are kernel matrices, giving a normalized form of similarity only to the K nearest neighbors.
In addition to introducing a layer bias, we added convergence detection based on the first and second discrete derivatives, and parallelized each iteration for an m−fold performance increase when fusing m networks. Our weighted network fusion method is freely available at https://github.com/maxconway/SNFtool.
Results
Estimating optimal bacterial response in multidimensional phenotypic spaces
A common assumption in systems biology is that microorganisms tend to shape their metabolic network in order to maximize the growth rate (biomass). However, whether the biomass is the right objective for the analysis of the metabolism is still a matter of debate [33]. Furthermore, there is increasing evidence that bacteria have to cope with multiple, sometimes competing, objectives that must be fulfilled simultaneously [34]. It is also likely that evolution has shaped cells in order to reach an optimal tradeoff between all their objectives [35]. This suggests that a singleobjective approach, the maximization of the growth rate, may not be appropriate in many systems biology applications.
To enable a multiobjective view of the metabolism, we base our method on a multilevel linear program aimed towards classifying the bacterial response associated with external or growth media conditions in any multidimensional objective space chosen by the researcher. Our pipeline can be used to design environmental conditions that likely lead to a desired phenotype. Due to its realvalued kernel, it is well suited for applications involving various perturbations, such as different codon usage bias or continuous genetic manipulations.
When maximizing for acetate and formate production, E. coli shows greater variability and different outcomes in different conditions. Conversely, when maximizing for succinate and ethanol production, less conditions are able to ensure high succinate or ethanol production. Interestingly, many conditions able to ensure high ethanol production are anaerobic. In the acetatebiomass space, the two tradeoff curves intersect, due to the fact that in one anaerobic condition E. coli is able to produce 22.72 mmol h ^{−1} gDW ^{−1} of acetate while keeping a growth rate of 0.27 h ^{−1}. This configuration is not reached in aerobic conditions. We also note that the presence of oxygen strongly affects the bacterial response in the ethanolbiomass space (e.g. compared with the formatebiomass, acetatebiomass and succinatebiomass spaces).
The distinction between aerobic and anaerobic conditions is of great metabolic importance, and it is reliably recorded in the Colombos condition database. From Fig. 3, we note that under almost all circumstances, the aerobic Pareto front covers more “metabolic space”, since adding oxygen allows more possible metabolic configurations. Furthermore, a large number of conditions excrete none of each of the byproducts. This is also unsurprising, since excreting these byproducts effectively means excreting energy, which has obviously been selected against. We also note that, in the graphs showing biomass and byproduct excretion, for all byproducts except for succinate, there is a clear scarcity of conditions in the high byproduct, low biomass region (top left). We hypothesize that this is due to the fact that if the overall metabolism of the cell is slow, i.e. low biomass, it will be more difficult for the cell to spare resources to give high byproduct production while still remaining alive.
Interestingly, only a few conditions are able to produce simultaneously succinate and ethanol. In fact, the succinateethanol and biomassethanol subspaces are where the difference between aerobic and anaerobic condition is more evident. This is mainly due to the difficulty in producing ethanol and high amount of biomass in anaerobic condition (maximum 10.71 mmol h ^{−1} gDW ^{−1} of ethanol, with a biomass of 0.27 h ^{−1}). Indeed, the tradeoff in anaerobic condition is far from the trade off in aerobic condition especially in the lowbiomass and highbiomass intervals. This is not observed in the succinatebiomass phenotypic space, where the aerobic tradeoff is close to the aerobic tradeoff even at high biomass production rates.
Validation on a phenomics dataset of growth conditions
We report that our method consistently outperforms the methods for integration of gene expression data in FBA that employ a linear map between gene expression values and the multiplicative factors for the flux bounds (an approach used, e.g., in [38]). For instance, with the dataset used in this study, linearly increasing and decreasing bounds according to gene expression data only leads to a Pearson correlation of 0.13 (pvalue=0.67) and, surprisingly, to a negative Spearman correlation of −0.19 (pvalue=0.52) between predicted and measured values.
The inclusion of the underground reactions causes a relevant effect in the prediction of growth and production of additional objectives performed through fluxbalance analysis. This is due to the fact that underground reactions share several metabolites with the iJO1366 reactions. Some underground reactions also take part in pathways leading to biomass precursors.
Multiomic and multicondition network fusion
In the transcriptomic layer, a condition is represented by a gene expression array; conversely, in the phenotypic layer, it is represented by a flux rate array. Our phenotypic flux data and, to a lesser extent, Colombos expression data displayed high kurtosis, indicating that both layers display many outliers. We do not remove outliers for two main reasons: first, in some cases this would remove most of the data; second, this represents a biologically reasonable allornothing regulation, as one would see in a bistable system created by positive feedback in regulation (for instance [39]).
Reactions identified as important classifiers by the decisiontreebased prediction of clusters
Reaction name  Importance (rounded %)  Cluster 1 LB  Cluster 3 UB 

Biomass generation  17  1.009824  1.0079 
CHEBI:44800 production  17  0.001353165  0.001350586 
5deoxyribose production  17  0.0002332694  0.000232825 
pCresol production  17  0.0002251908  0.0002247617 
CHEBI:16490 production  17  2.019649e06  2.0158e06 
Other  15 
To detect those fluxes and expressions that are most indicative of the three configurations, we employed a regression technique based on recursive partitioning [41]. This suggested a number of decision trees with only two rules splitting on flux values, each able to assign the correct cluster in 97 % of cases. We validated these fluxes as genuinely important via bootstrapping: on repeated samples of 80 % of the data, the same fluxes were detected as important. The fluxes we found to be closely associated with clustering are reported in Table 2. Of these, we can immediately associate biomass generation and 5deoxyribose production with growth rate. The other exchange reactions may have specific relationships with growth rate configurations, or their detection may be due to a general correlation between rates of bacterial growth and excretion of byproducts. The orange and green bars in Fig. 5 show how effective these fluxes are in partitioning the data into clusters. Interestingly, two out of five predictors are underground reactions.
Discussion
Multilayer networks are increasingly being associated with microorganisms, with a rapidly growing number of applications in systems biology [42]. Bacterial metabolism can itself be thought of as a multilayer network. This is due to the fact that a microorganism builds complexity in time and space by evolutionary increasing its capability of performing computation through chemical reactions over time, increasing the so called “bacterial computational capability” [43], usually through late transfer of modules that have spread across bacterial species. Bacteria have therefore hierarchical regulation machineries for processing gene expression to RNA, RNA to proteins, proteins to posttranslational modified (or complexes of) proteins, and finally complexes of elements to structural organized macrostructures.
While singlelayer networks are able to model interactions between components, their strongest assumption, i.e. all the interactions between two nodes can be mapped as a single link, does not always hold in biology. Interactions between the same nodes can occur at different levels and in different settings [44, 45]. For instance, in case of a compendium of environmental growth conditions for a microorganism, the bacterial response is measured at different omic levels. In this study, we analyzed the transcriptomic and fluxomic levels.
Conclusion
When we are interested in classifying similarity among growth conditions, and in finding communities structure within their network, we must account for the different omic levels where the response to these growth conditions can be measured.
After building a genomescale reconstruction of Escherichia coli that takes into account underground metabolism, we map a set of growth conditions to the biomassacetateformate and biomasssuccinateethanol spaces of flux rates. This allows us to estimate the metabolic potential of E. coli, and to predict regions of the fluxomic space where the bacterium operates in different conditions. As a result, each condition is represented by a gene expression profile in the transcriptomic layer, and by a flux profile in the fluxomic layer. In both layers, a network of conditions is built taking into account similarity between gene expression profiles and predicted flux rates respectively. We then fuse the resulting multiplex network into a singlelayer network, by introducing a bias to represent the weight put on the phenotype layer rather than on the gene expression layer. This weight should be related to the accuracy of the genomescale metabolic reconstruction. For instance, for large and exhaustive models (e.g. the E. coli adopted in this study), the information we gather after taking into account the predicted phenotype is likely to be more detailed than merely considering the gene expression levels.
Our multiomic network fusion method builds on SNF to allow us to: (i) specify a bias between layers; (ii) calculate a similarity matrix in a robust fashion despite extremely high kurtosis distributions; and (iii) attribute the clusters in the fused similarity matrix to a subset of dominant reactions.
Integrating multiplex networks of omic levels in combination with metabolic models serves as an indicator of the behavior of the phenotype across the phasespace of all possible gene expression profiles, and their associated metabolic networks. This paves the way for the highlevel topological understanding of multiomic and multicondition models. Assigning a specific metabolic model to each experimental condition proves useful in a number of applications. For instance, using the phasespace of conditions, one can predict and investigate metabolic regions in which the bacterium can work under many experimental settings, as well as unapproachable regimes, or regions in which it can grow only in specific conditions. Furthermore, given topological information on the space of experimental conditions, one can infer the position of nontested or incomplete expression profiles.
Multiomic data modelling will also provide insights into the complexity of regulatory metabolic mechanisms. However, such complexity could be meaningfully analyzed and used in metabolic engineering only through comparative analysis of the organismal response to thousands of environmental conditions. It is noteworthy that this approach resembles the deep phenotyping techniques used in medicine.
Availability
The datasets supporting the conclusions of this article are included within the article (and its additional files). Our weighted network fusion method for genomescale models is freely available at https://github.com/maxconway/SNFtool.
Abbreviations
 SNF:

Similarity network fusion
 FBA:

Flux balance analysis
 GPR:

Geneproteinreaction
 LB:

Lower bound
 UB:

Upper bound
Declarations
Acknowledgements
We would like to thank Naruemon Pratanwanich for initial discussion regarding this manuscript.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Chalise P, Koestler DC, Bimali M, Yu Q, Fridley BL. Integrative clustering methods for highdimensional molecular data. Transl Cancer Res. 2014; 3(3):202.PubMedPubMed CentralGoogle Scholar
 Chindelevitch L, Trigg J, Regev A, Berger B. An exact arithmetic toolbox for a consistent and reproducible structural analysis of metabolic network models. Nat Commun. 2014; 5:4893.View ArticlePubMedPubMed CentralGoogle Scholar
 Saha R, Chowdhury A, Maranas CD. Recent advances in the reconstruction of metabolic models and integration of omics data. Curr Opin Biotechnol. 2014; 29:39–45.View ArticlePubMedGoogle Scholar
 Bordbar A, Monk JM, King ZA, Palsson BO. Constraintbased models predict metabolic and associated cellular functions. Nat Rev Genet. 2014; 15(2):107–20.View ArticlePubMedGoogle Scholar
 Faith JJ, Hayete B, Thaden JT, Mogno I, Wierzbowski J, Cottarel G, et al.Largescale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biol. 2007; 5(1):8.View ArticleGoogle Scholar
 Shen R, Olshen AB, Ladanyi M. Integrative clustering of multiple genomic data types using a joint latent variable model with application to breast and lung cancer subtype analysis. Bioinformatics. 2009; 25(22):2906–12.View ArticlePubMedPubMed CentralGoogle Scholar
 Wang B, Mezlini AM, Demir F, Fiume M, Tu Z, Brudno M, et al.Similarity network fusion for aggregating data types on a genomic scale. Nat Methods. 2014; 11(3):333–7.View ArticlePubMedGoogle Scholar
 Orth JD, Conrad TM, Na J, Lerman JA, Nam H, Feist AM, et al.A comprehensive genomescale reconstruction of Escherichia coli metabolism. Mol Syst Biol. 2011; 7(1):535.View ArticlePubMedPubMed CentralGoogle Scholar
 Notebaart RA, Szappanos B, Kintses B, Pál F, Györkei Á, Bogos B, et al.Networklevel architecture and the evolutionary potential of underground metabolism. Proc Natl Acad Sci. 2014; 111(32):11762–7.View ArticlePubMedPubMed CentralGoogle Scholar
 Guzmán GI, Utrilla J, Nurk S, Brunk E, Monk JM, Ebrahim A, et al.Modeldriven discovery of underground metabolic functions in escherichia coli. Proc Natl Acad Sci. 2015; 112(3):929–34.View ArticlePubMedPubMed CentralGoogle Scholar
 Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, et al.Gapped BLAST and PSIBLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997; 25(17):3389–402.View ArticlePubMedPubMed CentralGoogle Scholar
 Palsson BØ. Systems biology: Constraintbased reconstruction and analysis: Cambridge University Press; 2015.
 Machado D, Herrgård M. Systematic evaluation of methods for integration of transcriptomic data into constraintbased models of metabolism. PLoS Comput Biol. 2014; 10:1003580.View ArticleGoogle Scholar
 Angione C, Costanza J, Carapezza G, Lió P, Nicosia G. A design automation framework for computational bioenergetics in biological networks. Mol BioSyst. 2013; 9(10):2554–64.View ArticlePubMedGoogle Scholar
 Angione C, Lió P. Predictive analytics of environmental adaptability in multiomic network models. Sci Rep. 2015; 5:15147.View ArticlePubMedPubMed CentralGoogle Scholar
 Angione C, Pratanwanich N, Lió P. A hybrid of metabolic flux analysis and bayesian factor modeling for multiomics temporal pathway activation. ACS Synth Biol. 2015; 4(8):880–9.View ArticlePubMedGoogle Scholar
 Durot M, Bourguignon PY, Schachter V. Genomescale models of bacterial metabolism: reconstruction and applications. FEMS Microbiol Rev. 2009; 33(1):164–90.View ArticlePubMedGoogle Scholar
 Mar JC, Matigian NA, MackaySim A, Mellick GD, Sue CM, Silburn PA, et al.Variance of gene expression identifies altered network constraints in neurological disease. PLoS Genet. 2011; 7(8):1002207.View ArticleGoogle Scholar
 Pál C, Papp B, Hurst LD. Highly expressed genes in yeast evolve slowly. Genetics. 2001; 158(2):927–31.PubMedPubMed CentralGoogle Scholar
 Fraser HB, Hirsh AE, Steinmetz LM, Scharfe C, Feldman MW. Evolutionary rate in the protein interaction network. Science. 2002; 296(5568):750–2.View ArticlePubMedGoogle Scholar
 Paltanea M, Tabirca S, Scheiber E, Tangney M. Logarithmic growth in biological processes. In: Computer Modelling and Simulation (UKSim), 2010 12th International Conference On. IEEE: 2010. p. 116–21.
 Guimaraes JC, Rocha M, Arkin AP. Transcript level and sequence determinants of protein abundance and noise in Escherichia coli. Nucleic Acids Res. 2014; 42(8):4791–9.View ArticlePubMedPubMed CentralGoogle Scholar
 Macfarlane R. An enzyme cascade in the blood clotting mechanism, and its function as a biochemical amplifier: Nature Publishing Group; 1964.
 Nishida E, Gotoh Y. The map kinase cascade is essential for diverse signal transduction pathways. Trends Biochem Sci. 1993; 18(4):128–31.View ArticlePubMedGoogle Scholar
 Firczuk H, Kannambath S, Pahle J, Claydon A, Beynon R, Duncan J, et al.An in vivo control map for the eukaryotic mRNA translation machinery. Mol Syst Biol. 2013; 9(1):635.View ArticlePubMedPubMed CentralGoogle Scholar
 Shimizu K. Metabolic flux analysis based on 13clabeling experiments and integration of the information with gene and protein expression patterns. Adv Biochem Eng Biotechnol. 2004; 91:1.PubMedGoogle Scholar
 Peng L, Shimizu K. Global metabolic regulation analysis for Escherichia coli k12 based on protein expression by 2dimensional electrophoresis and enzyme activity measurement. Appl Microbiol Biotechnol. 2003; 61(2):163–78.View ArticlePubMedGoogle Scholar
 Angione C, Carapezza G, Costanza J, Lió P, Nicosia G. Pareto optimality in organelle energy metabolism analysis. Comput Biol Bioinformatics IEEE/ACM Trans. 2013; 10(4):1032–44.View ArticleGoogle Scholar
 Estrada E, GómezGardeñes J. Communicability reveals a transition to coordinated behavior in multiplex networks. Phys Rev E. 2014; 89(4):042819.View ArticleGoogle Scholar
 Nicosia V, Bianconi G, Latora V, Barthelemy M. Growing multiplex networks. Phys Rev Lett. 2013; 111(5):058701.View ArticlePubMedGoogle Scholar
 Arias CF, Catalán P, Manrubia S, Cuesta JA. toyLIFE: a computational framework to study the multilevel organisation of the genotypephenotype map. Sci Rep. 2014; 4:7549. [doi:10.1038/srep07549].View ArticlePubMedPubMed CentralGoogle Scholar
 Fong SS, Joyce AR, Palsson BØ. Parallel adaptive evolution cultures of Escherichia coli lead to convergent growth phenotypes with different gene expression states. Genome Res. 2005; 15(10):1365–72.View ArticlePubMedPubMed CentralGoogle Scholar
 Birch EW, Udell M, Covert MW. Incorporation of flexible objectives and timelinked simulation with flux balance analysis. J Theor Biol. 2014; 345:12–21.View ArticlePubMedGoogle Scholar
 Schuetz R, Zamboni N, Zampieri M, Heinemann M, Sauer U. Multidimensional optimality of microbial metabolism. Science. 2012; 336(6081):601–4.View ArticlePubMedGoogle Scholar
 Zakrzewski P, Medema MH, Gevorgyan A, Kierzek AM, Breitling R, Takano E. Multimeteval: comparative and multiobjective analysis of genomescale metabolic models. PloS ONE. 2012; 7(12):51511.View ArticleGoogle Scholar
 Meysman P, Sonego P, Bianco L, Fu Q, LedezmaTejeida D, GamaCastro S, et al.Colombos v2. 0: an ever expanding collection of bacterial expression compendia. Nucleic Acids Res. 2014; 42(D1):649–53.View ArticleGoogle Scholar
 Hui S, Silverman JM, Chen SS, Erickson DW, Basan M, Wang J, et al.Quantitative proteomic analysis reveals a simple strategy of global resource allocation in bacteria. Mol Syst Biol. 2015; 11(2):784.View ArticlePubMedPubMed CentralGoogle Scholar
 Brandes A, Lun DS, Ip K, Zucker J, Colijn C, Weiner B, et al.Inferring carbon sources from gene expression profiles using metabolic flux models. PloS ONE. 2012; 7(5):36947.View ArticleGoogle Scholar
 Ferrell JE. Feedback regulation of opposing enzymes generates robust, allornone bistable responses. Curr Biol. 2008; 18(6):244–5.View ArticleGoogle Scholar
 Angione C, Conway M, Lió P. Spectral clustering performed on the fused network. https://github.com/maxconway/supplementarydata/blob/master/BMCbioinformatics2016/Fused_similarities.csv.
 Therneau T, Atkinson B, Ripley B. Rpart: Recursive Partitioning and Regression Trees. 2015. R package version 4.19. http://CRAN.Rproject.org/package=rpart.
 Castellani G, Intrator N, Remondini D. Systems biology and brain activity in neuronal pathways by smart device and advanced signal processing. Front Genet. 2014; 5:253.View ArticlePubMedPubMed CentralGoogle Scholar
 Angione C, Costanza J, Carapezza G, Lió P, Nicosia G. Analysis and design of molecular machines. Theor Comput Sci. 2015; 599:102–17.View ArticleGoogle Scholar
 Boccaletti S, Bianconi G, Criado R, Del Genio C, GómezGardeñes J, Romance M, et al.The structure and dynamics of multilayer networks. Phys Rep. 2014; 544(1):1–122.View ArticleGoogle Scholar
 Moni MA, Liò P. Networkbased analysis of comorbidities risk during an infection: Sars and hiv case studies. BMC Bioinformatics. 2014; 15(1):333.View ArticlePubMedPubMed CentralGoogle Scholar