article: Nonalcoholic fatty liver disease (NAFLD) is associated with metabolic syndromes spanning a wide spectrum of diseases, from simple steatosis to the more complex nonalcoholic steatohepatitis. To identify the deregulation that occurs in metabolic processes at the molecular level that give rise to these various NAFLD phenotypes, algorithms such as pathway enrichment analysis (PEA) can be used. These analyses require the use of predefined pathway maps, which are composed of reactions describing metabolic processes/subsystems. Unfortunately, the annotation of the metabolic subsystems can differ depending on the pathway database used, making these approaches subject to biases associated with different pathway annotations, and these methods cannot capture the balancing of cofactors and byproducts through the complex nature and interactions of genome-scale metabolic networks (GEMs). Here, we introduce a framework entitled Minimum Network Enrichment Analysis (MiNEA) that is applied to GEMs to generate all possible alternative minimal networks (MiNs), which are possible and feasible networks composed of all the reactions pertaining to various metabolic subsystems that can synthesize a target metabolite. We applied MiNEA to investigate deregulated MiNs and to identify key regulators in different NAFLD phenotypes, such as a fatty liver and liver inflammation, in both humans and mice by integrating condition-specific transcriptomics data from liver samples. We identified key deregulations in the synthesis of cholesteryl esters, cholesterol, and hexadecanoate in both humans and mice, and we found that key regulators of the hydrogen peroxide synthesis network were regulated differently in humans and mice. We further identified which MiNs demonstrate the general and specific characteristics of the different NAFLD phenotypes. MiNEA is applicable to any GEM and to any desired target metabolite, making MiNEA flexible enough to study condition-specific metabolism for any given disease or organism.
  Nonalcoholic fatty liver disease (NAFLD) is the most common cause of liver disease in western countries [1], affecting an estimated 25% to 45% of the general US population, though this percentage is much higher among individuals suffering from obesity and diabetes [1,2]. There are two categories of NAFLD: steatosis, which is an accumulation of fat in the liver, and nonalcoholic steatohepatitis (NASH), which consists of the additional presence of liver inflammation and hepatocellular injury with or without fibrosis [2]. Gene deregulation occurs when the cell no longer maintains precise control over certain genes, and therefore expressed proteins, and is associated with various diseases such as NAFLD. Determining which genes are deregulated in disease and how these genes are deregulated can increase the understanding of any disease as well as provide new pathways for therapeutic treatments. To study gene deregulation, pathway enrichment analysis (PEA) uses maps of metabolic processes and subsystem reactions to determine if a given list of genes (gene set) is associated with a certain biochemical pathway, shedding important conceptual insight into gene deregulation. Many PEA algorithms, such as ConsensusPathDB [3] and Piano [3,4], have been developed to identify biological processes based on gene sets, while other algorithms, such as IMPaLA [5], are based on both gene and metabolite sets. While these algorithms can successfully provide insights about various disease phenotypes, they use predefined metabolic maps that are limited to our current knowledge of these hugely intricate pathways and can differ depending on the pathway database used, meaning they could miss information about the wide range of complex metabolic interactions. These approaches are also subject to biases associated with the different pathway annotations, and they cannot capture the balancing of cofactors and byproducts through the complex nature and interactions of genome-scale metabolic networks (GEMs). To overcome this, a method has been developed based on elementary flux modes (EFMs), which are non-decomposable flux distributions in metabolic networks [6]. These flux distributions indicate a conceptual description of metabolic pathways. This method computes EFMs from a given GEM and uses tissue-specific transcriptomics to identify a subset of tissue-specific EFMs. However, the enumeration of all EFMs in a GEM can be computationally intractable, and EFMs are not necessarily specific to a target metabolite synthesis, meaning that targeted pathways are not always able to elucidated using EFMs. Modeling methods of biological networks, such as metabolic, signaling, and protein-protein interaction can be used as a framework to study biological functions of the liver for the prevention and treatment of liver-associated disease [7]. For example, one can use graph-based methods (GBM) that require network topologies and gene expression profiles as inputs to extract deregulated subnetworks that occur between two conditions [8–10]. These methods have been applied to protein-protein networks for prostate cancer study, and they are potentially applicable to metabolic networks. GBM methods use graph-theoretic properties on network topologies but miss information about additional constraints, such as mass balance [11] and thermodynamics [12]. Ideally, the results for these studies would include a set of mass balanced subnetworks that could be used to understand the carbon, energy, and redox flows from precursor metabolites to target metabolites and complex metabolic tasks. Here, we propose a method called Minimal Network Enrichment Analysis (MiNEA) that compares two conditions using transcriptomics, proteomics, and metabolomics to identify deregulated minimal networks (MiNs), which are reactions pertaining to various metabolic subsystems that can synthesize a target metabolite. The MiNEA algorithm works by formulating metabolic tasks (MTs) to mimic the various NAFLD phenotypes, such as lipid droplet formation, lipoapoptosis, liver inflammation, and oxidative stress. These MTs include the known set of metabolites necessary to form the modeled phenotype. The cell can use different pathways to form the same metabolites and create the same phenotypes, and the purpose of MiNEA is to identify deregulated alternative routes for a given MT between two conditions, such as a control and treatment. Additionally, MiNEA uses these MTs to compute alternative thermodynamically feasible MiNs applying thermodynamic constraints [12,13] to a mouse GEM [14]. We used mouse and human liver sample expression data [1,15] and identified deregulated metabolic processes in mice and humans that potentially lead to the NAFLD phenotypes. MiNEA identified an upregulation in the oxidative stress synthesis network in mice, but this network was found to be downregulated in humans. The cholesterol and triacylglycerol synthesis networks were deregulated in humans only, while the ceramide synthesis network was only deregulated in mice. We found downregulated reactions in the synthesis network for cholesteryl esters, cholesterol, and alanine in both humans and mice. We further identified downregulated reactions in the superoxide anion (SOA) synthesis network in humans, specifically in NASH as opposed to steatosis, while upregulation was found in the SOA synthesis network in mice. This perturbation through the SOA synthesis network in NASH suggested an unbalanced ceramide synthesis, and studies have shown that the ceramide is a key regulator of apoptosis and promotes fibrosis in the hepatic steatosis model [16,17]. MiNEA can generate MiNs for any target metabolite, e. g. a metabolite produced under specific phenotypes or a biomass building block [18] that is needed for cell growth. Additionally, MiNEA can integrate condition- and context-specific omics data to understand the deregulated phenotypes associated with a set of differentially expressed genes. These characteristics make MiNEA a versatile tool for exploring and understanding different metabolic phenotypes.
  Nonalcoholic fatty liver disease (NAFLD) has been defined as a metabolic disease associated with insulin-resistance syndrome [19]. To study the differences seen in NAFLD phenotypes between mouse and human manifestations, we collected human expression data from the three diagnosis groups, normal (N), steatosis (S), and nonalcoholic steatohepatitis (NS), and mouse expression data from control and DDC-supplemented diet conditions for three genetically different mouse strains, AJ, B6, and PWD. DDC-supplemented diet reproduces steatosis and NASH phenotypes (see Materials and Methods for detail). We refer to these data as human expression data and mouse expression data throughout this section and integrate it into a mouse model iMM1415 [14] that is constructed based on a human GEM Recon1 [20]. Sigurdsson and colleagues found that the mammalian organism with the highest number of genes homologous to Recon 1 genes was the mouse (Mus musculus) (1,415 genes, 97%). We compared iMM1415 and Recon1 and found that the iMM1415 shares 98% of reactions with Recon1 and in the remaining 2% reactions more than 75% (more than 1. 5% of total reactions) were associated with the extracellular transport mechanism. This suggests that Recon1 and iMM1415 have a very similar metabolism. Computed minimal networks do not change between Recon1 and IMM1415 due to less variability between the two models. Thus, we use iMM1415 as a generic model. If we impose tissue specific constraints or organism specific external medium, we expect some changes in the minimal networks, mainly in the number of alternative networks. However, such constraints involve a lot of uncertainty and we chose to use a generic model and our results here can be viewed as an upper bound on the possible solutions. We analyzed both human and mouse MTs using the generic model (iMM1415) with the assumption that both types of mammalian cells have a similar metabolism (See S1 Text). In the following sections, we analyzed the deregulation of minimal networks in core reactions and enriched minimal networks based on the deregulation (up- or downregulated) genes. GEMs represent an entire cellular metabolic network in the form of mathematical constraints [21]. These network reconstructions have grown rapidly in the last decades, and now many GEMs for different organisms are available [22]. GEMs and condition-specific experimental data can be employed to generate and test hypotheses using the Minimal Network Enrichment Analysis (MiNEA) framework that has been developed in this study. Compared to a method that optimizes only some part of GEMs to find subnetworks [23], MiNEA optimizes whole GEMs to find minimal networks. To examine the data that can be derived from MiNEA in terms of NAFLD, we first wanted to compute MiNs for given MTs (see Materials and Methods) that were significant phenotypes in NAFLD. Instead of representing a main linear route between precursors and targets of MTs, a MiN is composed of many subsystems integrated into subnetworks that must be active, meaning that reactions of subnetworks carry flux, to fulfill the target MT. The MiNEA algorithm facilitates enumeration of alternative MiNs and provides more flexibility to the analysis of different metabolic phenotypes and their respective environmental and genetic perturbations. The use of MTs for analysis with MiNEA allows this method to be more easily generalized and applied to the study of other metabolic phenotypes and diseases. A summary of the MiNs calculated for MTs (see Materials and Methods) are shown in Table 1. The shortest MiNs were for the synthesis of hydrogen peroxide (H2O2) (network size = 37), and the longest were for the synthesis of cholesteryl ester (network size = 131; Table 1). We found the greatest number of possible alternatives for the hexadecanoate (HDCA) synthesis and the least number possible for the phosphatidylserine (PS) synthesis (Table 1), suggesting that the HDCA synthesis is more flexible and the PS synthesis less flexible towards alternative formation compared to the rest of the MTs. Reactions that overlapped between all alternative MiNs were called high-frequency reactions (HFRs), and the percentage of HFRs shows how similar or divergent a MiN is compared to the other MiNs. The percentage of HFRs was between 62% to 91% across the MTs (Table 1). In the superoxide anion (SOA) and the hydrogen peroxide (H2O2) synthesis, 62% and 89% of the reactions were HFRs, respectively, which suggests that the alternative MiNs for the SOA synthesis were more divergent compared to the H2O2 synthesis. As an expanded example, a MiN from the H2O2 synthesis comprises many reactions from various metabolic subsystems (Fig 1). Each MiN represents a group of active reactions within multiple metabolic subsystems/pathways that are required for a MT, while a pathway is a group of annotated reactions. In the example in Fig 1, most of the active reactions for the H2O2 synthesis are from the pentose phosphate pathway (PPP) and the glycolysis/gluconeogenesis pathway. PEA differs in that it identifies a marked deregulation in a pathway as compared to MiNEA, which identifies deregulation in a MiN that could include multiple pathways. We further analyzed the HFRs within each phenotype and found that the most represented metabolic pathway in the apoptosis and inflammation (AI) phenotype was related to amino acid metabolism (# phenotypic HFRs = 36; Table 1; Fig 2). Interestingly, while the lipid droplet (LD) phenotype has many associated MTs, it shares 27 reactions in which 5 reactions are from the glycolysis/gluconeogenesis pathway (Fig 2). These shared reactions can be constitutive candidates for the LD phenotypes. The HFRs forming the subsystems in Fig 2 are potential candidates for the key regulators of the NAFLD phenotypes. The HFRs within each MT create a topology of reaction hubs in metabolic networks that can be similar to protein hubs in protein-protein interaction networks [24,25]. The misregulation of these hubs in disease can explain the origin of disease phenotypes. Therefore, we performed an analysis of HFR deregulation to study the NAFLD-specific phenotypes. We identified the up- and downregulated reactions from mouse and human expression data (described in S1 Text) and analyzed the HFRs for which the associated genes were found to be deregulated. In humans, the number of downregulated HFRs was higher compared to upregulated HFRs for the NS vs N and NS vs S across all phenotypes (Table 2). In mice, however, the NS vs N, NS vs S, and S vs N comparisons showed a higher number of upregulated HFRs compared to downregulated ones for the oxidative stress (OS) phenotype, but for the LD and AI phenotypes, we observed similar patterns of up- and downregulated HFRs to humans (Table 2). We further analyzed the deregulation in reaction hubs for each MT to identify the most perturbed MTs in NS and S and the differences between NS and S. In humans, the percentage of downregulated HFRs was higher than upregulated ones for the synthesis of all metabolites in the NS vs N and NS vs S states, while the percentage of upregulated HFRs was higher than downregulated ones for the S vs N state (Fig 3). There was a higher percentage of downregulated HFRs for the cholesterol (CHOL) and cholesteryl ester (CHOL_ES) synthesis networks in NS vs S and for the superoxide anion (SOA) and H2O2 synthesis networks in NS vs S (Fig 3). This means that downregulated HFRs have a higher impact on NASH phenotypes. Since the mouse strains AJ and PWD show phenotypes NS and S, respectively, and the mouse strain B6 shows low S and NS phenotypes after feeding with the DDC diet (Materials and Methods), we analyzed B6 separately. For mice, in NS vs N (the comparison associated to AJ) and B6 DDC vs N comparisons, we found that the percentage of downregulated HFRs in the SOA and H2O2 synthesis networks was higher than for upregulated ones (Fig 3). This indicates synthesis networks of oxidative stress were found to be perturbed for AJ and B6 mice. Additionally, for the NS vs S in mice, an elevated percentage of upregulated HFRs was identified for the synthesis of the triacylglycerol (TAG), phosphatidylinositol (PAIL), diacylglycerol (DAG), and phosphatidic acid (PA) metabolites; a high percentage of downregulated HFRs was found for the alanine (ALA) and glutamine (GLU) synthesis networks (Fig 3). This is consistent with the observation of lipid droplet perturbation in NASH and steatosis [26], and it further suggests that the synthesis of TAG, PAIL, DAG, and PA are important for lipid droplet formation. For the synthesis of ROS, such as H2O2, the human NS vs N case presented with a high percentage of downregulated HFRs, while the same case in mice presented with a high percentage of upregulated HFRs (Fig 3). The reactions of the PPP associated with genes TKT1, TKT2, GND, and TALA were identified as downregulated HFRs in the NS vs N human. However, for the NS vs N mouse and B6 DDC vs N these reactions were identified as upregulated. PPP is known as a source of nicotinamide adenine dinucleotide phosphate (NADPH) that prevents oxidative stress. Thus, the different patterns of regulation observed for these reactions can unbalance both NADPH concentration and H2O2 synthesis. Interestingly it has been reported that the unbalance in NADPH production through PPP by either an under or over-production of NADPH may induce oxidative stress [27,28]. Interestingly, Mardinoglu and colleagues [29] found an alteration in the hepatic expression levels of associated enzymes of the de novo glutathione synthesis and reported that redox imbalance is associated to metabolic dysfunction and the development of NAFLD. The acetyl CoA C acetyltransferase (ACACT1rm) and phosphoglycerate mutase reactions of the ceramide synthesis networks were identified as downregulated HFRs for the NS vs S and NS vs N in humans along with the NS vs N in mice. Thus, perturbation in the genes associated with these reactions can affect the ceramide level in NASH in both mice and humans. In line with this finding, a study suggested that deregulated ceramide production promotes liver injury and the development of NASH through disruption of endoplasmic calcium homeostasis, as well as through the inhibition of autophagy [30]. Furthermore, another study suggests that the overexpression of acid ceramidase decreases hepatic ceramide levels which is a major contributor of metabolic disease, such as obesity and hepatic steatosis [31]. The integration of metabolomics can change the size and structure of the MiNs for MTs by eliminating thermodynamically infeasible reaction directionalities. In mice, there are available metabolic profiles for the N, NS, and B6 cases [1]. We integrated these profiles into the iMM1415 [12,14] and found that, in this case, there are not marked changes in the MiNs and the result in the HFRs, which means we did not find a marked impact on the synthesis networks of NAFLD phenotypes. A further sensitivity analysis can guide future metabolomics studies to target metabolites that have a marked thermodynamic impact on a metabolic network [35]. Thus, the information content of such metabolite sets can be influential and more relevant to the phenotypes under study. We identified all possible minimal networks of given MTs separately and enforced at least 80% of the maximum synthesis rate. However, the" normal physiology" , such as activation of normal liver functions can employ a weighted combination of the networks derived by the alternatives in the set of MTs that are characteristic of the" normal physiology" . The exact" weighted combination" of MTs is rather uncertain for all MTs simultaneously, as it could create more artifacts. Therefore, we suggest generating minimal networks of given MTs separately. For a given physiology if the weights of MTs are known, then MiNEA can be modified and extended as a path forward for capturing the physiology. Gene protein reaction (GPR) rules should be able to resolve the regulation of reaction’s rates, but only if the GPR is well curated and validated with experimental data during the reconstruction process. Based on GPR rules, if a reaction is associated with a mixture of upregulated and downregulated genes, we did not consider the reaction as regulated. Here, we chose a rather conservative threshold to resolve these uncertainties, and provide an upper bound with reduced noise due to uncertainty. However, the MiNEA is flexible and can include such improved GPR associations. We should acknowledge that gene regulation based on mRNA levels is not sufficient to determine reaction rates (fluxes) because post-translational modifications, protein translation alone, and kinetics of the enzymes are really important in determining the mapping between mRNA levels and reaction rates. Depending on the availability of proteomics (and phosphoproteomics) data, and even kinetics (such V_max), they can be integrated within the MiNEA method. We applied MiNEA for the study of deregulated metabolic processes rather signaling and regulatory processes because metabolism is better characterized, as is shown with the increased availability of GEMS for many organisms [36]. Despite the challenges in reconstructing constraint-based signaling networks, recently such reconstructions have started to become available. Signaling networks for the toll-like receptor (TLR) and epidermal growth factor receptor (EGFR) are now available, and MiNEA can be extended to include these networks [37,38], meaning that MiNEA could easily be applied to the study of deregulation in signaling networks.
  The microarray gene expression data pertaining to the three diagnostic groups (normal [N], Steatosis [S], and NASH [NS]) were collected from the ArrayExpress public repository for microarray data under the accession number E-MEXP-3291 [15]. We performed an analysis of the differentially expressed genes (DEGS) and a pairwise comparison between diagnostic groups. To control the false discovery rate at level of 0. 05, multiple hypothesis testing was used [39]. NASH phenotypes can be reproduced in mouse models by treatment of chronic intoxication of with 3,5-diethoxycarbonyl-1,4-dihydrocollidine (DDC) -supplemented diet [1]. Three genetically different mouse strains AJ, B6, and PWD were fed with a DDC-supplemented diet. Under DDC-supplemented diet, the steatosis and NASH phenotypes were the most obvious in the PWD and AJ strains, respectively, while for the B6 strain, low levels of the steatosis and NASH phenotypes were observed [1]. The B6 mouse strain showed less induction of the NASH phenotype than the AJ strain [1]. The normal (N) feeding group is associated with mice fed the control diet. All mice strains, AJ, B6 and PWD, were fed this diet and are labeled as N≈AJ, N≈PWD or N≈B6 to indicate the control diet and mouse strain. After feeding the DDC-supplemented diet, AJ and PWD showed high NASH phenotypes and high steatosis phenotypes, respectively [1]. Thus, under the DDC-supplemented diet, the symbol for steatosis (S) is associated with the PWD mouse strain and the symbol for NASH (NS) is associated the AJ mouse strain (with DDC-supplemented diet: NS≈AJ, S≈PWD). Thus, NS vs N represents the AJ mouse strain treated with the DDC diet vs the AJ mouse strain fed the control diet, S vs N represents the PWD mouse strain treated with the DDC diet vs the PWD mouse strain fed the control diet, and NS vs S represents the AJ mouse strain treated with the DDC diet vs the PWD mouse strain fed the DDC diet. Metabolite concentration data and RNA-seq gene expression data from mouse liver samples of the three mouse strains (AJ, B6, and PWD) measured under the control and DDC-supplemented diet conditions were collected from the work of Pandey et al. [1]. The R-package “edgeR” [40] was used to identify the differentially expressed genes (DEGs) using three biological replicates of each mouse strain for both the control and the DDC-treated conditions. Strain-wise identification of the DEGs between the DDC-treated and control states was performed separately for each mouse strain. We used the Benjamini-Hochberg procedure implemented in edgeR to control the false discovery rate to a level of 0. 05 [39]. Steatosis can occur due to the accumulation of lipid droplets in the liver, which can subsequently lead to NASH upon liver inflammation [41]. Lipid droplets are composed of many lipid metabolites [42,43], which are summarized in Table 6. Reactive oxygen species (ROS), such as the superoxide anion, can also damage hepatic membranes and play an important role in the development of NASH [44]. To study these various phenotypes and their role in the various forms of liver disease, the metabolites associated with their various phenotypes were examined. Table 6 summarizes the key metabolites associated with the major liver disease phenotypes: lipid droplets, liver inflammation, apoptosis, and oxidative stress. We investigate the synthesis networks of these metabolites using MiNEA, where we called the synthesis of a metabolite as a metabolic task (MT). The cell can use various different pathways depending on the current state and conditions of the environment to fulfill its immediate metabolic needs. In order for MiNEA to identify deregulated alternative routes for a given MT between two conditions (Fig 4), the inputs required are a metabolic model, a list of MTs, and gene or protein expression data. In this work, we used a mouse GEM (iMM1415) reconstructed by Sigurdsson and colleagues [14], a list of MTs which are the synthesis of the metabolites shown in Table 6, and gene expression data of mouse and human liver samples [1,15]. The MiNEA calculations begin by applying thermodynamic constraints to the model as described previously [12] to eliminate thermodynamically infeasible reactions, meaning that reactions from the metabolic network can carry fluxes only if thermodynamics allows [13] (Fig 4; step 1). Then, it enumerates all the thermodynamically feasible minimal-size networks (MiNs) that are active for the given list of provided MTs (Fig 4; step 2), and all of these MiNs are composed of reactions that carry non-zero flux. For this step, the following mixed-integer linear programming (MILP) problem is applied with the objective of minimizing the number of reactions that carry flux or maximizing the number of reactions that cannot carry flux, while enforcing that the network should be able to synthesize the metabolites listed in Table 1: Maximize∑i=1RGEMzrxn, i S. v=0 0≤vi≤vmax, i*Zrxn, i vMT≥c*VMT, max RGEM is the number of reactions in a GEM and zrxn, i is a binary variable associated with reaction i. Reversible reactions were split into two forward reactions. VMT, max is the maximum yield to produce a metabolite associated with a MT. The parameter c is generates MiNs that allow for flexibility on the yield from the MT. We chose c = 0. 8 to allow for a yield of at least 80% of the maximum yield for the associated MT. All alternative MiNs that are the minimum size (msize) for the synthesis of the metabolites in Table 1 are enumerated as described by Figueiredo and colleagues [48] (Fig 4; step 3). One can enumerate many alternative networks larger than msize, such as msize+1 and msize+2. Each MiN is a subnetwork comprising a set of reactions, metabolites, and reaction-associated genes through gene-protein-reaction (GPR) association. In step 4, the deregulated genes and proteins that differ between the two conditions under study are identified (Fig 4). Finally, a hypergeometric test is performed on the sets of deregulated genes or deregulated reactions to identify the deregulated MiNs (Fig 4; Step 5). The significance of the MiN based on the deregulated genes was calculated using the hypergeometric probability density function (P), P (k, N, K, n) = (Kk) (N−Kn−k) (Nn), (1) where k and K are the numbers of deregulated genes and the number of total genes in a given MiN, respectively, and n and N are the total number of deregulated genes and the total number of genes of a metabolic network, respectively. In this context, the term “deregulated genes” is used for both up- or downregulated genes. A reaction can have three states: upregulated, downregulated, or unregulated. The regulation of a reaction was determined based on its associated differentially expressed genes. According to this metric, a reaction was identified as upregulated or downregulated if the corresponding genes are only upregulated or only downregulated, respectively. A reaction that is associated with a mixture of up- and downregulated genes is not characterized as regulated due to the inconsistency of gene expression. Up- and downregulated MiNs, which could contain various combination of up- or downregulated reactions, were identified in Step 4 of MiNEA (Fig 1) based on the total number of up- or downregulated reactions, respectively. As expected, up- and downregulated MiNs comprise markedly high numbers of up- and downregulated reactions, respectively. The significance of a MiN based on upregulated or downregulated reactions was calculated using multivariate Fisher' s hypergeometric distribution. This method has been previously used for the selection of tissue-specific elementary modes using gene expression data [6]. To identify significantly upregulated MiNs in a given set of MiNs, we selected those MiNs that contained an elevated number of upregulated reactions and as few as possible downregulated ones, while for the identification of significantly downregulated MiNs, we selected MiNs with an elevated number of downregulated reactions and as few as possible upregulated ones. Assume R and T are the total numbers of reactions in a GEM and a MiN, respectively, which can be decomposed as follows: R=Rup+Rdown+RnoT=Tup+Tdown+Tno (2) Rup, Rdown, and Rno represent the number of upregulated, downregulated, and unregulated reactions in a GEM, respectively. In the context of a MiN, T is the total number of reactions, and Tup, Tdown, and Tno are the number of upregulated, downregulated, and unregulated reactions, respectively. To consider a MiN as upregulated, we need to ensure that the pair (Tup, Tdown) in the MiN of T reactions cannot arise by chance in the context of the whole network. To obtain a better upregulation by chance, the p value was computed using Eq (3). Note that an equal or better outcome than the pair (Tup, Tdown) satisfies two conditions: (i) the number of downregulated reactions is smaller than or equal to Tdown, and (ii) the number of upregulated reactions is greater than or equal to Tup, whereas the number of reactions in the MiN remains unchanged. We can compute the p value for downregulated reactions the same way as for the upregulated reaction simply by changing up with down and down with up in the above equation. We calculated the deregulated percentage (DRP) for each MiN, which indicates the percentage of up- or downregulated genes in a MiN. For example, if a given MiN comprises 20 genes, and 5 genes of the MiN are deregulated, the then DRP is equal to 0. 4. Since the active reactions in a MiN were classified as either up- or downregulated, we calculated both upregulated percentage (UPR) and downregulated percentage (DnRP) for each MiN. To extend the degree of confidence of the results from MiNEA, we introduced the alternative minimal network frequency (AMiNF) metric that determines the significantly deregulated percentage of MiNs for a given MT. If all the alternative MiNs for a MT are significantly deregulated, AMiNF takes the value 1. The AMiNF is 0 when none of the MiNs for a MT are significantly deregulated. MiNEA was implemented using Matlab r2016a (The Mathworks, Natick, MA, USA), and MILP problems were solved using CPLEX solver (ILOG, Sunnyvale, CA, USA). The MILP gaps for all problems were converged to less than 0. 05% in less than 2,400 s. We also plan to make the MiNEA algorithm available as a tool for distribution to the community.
summary: This work aims to introduce a network-based enrichment analysis using metabolic networks and transcriptomics data. Previous pathways/subsystems enrichment methods use predefined gene annotations of metabolic processes and gene annotations can differ based on different resources and can produce bias in pathway definitions. Thus, we introduce a framework, Minimum Network Enrichment Analysis (MiNEA), which first finds all possible minimal-size networks for a given metabolic process/task and then identifies deregulated minimal networks using deregulated genes between two conditions. MiNEA also identifies the deregulation in key reactions that are overlapped across all possible minimal-size networks. We applied MiNEA to identify deregulated metabolic tasks and their synthesis networks in the steatosis and the nonalcoholic steatohepatitis (NASH) diseases using a metabolic network and transcriptomics data of mouse and human liver samples. We identified key regulators of NASH for the synthesis networks of hydrogen peroxide and ceramide in both humans and mice. We also identified opposite deregulation in NASH for the phosphatidylserine synthesis network between humans and mice. MiNEA finds synthesis networks for a given target metabolite and due to this it is flexible to study deregulation in different phenotypes. MiNEA can be widely applicable for studying context-specific metabolism for any organism because the metabolic networks and context-specific gene expression data are now available for many organisms.
keywords: steatosis
  anatomical pathology
  medicine and health sciences
  pathology and laboratory medicine
  metabolic networks
  diet
  cytopathology
  liver diseases
  animal models
  hemorrhagic fever with renal syndrome
  model organisms
  nutrition
  metabolites
  network analysis
  experimental organism systems
  gastroenterology and hepatology
  research and analysis methods
  infectious diseases
  computer and information sciences
  animal studies
  gene expression
  mouse models
  biochemistry
  genetics
  biology and life sciences
  viral diseases
  fatty liver
  metabolism
title: Investigating the deregulation of metabolic tasks via Minimum Network Enrichment Analysis (MiNEA) as applied to nonalcoholic fatty liver disease using mouse and human omics data