# THE IMPACT OF SYSTEMS MEDICINE ON HUMAN HEALTH AND DISEASE

EDITED BY: Adil Mardinoglu and Jens Nielsen PUBLISHED IN: Frontiers in Physiology and Frontiers in Neuroscience

#### *Frontiers Copyright Statement*

*© Copyright 2007-2017 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.*

*The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.*

*Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.*

*Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.*

*As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.*

*All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.*

ISSN 1664-8714 ISBN 978-2-88945-140-1 DOI 10.3389/978-2-88945-140-1

## About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

# Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

# Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

# What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# **THE IMPACT OF SYSTEMS MEDICINE ON HUMAN HEALTH AND DISEASE**

Topic Editors:

**Adil Mardinoglu,** Chalmers University of Technology & KTH - Royal Institute of Technology, Sweden **Jens Nielsen,** Chalmers University of Technology & KTH - Royal Institute of Technology, Sweden

Complex disorders including obesity, diabetes, fatty liver disease, cardiovascular disease and cancer are results from a combination of genetic, environmental and lifestyle factors. The prevalence of such disorders has increased dramatically in the last two decades and there is an urgent need for the development of new prognostic tools for the treatment of such diseases. However, this requires a deep understanding of the underlying molecular mechanisms involved in the occurrence of the diseases. With the advances in high throughput technologies, biological components of cells can be measured with a very high resolution and these data can be used for investigating whole systems properties using a network-based approach. Systems medicine provides an integrative platform for studying the interactions between the biological components of the cell using a holistic approach and generating mechanistic explanations for the emergent systems properties. This inter-disciplinary field of study allows for understanding biological processes of cells in health and disease states, gaining new insights into what drives the appearance of the disease and finally identifying proteins and metabolites implicated in human disease. Systems medicine utilizes mathematical approaches to generate models which can be employed for designing new sets of experiments and for mapping the response of the system to perturbations quantitatively. These models as well as the developed tools can accelerate the emergence of personalized medicine which can transform the practice of medicine and offer better targets for drug development with minimum side effects.

**Citation:** Mardinoglu, A., Nielsen, J., eds. (2017). The Impact of Systems Medicine on Human Health and Disease. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-140-1

# Table of Contents


# Editorial: The Impact of Systems Medicine on Human Health and Disease

Adil Mardinoglu1, 2 \* and Jens Nielsen1, 2 \*

*<sup>1</sup> Department of Biology and Biological Engineering, Chalmers University of Technology, Gothenburg, Sweden, <sup>2</sup> Science for Life Laboratory, KTH–Royal Institute of Technology, Stockholm, Sweden*

Keywords: systems medicine, metabolic networks and pathways, systems biology, networks medicine, biological networks, metabolic diseases

**Editorial on the Research Topic**

### **The Impact of Systems Medicine on Human Health and Disease**

Complex disorders including obesity, diabetes, fatty liver diseases, cardiovascular disease (CVD), and cancer are results from a combination of genetic, environmental, and lifestyle factors. The prevalence of such disorders has increased dramatically in the last two decades and there is an urgent need for the development of new prognostic tools for the treatment of such diseases. However, this requires a deep understanding of the underlying molecular mechanisms involved in the occurrence of the diseases. With the advances in high throughput technologies, biological components of cells can be measured with a very high resolution and these data can be used for investigating whole systems properties using a network-based approach. Systems medicine provides an integrative platform for studying the interactions between the biological components of the cell using a holistic approach and generating mechanistic explanations for the emergent systems properties (Mardinoglu and Nielsen, 2012). This inter-disciplinary field of study allows for understanding biological processes of cells in health and disease states, gaining new insights into what drives the appearance of the disease, and finally identifying proteins and metabolites implicated in human disease. Systems medicine utilizes mathematical approaches to generate models which can be employed for designing new sets of experiments and for mapping the response of the system to perturbations quantitatively. These models as well as the developed tools can accelerate the emergence of personalized medicine which can transform the practice of medicine and offer better targets for drug development with minimum side effects.

Genome Scale Metabolic Models (GEMs) are important tools in systems biology and employed for simulating the metabolism of cells/tissues in different states (Mardinoglu et al., 2013). GEM includes all known metabolic reactions and associated protein coding genes in a particular cell or tissue and the reconstruction of GEMs is costly and time consuming process. Pacheco et al. benchmarked the algorithms that have been used in the automated reconstruction of contextspecific GEMs. This may allow researchers to identify limitations of each method and help them to increase the quality of context-specific reconstruction algorithms and as well as the generated models. Zhang and Hua reviewed the potential applications of GEMs in systems medicine and industrial biotechnology. The authors described the key concepts and assumptions used in the application of the GEMs and provided detailed explanation about the recent applications which may promote the use of GEMs by biologists.

GEMs have been widely used in the discovery of biomarkers and drug targets that can be used in the development of efficient treatment strategies (Mardinoglu and Nielsen, 2015). These integrative models can also be used in the stratification of the patients as well as the identification of the subjects

Edited and reviewed by: *Raina Robeva, Sweet Briar College, USA*

> \*Correspondence: *Adil Mardinoglu adilm@scilifelab.se Jens Nielsen nielsenj@chalmers.se*

#### Specialty section:

*This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology*

Received: *10 October 2016* Accepted: *03 November 2016* Published: *24 November 2016*

#### Citation:

*Mardinoglu A and Nielsen J (2016) Editorial: The Impact of Systems Medicine on Human Health and Disease. Front. Physiol. 7:552. doi: 10.3389/fphys.2016.00552* at risk of developing complex diseases. CVD includes all the diseases of the heart and circulation and continues to constitute the leading cause of death globally. Therefore, there is an urgent need to develop tools and methods for identifying individuals at risk of developing CVD. Björnson et al. reviewed the current CVD risk scores and discussed how systems medicine could improve the identification of risk and maximize personalized treatment benefit.

Tumor cells alter their metabolism to maintain their growth and proliferation. Ghaffari et al. reviewed the metabolic alterations that are known to occur in cancer and underlined the use of genome-scale metabolic modeling approach in perceiving a system level perspective of cancer metabolism which can be used in discovery of biomarkers and drug targets. Özcan and Çakır have reconstructed glioblastoma multiforme (GBM) specific GEMs based on gene expression data and used their context-specific model to predict metabolic flux distributions in the brain tumor cells. The authors identified the GBM specific metabolic alterations and provided a comprehensive coverage of tumor metabolism. Moreover, flux distributions in glycolysis, glutaminolysis, TCA cycle, and lipid metabolism were predicted and validated by additional computational analysis and literature information.

Zhang et al. evaluated the potential value of current major systems biology-based approaches, including genome wide association studies, gene regulatory networks, and protein– protein interaction networks and GEMs. The authors discussed how integrative analysis of personal multi-omics data may provide increased understanding of personal genotype– phenotype relationships. Moreover, the potential benefit of integrating different type of networks in increased understanding of the interactome is discussed. Recently, metabolic networks have been integrated with other biological networks and these networks have been used in the analysis if omics data obtained from different clinical conditions (Lee et al., 2016a,b). In the same context, Brown reviewed the challenges in the generation of SuperModels where computational models were integrated with various data types and analytics. The author has also discussed how these SuperModels may assist in the development of personalized and precision medicine.

# REFERENCES


Stable isotope assisted metabolomics techniques have been used in systems medicine for metabolic flux analysis and pathway discovery. It is essential to understand the pathophysiology of dyslipoproteinemia in humans since it has major role in the progression of complex diseases including CVD and diabetes. The use of tracers labeled with stable isotopes and mathematical modeling may provide a powerful tool for probing lipid and lipoprotein kinetics in vivo. Adiels et al. reviewed the recent kinetic studies and discussed how these studies may improve our understanding of impaired human lipoprotein metabolism. Even though, stable isotope labeling analyses have been highly targeted, in recent years, tools for the global non-targeted detection, quantification, and computational analysis of isotopic enrichment have become available. Weindl et al. discussed how such non-targeted stable isotope labeling analyses can be applied for systems medicine applications.

Finally, Shmelkov et al. described an approach to improve in silico identification of a comprehensive ensemble of targets for any drug weighted by the expression of those receptors in relevant tissues. Their approach may increase the sensitivity of target detection and allow for systematic integration of bioactivity/docking scores between drugs/compounds and proteins with the expression patterns of those proteins in human tissues.

The work presented in this Research Topic addresses many of the recent progress in the applications of the GEMs is medical applications. The integration of GEMs with other biological networks and generation of whole cell models may foster the development of personalized medicine which may increase benefits and reduce risks for patients.

# AUTHOR CONTRIBUTIONS

AM and JN wrote the editorial together and approved its final version.

# ACKNOWLEDGMENTS

This work was supported by grants from the Knut and Alice Wallenberg Foundation.

Mardinoglu, A., and Nielsen, J. (2015). New paradigms for metabolic modeling of human cells. Curr. Opin. Biotechnol. 34, 91–97. doi: 10.1016/j.copbio.2014. 12.013

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Mardinoglu and Nielsen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Benchmarking Procedures for High-Throughput Context Specific Reconstruction Algorithms

Maria P. Pacheco<sup>1</sup> , Thomas Pfau1, 2 and Thomas Sauter <sup>1</sup> \*

*<sup>1</sup> Systems Biology Group, Life Sciences Research Unit, University of Luxembourg, Luxembourg, Luxembourg, <sup>2</sup> Department of Physics, Institute of Complex Systems and Mathematical Biology, University of Aberdeen, Aberdeen, UK*

Recent progress in high-throughput data acquisition has shifted the focus from data generation to processing and understanding of how to integrate collected information. Context specific reconstruction based on generic genome scale models like ReconX or HMR has the potential to become a diagnostic and treatment tool tailored to the analysis of specific individuals. The respective computational algorithms require a high level of predictive power, robustness and sensitivity. Although multiple context specific reconstruction algorithms were published in the last 10 years, only a fraction of them is suitable for model building based on human high-throughput data. Beside other reasons, this might be due to problems arising from the limitation to only one metabolic target function or arbitrary thresholding. This review describes and analyses common validation methods used for testing model building algorithms. Two major methods can be distinguished: consistency testing and comparison based testing. The first is concerned with robustness against noise, e.g., missing data due to the impossibility to distinguish between the signal and the background of non-specific binding of probes in a microarray experiment, and whether distinct sets of input expressed genes corresponding to i.e., different tissues yield distinct models. The latter covers methods comparing sets of functionalities, comparison with existing networks or additional databases. We test those methods on several available algorithms and deduce properties of these algorithms that can be compared with future developments. The set of tests performed, can therefore serve as a benchmarking procedure for future algorithms.

Keywords: metabolic networks and pathways, metabolic reconstruction, constraint-based modeling, tissue specific networks, benchmarking, validation

# 1. INTRODUCTION

Metabolic network reconstructions become ever more complicated and complete with reconstructions like Recon2 (Thiele et al., 2013) or HMR (Mardinoglu et al., 2014) containing more than 7000 reactions. While these reconstructions are a great tool for the analysis of the potential capabilities of an organism, one challenge faced by many researchers is that different cell types in multicellular organisms exhibit diverse functionality and the global generic network is too flexible. This issue has been addressed in two ways, by manually generating tissue specific models (Gille et al., 2010; Quek et al., 2014) or by creating algorithms for automatic reconstructions

#### Edited by:

*Adil Mardinoglu, Chalmers University of Technology, Sweden*

#### Reviewed by:

*Osbaldo Resendis-Antonio, Instituto Nacional de Medicina Genomica, Mexico Markus J. Herrgard, Technical University of Denmark, Denmark*

> \*Correspondence: *Thomas Sauter thomas.sauter@uni.lu*

#### Specialty section:

*This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology*

Received: *23 September 2015* Accepted: *14 December 2015* Published: *22 January 2016*

#### Citation:

*Pacheco MP, Pfau T and Sauter T (2016) Benchmarking Procedures for High-Throughput Context Specific Reconstruction Algorithms. Front. Physiol. 6:410. doi: 10.3389/fphys.2015.00410* Pacheco et al. Benchmarking Context Specific Reconstructions Methods

(Becker and Palsson, 2008; Jerby et al., 2010; Zur et al., 2010; Agren et al., 2012; Wang et al., 2012; Vlassis et al., 2014; Yizhak et al., 2014; Robaina Estévez and Nikoloski, 2015). Ryu et al. (2015) and Robaina Estévez and Nikoloski (2014) recently reviewed this field and give a good overview of the available reconstructions and point to many algorithms used in this context. While Ryu et al. (2015) are more concerned with the the state of the reconstructions, Robaina Estévez and Nikoloski (2014) focused on the applicability and properties of the available algorithms. With that many methods available, the method selection is difficult, and it is an enormous effort to try and distinguish which network, of a set of generated networks is best. Quality assessment is therefore essential but the methods used to evaluate the currently available algorithms are very diverse and it is difficult to compare them with each other. There are several approaches for validation which can essentially be split into two different categories: Consistency testing and Comparison based testing. The first is concerned with robustness against noise, e.g., missing data, and whether distinct sets of input data yield distinct models. The second commonly aims at validating the resulting model against other models or against additional data. Comparison tends to be the more common approach so far, while consistency is often ignored. This leads to the problem that algorithms are often prone to be over-specific to the comparison dataset (e.g., parameters like expression thresholds or weights working well for only one specific tissue). While comparison methods validate the reconstructed model, they are however not validating the consistency. Thus, it is possible that small differences in the input dataset can lead to vastly different networks, or even very diverse datasets yield the same models. The latter is particularly true if e.g., a biomass function is set as objective function, since it will lead to the inclusion of a multitude of reactions, which might not be necessary if a specific tissue is supplied with some metabolites by other tissues. To investigate the quality of automatically reconstructed networks it is therefore necessary to rigorously test them. In the following paragraphs, we describe multiple methods that were used in the past. **Table 1** also gives an overview of these approaches, and details which concept was used for validation of which algorithm.

# 1.1. Methods for Testing Algorithmic Consistency

The idea of consistency testing covers two major aspects: Robustness of the method and its capacity to distinguish slightly different contexts.

If feasible, random cross validation of the resulting models for a given set of input data can help to determine the robustness of the method with respect to noisy data (Vlassis et al., 2014). Left-out cross-validation allows identifying the reactions that if left-out from the input set would nevertheless be included (or excluded for inactive reactions) in the output model as their inclusion is supported by other reactions of the input set (Pacheco et al., 2015). The robustness of algorithms against noise can also be assessed by adding noise to the expression data i.e., by using a weighted combination of real and random data (Machado and Herrgård, 2014). The main issue of using random and left-out cross validation with most of the TABLE 1 | Overview of methods used for validation of automated tissue specific reconstruction algorithms.


current algorithms is that running times of several hours makes decent cross-validation with hundreds of test and validation sets infeasible. While small cross validation runs (e.g., when multiple sources of input data are available and only some sets are considered, Jerby et al., 2010) can give an indication of robustness, they cannot replace random sampling runs, which reflect noisy data much better.

To test the diversity of generated networks, many algorithms are employed to generate multiple networks and those networks are then investigated for dissimilarity (Becker and Palsson, 2008; Wang et al., 2012; Agren et al., 2014; Pacheco et al., 2015; Uhlén et al., 2015). If networks of similar cell types group together in a clustering and networks of divergent cell types are further apart, this indicates that the method does indeed generate specific networks. While it is desirable to obtain distinct networks for distinct tissues, the optimal method should not be too sensitive to small changes in the input data. Otherwise the resulting networks are prone to overfitting to the provided input data.

# 1.2. Methods for Comparison Based Testing

Comparison based testing is commonly employed to show the advantages of the presented algorithm compared over previous algorithms or to show the quality of the reconstructed network based on additional, formerly unknown, data. While the former has been employed for the validation of some algorithms (Wang et al., 2012; Vlassis et al., 2014; Robaina Estévez and Nikoloski, 2015), and becomes more important with an increasing number of available methods, it has also recently been used to compare multiple methods systematically (Machado and Herrgård, 2014; Robaina Estévez and Nikoloski, 2014). In the review by Machado and Herrgård (2014) 8 different methodologies (including GIMME, Becker and Palsson, 2008, iMAT, Zur et al., 2010 and a method by Lee et al., 2012) where tested on an independent dataset. However, their focus was on comparing the quality of flux value predictions, i.e., flux bounds specific to a condition in Escherichia coli and yeast, and not the reconstruction of tissue specific networks, i.e., the extraction of an active subnetwork.

# 1.2.1. Comparison Against Manually Curated Networks

Comparison to a manually curated tissue was employed by Agren et al. (2012) for the INIT algorithm, when they compared their automatically generated liver reconstruction to HepatoNet. However, they were restricted to a comparison on the gene level, since the source network used by INIT was the HMR database (Mardinoglu et al., 2013), while HepatoNet used its own identifiers. As they mention the difference between the reconstructed and manually curated models was partially due to absence of genes from HMR that were present in HepatoNet. Simultaneously, it is likely that the curators of HepatoNet lacked information on some of the genes present in HMR. Thus, to validate a methodology it is necessary for both the "reference" network and the source network to be compatible.

## 1.2.2. Comparison Against Additional Datasets and Databases

Similarly, many methods compare the resulting reconstructions to additional databases that contain tissue localization data (like BRENDA, Schomburg et al., 2013, HPAm Uhlén et al., 2015 or the Gene Expression Omnibus, Barrett et al., 2013), which was performed for multiple reconstruction methods (Shlomi et al., 2008; Wang et al., 2012; Robaina Estévez and Nikoloski, 2015). The common approach is to check for matches of either genes or proteins that the algorithm assigned to the tissue. This validation (and the results) are however highly dependent on whether the reconstruction method aims at creating a consistent network, or whether it allows inconsistent reactions to be part of the reconstruction. The latter will very likely increase the amount of correctly assigned genes, as enzymatic activities that cannot carry flux in the source reconstruction, would otherwise be excluded. In addition, when extracting reactions from a source network, the associated gene-protein reaction relations are commonly not altered. Thus, genes, which are inactive in a specific tissue show up as assigned to the tissue. Removing them however, could potentially be problematic if the tissue does express the removed gene under a specific condition. In this instance the tissue reconstruction would no longer contain information about this fact, and would indicate wrong potentials of the tissue. Another method that could be used as an assessment for predictive quality of an algorithm was performed by Folger et al. (2011) and subsequently by Pacheco et al. (2015). They used gene silencing data from an shRNA screen and compared it with gene essentiality predictions from a flux balance analysis (FBA) screen. The cancer network generated in this work showed an enrichment of essential genes in the genes indicated in the shRNA screen. In Pacheco et al. (2015), the list of essential genes predicted by FASTCORMICS was further compared to essential genes predicted by PRIME, MBA, mCADRE, and GIMME. Likewise bibliographic approaches have been employed to determine the agreement of reactions belonging to a certain subsystem in the reconstructed network and those subsystems being mentioned in connection with the reconstructed tissue in the literature (Shlomi et al., 2008).

To assess the predictive capability of the Model Building Algorithm (MBA), Jerby et al. (2010) used flux data from a study performed in primary rat hepatocytes and compared the ability of the source reconstruction and the generated reconstruction to predict internal fluxes given the exchange fluxes (and vice versa). This allowed them to assess whether the tissue specific network was indeed performing better in estimating the internal fluxes than the generic reconstruction (in this instance Recon1). They could show that indeed the tissue specific network had a better capability to capture the actual fluxes than the generic reconstruction. This concept was also used by Machado and Herrgård (2014) in their assessment of multiple methods for network contextualization. However, while contextualization commonly aims at altering flux bounds, which leads to a good comparability of flux measurements with predictions, tissue specific reconstruction is aiming at determining the network available in a given tissue. This means that bounds from the underlying source reconstruction are used and these are often unsuitable for the tissue of interest. But as shown by Jerby et al. (2010), even the pure network structure alteration can already improve the agreement between network fluxes and measured data, at least on a qualitative level.

A method developed by Shlomi et al. (2009) to compare the resulting network for the effects of inborn errors of metabolism (IEM) is also often used in model quality assessment. The concept is, briefly, to analyse flux ranges of the exchange reactions of the created network and compare them with clinical indications of increased or decreased metabolite levels. This concept has also been used for assessment of Recon2 (Thiele et al., 2013) who investigated a diverse set of IEMs and could show their effect even on the level of a generic reconstruction. Similarly, the authors of PRIME (Yizhak et al., 2014) used experimentally measured uptake and excretion rates and compared them to the secretion rates determined by the models their algorithm generated. While the former approach is commonly used to provide a qualitative assessment of increase or decrease in production potential, the latter results in a quantitative comparison. However, it requires the availability of uptake and secretion rates, which are commonly only available for cell lines and could be largely different in real tissues.

Another common approach to investigate the quality of reconstructions is the comparison with lists of metabolic functions. This approach is both used to validate automated reconstructions (Jerby et al., 2010; Wang et al., 2012) as well as manual reconstructions (Gille et al., 2010). The aim is to establish whether the reconstruction supports the current knowledge of the target tissue (e.g., a liver reconstruction should support the conversion of ammonia to urea), and to show that there are no structural issues in the reconstructed network (e.g., free regeneration of ATP or reductants).

# 1.3. A Benchmark for Testing Tissue Specific Reconstruction Algorithms

In this paper we present a potential benchmark that is using several of the mentioned methodologies to assess the consistency and quality of reconstructed networks and tested it with several of the available algorithms.

There are however multiple obstacles, when defining a benchmark for contextualization algorithms. There is no such thing as a "perfect" measurement, which will always leave us with noisy data to incorporate. Furthermore, we do not yet have a contextualized model that perfectly reflects a given context which could be used as a target model. In addition, the global reconstructions are not yet complete, and will likely never be and finally, there is a wide variety of data that can be used to contextualize models. Thus, to define a benchmark we will address these questions by generating networks which we define as reference networks for out testing.

The actual benchmark is preceded by a characterization of the algorithms, in which the similarity level of the context-specific reconstructions obtained with real and artificial input data is assessed. In the latter test, artificial models of different sizes were built and 50, 60, 70, 80, and 90% of the reactions of these networks were used as input for the tested algorithms. The capacity of the algorithm to distinguish between different models was compared for the different percentages of input data.

In the actual benchmark, the confidence level of the reactions included in the context-specific reconstructions using real data was assessed by matching z-scores obtained by the Barcode method (McCall et al., 2011) that basically indicate the difference in intensity between the measured intensity and the intensity distribution observed in an unexpressed state and through a comparison against the confidence score at the proteomic level of the Human Protein Atlas (Uhlén et al., 2015). In a second comparison, artificial models were built and 50, 60, 70, 80, and 90% of the reactions of these networks were used as input for the tested algorithms and the output models were then compared to the complete input model. The context-specific networks obtained with the real data were also tested for the functionalities established by Gille et al. (2010).

# 2. MATERIALS AND METHODS

# 2.1. Models Used for Benchmarking

There are currently two competing global reconstructions for humans available: Recon2 (Thiele et al., 2013) and HMR2 (Mardinoglu et al., 2013). To be able to test multiple validation techniques, we needed to select one of those reconstructions as the source network used by the tested algorithms. We decided to employ Recon2, as we used functionalities originating from HepatoNet (Gille et al., 2010), a model based on Recon1 (Duarte et al., 2007) and largely incorporated into Recon2. However, we still had to modify Recon2 to allow the algorithms to fully reconstruct HepatoNet (the procedure can be found in **Supplementary File 1**). HepatoNet was also adapted to match reactions and metabolites

#### TABLE 2 | Algorithms available for tissue specific metabolic network reconstruction.


*Most methods can use expression data as input but there are some that need additional inputs.*

with Recon2. This modified Recon2 was used as source model for all runs.

In addition to HepatoNet as a comparison model for real data, we constructed ten artificial sub-networks from Recon2. Those networks were generated to be approximately equally spaced in a range between 1000 and 3500 reactions. They were generated by randomly removing up to 4500 reactions from our Recon2 version and determining the consistent part of the remaining model. The first model within ±50 reactions of equally spaced points in the interval [1000, 3500] was selected as representative for this point. The models and model sizes can be found in **Supplementary File 5**.

# 2.2. Characterization of the Algorithms

There are many algorithms available for tissue-specific metabolic network reconstructions (see **Table 2**). In this section we will detail the algorithms used in our study and give reasons, why others were excluded.

In order to test the algorithms with real data, liver models were built by the tested algorithms using as input 22 arrays from different datasets downloaded from the Gene Expression Omnibus (GEO; Edgar et al., 2002) database (**Supplementary File 2**). The same data was also used for the cross-validation assays.

### 2.2.1. GIMME (Becker and Palsson, 2008) and iMAT (Zur et al., 2010)

For the benchmarking of the GIMME (Becker and Palsson, 2008) and the iMAT (Zur et al., 2010) algorithms, the implementation provided by the COBRA toolbox (Schellenberger et al., 2011) was used with an expression threshold corresponding to the 75th percentile. The proceedExp option was set to 1 as the data was preprocessed. For GIMME, the biomass objective coefficient was set to 10−<sup>4</sup> .

# 2.2.2. INIT (Agren et al., 2012)

In the original paper, INIT (Agren et al., 2012) assigns weights to the genes associated to the input model that were computed by dividing the gene expression in the tissue of interest by the average expression across all tissues. As for the first experiment, only liver arrays were available, z-scores obtained by the Barcode discretization method (Zilliox and Irizarry, 2007; McCall et al., 2011), were used as weights (see below).

## 2.2.3. RegrEx (Robaina Estévez and Nikoloski, 2015)

The RegrEx implementation in the supplementary files of Robaina Estévez and Nikoloski (2015) was used. This algorithm has previously only been used with RNA-seq data and therefore no established discretization method exist for microarray data. In order to allow a comparison with the others methods, the intensity values after frma normalization and the standard variation were directly mapped to the reactions of the model using the Gene-Protein-Reaction rules (GPR). For reactions that are not associated to any gene, the expression and the standard deviation were set to 0 and 1000, respectively.

# 2.2.4. Akesson (Åkesson et al., 2004)

For this algorithm, the data was normalized with the frma normalization method and then discretized with Barcode. Genes with z-scores below 0 in 90% of the arrays, were considered inactive and the bounds of the associated reactions, taking into account the Gene-Protein-Reaction rules (GPR), were set to 0. FASTCC (Vlassis et al., 2014) was then run to remove reactions that are unable to carry a flux.

# 2.2.5. FASTCORE z-score

For FASTCORE z-score, the expression data was normalized with the frma method and discretized using Barcode. Barcode uses previous knowledge on the intensity distribution across thousands of arrays to calculate for each probe set of the analysed array the number of standard deviations to the median of the intensity distribution for the same probe set in an unexpressed state. Genes with a z-score above 5 in 90% of arrays are considered as expressed and mapped to the reactions according to the Gene-Protein-Reaction rules (GPR) to obtain a core set that is fed into FASTCORE (Vlassis et al., 2014).

# 2.2.6. FASTCORMICS (Pacheco et al., 2015)

The expression values were first normalized with frma, converted into z-scores using Barcode (McCall et al., 2011) and further discretized using an expression threshold of 5 z-scores and an unexpression threshold of 0 z-score. Genes with 90% of the arrays above the expression threshold are assigned a score of 1 while those below the unexpression threshold are assigned a score of −1. All other genes are associated with a discretization score of 0. These scores are then mapped onto the model using the Gene-Protein-Reactions rules to obtain lists of core and unexpressed reactions. Unexpressed reactions are excluded from the model.

The FASTCORMICS workflow allows the inclusion of a medium composition, which was not used in the tests, as the aim was to provide the same information to all algorithms. A modified version of FASTCORE is then run that maximizes the inclusion of core reactions while penalizing the entry of non core reactions. Note that transporter reactions are excluded from the core set but are not penalized.

# 2.2.7. Context-Specific Reconstruction Algorithm that were not Tested

PRIME and tINIT were not included in the tests as they require, in addition to expression data, growth rates for PRIME and information on tissue functionalities for tINIT. Determination of growth rates in multicellular organisms is restricted to cell lines or cancerous cells, as most other cell types are finally differentiated and therefore no longer divide. Since growth rates are an essential part of PRIME it was excluded from the tests. While functionalities are available for some metabolically very active tissues (like kidney and liver), they are often not available for others. Since we wanted to test a wide range of potential tissues, we decided not to employ functionalities in our input set. Therefore, tINIT would be reduced to INIT as the remaining functionality is the same. Since we wanted to focus on gene expression data, which is currently the most readily available type of data, we did not add metabolomic information into our screens. GIM3E would need this type of information and was therefore not tested. Finally, MBA, Lee and mCADRE took more than 5 days for a single run on 2 cores of our cluster and where therefore not included.

# 2.2.8. Similarity of the Context-Specific Models and Algorithm-Related Bias

The similarity level between the context-specific models built by the tested algorithms was assessed by computing the Jaccard index between each pair of models. The matrix containing the Jaccard indices was then clustered using Euclidian distance. Further, for each context-specific model, the number of reactions found by only 1, 2 up to all of the methods was computed and represented as a stacked boxplot. The colored areas represent the different models built by the tested algorithms and for each bin the colored area is proportional to the number of shared reactions.

# 2.2.9. Sensitivity and Robustness Testing Using Artifical Data

While there are methods that take continuous expression measurements into account (Colijn et al., 2009; Lee et al., 2012, and reviewed in Machado and Herrgård, 2014), other methods require the user to define sets of reactions that are present (FASTCORE, MBA) or perform some form of discretization to determine the presence or absence of a gene or a reaction (Akesson, GIMME, iMAT, FASTCORMICS). The latter types of methods, using some form of presence/absence calls can be more rigorously tested for robustness, as a target model can be used to provide the present and absent genes/reactions.

We also tested these algorithms using the artificially created networks. The test was performed as follows:

The potential available information was defined as the sets of reactions present in each submodel and absent from each submodel. Based on this data different percentages of input information (50%, 60%, 70%, 80%, and 90%) were provided to the algorithms. The same random samples were provided to the tested algorithms to allow a further comparison between the algorithms (generating a total of 5000 models for each algorithm). To be able to use reaction data, we modified the implementation of the GIMME algorithm to allow the direct provision of the ExpressedRxns and UnExpressedRxns fields. The model similarities were assessed by calculating the Jaccard index between each pair of models generated for input sets from different target models. In addition, the internal distances of all models generated for one target model were calculated (a total of 50,000 comparisons per algorithm). Furthermore, the corresponding models for each algorithm and each tested input percentage were compared, to obtain the inter-algorithm distance.

#### 2.2.10. Robustness Testing Using Real Data

For the cross-validation, 20% of the reactions were removed from the core set and transferred to the validation set. The number of these reactions that were included in the output model was determined and a hypergeometic test was computed. The process was repeated 100 times randomizing at each iteration the core set to form different validation sets. For algorithms that take continuous data as input, the cross-validation assay was adapted as follows: 20% of the gene-associated reactions were removed from the input set by setting the expression to 0 and the standard deviation to 1000 for RegrEX and the rxnsScores to 0 for INIT. But only reactions considered to be expressed with a high confidence level formed the validation set i.e., for INIT only reaction with z-scores above 5 and with expression value above 10 for RegrEX. For Akesson the validation set was composed of inactive reactions. The results for Akesson have to be taken with care as the validation set is only composed of 4 reactions. This is due to Barcode only indicating very few genes as absent, which led to only about 40 reactions being removed from Recon2.

# 2.3. Benchmarking with Real Data

### 2.3.1. Confidence Level of the Reactions

The z-scores computed by Barcode give the number of standard deviations of a gene expression level above the mean of the same genes in an unexpressed state. The z-scores of the genes were mapped to the reactions of Recon2 (Thiele et al., 2013), HepatoNet (Gille et al., 2010) and to the context-specific models built by the different workflows using the Gene Protein Rules (GPR). In the same way, the confidence levels assigned by the Human Protein Atlas (HPA) to the proteins of the database were mapped to the reactions of the different context-specific models.

### 2.3.2. Comparison Between Different Tissue Models

The ability of the algorithm to capture metabolic variations among tissues was tested using the GSE2361 dataset (Ge et al., 2005) downloaded from Gene Expression Omnibus (GEO) that contains 36 types of normal human tissues. Twenty-one of the 36 tissues matched tissues in the Human Protein Atlas. The confidence levels of the proteins in the different tissues were first matched to the modified version of Recon2 to determine if proteins with high and medium confidence level are ubiquitously expressed or expressed in a more tissue specific manner. Then the confidence levels were matched to the corresponding contextspecific models to verify if the variation observed among the tissue context-specific models matched the one observed in the Human Protein Database.

To further access the quality of the reconstructed models, the fraction of reactions of the Recon2 pathways that are active in the output models were computed. The obtained matrix was then clustered in function of the Euclidean distance (see Supplementary Figure 6).

# 2.4. Benchmarking with Artificial Data

The runs using artificial data, performed for sensitivity and robustness analysis, were also used to provide an additional benchmarking measurement for the algorithms. Sensitivity, specificity and false discovery rate were calculated by comparison of the reconstructed networks with the respective target network. The artificial nature of these networks allowed us a complete knowledge of the actual target thus making these calculations possible.

# 2.5. Network Functionality Testing

Function testing is commonly achieved, by defining a set of metabolites that are available and can be excreted and requiring other metabolites to be produced/consumed or a reaction to be able to carry flux. The input and output can either be cast into a linear problem by adding importers and exporters or by relaxing the steady state requirement for the imported and exported metabolites. Gille et al. (2010) used the latter definition and we adapted this approach using the following modification of the standard FBA approach:

$$\begin{aligned} \min & \quad \sum \boldsymbol{\nu}\_{i}^{+} + \boldsymbol{\nu}\_{i}^{-} \\ \text{s.t. } & b\_{l} \leq \boldsymbol{S}' \ast \boldsymbol{\nu}' \leq b\_{u} \\ & \boldsymbol{0} \leq \boldsymbol{\nu}\_{i}^{+} \leq \boldsymbol{u}b\_{i} \quad \forall i \in internal\ reactions \\ & \boldsymbol{0} \leq \boldsymbol{\nu}\_{i}^{-} \leq -lb\_{i} \quad \forall i \in internal\ reactions \\ & \boldsymbol{\nu}\_{i}^{+} - \boldsymbol{\nu}\_{i}^{-} = \boldsymbol{0} \quad \forall i \in exchange\ reactions \end{aligned}$$

$$\text{with } S' = [\text{S, } -\text{S}] \\ \text{and } \nu' = \begin{bmatrix} \nu^+ \\ \nu^- \end{bmatrix}$$

$$b\_{l,i} = \begin{cases} -10000 & \forall i \in improved\text{ }metadata(-) = \text{ } \\ -1 & \forall i \in produced\text{ }objects\text{ }i \text{ }s \text{ } \\ 1 & \forall i \in consumed\text{ }object\text{ }i \text{ }s \text{ } \\ 0 & \text{ }else \end{cases}$$

and bu,<sup>i</sup> = 10000 ∀i ∈ exported metabolites(+/ =) −1 ∀i ∈ produced objectives(+) 1 ∀i ∈ consumed objectives(−) 0 else

The test is considered to be successful if there is a non zero value for all evaluators when calculating S ′ · v ′ .

# 2.6. Computational Resources

Except for RegrEx, all runs using the liver data were performed on two cores of a 2.26Ghz Xeon L5640 processor on the HPC system of the University of Luxembourg (Varrette et al., 2014) to achieve comparable running times. Tissue comparison runs and artificial simulation runs were performed on the same cluster but not limited to specific node types.

# 3. RESULTS

# 3.1. Characterization of the Algorithms

## 3.1.1. Similarity of the Context-Specific Models and Algorithm-Related Bias

The aim of this characterization step is to categorize the algorithms based on the similarity of their output models in order to gain insight into algorithm-related bias, requirements of the algorithms i.e., thresholds and more importantly when to use which algorithms. In an ideal case, one would expect that when fed with the same input data, the different algorithms would produce similar networks. But when comparing the contextspecific liver models generated with the different algorithms and HepatoNet, only 530 reactions were found in all networks and 77 reactions of our version of Recon2 were inactive in all context-specific models and HepatoNet. The 530 reactions were found among 54 different subsystems, including reactions belonging to pathways expected in all tissues like i.e., the Krebs cycle, glycolysis/gluconeogenesis, but also pathways that were described to take place mainly in the liver, like i.e., bile acid synthesis (Rosenthal and Glew, 2009; Wang et al., 2012) or some reactions of the vitamin B6 pathway (pyridoxamine kinase, pyridoxamine 5′ -phosphate oxidase and pyridoxamine 5 ′ -phosphate oxidase; Merrill et al., 1984). This huge variability is due to workflow-related bias and to different strategies and aims of the algorithms. FASTCORE (Vlassis et al., 2014), expects as input a set of reactions with a high confidence level which are assumed to be active in the context of interest and therefore all core reactions are included in the output model (**Table 3**). In contrast, FASTCORMICS (Pacheco et al., 2015) only includes a core reaction if it does not require the activation of reactions with low z-scores. The main objective of GIMME (Becker and Palsson, 2008) is to build a model by maximizing a biological function. The input expression data is used to identify, which reactions are not required for the objective and can function therefore be removed from the model due to low expression values (**Table 3**). iMAT (Zur et al., 2010; Lee et al., 2012) and RegrEx (Robaina Estévez and Nikoloski, 2015) maximize the consistency between the flux and the expression discarding reactions that have high expression values if necessary, which might be problematic if reactions have to be included in the model like i.e., the biomass function. INIT (Agren et al., 2012) uses weighted activity indicators as objective, with those having stronger evidence being weighted higher. Whereas the Akesson's (Åkesson et al., 2004) algorithm aims to eliminate non expressed reactions.

The models, when clustered in function of the Jaccard Similarity Index (**Figure 1**), form 2 branches and an outlier: HepatoNet. The first cluster is composed of algorithms that TABLE 3 | Models numerics: Size, number of input reactions with high expression, respectively z-score levels, fractions of input reactions set included in the output models, number of genes-associated reactions in the model and running time.


\**Note that RegrEx was run on a different computer with an Intel(R)Xeon(R)CPU E3 1241-v3 @ 3.50 GHz processor.*

take as input continuous data and attempt to maximize the consistency between the data and the Akesson algorithm that eliminates inactive reactions. The second cluster is composed of algorithms that discretize the data in expressed and nonexpressed genes. Among this cluster, a second subdivision is observed between the algorithms that used z-score converted data (i.e., FASTCORE z-score and FASTCORMICS) and the ones that use normalized data without further transformation.

Overall the highest similarity level are found between FASTCORE z-score and FASTCORMICS with a score of 85% of similarity followed by iMAT and GIMME with 77% of similarity. The lowest similarity level is found between FASTCORMICS and HepatoNet with only 26% of overlap. The largest overlap between HepatoNet and context-specific reconstructions is found for INIT with 43% of similarity. Note that the INIT model although having as input Barcode discretized data does not cluster with FASTCORE z-score or with the FASTCORMICS models but with RegrEx, suggesting that the choice to consider continuous data rather than defined core set has a larger impact on the output models.

As the algorithms were fed with the same input data, reactions that are predicted by one or only few algorithms are more likely to be algorithm-related bias (**Figure 2**). The Akesson model that contains 98.56% of the input model includes the largest number of reactions (201) that are absent in the others models.

The reactions included in the FASTCORE, FASTCORMICS, iMAT and GIMME models are for 97%, 98%, 96%, respectively 89% supported by at least 3 other algorithms and display a similar profile shifted to the right. HepatoNet, INIT and the Akesson's model share 92%, 83%, respectively 91% with 3 other algorithms and have different profiles from the algorithms of the first group composed of algorithm that include a discretization step.

In summary, discretization-based algorithms show the highest similarity level and therefore the lowest number of reactions due to potential algorithm-related bias.

# 3.1.2. Sensitivity and Robustness Testing Using Artifical Data

Since we noticed that there are two sets of algorithms among the discretizing algorithms, we decided to further test their properties

with artificial networks by comparing resulting models from multiple runs for different models and levels of completeness of input data.

**Figure 3** provides the average similarities for all models reconstructed for each target model at different available information percentages (A full set of mean similarities for each percentage and each artificial model along with the data for the plots is provided in **Supplementary Files 1**, **3**). Each square represents the mean Jaccard index of the all combinations of networks generated for different input networks [e.g., (1,2) is the average similarity of all networks generated for models 1 to all networks generated for model 2]. The diagonal represents the internal similarity of all networks generated for one model. When 90% of the data is available, all the algorithms are able to distinguish variation between the different models. But with a less complete data set, inclusive algorithms lose in specificity and therefore also progressively lose the capacity to distinguish between different models. Further with 30 and 50% reactions missing, it would be expected that the algorithms get less robust, but Akesson and GIMME only show a modest decrease of robustness (as shown in the diagonal). A similar behavior for the GIMME algorithm was also described by Machado and Herrgård (2014) in a experiment where noise was progressively added to the input data to finally obtain a random input dataset. GIMME showed the same average error in prediction for the random and original data (Machado and Herrgård, 2014), suggesting that due to the optimization of the biomass function, the expression data has a reduced impact on the model building. Comparing the models resulting from runs with different completeness

of input data illustrates that the methods tend to converge on more complete data sets, with the Akesson approach and GIMME being more inclusive and the FASTCORE family being more exclusive (see **Figure 4** and **Supplementary File 4**). While initially, with incomplete data, the methods are distinguishable by the networks generated, this difference becomes smaller with additional knowledge.

### 3.1.3. Robustness Testing Using Real Data

In order to further evaluate the confidence level of the reactions included in the different context-specific models a 5-fold crossvalidation was performed. The experiment was repeated 100 times with a different validation set. GIMME, iMAT, and FASTCORMICS show the highest robustness, followed by FASTCORE and FASTCORE z-score (see **Table 4**). Algorithms that maximize the consistency between the data and the flux, e.g., INIT and RegrEx, are less robust with insignificant pvalue. For Akesson no hyper-geometric test was performed as the validation set was too small to obtain a reliable p-value. Note that for context-specific reconstruction algorithms a tradeoff has to be found between robustness and the capacity to capture differences between similar contexts. For this reason, a too high robustness might not be desirable as it would imply that the algorithm might lose in resolution power, i.e., the ability to distinguish between different sets of input data. Therefore, it is also advisable to not test for robustness without testing the resolution power.

# 3.2. Benchmarking with Real Data 3.2.1. Confidence Level of the Reactions Included in the Different Models

As shown by the previous similarity test, there are several alternative approaches to build context-specific models. To assess the confidence level of a reconstruction, one can quantify the confidence level of the reactions included by each algorithm. Context-specific algorithms assume that the higher the reactions associated expression levels, the more likely the reactions to be active. Following this logic, contextspecific reconstructions should be enriched for higher expression levels. As the background level is non negligible and highly dependent on the probes, we corrected for probe effect using the Barcode method. The z-scores computed by Barcode translate the number of standard deviations to the intensity distribution of the same genes in an unexpressed state. The z-scores of the genes mapped to the reactions of Recon2 (Thiele et al., 2013), HepatoNet (Gille et al., 2010) and to the context-specific models


built by the different algorithms show that the distribution of the z-scores are for most models shifted, as expected, toward higher z-scores values with a significant p-value for all contextspecific models except RegrEX (Robaina Estévez and Nikoloski, 2015). Algorithms that use a discretization method show a larger shift to the right than algorithms that maximize the consistency between the flux and the data. Within this group the FASTCORMICS (Pacheco et al., 2015) shows the most significant shift toward the highest z-score values followed by FASTCORE z-score, GIMME (Becker and Palsson, 2008), and iMAT (Zur et al., 2010) (**Figure 5** and **Table 5**). Surprisingly, the consistent version of HepatoNet (Gille et al., 2010) is associated to slightly higher z-scores than Recon2 (Thiele et al., 2013) but significantly

(Fastcore/FASTCORMICS) algorithms is clearly visible.

lower than most discretization based automated context-specific reconstructions.

Further, unlike their competitors, all the discretization-based context-specific reconstructions show an enrichment of genes with a high and medium confidence scores to be expressed at the protein level (Uhlén et al., 2015). A stronger enrichment is observed for FASTCORE z-score and FASTCORMICS with 46 and 50% of the gene associated reactions having a high or medium confidence level **Table 6**, respectively. GIMME and iMAT include 28 and 30% reaction with high or medium confidence levels, respectively. Again surprisingly, HepatoNet does not show an enrichment for high and medium confidence levels.

FASTCORMICS) are enriched for higher z-score values.

TABLE 5 | Comparison between the z-score distribution associated to the models build by the different methods.


*A low p-value indicates that the z-score distribution of Model 2 is shifted toward higher values (to the right) compared to Model 1 (Kolmogorov-Smirnov test).*

In summary, dicretization-based algorithms include reactions with a higher confidence level at the transcriptomic and proteomic level than their competitors.

### 3.2.2. Comparison Between Different Tissue Models

The aim of a context-specific algorithm, as indicated by the name, is to build models that capture the metabolism of a cell for a given context and therefore these algorithms have to be able to capture variations in the metabolism of different tissues. To pass the following test, context-specific algorithms not only have to be sensitive (or to have a high resolution power) in order capture metabolic difference between tissues, but the reconstructions for different tissues have to be enriched for high or medium confidence levels based on HPA. The last criteria allows to identify algorithms that build different models based on noise or algorithm-related bias. In order to assess the variation among tissues in HPA, the genes with high, medium and low confidence levels for 48 different tissues were mapped to the input model Recon2, showing that very few reactions have a high or medium confidence level in all tissues. In summary, most reactions with high and medium confidence scores have a more tissue-specific expression (**Figure 6**).

A similar experiment was performed with context-specific reconstructions built by the tested algorithms, in which the number of algorithms that shared a reactions was assessed (see **Figure 7**). For RegrEX, INIT and Akesson models, the majority of reactions are found in all tissues. For GIMME, most reactions are either tissue-specific or present in all the tissues. In contrast, the models built by the members of the FASTCORE family show a distribution similar to the that obtained in **Figure 6**, for HPA. For iMAT only 8 models could be obtained as the computational demands for the reconstructions of the others tissues surpasses the number of core available and the maximal running of 5 days. When looking at the confidence levels associated with the 21 different tissue-specific models, FASTCORE z-score and FASTCORMICS show in 20 out of 21 the highest percentage of reactions with a high or medium confidence level (see **Figure 8**). The size of the different tissue metabolic models built by the tested algorithm can be found in the **Supplementary File 6**).

The quality of the tissue-specific models built by the different algorithm were accessed by focusing on selected pathways known to have a more tissue-specific expression, namely bile acid synthesis and heme synthesis. The bile acid synthesis occurs in liver, although one or the other enzyme of the pathways might occasionally be expressed by other tissues (Rosenthal and Glew, 2009; Wang et al., 2012). As expected the FASTCORE family, GIMME and iMAT predicted that the highest fraction of active reactions are found in the liver followed by the foetal liver for the FASTCORE family members and iMAT and by placenta and foetal liver for GIMME. Whereas, the INIT models of skin, bone marrow, corpus, thalamus, pituitary gland and foetal liver had a higher fraction of active reactions than the liver model. Thirteen out of thirty-six of the tested Akesson models predicted 90% and more reactions of the bile acid pathway as active. RegrEX predicted a slightly higher fraction in the thalamus than in the liver and a comparable fraction in the ovary, the foetal brain and the corpus (**Supplementary Files 1**, **6**).

The heme synthesis that occurs mainly in the developing erythrocytes and in the liver (Ajioka et al., 2006), was given as 100% active by the FASTCORE family and completely inactive by GIMME and iMAT in the liver. But these two algorithms predicted the pathway to be active in other tissues. As a matter fact, all the algorithms predicted the pathway to be active in others tissues than the liver. INIT, RegrEX and Akesson included this pathway in 20, 22, and all tested 36 tissues, respectively. Fewer models of the FASTCORE family contained reactions of this pathway: uterus and tyroid for FASTCORMICS and spleen, placenta, uterus, thyroid, skin, bone marrow, amygdala, lung and foetal liver for FASTCORE.


TABLE 6 | Number, percentage of gene-associated reactions and percentage of reactions of each context-specific reconstruction that have a high, medium and low confidence score to be expressed at the protein level.

*An enrichment in high and medium confidence level is observed for discretization-based algorithms (GIMME, iMAT, FASTCORE z-score and FASTCORMICS).*

# 3.3. Benchmarking with Artificial Data

To further evaluate the quality of the algorithms, we also used the artificial data (see Section 3.1.2) to benchmark the algorithms. Comparing the resulting models with the target models, we again see that for more complete input sets, the model quality tends to become more similar (see **Figure 9**). It is interesting to note that the false discovery rate (FDR) of FASTCORE for higher percentages is similar to those of the inclusive models, while FASTCORMICS achieves a better FDR. This indicates alternative routes to activate reactions. In general, there is again the tradeoff between adding too much or too little. It is however interesting that the exclusive algorithms tend to miss targets and their sensitivity is independent on the size of the target model while this is different on inclusive algorithms. Exclusive algorithms show a better FDR than inclusive algorithms. Further, for smaller target models, the loss in precision of inclusive algorithms (1-FDR) is more pronounced for 50 and 70% of the input data, as the inclusive algorithms tend to overestimate the actual model. Similar to the previous experiment, it would be expected that the sensitivity (robustness) would decrease with an increased percentage of missing data. But the inclusive algorithms show an invariant sensitivity in function of the available data suggesting that the expression data has reduced impact on the model building. The specificity for the exclusive algorithms is independent of the target model size and are less affected by the increased missing data than the inclusive algorithms. The sizes of the different reconstructed models also indicates the trend for convergence, and a figure showing the converging sizes is provided in **Supplementary File 1**.

# 3.4. Functionality Testing

Functional testing allows us to assess which known functions of a specific tissue are captured by a reconstruction. We used the set of functions defined in HepatoNet and formalized in Section 2.5 for the liver and tested them on all reconstructed networks. We noticed that the success rate of HepatoNet and the generic

with a high or medium confidence level that are shared between 1, 2, 3, up to 48 tissues of the Human Protein Atlas. Reactions with a high confidence level tend to have a tissue-specific expression.

reconstruction Recon2 are comparable with 244 vs. 247 of 310 network tasks and 109 vs. 98 of 123 physiological tasks for Recon2 and HepatoNet, respectively. The discrepancy with the original publication is likely due to alternative solutions and we noticed that HepatoNet allows free production of NADH and thereby ATP (see Table 2 in **Supplementary File 1**). The discrepancy between the consistent and inconsistent HepatoNet is due to the formulation of the functionalities, which do not require exchange reactions but modify the b vector, thus generating implicit importers and exporters and allowing inconsistent parts of the network to carry flux. We also noticed an important issue with functional testing: For random models, the larger the models, the higher the functionality score (with R <sup>2</sup> = 0.869 and 0.915 for network and physiological functions, respectively). To illustrate this issue, we generated 400 random networks by removing a random number of up to 2000 reactions from the consistent part of Recon2 and subsequently removing all reactions which could no longer carry any flux. We then tested all network and physiological functions on these networks. The results can be seen in **Figure 10**, for both the network and physiological tests.

Blue circles represent the random networks; the consistent HepatoNet and the original HepatoNet are displayed in orange, and show a strong enrichment in functionalities. The higher number of functionalities covered in HepatoNet stems from several reactions which are inconsistent, but can be used in a functional testing as described above. We also marked the models generated using the GEO dataset for liver, which score similar to equally sized random models. One of the main reasons for the strong correlation between model size and successful tests

computational complexity of iMAT it was only possible to generate 14 out of 36 tissue models.

model/input data. The model sizes are: Model 1: 961, Model 4: 1876, Model 7:2629, Model 10: 3455 While the quality of the FASTCORE models is independent of the target model size, the inclusive approaches tend to largely overestimate smaller models, when insufficient data is available. A plot with all models can be found in Supplementary File 1.

is the amount of "positive" testing. Many tests are concerned with some type of biosynthesis or degradation and a larger model is more likely to be able to fulfil these requirements than a smaller model. But even using e.g., the biomass function (like GIMME) as part of the input, the models do not get significantly better than a random model on expression data for liver. None of the algorithms tested achieves high scores in the functionality test and several algorithms are on the lower end of the random network reference. A plot showing the tests passed by the different algorithms is supplied in **Supplementary File 7**. tINIT could potentially surpass most other algorithms on this test, as it includes functionality information in its reconstruction routine. However, the formulation of tINIT functions is again slightly different from the formulation in HepatoNet and thus not directly compatible.

# 4. DISCUSSION

The primary aim of this work was to review and discuss the existing validation methods and to propose a unified benchmark for the assessment of context-specific reconstruction algorithms. This benchmark will help to identify potential deficiencies of existing and new algorithms and by such increase the quality of context-specific reconstruction algorithms and the models they generate. Although the tested algorithms were validated by their authors in order to be published, the validation methods applied are often incomplete, e.g., a particular aspect of the output model fitting the context of the paper is tested like the ability to produce lactate from glucose in cancer models, leaving other pathways unconsidered. Further, discretization thresholds and other free parameters of the algorithms are likely to be set to optimally fit a particular dataset. Thus, when used in another context the algorithm might perform worse than expected from the original publication. The need of a unified benchmark is nicely illustrated by **Figure 1** which shows that despite being fed with the same inputs, the output models vary considerably from each other e.g., the output models of RegExp and FASTCORE that share only around 30% of the reactions.

Part of the variance between the output models is due to different aims and philosophies of the tested algorithms but also due to algorithm-related bias. The second aim of this work was to demonstrate to the users that the context-specific reconstruction algorithms are not equivalent and that the choice of the algorithm and selection of parameter settings for the algorithms have to be performed with care respecting the philosophy of the tested algorithm. For example, GIMME maximizes a chosen biological function and when using GIMME the user assumes that the metabolism of a cell is aimed at the fulfilment of this function. While this biological function can be assumed to be growth for many microorganisms or cancer cells, it is likely to be more complex for multicellular organisms, where multiple "objectives" have to be balanced. In the same way, FASTCORE takes as input core reactions that are always included in the output model and therefore a higher threshold corresponding to a higher confidence level should be set when using FASTCORE.

Although the parameters were set according to the original papers, we are aware that some of the tested algorithms might perform better with a different parameter setting. We decided nevertheless when possible not to change the original parameter settings of the algorithm. First, because the main objective of this paper is not to assess existing algorithms but to propose a benchmark to validate context-specific algorithms. Second the finding of the optimal parameter setting is a computational demanding processes that would require i.e., crossvalidations or other criteria that are not always available. Finding the optimal parameter setting is beyond the scope of a benchmark and rises other questions like overfitting to the data. Third, algorithms should be sufficiently robust to be applied to other datasets with the optimal settings as defined by the authors. As a general principle, in order to avoid overfitting, the parameter estimation should not be performed on the same data than the one used for model generation. We therefore encourage the authors and the users of these algorithms to test them with others parameter settings that might be more appropriate.

The benchmark that we suggest and for which we provide the scripts (http://systemsbiology.uni.lu/software) is based on several criteras:

First of all the algorithms have to produce models of high quality that include genes or reactions that are supported by some evidence to be expressed in the context of interest. This aspect was assessed in the workflow by mapping Barcode z-scored gene information and confidence levels established by the Human Protein Atlas to the models. Context-specific reconstruction that extract sub- networks composed only of active reactions in the context of interest from a general reconstruction tend to produce output models that are enriched for genes with high z-scores and a high confidence level to be expressed at the protein level. Indeed although the activity does not correlate perfectly with expression intensities, it was shown that algorithms that exclude reactions with low expression values show a better predictive power than the generic models from which they were extracted. Both tests show that algorithms that perform a discretization of the input data perform better in these tests than algorithms that maximize the consistency between flux values and the data.

We noticed that within the discretizing algorithms, there are two conceptually distinct approaches when considering unsupported reactions. An inclusive concept which considers unknown data as present and an exclusive concept that considers unknown data as absent. Inclusive concepts tend to produce larger networks and score lower, when comparing the networks to additional data, while exclusive concepts tend to produce smaller networks and score higher.

This can be considered as algorithm related bias and it is likely that when multiple algorithms are supplied with the same inputs, reactions that are found by only one or only few algorithms are more likely to be due to algorithm-related bias. Algorithm related bias is not negligible as shown by the huge variability of liver reconstructions with e.g., up to 30% of the reactions being different between the FASTCORE and RegrExp algorithm (**Figure 1**).

Further, algorithms have to be robust to noise but nevertheless be precise enough to capture the variations in the metabolism of a cell in different contexts i.e., different cell types, different states e.g., healthy vs. disease and eventually between different patients. These two criteria were tested using both experimental and artificial data. Algorithms like GIMME are performing extremely well in the cross-validation assay but score low in the tissue comparison test, as GIMME produces quite similar reconstructions for the different tissues tested. The algorithms using an inclusive concept tend to be more robust to noisy data but have a reduced resolution power. In contrast, exclusive algorithm are less robust as they tend to only recover reactions that are supported by the input data or reactions that are needed to obtain a consistent model, which allow a greater resolution power. Therefore, among the tested algorithms, the FASTCORE family capture best the variation between the different tissues. Further, the confidence level of the reactions included in the 21 tissue models showed that the variability captured by the FASTCORE family models, was not due to noise or algorithm related bias. In the same aspect, the artificial model test gave some interesting insight into the quality of the reconstruction algorithms. While both groups of algorithms, including and excluding, generated about the same model when perfect information was available, they start to diverge at lower amounts of available data. In particular, with less information available the exclusive algorithms underestimate the target network and the including ones overestimate it. While this is to be expected it indicates that the use of two algorithms can give a good approximation of the quality of the available input data and completeness of the reconstruction. If both types of algorithms (inclusive and exclusive) do diverge substantially, it is likely that a relevant amount of input information is missing and that the "true" model is somewhere in between. Similarly, if the models are almost identical, it is likely that the input information and the reconstruction quality is high. GIMME will always include the objective function and all reactions necessary for this function to carry flux. Therefore, those reactions might influence the network size considerably. One advantage of an exclusive concept in this respect, is that its variability is less target model dependent than an inclusive approach. For smaller models, the FDR for inclusive models tends to rise much more rapidly with a more incomplete input data set than for larger models. As we commonly are unaware of the actual size of the target network, this might cause problems when using inclusive approaches.

Another important aspect is the computational demand. To determine the processing time we decided when possible not to change the solver used in the original paper as we noticed that algorithms like e.g., RegrEX are sensitive to the used solver, with gurobi finding an initial solution guess faster than e.g., cplex and thus the result returned by cplex being unusable for the algorithm. The range of computational times is however substantial, with fast algorithms running in seconds to minutes and others taking hours or even days. One of the greatest advantages of faster algorithms, is their capability to be more thoroughly evaluated using cross-validation techniques, which is infeasible for an algorithm running several days. We also observed an issue when running the INIT algorithm. For unknown reasons, the algorithm consistently stopped after 10 h of computation. In particular, the resulting models were odd at best, as they should be close to the models generated by FASTCORE, and in the artificial test, should be optimal on optimal inputs. However, the artificial test was far from optimal, and we assume that the solver does terminate computation at some point.

Finally, we also assessed the capacity of the contextspecific reconstruction to pass the functional test as established in Gille et al. (2010). We found that no algorithm outperforms random models, but that a fitted model can indeed show higher scores without adding more reactions, as seen in **Figure 10**. Unfortunately, obtaining functional data is a very time consuming process and necessitates intensive literature research every time a new tissue model is created. The failure of the tested algorithms in the functional test is mainly due to the high number of non-gene associated reactions in the generic input model (one third of Recon2) and due to the reactions associated to genes with low expression levels. The tested algorithms extract a sub-network from the input model that includes all or most reactions associated with high expressions levels (core) and few reactions with low expression levels (non-core) in order to obtain a consistent model. A slightly different core reactions set, can cause the core reactions to be connected in a different way and as a result the model displays different functionalities. As the choice of the non-core reactions is to a large extent not guided by the data, the obtained functions are random as shown by the functionality test. Interestingly, the reactions found in HepatoNet do have weak evidence when compared to HPA or z-scores, which partially provides another explanation for the inability of the tested algorithms to recover these activities. This however indicates that the general reconstruction currently used lacks either the correct gene-protein-reaction associations for several reactions necessary for the functionalities in liver, that there are alternative pathways missing in the reconstruction and the reactions used in HepatoNet are not the "true" reactions, that the functions are incorrectly assumed to be available in liver or that the functionality lacks information about the consumed cofactors. Indeed, as all the exchange reactions are closed, some reactions might not carry a flux as the associated cofactor cannot be regenerated. This would also explain why bigger models accumulate more functions. The larger the models, the higher the likelihood of internal loops that could allow a regeneration of cofactors. Further it might also indicate that transcriptomics alone might not be sufficient to build functionally correct models. Information on the uptake and excreted metabolite added to the input reactions set would probably increase the score of most algorithms. We did nevertheless not include this type of information in the input data as the latter is not available for in vivo tissues. While presence of importers and exporters does not influence the functional tests, they are however highly influenced by the availability of internal transporters.

Assuming that the defined functions are indeed present in liver, this would indicate the importance of algorithms like tINIT which do take these functionalities into account and which could, given the right reference network, indicate potential missing links in the current reconstructions. tINIT is nevertheless not able to capture metabolic differences between different tissue as shown in Uhlén et al. (2015), calling for a new generation of algorithms that capture metabolic variation and that are able to take as input functionalities. Note here that algorithms like PRIME that do not extract a subnetwork to obtain a contextspecific model, but modifies the bounds of the reactions of the input model, will have regardless of the modeled celltype or context the same functionalities as the input model. Therefore, PRIME would score as high as the generic Recon2 in a qualitative test. Nevertheless, the approach used by PRIME is extremely dependant on the accuracy of the growth measurement and biomass formulation, leading to a very variable quality of the flux prediction (Yikzah et al, 2014). In a quantitative test aiming to predict the production rate of lactate by cancer cells, PRIME showed a lower correlation to the experimental data than FASTCORMICS (Pacheco et al., 2015). This suggests that building context-specific algorithms with the discretizationbased algorithms and then constraining the uptakes rates of several key amino-acids and glucose as performed in Pacheco et al. (2015) seems to be favorable. Further, as discussed in the main text, there is no unique function to which the metabolism of a non-cancerous pluricellular cell could be reduced and sofar is limited to handle one metabolic function.

In general, we would recommend to assess the quality of an algorithm based on a combination of functional tests for a reconstructed tissue always in comparison to random networks, confirmation using an independent source of information (e.g., proteomics data, when only using expression data for the reconstruction), and an assessment of algorithmic properties, like dependence on target or input model size and dependence on input data quality. For the latter we would suggest using artificial networks to provide a complete knowledge on the expected outcome.

# AUTHOR CONTRIBUTIONS

MP, TP, and TS designed the study. MP and TP implemented the validation methods and performed the calculations. MP, TP, and TS wrote the manuscript.

# ACKNOWLEDGMENTS

TS and TP are funded by the Life Science Research Unit, University of Luxembourg. MP was supported by a fellowship from the Fond National de la Recherche Luxembourg (AFR 6041230). We would like to thank the HPC team of the University of Luxembourg for their support.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys. 2015.00410

Supplementary File 1 | Supplementary tables and figures. Modifications of Recon to allow the reconstruction of hepatonet. Example of free NADPH production in hepatonet.

Supplementary File 2 | List of 22 microarrays used to build the liver models.

# REFERENCES


Supplementary File 3 | Data underlying the similarity plots of Figure 3.

Supplementary File 4 | Data underlying the similarity plots of Figure 4.

Supplementary File 5 | Artificial Models and used Recon and HepatoNet Models in SBML format.

Supplementary File 6 | The sizes of the models of the different tissues built by the tested algorithms and the faction of active reactions in each pathway for the different tissues.

Supplementary File 7 | Successes of the models in the functional tests based on the tests from HepatoNet.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Pacheco, Pfau and Sauter. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Applications of Genome-Scale Metabolic Models in Biotechnology and Systems Medicine

Cheng Zhang<sup>1</sup> and Qiang Hua1, 2 \*

*<sup>1</sup> State Key Laboratory of Bioreactor Engineering, East China University of Science and Technology, Shanghai, China, <sup>2</sup> Shanghai Collaborative Innovation Center for Biomanufacturing Technology, Shanghai, China*

Genome-scale metabolic models (GEMs) have become a popular tool for systems biology, and they have been used in many fields such as industrial biotechnology and systems medicine. Since more and more studies are being conducted using GEMs, they have recently received considerable attention. In this review, we introduce the basic concept of GEMs and provide an overview of their applications in biotechnology, systems medicine, and some other fields. In addition, we describe the general principle of the applications and analyses built on GEMs. The purpose of this review is to introduce the application of GEMs in biological analysis and to promote its wider use by biologists.

Edited by:

*Xiaogang Wu, Institute for Systems Biology, USA*

#### Reviewed by:

*Stephen Fong, Virginia Commonwealth University, USA Kazuyuki Shimizu, Keio University, Japan Hiroshi Shimizu, Osaka University, Japan*

> \*Correspondence: *Qiang Hua qhua@ecust.edu.cn*

#### Specialty section:

*This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology*

Received: *29 October 2015* Accepted: *15 December 2015* Published: *07 January 2016*

#### Citation:

*Zhang C and Hua Q (2016) Applications of Genome-Scale Metabolic Models in Biotechnology and Systems Medicine. Front. Physiol. 6:413. doi: 10.3389/fphys.2015.00413* Keywords: genome-scale metabolic models, systems biology, metabolic capability analysis, in silico metabolic engineering, systems medicine

# INTRODUCTION

Genome-scale metabolic models (GEMs) are reconstructions of the metabolic networks of many kinds of cells, including those of microorganisms, plants, and mammals. In some cases, GEMs could represent the whole tissue or body of a multicellular organism. In these metabolic networks, the gene-protein-reaction (GPR) relationships are annotated. In addition, all the reactions in GEMs are mass- and energy-balanced, ensuring stoichiometric balance. Thus, GEMs enable researchers to conduct system-level metabolic response analysis and flux simulation, which is not possible using general metabolic pathway databases such as KEGG. Furthermore, since GPR relationships are included in GEMs, other omics data such as transcriptomic and proteomic data could be systematically integrated into GEMs. Thus, GEM-based multi-omic analyses are more informative with stoichiometric balance and could possibly provide deeper biological insights.

In the past 15 years, GEMs have garnered considerable research attention. In 2000, the first GEM, a model of Escherichia coli MG1655, was reported (Edwards and Palsson, 2000). A few years later, a yeast GEM was published (Doerks et al., 2002), thus initiating a new era for systems biology. In the beginning, researchers tried to use GEM-based in silico simulations to guide the rational design of industrial microorganisms (hereafter referred to as in silico metabolic engineering). In 2003, a method called OptKnock (Burgard et al., 2003) was published and it employed a bi-level optimization program to search for reaction knockout targets that would yield overproduction of a desired biochemical while maintaining optimal growth. Following that, a series of in silico metabolic engineering methods were developed for various gene manipulations other than knockout (Pharkya et al., 2004; Pharkya and Maranas, 2006; Choi et al., 2010; Ranganathan et al., 2010; Park et al., 2012; Chowdhury et al., 2014; Mahalik et al., 2014), leading to a marked expansion in the usage of GEMs. Furthermore, many of the in silico metabolic engineering methods were experimentally validated (Fong et al., 2005; Izallalen et al., 2008; Asadollahi et al., 2009; Brochado et al., 2010; Choi et al., 2010; Yim et al., 2011; Xu et al., 2011; Park et al., 2012; Ranganathan et al., 2012; Otero et al., 2013; Kim et al., 2014), which showed the power of GEMbased applications. With the development of systems biology, GEMs were also used as scaffolds for systematic integration of omics data because GEMs could be used to reconstruct the relationship among genes, enzymes, and metabolism. Numerous algorithms have been developed to integrate various types of omics data such as thermodynamics (Henry et al., 2007), transcriptomics/proteomics (Becker and Palsson, 2008; Colijn et al., 2009; Zur et al., 2010), fluxomics (Wiback et al., 2004), and metabolomics (Cakir et al., 2006). In return, the integration of omics data could improve the prediction of GEMs. More recently, GEM has been applied to systems medicine. Since the reconstruction of the first global GEM for humans, Recon 1, which was established in 2007 (Duarte et al., 2007), researchers have started to explore the possibility of clinical applications of GEMs and have reported several successful cases (Agren et al., 2014; Gatto et al., 2014; Jerby-Arnon et al., 2014). In fact, GEMs and their applications have received considerable attention recently.

Although GEMs are becoming increasingly popular, they are not easy to understand or use by non-experts. The complex code and script usually used for GEM-based computational applications and analyses are not readily available to the community of biologists, greatly hampering the wide usage of GEMs. In this review, we describe the key concepts and assumptions of GEMs. In addition, we describe the general principle of the applications and analyses built on GEMs. The information presented here is expected to promote the spread of GEM usage by biologists.

# BASIC CONCEPT OF GEMs

As mentioned above, GEMs are metabolic networks. **Figure 1A** shows a partly visualized glycolysis pathway in a GEM of E. coli, and within this part, we can see that metabolites are linked with each other by reactions, which are associated with enzymes, which are encoded by genes. It should be noted that the stoichiometric coefficient in metabolic reactions in **Figure 1A** (as shown in **Figure 1B**) could not be visualized in a graph. Therefore, GEMs employ a stoichiometric matrix (S matrix) to represent all the coefficients in metabolic reactions (**Figure 1C**). In the S matrix, the ijth element represents the stoichiometric coefficient of the ith metabolite in the jth reaction in the GEM. If the coefficient is positive, the metabolite is produced; otherwise, it's consumed. In addition, the GPR relationships in GEMs are simplified into a two-dimensional binary matrix showing the association between genes and reactions (**Figure 1D**), in which the ijth element is one if the ith reaction is associated with the jth gene, and it's zero if they aren't associated.

GEMs have several notable features: (1) They are collections of existing knowledge of the metabolism of a specific organism, and in most GEM-based applications, it's assumed that the metabolic network is complete, with very few exceptions, such as for gap finding and gap filling (Latendresse, 2014). (2) They are stoichiometric-balanced networks, which means mass as well as energy balance, reduction, and proton balance are well considered. (3) GPR relationships are annotated in GEMs, but the interactions are not quantitatively described. (4) Even though GEMs describe the metabolism, concentrations of metabolites are not directly included and flux balance analysis (FBA; Orth et al., 2010) is employed for flux simulations, which assumes that there is no (unexpected) accumulation of metabolites within GEMs.

# USING GEMs FOR ESSENTIALITY AND SYNTHETIC LETHALITY ANALYSIS

As mentioned above, since GEMs are complete metabolic networks, they can be used for gene/reaction essentiality analysis (EA; Edwards and Palsson, 2000). In general, EA identifies all essential genes or reactions whose knockout will disable a specific biological function through FBA. EA could be easily implemented in silico using GEMs by enumerating all single gene/reaction knockouts and testing whether their biological objectives are still functioning. In addition, synthetic lethality analysis (SLA), which scans for combinatory knockouts of multiple reactions/genes that lead to blocking of the target biological function, could also be implemented in a similar way. And recently, several methods have been developed to perform advanced SLA efficiently (Suthers et al., 2009; von Kamp and Klamt, 2014; Pratapa et al., 2015; Zhang et al., 2015).

It's generally believed that gene/reaction EA could be performed by topologic analysis of the metabolic network. However, since the stoichiometric coefficients are absent in topologic metabolic networks, they're less accurate. For example, **Figure 2** shows the topologic network of the toy model from **Figure 1**. Based on its topologic properties, this metabolic work can use D-glucose-6-phosphate, NAD, and phosphate as substrates and produce 3-phospho-D-glycerate, NADH, and a proton. However, this pathway always consumes more ADP than it produces, and produces more ATP than it consumes. Therefore, this pathway will be blocked without ADP supplementation and this finding was not possible by topologic analysis.

Essentially, if a GEM is well established, its EA and SLA results could be very accurate. For example, in the most used E. coli and S. cerevisiae GEMs, around 90% of the predicted essential genes have been validated in vivo (true-negative; Feist et al., 2007; Heavner et al., 2013). This is within expectation, because if a function is blocked in silico, it's very unlikely that there could be a complimentary solution in vivo to recover it. The explanation for the very few false-negative predictions (negative growth in silico and positive growth in vivo) is that there's a knowledge gap, such as unknown enzyme or unknown function of an existing enzyme, which leads to the underestimation of the capability of the GEM. On the other hand, even if the GEMs are 100% complete, there may still be false-positive predictions since the missing information of regulation and protein (enzyme) efficiency could lead to extra constraints to GEMs, thereby rendering a nonessential reaction/gene in silico essential in vivo. It's worth mentioning that, after a certain period of adaptive

respectively. G6P, D-glucose-6-phosphate; F6P, D-fructose-6-phosphate; FDP, D-fructose-1-6-bisphosphate; G3P, glyceraldehyde-3-phosphate; 13DPG, 3-phospho-D-glyceroyl-phosphate; 3PG, 3-phospho-D-glycerate; and Pi, phosphate.

evolution, a false-positive knockout could become nonessential in vivo again (Patil et al., 2005).EA and SLA have mainly been used to validate newly constructed GEMs and in recent years, EA and SLA were applied to study of systems medicine (see Section Using GEMs in Studies of Systems Medicine).

# USING GEMs AS SCAFFOLDS FOR MULTI OMICS DATA INTEGRATION AND INTERPRETATION

Recently, increasing volumes of transcriptomic, proteomic, and metabolomics data are becoming publically available, and it's believed that GEMs are good scaffolds to make use of these multi omics data. In GEMs, omics data could be quantitatively integrated as constraints for metabolic fluxes, thereby allowing systematic and quantitative evaluation of these data, which was not possible using traditional metabolic networks. This is the most significant advantage of using GEMs as scaffolds.

Although, GEMs are metabolic networks, the most used omic data for GEMs are transcriptomic and proteomic. This is because the technic is really advancing in the field and makes large number of high quality transcriptomic and proteomic data available. However, since the GPR relationships are qualitative in GEMs (**Figure 1C**), one needs to make assumptions to define the quantitative relationship between gene/protein expression and metabolic fluxes when integrating transcriptomic or proteomic data into GEMs. This is problematic because the complicated relation between fluxes and expression level of genes and enzymes in vivo are unlikely to be captured by a general assumption (MacHado and Herrgård, 2014). On the other hand, there're many well-defined approaches to integrate fluxomics and metabolomics, data (Khodayari et al., 2014; Martín et al., 2015; Miskovic et al., 2015). However, it's very difficult (if not impossible) to get genome scale data of them. Hence, we suggest that even though omics data are integrated, one should be skeptical about the quantitative results of simulations or predictions from GEMs.

Nonetheless, we believe that it is better to qualitatively interpret the omics data using GEMs. For instance, it would be much more reliable to use omics data to determine the presence or absence of reactions and to construct high-quality and specific GEMs (Zur et al., 2010; Agren et al., 2012, 2014; Mardinoglu et al., 2013; Yizhak et al., 2014). In addition, many researchers started to integrate significance of differential expression of genes with GEMs rather than their quantitative expression to interpret the biological information behind omic data (Patil and Nielsen, 2005; Cakir et al., 2006; Jensen and Papin, 2011; Fang et al., 2012; Navid and Almaas, 2012). Moreover, qualitative interpretation of omics data with GEMs have recently been applied to systems medicine (see Section Using GEMs in Studies of Systems Medicine). These studies demonstrated the usefulness of GEMs as scaffolds. In short, we suggest that GEMs are powerful platforms for integration of omics data for gaining biological insights rather than quantitative results.

# USING GEMs FOR IN SILICO METABOLIC ENGINEERING

Using GEMs for in silico metabolic engineering has been a widely discussed topic for years. It's generally believed that GEM-based methods could predict gene modification strategies for overproduction of desired biochemicals and thus, accelerate the overall metabolic engineering process. In the last decade, various kinds of in silico metabolic engineering methods had been developed and many of them were applied experimentally (Kim et al., 2015; Long et al., 2015; MacHado and Herrgård, 2015).

Although in silico metabolic engineering methods seemed quite different from each other, they follow a similar procedure: (1) they define what a desired strain is and (2) identify approaches that push the wild-type strain to become the desired one. So far, a variety of approaches were used in in silico metabolic engineering, such as reaction/gene knock-out (Burgard et al., 2003; Patil et al., 2005; Kim et al., 2012; Ren et al., 2013; Ruckerbauer et al., 2014; Zhang et al., 2015), overexpression/suppression (Pharkya and Maranas, 2006; Choi et al., 2010; Ranganathan et al., 2010; Park et al., 2012; Chowdhury et al., 2014), foreign pathway knock-in (Pharkya et al., 2004), and swapping the co-factor for a target enzyme (NADH to NADPH or vice versa; King and Feist, 2013). However, the methods for knock-out identification are the majority since a knockout is much easier to define in silico than up-/downregulation of genes as mentioned before. On the other hand, different methods could have independent definition of desired strains. For instance, some of the methods define the desired strain by simply setting thresholds for growth and production, respectively, and others could define the desired strain following some biological assumptions (Edwards et al., 2001; Segrè et al., 2002).

Interestingly, methods pursuing different type of desired strains could all lead to experimentally valid strategies for metabolic engineering (Fong et al., 2005; Trinh et al., 2008; Fowler et al., 2009; Choi et al., 2010; Yim et al., 2011; Ng et al., 2012; Nocon et al., 2014), but the production of target products predicted in silico seldom achieved in vivo. The explanation to this is complicated, and could come from both the computational and experimental side. However, one of the key reasons should be that GEM with only metabolic network is not enough to quantitatively predict the behavior of strains in vivo. In conclusion, we suggested that all kinds of in silico metabolic engineering methods are instructive, but it's better to use them for gaining information rather than to develop exact strategies.

# USING GEMs IN STUDIES OF SYSTEMS MEDICINE

Using GEMs for systems medicine studies have recently been highlighted (Mardinoglu and Nielsen, 2015; Yizhak et al., 2015). GEMs simulate the human metabolism in a holistic way, and this greatly advances systems medicine studies by enabling systematic evaluation of metabolic feature of human disease. Great efforts had been made in reconstructing GEMs of human, and there're now several publically available generic human metabolic networks such as Recon 1, Recon 2, EHMN, and HMR (Duarte et al., 2007; Ma et al., 2007; Agren et al., 2012; Thiele et al., 2013). In addition, since the technology is advancing, tissue specific or cell specific genomic, proteomic and transcriptomic data are becoming available (Cancer Genome Atlas Research Network, 2008; Uhlén et al., 2015). These led to rapid development in reconstruction of high quality tissue or cancer specific GEMs (Zur et al., 2010; Agren et al., 2012, 2014; Mardinoglu et al., 2013) and, therefore, enabled more confident interpretation of metabolism of diseases.

For instance, cancer specific GEMs together with EA and SLA analysis were recently used for identification of oncogenes/metabolites and biomarkers for diagnosing specific cancer (Agren et al., 2014; Jerby-Arnon et al., 2014; Gatto et al., 2015; Gatto and Nielsen, 2015). Since this procedure mainly uses the true-negative part of EA and SLA, the analysis could be highly reliable. For example, (Agren et al., 2014) identified 101 drug targets for liver cancer treatment; and 83 of them are currently in use or have shown strong correlation with cancer progression. In addition, together with multi-omic data, GEMs were used to find the mechanistic explanation of various diseases. By interpreting clinical omic data with GEMs, the mechanistic understanding of non-alcoholic fat liver disease and type two diabetes were reported (Mardinoglu et al., 2014; Väremo et al., 2015). Moreover, GEMs were also used to explore the effect of microbiota (Ji and Nielsen, 2015). By simulate and predict the interaction of gut microbiota and their effect on hosts, several recent studies revealed that microbiota modulate the amino acid and glutathione metabolism of their host (Shoaie et al., 2013, 2015; Mardinoglu et al., 2015). These exciting studies exhibited the great potential of GEMs in the field of systems medicine, and hopefully there would be much more excellent works coming out.

# DISCUSSION

GEMs are very useful platforms and tools for systems biology, but they're still very young compared to traditional ones. Fluxes of reactions could be quantitatively simulated using GEMs, although caution should be exercised before drawing conclusions based on simulated fluxes owing to the huge solution space of GEMs (Reed, 2012). Although solution space could be reduced by adding constraints through integration of omics data, it would be better to gain biological insights by qualitative interpretation of omics data rather than quantitative fluxes.

In order to achieve accurate quantitative prediction, the scope of GEMs should be expanded. The establishment of ME-models set a good example for this (Thiele et al., 2009, 2012). In ME-models, the interaction of genes (mRNA), enzymes, and metabolic fluxes are quantitatively expressed, enabling proper integration of transcriptomic and proteomic data. However, it is still difficult to integrate metabolomics data into ME-models. A potential option to integrate metabolite concentration into GEMs is cybernetic modeling. However, to date, there has been no study on genome-scale cybernetic modeling because there are too many parameters to simulate, making it computationally infeasible.

In general, no model is perfect. Genome-scale modeling methods are still under development and have several drawbacks. In addition, it has been recently reported that many published GEMs are of low qualities (Chindelevitch et al., 2014; Ravikrishnan and Raman, 2015). Therefore, they should be used with caution. As concluded in this review, GEMs are more suitable for qualitative applications at this stage, such as EA and SLA analysis. When using GEMs for quantitative applications such as in silico metabolic engineering, one should be aware of the key assumption behind the method and take the results as instructions. However, it should also be noted that, GEMs are open platforms and have great potential in a wide array of applications. Currently, GEMs are used for simulating the interactions between multiple organisms, multiple tissues (Bordbar et al., 2011), and even between microbiota and human tissues. On the other hand, EA and SLA were developed years ago, but they were not used in the discovery of anti-cancer drugs until recent years. These are good examples of how to explore novel applications based on classical methods. Thus, in future, GEMs can be expected to be more widely used in biotechnology, bioengineering, and many other fields.

# AUTHOR CONTRIBUTIONS

CZ conducted the writing of this paper, and QH modified and edited it.

# ACKNOWLEDGMENTS

This work was funded by National Basic Research Program of China (973 Program) (2012CB721101), and China Scholarship Council.

# REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Zhang and Hua. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Personalized Cardiovascular Disease Prediction and Treatment—A Review of Existing Strategies and Novel Systems Medicine Tools

#### Elias Björnson1, 2, Jan Borén<sup>2</sup> and Adil Mardinoglu1, 3 \*

<sup>1</sup> Department of Biology and Biological Engineering, Chalmers University of Technology, Gothenburg, Sweden, <sup>2</sup> Department of Molecular and Clinical Medicine/Wallenberg Laboratory, University of Gothenburg, Gothenburg, Sweden, <sup>3</sup> Science for Life Laboratory, KTH – Royal Institute of Technology, Stockholm, Sweden

Cardiovascular disease (CVD) continues to constitute the leading cause of death globally. CVD risk stratification is an essential tool to sort through heterogeneous populations and identify individuals at risk of developing CVD. However, applications of current risk scores have recently been shown to result in considerable misclassification of high-risk subjects. In addition, despite long standing beneficial effects in secondary prevention, current CVD medications have in a primary prevention setting shown modest benefit in terms of increasing life expectancy. A systems biology approach to CVD risk stratification may be employed for improving risk-estimating algorithms through addition of high-throughput derived omics biomarkers. In addition, modeling of personalized benefit-of-treatment may help in guiding choice of intervention. In the area of medicine, realizing that CVD involves perturbations of large complex biological networks, future directions in drug development may involve moving away from a reductionist approach toward a system level approach. Here, we review current CVD risk scores and explore how novel algorithms could help to improve the identification of risk and maximize personalized treatment benefit. We also discuss possible future directions in the development of effective treatment strategies for CVD through the use of genome-scale metabolic models (GEMs) as well as other biological network-based approaches.

Keywords: patient stratification, risk estimation, metabolism, systems medicine, systems biology, network medicine

# INTRODUCTION

Cardiovascular disease (CVD), specifically ischemic heart disease and stroke, remains to be the world leading cause of death by a considerable margin (World Health Organization, 2012). It also remains a challenge to accurately predict who is going to develop CVD. For this purpose, several CVD risk-estimating algorithms including the Framingham risk score (Wilson et al., 1998), Reynolds risk score (Ridker et al., 2007), Pan European score (SCORE; Conroy et al., 2003), ASSIGN Scottish algorithm (Woodward et al., 2007), and QRISK2 UK algorithm (Hippisley-Cox et al., 2008) have been developed (Simmonds and Wald, 2012). The purpose of these algorithms are, by considering traditional risk factors for CVD such as age, BMI, smoking status, and blood lipid parameters (**Table 1**), to estimate the 10-year risk of a CVD-event so that preventative measures

### Edited by:

Jiarui Wu, Shanghai Institutes for Biological Sciences, China

#### Reviewed by:

Ranjan K. Dash, Medical College of Wisconsin, USA Marie Csete, Marie Csete Consulting, USA

> \*Correspondence: Adil Mardinoglu adilm@scilifelab.se

#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology

Received: 03 September 2015 Accepted: 06 January 2016 Published: 26 January 2016

#### Citation:

Björnson E, Borén J and Mardinoglu A (2016) Personalized Cardiovascular Disease Prediction and Treatment—A Review of Existing Strategies and Novel Systems Medicine Tools. Front. Physiol. 7:2. doi: 10.3389/fphys.2016.00002



An "X" marks the inclusion of a parameter in the risk score in question.

can be initiated for people who will benefit from this intervention. However, the current algorithms have been developed for population-based prediction of CVD and not for personalized prediction, making the task of predicting exactly who is going to develop CVD difficult. For this reason, even though drugs such as statins have shown tremendous benefit in secondary prevention, in a primary prevention setting the benefits have arguably been modest. Preventative intervention is likely beneficial in a subset of the population, hence accurate risk stratification is an essential tool to enable effective preventative treatment. Rapid and continuous efforts are needed to develop novel biomarkers for achieving high diagnostic accuracy to predict CVD.

Technical breakthroughs have enabled unprecedented progress in the field of omics (i.e., genomics, transcriptomics, proteomics, metabolomics, and lipidomics). Arguably, this should result in great potential in the field of biomarker discovery. Publications in the field of biomarker discovery have increased dramatically over the past two decades, however the increase in the number of clinically useful biomarkers have been meager (Drucker and Krapfenbauer, 2013). In the area of drug development, there is a need for new effective preventative drugs for CVD. But even the most effective drug must be given to the correct subjects. An important distinction must be made between accurate risk identification and accurate personalized prediction of treatment benefit. In a clinical setting, this means that the following two questions should be able to be answered by a CVD risk score as accurately as possible: (i) Will this patient develop CVD within a certain time period? (ii) What is the increase in life expectancy and disease-free years if this particular patient initiates this particular (drug-based or life style-based) intervention? In this review, we discuss the challenges associated with the current CVD risk-estimating algorithms as well as the potential of a systems biology approach to produce better risk scores as well as more effective CVD drugs.

# CURRENT CHALLENGES IN CVD RISK PREDICTION

The ultimate goal of a CVD risk-estimating algorithm is to accurately predict who and when someone is going to develop CVD. This ability should not be confused with the ability of an algorithm to predict how many out of a population will develop CVD during a certain time period. Thus, populationbased prediction is different from personalized prediction. In a study by van Staa et al. (2014) this question was addressed by following 1.8 million subjects for an average of 3.3 years. The three widely used risk prediction algorithms Framingham, ASSIGN, and QRISK2 were evaluated to see if the risk scores accurately predicted not only population-based risk but also personalized risk of CVD. To achieve this, the three risk scores were applied at each of the 1.8 million subjects and compared to a competing risk Cox proportional hazard (CRCPH) model. The study reported that the algorithms accurately predicted how many CVD events would occur in the population, and accurately predicted low-risk subjects. However, for high-risk subjects the three algorithms agreed modestly with the CRCPH model. What this study illustrates is that the Framingham, ASSIGN, and QRISK2 CVD risk scores accurately estimate population-based risks and do identify low risk subjects but the algorithms do not accurately predict who is going to develop CVD.

Predicting benefit from an intervention at a personalized level may be a very valuable tool in CVD treatment. Ferket et al. (2012) estimated how much personalized benefit is gained from statin therapy in a population of 2428 Dutch people. A microsimulation model was used to create a personalized calculator of gains in total and CVD-free life expectancy with statin therapy, and the results of the model for each person was compared with the CVD risk predicted by SCORE. The authors observed an average of 0.3 years of increased life expectancy and 0.7 years of increased CVD-free life expectancy gained from an average of 18.3 years of statin therapy. These gains from statin therapy was considered modest, especially considering that side effects were ignored by the model. Further on, statin therapy is currently encouraged with increasing age due to its correlation with higher CVD risk scores. However, importantly; due to competing risk of death from other diseases, it might not follow that increased 10-year risk of CVD implies larger benefit from statin therapy. For example, as stated in the paper "both a 55-year-old nonsmoking woman with a ten-year CVD mortality risk of 2% and a 65-year-old male smoker with a ten-year CVD mortality risk of 15% might both gain one year of CVD-free life expectancy with statin therapy." For the entire population, 25% with a low SCORE risk achieved equal or larger gains in CVD-free life expectancy than the median gain in participants with a high SCORE risk estimation. This distinction between risk of CVD and benefitof-treatment may appear subtle but is important. For secondary prevention, statin therapy have shown tremendous benefit, but what this study illustrates is the challenge of primary prevention treatment decision and that there exist a need for risk scores which also estimates personalized benefit of treatment.

# CURRENT CVD BIOMARKER DISCOVERY

With the recent advances in metabolomics technologies, hundreds to thousands of metabolites can be simultaneously detected in tissues and biofluids (e.g., blood and urine) to provide a snapshot of the current physiology. Metabolic signatures of obesity (Newgard et al., 2009), future insulin resistance, T2D (Wang et al., 2011), CVD (Shah et al., 2010; Magnusson et al., 2013), NAFLD, and different types of cancer (Ganti and Weiss, 2011; Tan et al., 2012; McDunn et al., 2013; Zeng et al., 2014) have been characterized for identification of associated risk factors as well as for discovery of novel biomarkers.

Branched chain amino acids (BCAAs), valine, leucine, and isoleucine as well as aromatic amino acids, tyrosine, and phenylalanine were discovered to predict the development of diabetes, which is strongly associated with CVD (Wang et al., 2011). Moreover, BCAAs together with the urea cycle metabolite levels in the plasma were used to predict the development of CVD (Shah et al., 2010). Magnusson et al. (2013) developed a method called diabetes-predictive amino acid (DM-AA) score using the metabolic signature of three amino acids (tyrosine, phenylalanine, and isoleucine) and showed that the plasma level of these amino acids correlated with intima-media thickness, plaque formation and exercise-induced myocardial ischaemia, which are three signs of CVD-related abnormalities. The authors also followed 4577 subjects for an average of 12 years, of which 253 suffered a CVD event. Compared to subjects with lowest quartile values of DM-AA score the odds ratio for CVD development were 1.27, 1.96, and 2.20 for quartile 2, 3, and 4, respectively.

Insulin resistance (IR) has been strongly linked to increased risk of CVD (Ginsberg, 2000), yet no measure of IR is included in the current risk-estimating algorithms (**Table 1**). The so called Quantose IR algorithm has been developed to estimate IR using metabolomics and lipidomics data (Cobb et al., 2013). Quantose IR is apart from the level of fasting insulin based on α-hydroxybutyrate and the two lipid species 1-linoleoylglycerophosphocholine and oleate. This algorithm is an example of a possible improvement in the evaluation of IR through the need of only a fasting blood test and it may increase the accuracy of the current CVD riskestimating algorithms; however, this has not been systematically evaluated.

Recently, three lipid species TAG(54:2), CE(16:1), and PE(36:5) were discovered as useful for improving the Framingham risk score in 685 subjects of the prospective population-based Bruneck cohort (Stegemann et al., 2014). Addition of another three lipid species and exclusion of HDLcholesterol and total cholesterol from the Framingham risk score resulted in an additional improvement. Framingham risk score has also been improved by adding the three microRNAs including miR-126, miR-223, and miR-197 as biomarkers of CVD (Zampetaki and Mayr, 2012). Moreover, Bolton et al. (2013) evaluated a panel of 27 single nucleotide polymorphisms (SNPs), discovered from genome-wide association studies, to predict the occurrence of coronary heart disease. Compared to a Cox proportional hazard model based on traditional risk factors, the addition of the SNP panel significantly improved the accuracy of the model. Hence, evident improvements upon the traditional risk scores estimated by the existing algorithms have already been achieved by omicsderived biomarkers of CVD. However, the gains are arguably modest.

# WHY HAVE SO FEW NEW BIOMARKERS BEEN DISCOVERED?

There exist a large discrepancy between the number of biomarker discovery publications and the number of new biomarker patents (Drucker and Krapfenbauer, 2013). For all diseases (not only CVD) only 1–2 new biomarkers were approved by the Food and Drug Administration each year in the US between 1995 and 2009 despite the enormous technical advances in the omics fields during the same period (Anderson, 2010). There are probably a number of reasons for this, including lack of standardized biomarker discovery pipeline, lack of good verification platform for large sample sets and lack of an underlying theory of biomarkers.

There are three categories in which newly discovered potential biomarkers fall into: chance, bias, and generalizability. The only category that may result in a potentially clinically useful biomarker is the latter. The risk of a false discovery increases with increasing number of measured parameters. Therefore, the current ability to measure hundreds to thousands of analytes in a single experiment will result in potential false discoveries. However, this problem can be remedied by commonly used statistical techniques and is therefore probably not the largest explanation to the lack of novel biomarkers.

The issue of bias is however not a problem to be overcome by statistical analysis techniques but is instead inherent in the experimental design. For example, when a biomarker study is commenced a study population is separated into a diseased group and a control group. However, when analyzing the characteristics of the groups, it might be discovered that the diseased group is also in average older and heavier than the control group. Is it then possible to say that a discovered biomarker is a biomarker of the disease, the age, or the weight? For this reason, the groups are often matched against each other to minimize known confounding factors. However, unknown confounding factors might still bias the study. The only remedy to this problem is randomization. Unfortunately, by definition, a biomarker discovery study can never be randomized thus making the risk of so called bias of inequality at baseline an inherent problem of biomarker discovery. How important this issue is and if it can explain the lack of accurate CVD biomarkers is currently unknown, but it is likely an important contributing factor. Bias can also be introduced if the samples from the different groups are treated differently throughout the analysis pipeline. It is therefore of vital importance that the handling and analysis of samples are conducted consistently. If there is bias of inequality at baseline between two groups, there is a risk that a measured parameter will correlate with an unknown confounding factor and not with the disease. The risk to have any discovery due to bias thus increases both with the number of confounding factors and with the number of parameters analyzed. To overcome this problem it might (paradoxically to the field of omics) be desirable to measure as few parameters as possible. Thus, one way of achieving maximum chance of detecting true biomarkers is to have a biomarker theory. An underlying theory would be able to a priori point to what should be measured, thus limiting the need to measure lots of parameters.

As an alternative to the search for a single biomarker of CVD, another approach is to use a panel of biomarkers. If such a panel is to be highly sensitive and highly specific it requires that the individual biomarkers are so called orthogonal against each other. This means that every biomarker adds diagnostic value to the panel rather than just co-vary with other markers. Recent technologies such as protein multiplex platforms do invoke hope that effective biomarker-panels of CVD could be created and used in the clinic.

# NOVEL TOOLS IN SYSTEMS MEDICINE

Genome Scale Metabolic Models (GEMs) are employed for simulating the metabolism of cells/tissues. When generating a GEM, all known metabolic reactions in a particular cell or tissue are integrated into one network topology. Once the model has been constructed, it can be used in conjunction with flux balance analysis which allows for in silico metabolic simulation of the cell or tissue type in question (Mardinoglu and Nielsen, 2012, 2015; Mardinoglu et al., 2013; O'Brien et al., 2015; Yizhak et al., 2015). GEMs in combination with transcriptomics, proteomics, metabolomics, or lipidomics data have the potential to identify perturbed metabolic subnetworks in silico (Agren et al., 2012, 2014; Shoaie et al., 2013, 2015; Yizhak et al., 2013, 2014a,b; Galhardo et al., 2014; Mardinoglu et al., 2014; Gatto et al., 2015; Ghaffari et al., 2015; Varemo et al., 2015; Zhang et al., 2015). GEMs constitute a possible powerful tool in the area of human complex disease since it enables the potential of pathophysiological understanding of a disease (Ryu et al., 2015).

Another interesting tool in systems medicine is protein– protein interaction (PPI) networks (Rolland et al., 2014). PPI networks has the potential to provide useful information in CVD, since each protein is placed in a larger network context and thus alterations in proteins in the diseased state can be compiled and translated into meaningful biological tasks. For example, if 100 different proteins are shown to be altered in the blood macrophages or endothelial cells of people with CVD and 80 of them happen to be highly connected, shown by a PPI, then that part of the network and the related metabolic function could be concluded to be perturbed in the diseased state. Further on, if a few of the proteins are shown to interact with lots of the other disease-related proteins, these highly connected proteins might be central to the disease progress itself. Thus, PPIs could identify central hubs in the disease-network, hubs that might provide pathophysiological understanding and be suitable as drug targets.

As mentioned, an a priori theory of biomarkers could aid in biomarker discovery. A theory of biomarkers could be created through the use of GEMs and PPI networks. A hypothetical example for use of GEMs in CVD would be to model the metabolism of cell types in the blood, for example macrophages, endothelial cells or myocardial cells. If this would be done, predictions about the metabolism of these cells and possible metabolic alterations in CVD could be enabled. Specifically, if GEMs would provide a mechanistic understanding of for example macrophages and their possible metabolic alterations in CVD, a limited set of plausible biomarkers (proteins or metabolites) could be selected and measured independently in a biomarker discovery study. This approach, coupled with stringent experimental biomarker discovery design would limit the risk of bias and could increase the chance of discovering clinically useful biomarkers.

A concrete example for using GEMs which could be relevant to CVD involves macrophage activation. Since there is a link between inflammation and CVD and since macrophages play an important role in the build-up of atherosclerotic plaques, studying the metabolism of macrophages could aid in the understanding of CVD. Bordbar et al. (2012) used genome-scale metabolic modeling in combination with transcriptomics, proteomics, and metabolomics to reveal the metabolic features and modulators of macrophage activation. They identified metabolites which enhanced (glucose and arginine) and suppressed (tryptophan and vitamin D3) macrophage activation. These particular metabolites were previously known to be associated with immunoactivation but the mechanism was unknown. Such a mechanistic insight into what regulates macrophages could help in designing effective interventions. In this case, the plausible intervention would be to limit glucose and arginine intake and increase tryptophan and vitamin D3 intake to decrease the activity of the blood-macrophages. Probably, an intervention like this is not as straight-forward but it does provide a rational approach for the development of treatment strategies which could be tested empirically.

Heart performance is naturally relevant for cardiovascular health and is plausibly affected by the heart's energy metabolism. Little is however known about the energy metabolism of the heart in humans in vivo, during varying nutrient conditions and pathological conditions such as heart failure and diabetes. In order to simulate cardiac performance Karlstädt et al. (2012) developed CardioNet—a GEM covering the metabolism of the cardiomyocyte. Simulations for different nutrient conditions were performed and the efficiency (how much ATP was produced compared to substrate and oxygen consumption) of the heart was evaluated. Differences were seen when comparing different combinations of substrates in terms of cardiac output. The authors observations suggested that high levels of the ketone body acetoacetate (which can be seen in for example diabetes) would decrease cardiac output and increase ROS production indicating possible decreased cardiac contractility. It is currently not known how e.g., diabetes could affect cardiac health. The study by Karlstädt et al. provide a possible pathophysiologic mechanisms of heart malfunction related to diabetes and more generally provide a framework for evaluating how varying oxygen and nutrient conditions could affect the heart.

In order to simulate the entire human cellular and tissue functions in a holistic approach, a whole cell/tissue model could be used. One example of a whole-cell model of the human pathogen Mycoplasma genitalium has been successfully developed and simulation of dynamic cellular states has been demonstrated (Karr et al., 2012). This holistic approach has not yet been employed on human cells but does show the potential use of such models. This process typically involves construction and employment of metabolic, regulatory, signaling, and PPI networks in conjunction with GEMs. The COBRA Toolbox (Schellenberger et al., 2011) and RAVEN Toolbox (Agren et al., 2013) which are valuable supports for researchers in genome-scale metabolic modeling should also be expanded to deal with simulation of these integrative models. Considering the 3675 protein coding genes (18% of the genome) in the generic human GEM HMR2 (Mardinoglu et al., 2014) and their interactions with other proteins in biological networks, such integrated computational models may provide further information about the relationship between the genotype and phenotype of CVD.

There are a number of hurdles to overcome for successful simulation of human metabolism in a biologically relevant matter. Reconstructing GEMs involves correctly defining, for each metabolic reaction, the stoichiometry, the substrate(s), the product(s), the enzyme(s), and the gene(s) which characterize that specific reaction. This information has to be correct for thousands of reactions. During the generation of the GEMs, the network often needs to be so called gap-filled in order for the network to be connected and complete. This gapfilling step is one source of potential errors in the model. Compartmentalization of the reactions is also a relevant issue, not least when constructing human GEMs. It is often not known where a reaction occurs in the cell and whether the substrate/product can be exported/imported into other compartments. Even though there is an extensive effort in defining the subcellular localization of proteins (Kampf et al., 2014; Uhlén et al., 2015), the complete draft information will not be available for another few years.

Another issue relevant for human cell specific metabolic models regards defining the environment. In microbial conditions, the growth media is very well-defined so that the possible uptake and secretion fluxes are also known. For human cells the environment is much less known, which can greatly affect the behavior of the model. For a GEM to simulate the function of a cell/tissue accurately a so called objective function needs to be defined. Usually, maximization of growth is used as an objective function for microorganisms. However, defining an objective function for human cells is not as straight-forward. For human cells, this could feasibly be very context specific, depending on for example regulation and signaling effects. Integrating regulatory and signaling networks with GEMs could therefore be important in order to capture biologically relevant behavior. This integration is however a challenging task due to increase in size of the networks. A GEM usually needs a pre-defined biomass equation. The biomass equation greatly influences the behavior of the model (directs the fluxes) and thus the model is very sensitive to the definition of the biomass equation. A number of issues has been raised on this topic and the genome-scale metabolic modeling community has responded successfully (Chindelevitch et al., 2015; Ebrahim et al., 2015).

Lastly, a model is often validated by its predictive ability, for example for a microorganism GEM to predict the growth rate and production rate of various substances. However, models are rarely shown to not be able to perform infeasible tasks. The unknowns in cell biology coupled with the degrees of freedom in the generated networks makes genome scale modeling challenging. However, several cancer related studies, testifying to the value of the genome scale modeling in portraying a networklevel view of the cancer metabolism and in discovery of novel drug targets and biomarkers have been recently reviewed (Yizhak et al., 2015) and a similar framework could plausibly be used for CVD.

In conclusion, placing high-dimensional omics data in a network context, whether through the use of GEMs, PPIs, or other networks (e.g., regulatory and signaling), may allow for an increased pathophysiologic understanding of CVD. In addition, GEMs together with other networks could provide a rational approach to biomarker discovery, limiting the risk of bias and increasing the chance of improving CVD risk scores (**Figure 1**). However, important limitations do currently exist regarding the biological relevance of human GEMs.

# NETWORK MEDICINE AND DRUG DEVELOPMENT

As stated, network-dependent analyses may allow for identification of metabolic perturbations in CVD. Biological networks have arguably evolved to be robust. For example, single blockade of 85–90% of all proteins in yeast do not result in any noticeable phenotypic alterations (Peters, 2013). Similarly, knock-out studies in mice suggest that only 10% of all potential drug target genes would have any effect as single targets (Peters, 2013). In the traditional reductionist approach to drug development, a disease modifying activity is reduced to a single

target. While this can be effective for certain diseases, it may not be enough for treatment of a complex disease such as CVD. CVD specifically could have multiple or complex causes which result in network-level perturbations. If this is the case, an alternative approach to CVD drug development would be identification of network-level perturbations and developing drugs that can affect the network rather than only a single protein.

The upcoming branch of network medicine or polypharmacology, integrates systems biology tools with pharmacology. Recently, a drug-target and a target–target interaction network was constructed to identify which targets of CVD drugs that possesses the most interconnectedness with drugs and other targets (Zheng et al., 2014). These targets have high probability of being important hubs in the CVD-related metabolic networks and thus interesting to treat with a multi-target compound. Subsequent virtual screening of compounds revealed several potential multi-target drug candidates and in vitro validation of five randomly selected candidate compounds revealed that four of them could indeed bind to these targets and thus possibly affect the CVDrelated metabolic network. However, this approach to drug discovery could perhaps also increase the risk of adverse effects precisely because the compounds in question are unspecific. Nevertheless, this method illustrates how a polypharmacological approach to CVD drug development could be conducted. If these types of methods of drug development will produce effective CVD-risk lowering interventions remains however to be seen.

Risk scores based on multi-biomarker panels might also aid in system-level drug development. If a potential drug affects a single target but does not affect a plethora of other biomarkers, this could provide an early indicator that the drug candidate might not prevent CVD. However, if multiple markers change after administration of a potential drug candidate, that might be indicative of reduced risk of CVD and a potentially successful drug. Risk scores based on multi-biomarker panels could of course similarly be used for evaluation of other types of interventions such as diet, and not only drug-based interventions. The field of polypharmacology is, albeit promising, still new. Future efforts in this area could hopefully result in the development of novel preventative CVD medications.

# CONCLUSIONS

Systems medicine uses omics data for reconstruction of cellular networks. High dimensional omics data is often not easy to directly translate into biological meaning. Therefore, the systems medicine approach could, by integrating different kinds of omics data and putting them in a network context, enable pathophysiological understanding of a disease in question. Systems medicine aims at identifying how the integrated network, rather than single genes or proteins, is altered in a diseased state. This approach allows for identification of perturbed subnetworks and may, apart from providing pathophysiologic understanding of the disease, also create a base to predict biomarkers and identify subnetworks as drug targets. This information could lead to more accurate CVD risk scores as well as more effective drugs/interventions.

In conclusion, it is important for each patient to understand his/her own risk of CVD as well as likely benefit of treatment to weigh against any potential side effects, thus there is a need for accurate personalized risk scores in conjunction with personalized prediction of treatment benefit. As illustrated, current risk-estimating algorithms can in this setting be improved upon. Accurate risk scores, more effective drugs and personalized estimation of benefit from treatment are three much needed tools in the area of CVD prevention. A systems medicine approach can hopefully provide value in all these areas.

# REFERENCES


# AUTHOR CONTRIBUTIONS

All three authors actively contributed in writing and editing of the manuscript.

# ACKNOWLEDGMENTS

This work was supported by the Bill and Melinda Gates Foundation, the Knut and Alice Wallenberg Foundation, European project FP7 METACARDIS (grant agreement HEALTH-F4-2012-305312/METACARDIS), Novo Nordisk A/S and the Innovative Medicines Initiative Joint Undertaking under EMIF grant agreement n◦ 115372.

disease-relevant nodes of the human metabolic network. Nucleic Acids Res. 42, 1474–1496. doi: 10.1093/nar/gkt989


differentiates obese and lean humans and contributes to insulin resistance. Cell Metab. 9, 311–326. doi: 10.1016/j.cmet.2009.02.002


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Björnson, Borén and Mardinoglu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Cancer Metabolism: A Modeling Perspective

#### Pouyan Ghaffari <sup>1</sup> , Adil Mardinoglu1, 2 and Jens Nielsen1, 2, <sup>3</sup> \*

*<sup>1</sup> Department of Biology and Biological Engineering, Chalmers University of Technology, Gothenburg, Sweden, <sup>2</sup> Science for Life Laboratory, KTH - Royal Institute of Technology, Stockholm, Sweden, <sup>3</sup> Novo Nordisk Foundation Center for Biosustainability, Technical University of Denmark, Hørsholm, Denmark*

Tumor cells alter their metabolism to maintain unregulated cellular proliferation and survival, but this transformation leaves them reliant on constant supply of nutrients and energy. In addition to the widely studied dysregulated glucose metabolism to fuel tumor cell growth, accumulating evidences suggest that utilization of amino acids and lipids contributes significantly to cancer cell metabolism. Also recent progresses in our understanding of carcinogenesis have revealed that cancer is a complex disease and cannot be understood through simple investigation of genetic mutations of cancerous cells. Cancer cells present in complex tumor tissues communicate with the surrounding microenvironment and develop traits which promote their growth, survival, and metastasis. Decoding the full scope and targeting dysregulated metabolic pathways that support neoplastic transformations and their preservation requires both the advancement of experimental technologies for more comprehensive measurement of omics as well as the advancement of robust computational methods for accurate analysis of the generated data. Here, we review cancer-associated reprogramming of metabolism and highlight the capability of genome-scale metabolic modeling approaches in perceiving a system-level perspective of cancer metabolism and in detecting novel selective drug targets.

Keywords: cancer metabolism, metabolic modeling, systems biology, systems medicine, genome scale metabolic reconstruction, metabolic networks and pathways, tumor metabolism

# INTRODUCTION

The past decades has seen a dramatic expansion in investigations on mechanism of cancer related metabolic adaptations, and this has resulted in accumulated evidences suggesting considerable association between several pathways in human metabolism and malignant transformation (Vander Heiden et al., 2009; Cairns et al., 2011; Schulze and Harris, 2012). Recently, the state of "deregulated cellular metabolism" was added by Hanahan and Weinberg as one of the hallmarks of cancer (Hanahan and Weinberg, 2000), reflecting the overall consensus around the idea of altered cellular metabolism through neoplastic progression (Hanahan and Weinberg, 2011). Consistent with this approach, investigation of cancer-associated metabolic alterations attracted considerable effort and resulted in successful identification of several selective metabolic targets that started to enter clinical trials (Vander Heiden, 2011; Galluzzi et al., 2013). Activation of oncogenes and deactivation of tumor suppressor genes, have been linked to cancer-associated metabolic reprogramming (Levine and Puzio-Kuter, 2010). Also, accumulation of metabolites such as 2-hydroxyglutarate, fumarate, and succinate has been proposed to drive oncogenesis, due to disruption in enzymatic activity of

### Edited by:

*Radhakrishnan Nagarajan, University of Kentucky, USA*

#### Reviewed by:

*Osbaldo Resendis-Antonio, Instituto Nacional de Medicina Genomica, Mexico Preetam Ghosh, Virginia Commonwealth University, USA*

> \*Correspondence: *Jens Nielsen*

*nielsenj@chalmers.se*

#### Specialty section:

*This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology*

Received: *23 September 2015* Accepted: *24 November 2015* Published: *16 December 2015*

#### Citation:

*Ghaffari P, Mardinoglu A and Nielsen J (2015) Cancer Metabolism: A Modeling Perspective. Front. Physiol. 6:382. doi: 10.3389/fphys.2015.00382*

isocitrate dehydrogenase (IDH), fumarate hydratase (FH), and succinate dehydrogenase (SDH), respectively (Isaacs et al., 2005; Selak et al., 2005; Dang et al., 2010). Additionally, some agents that conventionally have been used in cancer treatment–such as gemcitabine, 5-fluourouracil and methotrexate—are in fact metabolic enzyme inhibitors (Chabner and Roberts, 2005). In line with these findings, accumulated evidences from recent epidemiological studies supported influence of whole-body metabolism on tumor progression and drug response (Vander Heiden et al., 2009; Vander Heiden, 2011; Galluzzi et al., 2013).

High-throughput omics technologies experienced a dramatic breakthrough during the last decade enabling simultaneous measurement of interacting molecular components in context of complex cellular structure and characterization of transformed cellular processes at genome-scale. Accumulation of DNA, microRNA and protein expression measurements, along with more detailed metabolomics data, has revealed a broader view on cancer-associated metabolic shifts. The more advanced understanding of molecular and genetic events underlying neoplastic transformation acquired, the more complex portrait of tumor metabolism has emerged. Today, analyzing multilayer data generated via high-throughput technologies and integrating them into a descriptive unified model is a nontrivial challenge. Computational methods, capable of processing and integrating multi-dimensional data, have been employed in investigation of metabolic states in health and disease, and in identification of new selective targets and biomarkers. Here, we will review current knowledge in cancer metabolism, underlining capabilities of genome-wide metabolic models in opening new therapeutic windows.

# FUELS FOR CANCER CELLS

One of the best known features of most tumor cells is utilizing high amounts of glucose and metabolizing it differently from normal cells, converting pyruvate derived from glucose to lactate in the cytosol and secreting it rather than oxidizing pyruvate in mitochondria. Normal cells principally elevate conversion of glucose to lactate in hypoxic conditions, whereas this phenotype is common in transformed cancer cells even when oxygen is abundant, a phenotype first described by Otto Warburg more than fifty years ago and referred to as "The Warburg effect" or aerobic glycolysis (Warburg, 1956). Glycolysis can produce ATP faster but far less efficient than oxidative phosphorylation. This shift makes tumor cells to be dependent on high rate of glucose consumption to meet their biosynthesis, energy, and redox requirements (Semenza et al., 2001; Cairns et al., 2011). This metabolic phenotype of cancer cells, dramatic increase in glucose uptake, has been studied extensively and used clinically in cancer diagnosis to visualize tumors by 2-deoxy-2-[fluorine-18] fluoro-D-glucose positron emission tomography (18F-FDG PET; Som et al., 1980; Kelloff et al., 2005).

Elevated glycolytic flux promotes shunting of compounds into branched metabolic pathways to synthesis macromolecules needed for proliferation (Cairns et al., 2011). Downstream flux of glycolysis intermediates into oxidative and non-oxidative arms of the pentose phosphate pathway (PPP) produces ribose-5-phosphate and NADPH, two essential components for tumor cell growth. Ribose-5-phosphate is a precursor for nucleotide synthesis and NADPH is required to handle redox stress (Lunt and Vander Heiden, 2011; Dang, 2012). Another branched pathway from glycolysis is serine biosynthesis, which is important for nucleotide, amino acid, and lipid biosynthesis (Lunt and Vander Heiden, 2011). Sustained proliferation of some melanoma and breast cancer cells has been associated to amplification of phosphoglycerate dehydrogenase (PHGDH), the enzyme catalyzing the first step of serine biosynthesis (Locasale et al., 2011; Possemato et al., 2011). Furthermore, the metastasis of breast cancer cell has been associated with the up-regulation of the serine biosynthesis pathway (Possemato et al., 2011; **Figure 1**).

Glutamine plays a key role in sustaining rapid cell proliferation (Jain et al., 2012). Glutamine is the most abundant amino acid in culture media (Eagle et al., 1956) and plasma (Stein and Moore, 1954) and seems to be in excess relative to required amounts for protein and nucleotide synthesis (DeBerardinis et al., 2007). Tumor cells experiencing aerobic glycolysis need glutamine carbon to replenish TCA cycle intermediates and sustain enhanced biosynthetic metabolism to support cell proliferation (DeBerardinis et al., 2007; Lunt and Vander Heiden, 2011; Mullen et al., 2012). Glutamine is also a primary source of nitrogen for the cell (DeBerardinis et al., 2007). Under hypoxic conditions, glutamine can undergo reductive deamination to generate alpha-ketoglutrate and consequently oxaloacetate, pyruvate, and acetyl-CoA to sustain anabolic metabolism (Mullen et al., 2012). Although most cell culture-based studies of cancer metabolism have focused on the utilization and fate of glucose and glutamine, tumor cells in vivo have access to other sources of nutrients, like amino acids, to sustain the elevated proliferation rate. Measuring consumption and release profiles of metabolites from the NCI-60 panel of cell lines identified high correlation between glycine consumption and cancer cells proliferation rates (Jain et al., 2012). Exogenous serine uptake rate increases dramatically in tumor cells and deprivation of serine acts as a trigger to activate serine synthesis pathway and rapid inhibition of aerobic glycolysis, which results in an increased flux to the TCA cycle (Maddocks et al., 2012). Glycine and serine can be inter-converted by serine hydroxymethyltransferase (SHMT) and be used for one-carbon metabolism and nucleotide synthesis (Labuschagne et al., 2014; Boroughs and DeBerardinis, 2015). The directionality of this inter-conversion has critical effect on cancer cell proliferation. Exogenous serine can be used both for protein biosynthesis and it can be converted to glycine and one-carbon units needed for de novo nucleotide biosynthesis, whereas exogenous glycine cannot compensate for nucleotide synthesis (Labuschagne et al., 2014). These findings may reflect the fact that tumor cell proliferation is supported by serine rather than glycine. Glycolysis also plays an essential role for nucleotide biosynthesis (Lunt et al., 2015) and understanding relative consumption rate of exogenous amino acids compared to glucose-derived serine and glycine in transformed cells will be important.

In addition to glutamine, serine and glycine, other amino acids may also contribute to cancer cell proliferation. Branched-chain

FIGURE 1 | Overview of cancer-associated metabolic pathways. The main metabolic pathways that contributes to malignancy and offer potential drug targets are illustrated. Metabolic enzymes that have been associated with tumor initiation and growth are marked in red. GLUT1, glucose transporter; HK, hexokinase; 6PGD, 6-Phosphogluconate dehydrogenase; PFKFB3, 6-phosphofructo-2-kinase; GAPDH, Glyceraldehyde 3-phosphate dehydrogenase; PHGDH, phosphoglycerate dehydrogenase; PGAM1, Phosphoglycerate mutase 1; PKM, Pyruvate kinase; LDHA, lactate dehydrogenase A; MCT, Monocarboxylate transporter; PDH, Pyruvate dehydrogenase; CPT1, Carnitine palmitoyltransferase I; FASN, fatty acid synthase; RNR, ribonucleotide reductase; FH, fumarate hydratase; SDH, succinate dehydrogenase; IDH, isocitrate dehydrogenase; GLUD, glutamate dehydrogenase; GLS1, glutaminase 1; ASCT2, Amino-acid transporter 2; ACLY, ATP citrate lyase; ACC, acetyl-CoA carboxylase; ACSS2, Acetyl CoA synthetase2; DHFR, DHF reductase; TYMS, thymidylate synthase; HMGCR, HMG-CoA reductase; CK, choline kinase; G6P, glucose-6-phosphate; F6P, fructose-6-phosphate; PRPP, 5-phospho-alpha-D-ribose 1-diphosphate; IMP, inosine monophosphate; UMP, uridine monophosphate; dNTP, deoxynucleotide triphosphate; G3P, glyceraldehyde 3-phosphate; 3PG, 3-phosphoglycerate; 2PG, 2-phosphoglycerate; PEP, phosphoenolpyruvate; OAA, Oxaloacetate; AKG, α-ketoglutarate; HMG-CoA, 3-hydroxy-3-methyl-glutaryl coenzyme A; THF, tetrahydrofolate; 5,10 mTHF, 5,10-methylene tetrahydrofolate; DHF, dihydrofolate.

amino acids (BCAAs) are abundant amino acids in plasma (Stein and Moore, 1954; Meister, 1965), and growth of wild type IDH glioma, subgroup of brain tumors with poorest clinical treatment, is highly associated with expression of branched-chain amino acid transaminase 1 (BCAT1; Yan et al., 2009; Tönjes et al., 2013). Metabolomics profiling of patient-derived glioma samples also suggested correlation between increasing tumor grade and cysteine catabolism (Prabhu et al., 2014). Although, the current state of investigations suggest that amino acids primarily are used for protein synthesis in proliferating cells (Dolfi et al., 2013; Zhang et al., 2014), whereas catabolism of amino acids might be more important to generate ATP and maintain cellular redox state in nutrients limited condition.

In addition to the metabolism of carbohydrates and amino acids, lipids can also be used as an important fuel to supplement cancer cells proliferation requirements. Uptake of lipoproteins and free fatty acids (FFAs) from the bloodstream is the main source of satisfying lipid requirement in adult mammalian tissues. Although, fatty acids biosynthesis is limited to a subgroup of tissues, including adipose, liver and breast, reactivation of lipid synthesis is commonly observed in tumor cells with different sites of origin (Menendez and Lupu, 2007; Abramson, 2011). In vitro, glucose supplies significant amount of carbon needed for the de novo synthesis of lipids, however, in hypoxic condition or expression of oncogenic Ras, phospholipids uptake can contribute to lipid pools, compensating decreased flux of glycolytic carbon through pyruvate dehydrogenase (Kamphorst et al., 2013). Hypoxia may also change the fate of glutamine entering mitochondria by preferring reductive carboxylation to support tumor cell growth by activating fatty acid biosynthesis (Metallo et al., 2012; Mullen et al., 2012). Recently, regardless of relatively low levels in serum and culture media, acetate was proposed as an important carbon source for fatty acid biosynthesis and mitochondrial metabolism in hypoxic or highly glycolytic tumors. Acetyl-CoA synthetase enzyme (ACSS1; Björnson et al., 2015) and (ACSS2; Comerford et al., 2014; Schug et al., 2015) plays a key role in cancer cell proliferation, by capturing acetate and converting it to acetyl-CoA.

# GENOME-SCALE MODELING OF CANCER METABOLISM

Biological systems are complex interactive networks with interconnected set of components (e.g., metabolites, proteins, nucleic acids) and the function of the many different pathways connecting these components is highly regulated. Mutations may cause dysfunction of some of these regulatory or functional pathways, and this may lead to the emergence of dysfunctional phenotypes. Uncovering how these biological systems orchestrate their activities to support specific phenotypic transformation, e.g., normal to cancer, is a major challenge in medical science, and gaining insight into the mechanisms underlying these transformation may enable improved disease diagnostic, prognostic, and treatment strategies (Hyduke et al., 2013; Resendis-Antonio et al., 2015). Recent technological breakthroughs in high-throughput omics techniques and next generation sequencing (NGS) methods has transformed biomedicine into a data-rich discipline capable of providing deeper insights into phenotypic states of cells by simultaneous measurement of a large number of cellular components. However, analyzing very large sets of omics data with the objective to extract new biological knowledge is not a trivial effort. Despite this significant technological progress, we can still only measure fluxes in eukaryote cells for a limited number of reactions in central metabolism (Niklas et al., 2010). Systems biology approach in general, and genome scale metabolic models (GEMs) in particular, can facilitate this hurdle and bridge this gap by allowing multi-layer integration of omics data in the context of whole biological system (Mardinoglu and Nielsen, 2012; Mardinoglu et al., 2015; Yizhak et al., 2015). Furthermore, GEMs enable simulation of multi-species relationships in different metabolic states under dynamic environmental and genetic perturbations (Oberhardt et al., 2009; Shoaie et al., 2015; Zhang et al., 2015). In this context, GEMs are a powerful framework to analyze omics data in health and disease states and to investigate fundamental cellular mechanisms (Mardinoglu and Nielsen, 2015).

Reconstruction of a GEM is performed by assembling biochemical transformations, occurring within a specific cell or tissue into the metabolic network (Thiele and Palsson, 2010; Mardinoglu and Nielsen, 2015). In GEMs the constructed stoichiometry matrix of the network covers stoichiometric coefficients of the metabolic reactions supplemented by detailed mapping of the protein coding genes to their corresponding reactions. In general, the metabolic networks are assumed to be modeled in quasi-steady state, and reconstructed GEMs are analyzed with constraint-based modeling (CBM) techniques (**Figure 2**). CBM shapes feasible solution space by imposing physico-chemical constraints, including thermodynamics, mass balance, and minimum/maximum flux capacity boundaries. The generated models usually remain under-determined, with possible alternative flux distributions satisfying the constraints. Flux balance analysis (FBA) method is used to formulate the problem and select an optimal flux distribution by optimizing for an objective function such as maximum biomass yield or ATP production (Thiele and Palsson, 2010; Mardinoglu et al., 2013b; Simeonidis and Price, 2015; Yizhak et al., 2015).

To date, a number of generic GEMs of human metabolism including Recon1 (Duarte et al., 2007), Recon2 (Thiele et al., 2013), HMR (Mardinoglu et al., 2013a), and HMR2 (Mardinoglu et al., 2014) have been reconstructed (Duarte et al., 2007; Thiele et al., 2013; Mardinoglu et al., 2013a, 2014). These models represent an assembly of all reactions documented to take place in metabolism of human cells/tissues integrated with known genes catalyzing each reaction, and have been used to generate context based GEMs of healthy human cells/tissues as well as transformed cells. In recent years, availability of cancer related high-throughput omics data made it possible to map this data into the generic human GEMs and to reconstruct cancer-specific genome-scale models. Several methods aiming to acquire tissue-specific or condition-specific active metabolic networks from a generic model have been developed. One of the first attempts in this context was the Gene Inactivity Moderated by Metabolism and Expression (GIMME) algorithm, which uses mRNA expression data as input together with presumed metabolic objectives to develop the context-specific reconstructed models (Becker and Palsson, 2008). Shlomi et al. (2008) proposed a computational method to generate tissuespecific metabolic networks by integrating tissue-specific gene and protein expression data with generic human genome-scale metabolic network. Same authors developed another method to reconstruct tissue-specific genome-scale metabolic models by integrating transcriptomics, proteomics, metabolomics, and phenotypic data. They use their method to reconstruct a functional metabolic model of human hepatocytes (Jerby et al., 2010).

Folger et al. (2011) has used a core set of 197 highly expressed enzymes common to at least 90% of the cells in NCI-60 cell lines database to generate a generic genome-scale metabolic model of cancer. This model has been used to capture core metabolic functions shared by cancer cell lines and to identify potential anti-cancer drug targets (Folger et al., 2011). Moving forward, Agren et al. (2012) reconstructed active metabolic networks in genome-scale for 69 healthy and 16 cancer cell types based on protein abundance data. A comparative analysis of generated models allowed prediction of cancer-associated metabolic features with potential to be used in identification of novel drug targets (Agren et al., 2012). Similarly, Wang et al. (2012) reconstructed GEMs for 26 human cancer and

physico-chemical constraints. (D) FBA provides an optimum feasible flux distribution relevant to defined objective function and compatible with enforced constraints.

counterpart normal tissues, inferring tissue-specific metabolic models based on network topology and gene expression data. Pathway-level analysis of GEMs revealed eicosanoid metabolic pathway as a potential selective drug target (Wang et al., 2012). Jerby et al. developed a method to infer metabolic phenotypes by integrating transcriptomics and proteomics data derived from breast cancer patients into the GEMs. They identified a tradeoff linking decrease in tumor cells proliferation rate to evolved metastatic ability, as well as, to increased need for ROS detoxification (Jerby et al., 2012). Later, Goldstein et al. in an attempt to identify the role of p53 in regulating metabolic pathways, employed constraint-based modeling to characterize metabolic changes imposed by varying p53 status in human liver-derived tumor cells. Their results suggested that p53 may regulate glucose metabolism by preventing it to be shunted to growth promoting pathways such as glycolysis and the pentose phosphate pathway (PPP) and therefore inhibiting tumorigenesis (Goldstein et al., 2013). Feizi et al. used GEMs generated based on NCI-60 cell lines database to identify growthassociated metabolic sub-networks. They suggested critical role of concurrent synthesis and degradation of lipids in supplying energy for cell growth, and negative correlation of growthassociated sub-networks with colon cancer patients' survival (Feizi and Bordel, 2013). In addition to the cancer specific GEMs, Agren et al. (2014) reconstructed personalized cancer GEMs for 27 hepatocellular carcinoma (HCC) patients and used these models to predict 101 common and 46 patient-specific potential antimetabolites that could inhibit tumor growth in HCC patients. Recently, cell line-specific GEMs (CL-GEMs) were developed and successfully employed to identify cell-level metabolic phenotypes and selective drug targets (Yizhak et al., 2014; Ghaffari et al., 2015). We recently generated CL-GEMs for 11 human cancer cell lines with different site of origin and used them to investigated expressional heterogeneity of metabolic pathways across cancer cells. Furthermore, we used CL-GEMs to identify antimetabolites aiming to simultaneous inhibition of multiple growth supporting enzymes in proliferative cells. We predicted 60 common and 15 cell/cell group-specific potential antimetabolites and experimentally validated one of the identified anti-growth factors.

These studies demonstrate the successful employment of GEMs in deciphering metabolic foundations of neoplastic phenotypes and identification of new selective biomarkers and drug targets. Next generation of GEMs that integrate the cell metabolism with regulation and signaling are crucial but nontrivial challenge ahead, to tackle complexity of cancer cell metabolism. Up to now, a couple of methods have been developed to integrate the dynamic behavior of regulatory, signaling and metabolic networks in prokaryotes, and single-cell eukaryotes. The integrated FBA (iFBA) approach was used to generate an integrated model of central carbon metabolism in Escherichia coli by introducing ordinary differential equations (ODEs) and regulatory Boolean logic into FBA (Covert et al., 2008). Integrated dynamic FBA (idFBA) method proposed a FBA-based framework by incorporating metabolic, regulatory, and signaling processes through an integrated stoichiometric formalism, assuming fast reactions in quasi-steady state condition and introducing slow reactions in a time-delay manner (Lee et al., 2008). idFBA was applied to evaluate the phenotypic effects of environmental cues on Saccharomyces cerevisiae. However, large number of parameters required to be considered and accurate prediction of all the rate constants for the regulatory and signaling pathways is a long-standing challenge for this kind of approaches.

# CANCER METABOLIC HETEROGENEITY, TUMOR MICROENVIRONMENT, AND SITE-OF-ORIGIN

The progressive accumulation of knowledge of cancer biology has demonstrated that cancer is an extraordinarily heterogeneous and complex disease. Tumors evolve through a reiterative process of clonal extension, genetic diversification and selection within the highly dynamic and adaptive ecosystem of target tissue (Greaves and Maley, 2012). Cancer cells rewire their metabolism to satisfy demands of cell proliferation and survival within the constantly changing tumor microenvironment (Hanahan and Weinberg, 2011). Despite common tendency of cancer cells toward uniformity of metabolism to support growth stimulating adaptations, variable patterns of clonal architecture, genetic diversity, and environmental effects results in a heterogeneous metabolic signature of tumors (Greaves and Maley, 2012; Meacham and Morrison, 2013). Apparently, a single model of cancer metabolism cannot describe the diversity of metabolic alterations happening during tumor progression. Differences in tumor micro-environmental conditions and also nutrients availability may constrain the potential metabolic repertoire of the cancer cells. In vivo, nutrient availability varies for different cell types, inside and across tissues/organs. For example, the structure of liver enforces gradients of nutrients and oxygen accessibility through hepatocyte zones (Puchowicz et al., 1999). Considerable expressional variation has been observed in cell cycle markers, together with association between oxygen gradients and hypoxia-related gene expression in glioblastoma cancer cells (Patel et al., 2014). Differences in active state of pyruvate kinase have been associated with proliferating and nonproliferating cell populations within breast cancer, revealing an influence of glucose metabolism on tumor formation and growth (Israelsen et al., 2013; Iqbal et al., 2014; Wong et al., 2015). Comparing expressional patterns of metabolic genes from 22 different human tumor types demonstrated overall similarity of gene expression program between tumors and corresponding normal tissue (Hu et al., 2013). Reconstructed CL-GEMs for 11 human cancer cell lines have been employed to investigate different mRNA expression pattern of metabolic pathways, and this revealed expressional heterogeneity of metabolic pathways across cancer types which can be exploited to identify generic or cancer-specific metabolic targets (Ghaffari et al., 2015). Siteof-origin may define the effect of oncogenic drivers in different tissues, as Myc-driven lung and liver tumors display dissimilar phenotypes related to glutamine metabolism (Yuneva et al., 2012).

Tumor-host metabolic interactions also have strong effect on the cancer cell metabolic reprogramming. Cancer cells may attempt to induce growth favoring conditions by actively manipulating the tumor microenvironment which can directly influence tumor progression, metastasis, and redox status (Guillaumond et al., 2013; Shukla et al., 2014). Metabolic cooperation between intra-tumor cell populations can help cancer cells to handle spatial heterogeneity of environmental conditions and nutrients availabilities. It has been shown that lactate out-flux by hypoxic pancreatic ductal adenocarcinoma cells can fuel the growth of neighboring normoxic cancer cells (Sonveaux et al., 2008; Guillaumond et al., 2013). Better understanding of tumor metabolic heterogeneity and tumor-host interactions has a high potential to optimize therapeutic strategies for selectively targeting cancer. Moreover, identification of the secreted proteins or peptides from cancer cells which remodel the tumor microenvironment may allow for discovery of novel cancer biomarkers and therapeutic targets (Feizi et al., 2015).

# CONCLUDING REMARKS

Cancer cells are shown to experience characteristic changes in their metabolic programs, including increased uptake of glucose, enhanced rates of glutaminolysis and fatty acids synthesis, suggesting that metabolic shifts supports tumor cells growth and survival. Similarity some cancer-associated metabolic alterations are similar to the metabolic response of normal cells to growth-promoting signals, and this makes it difficult to separate neoplastic alterations clearly from the ones just reflecting increased cellular proliferation. However, different metabolic components target distinct oncogenic signaling pathways, and it is therefore important to elucidate the complex interaction between cellular metabolic network and oncogenic signaling network. Obviously, transformation from oxidative phosphorylation to aerobic glycolysis cannot simply explain the complete metabolic reprogramming event of tumor cells, and fairly little is known about different metabolic activities of cancer cells with diverse genetic and mutational background. Moreover, variances in metabolic profiles of tumors with same genetic lesion but different tissues of origin, suggests an open therapeutic window through investigation of complex metabolic interplay between tumor cells and stroma. It is becoming more evident that metabolic reprogramming in cancer cells is closely connected to hypoxia and ROS metabolism. Disrupting the balance between anti-oxidant production and increased biosynthetic activities in tumor cells seems to be a therapeutic opportunity.

Recent advances in high throughput omics technologies along with continuous progress in our understanding of cancer cell biology have portrayed a more complex and heterogeneous picture of tumor metabolism. Increased understanding of genetic and metabolic heterogeneity of cancer cells may open a new road toward development of new selective personalized diagnostic and therapeutic methods. GEMs, providing a mechanistic description of relationships between genes, metabolites and reactions within an interconnected functional metabolic network, have potential to integrate large-scale experimental datasets and extract knowledge out of their multi-dimensional complexity. Despite recent leap forward in employing GEMs to study cancer metabolism, more technical, and translational challenges are laying ahead. To date, the majority of models have been developed based on in vitro data that were produced under experimentally controlled conditions, but utilization of richer in vivo and clinical datasets, together with integration of cellular regulatory and signaling mechanism is critical for advancement of the field. Introduction of more metabolomics, fluxomics and growth related data along with other related omics data into the reconstruction and validation process of GEMs can help in increasing the accuracy and prediction power of these models. Furthermore, more and more data points toward the important role of tumor microenvironment in promoting plasticity and supporting cancer-associated metabolic adaptations, and modeling metabolic interaction between cancer and its environment can help in developing more effective therapeutic strategies. For this it may be necessary to integrate kinetics into stoichiometric genome-scale models.

In conclusion, metabolic reprogramming is crucial to support the uncontrolled proliferation, survival and migration of cancer cells, but at the same time renders tumors more vulnerable to metabolic perturbations. Identification of these metabolic dependencies at the genome-scale may provide an opportunity for optimized therapeutic intervention.

# REFERENCES


# AUTHOR CONTRIBUTIONS

PG, AM, and JN wrote the paper and all authors were involved in editing the paper.

# ACKNOWLEDGMENTS

This work was supported by grants from the Knut and Alice Wallenberg Foundation. JN acknowledges funding from the Novo Nordisk Foundation for research in his group.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Ghaffari, Mardinoglu and Nielsen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Reconstructed Metabolic Network Models Predict Flux-Level Metabolic Reprogramming in Glioblastoma

### Emrah Özcan and Tunahan Çakır\*

Computational Systems Biology Group, Department of Bioengineering, Gebze Technical University, Gebze, Turkey

Developments in genome scale metabolic modeling techniques and omics technologies have enabled the reconstruction of context-specific metabolic models. In this study, glioblastoma multiforme (GBM), one of the most common and aggressive malignant brain tumors, is investigated by mapping GBM gene expression data on the growth-implemented brain specific genome-scale metabolic network, and GBM-specific models are generated. The models are used to calculate metabolic flux distributions in the tumor cells. Metabolic phenotypes predicted by the GBM-specific metabolic models reconstructed in this work reflect the general metabolic reprogramming of GBM, reported both in in-vitro and in-vivo experiments. The computed flux profiles quantitatively predict that major sources of the acetyl-CoA and oxaloacetic acid pool used in TCA cycle are pyruvate dehydrogenase from glycolysis and anaplerotic flux from glutaminolysis, respectively. Also, our results, in accordance with recent studies, predict a contribution of oxidative phosphorylation to ATP pool via a slightly active TCA cycle in addition to the major contributor aerobic glycolysis. We verified our results by using different computational methods that incorporate transcriptome data with genome-scale models and by using different transcriptome datasets. Correct predictions of flux distributions in glycolysis, glutaminolysis, TCA cycle and lipid precursor metabolism validate the reconstructed models for further use in future to simulate more specific metabolic patterns for GBM.

Keywords: aerobic glycolysis, glutaminolysis, constraint-based models, omics data, tumor subtypes, GBM-specific metabolic model

# INTRODUCTION

Among malignant brain tumors, the most common one is glioblastoma (glioblastoma multiforme, GBM). It is also one of the most lethal cancer types, with a 5-year survival rate of only 3%, compared to an average of 30% for other types of brain tumors (Ostrom et al., 2013). This demands for a wellcharacterization of molecular mechanisms of glioblastoma cells to develop treatment strategies. Therefore, it is crucial to build computer models which can mimic major characteristics of the cancerous cells (Folger et al., 2011; Hadi and Marashi, 2014; Ghaffari et al., 2015; Yizhak et al., 2015). In glioblastoma, the most significant reprogramming occurs in the metabolic machinery of the cells. Major alterations associated with cancer metabolism such as Warburg effect (Shlomi et al., 2011) are also observed in glioblastoma (DeBerardinis et al., 2007; Wolf et al., 2010). Major ATP source is via aerobic glycolysis, although TCA cycle is still slightly active according to recent reports

#### Edited by:

Adil Mardinoglu, Chalmers University of Technology, Sweden

#### Reviewed by:

Qiang Hua, East China University of Science and Technology, China Thomas Pfau, University of Luxembourg, Luxembourg Hyun Uk Kim, Korea Advanced Institute of Science and Technology (KAIST)/ Danmarks Tekniske Universitet (DTU), South Korea

\*Correspondence:

Tunahan Çakır tcakir@gtu.edu.tr

#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Neuroscience

Received: 29 January 2016 Accepted: 26 March 2016 Published: 18 April 2016

#### Citation:

Özcan E and Çakır T (2016) Reconstructed Metabolic Network Models Predict Flux-Level Metabolic Reprogramming in Glioblastoma. Front. Neurosci. 10:156. doi: 10.3389/fnins.2016.00156 (DeBerardinis et al., 2007; Wolf et al., 2010; Ru et al., 2013). Another characteristics of metabolic remodeling of glioblastoma is the uptake of glutamine, which contributes to the replenishment of TCA cycle intermediates (DeBerardinis et al., 2007). Several other metabolic alterations occur in the metabolic flux patterns in glioblastoma, mostly due to increased flux toward lipid and nucleotide synthesis to sustain growth.

The systems approach to biology and medicine led to a number of computational approaches to study network-based alterations in cells in response to perturbations. The study of cancer metabolism via computational approaches has therefore shown a sharp increase in the last decade (Ghaffari et al., 2015; Yizhak et al., 2015). Genome-scale metabolic modeling is one of the highly preferred computational methods since it allows the investigation of the cellular flux state in genome-scale by only incorporating few constraints (Kim et al., 2012; Mardinoglu et al., 2013; Bordbar et al., 2014). A generic genome-scale metabolic model includes all potential biochemical reactions to be used by the associated organism. There are computational methods which process generic metabolic models by integrating with omics data such that condition specific metabolic models are reconstructed (Blazier and Papin, 2012; Saha et al., 2014). Such methods can be divided into two groups in terms of their algorithmic approach. The first group uses context-specific omics data directly to improve the prediction of metabolic flux distributions, such as E-Flux (Colijn et al., 2009), PROM (Chandrasekaran and Price, 2010), MADE (Jensen and Papin, 2011), tFBA (van Berlo et al., 2011), TEAM (Collins et al., 2012), and RELATCH (Kim and Reed, 2012). The second group processes the data to create context-specific models from generic metabolic models, such as GIMME (Becker and Palsson, 2008), iMAT (Shlomi et al., 2008), INIT (Agren et al., 2012), AdaM (Töpfer et al., 2012), mCADRE (Wang et al., 2012), and EXAMO (Rossell et al., 2013), which can later be used for flux calculation. It was shown that the two approaches have no clear superiority over each other (Machado and Herrgård, 2014).

There are studies which integrate omics data with the metabolic models in order to reconstruct context-specific cancer metabolic models. Human metabolic reconstruction Recon1 (Duarte et al., 2007) was used to reconstruct the first generic genome-scale model of cancer, aiming to capture main metabolic functions of many cancer types using cancer gene expression data (Folger et al., 2011). Agren et al. reconstructed a generic genomescale human metabolic model, to create genome-scale active metabolic networks for 69 different cell types including 16 cancer types using tissue specific proteome data (Agren et al., 2012). Recently published reviews (Ghaffari et al., 2015; Yizhak et al., 2015) survey the studies of cancer metabolism by reconstructed metabolic model approaches and discusses the challenges such approaches face.

In this study, metabolic alteration of glioblastoma was investigated using in-silico metabolic model reconstruction approach. The genome-scale brain metabolic model (Sertbas et al., 2014) reconstructed recently by our group was first modified by adding biomass growth reaction to reflect the tumor proliferation. Afterwards, the glioblastoma gene expression data from Gene Omnibus Database (Edgar et al., 2002) were integrated with the growth-implemented brain specific metabolic model to obtain GBM-specific metabolic models. The models predict major flux-level metabolic alterations and reprogramming associated with GBM, giving consistent results with both in-vitro and in-vivo studies.

# MATERIALS AND METHODS

# Genome-Scale Brain Metabolic Network for Brain Tumors

The genome-scale brain metabolic model iMS570 (Sertbas et al., 2014; Cakir, in press) reconstructed previously by our group possesses 630 metabolic reactions in and between astrocyte and neurons, which are controlled by 570 genes. iMS570 includes the fundamental pathways such as central carbon metabolism (glycolysis, pentose phosphate pathway, TCA cycle), lipid metabolism, nucleotide metabolism, amino acid metabolism (synthesis and catabolism), the well-known glutamate-glutamine cycle, other coupling reactions between astrocytes and neurons, and neurotransmitter metabolism. In total, 42 pathways are covered by the model. iMS570 does not have a growth reaction to simulate the proliferation of brain tumors since mammalian brain cells do not grow in non-tumor states. Therefore, an extended literature survey was performed to define a growth reaction for tumor proliferation in brain (See Supplementary File 1). The modified model which can grow in-silico thanks to the included biomass growth reaction is called iMS570<sup>g</sup> . The biomass composition was defined based on brain white matter since GBM is mostly observed in this tract as the parent tissue (Bohman et al., 2010; Omuro and DeAngelis, 2013; Cuddapah et al., 2014). Brain white matter has a high composition of lipid (54.9%) and protein (39.5%) (Brady et al., 2012). The percentages of glial and neuronal cells in the white matter were reported to be 94 and 6% respectively in a recent study by using a novel method based on tagging the DNA inside the nuclei with fluorescent proteins (Azevedo et al., 2009). These values were used to define the relative contributions of astrocyte and neuron cells to the biomass reaction in the model. Free amino acid composition for human brain reported by (Banay-Schwartz et al., 1992, 1993a,b) was used as the amino acid composition of the protein pool in the biomass reaction. (See Supplementary File 1 for details on biomass composition and contribution of cell types). One characteristic of GBM cells is the altered glutamine metabolism which manifests itself as decreased glutamine production, high glutamine uptake rate and glutaminolysis (Portais et al., 1993; DeBerardinis et al., 2007). Four extra reactions denoting the glutamine uptake and glutaminolysis metabolism were also included in iMS570<sup>g</sup> in order to cover the tumor-caused alterations in glutamine metabolism (See Supplementary File 1).

# GBM Transcriptome Datasets

Lee et al. used several published GBM transcriptome datasets in addition to their own study to investigate survival differences between GBM subtypes (Lee et al., 2008). The whole transcriptome dataset is stored in the public transcriptome database, GEO (Edgar et al., 2002), under GSE13041. The dataset covers gene expression data from different microarray platforms. We focused on a subset of the data from two different platforms (GPL96, Affymetrix Human Genome U133A Array and GPL570, Affymetrix Human Genome U133 Plus 2.0 Array) and analyzed them separately to document the effect of platform type on the results. A previous study defined three distinct subtypes of GBM tumors based on clustering analysis of transcriptome data (Phillips et al., 2006): Mesenchymal (Mes), ProNeural (PN), Proliferative (Pro). Here, PN type has a better prognosis, and has a more similar gene expression profile to normal brain and neurogenesis. The other two types have poor prognosis, and show resemblance to proliferative or mesenchymal-origin cells in terms of gene expression (Phillips et al., 2006). The data from GPL96 platform was analyzed by considering this classification, which was already implemented by the authors (Lee et al., 2008). Another dataset by Mangiola et al. (2013) was also used in this study. They investigated the relation between peritumoral tissue (brain adjacent to tumor) and GBM using gene expression profiles. Normal white matter was used as a control group. The transcriptome dataset is based on GPL96 platform, and it is stored in GEO database under GSE13276. The reason behind using another dataset was to test the effect of different datasets on the bioinformatic algorithms used in this study.

In total, five different GBM transcriptome datasets were formed for the purpose of this study: Three datasets of GBM subtypes for GPL96 platform of GSE13041, a dataset of GPL570 platform for GSE13041, and a dataset from GSE13276. All GBM samples used in our study were collected from tumor biopsies of GBM patients.

# Obtaining GBM-Specific Metabolic Models

iMS570<sup>g</sup> , the growth-implemented brain specific genome-scale metabolic network, was integrated with the GBM gene expression data mentioned in the previous section to generate contextspecific GBM metabolic models and metabolic flux distributions. Two alternative methods, GIMME (Becker and Palsson, 2008) and MADE (Jensen and Papin, 2011), were applied to generate GBM metabolic models and test the effect of different algorithms on the results (**Figure 1**). Friedmann-Morvinski et al. (2012) showed that GBM can originate not only from astrocytes but also from neurons. Therefore, GBM transcriptome data were mapped to both astrocytic and neuronal reactions in iMS570<sup>g</sup> in order to generate GBM metabolic models via GIMME and MADE. While the output of MADE is a context-specific flux distribution, the output of GIMME is a context-specific model which needs to be further processed to obtain a flux distribution.

### GIMME

GIMME (Gene Inactivity Moderated by Metabolism and Expression) algorithm uses binarized gene expression data and a genome scale metabolic network to generate a contextspecific reconstruction such that the highest consistency with the available data is ensured (Becker and Palsson, 2008). All five transcriptome datasets were used as experimental softconstraints to obtain corresponding GBM-specific metabolic models using GIMME algorithm. Transcriptome data were first binarized based on a specified threshold to obtain highly and lowly expressed genes. GIMME algorithm removes the reactions which correspond to gene expression levels below the specified threshold, and the algorithm adds a removed reaction back if the metabolic model cannot achieve the desired functionality. The desired functionality was used as biomass growth reaction in iMS570<sup>g</sup> . The threshold criteria for GIMME in this study was that the threshold must not be higher than the levels of some genes known to be upregulated in GBM. These genes are HK2, PKM2, GLS, ACLY, ACC, and FASN, which take roles in glucose, glutamine or lipid metabolisms (Wolf et al., 2010; Ru et al., 2013). Considering this criteria, we chose 1/2, 1/1, and 1/3 of "the mean of transcriptome data of related GBM dataset" as the thresholds in GIMME algorithm for GSE13041 (GPL96), GSE13041 (GPL570), and GSE13276 respectively. Five percent sensitivity analysis was applied for the chosen thresholds, and no significant change was observed in the calculated flux distributions (data not shown). The genes whose expression levels are higher than the threshold were assumed highly expressed and set to "1," and the genes whose expression levels are lower than the threshold were assumed lowly expressed and set to "0," to be used as an input to GIMME. GIMME functionality of COBRA (COnstraint-Based Reconstruction and Analysis; Schellenberger et al., 2011) Toolbox was used under MATLAB (MathworksInc., Natick, MA, USA) environment to run the algorithm. The number of the removed reactions by GIMME from iMS570<sup>g</sup> was 57, 55, and 54 for Mes, PN and Pro subtypes. For GPL570 based data 48 reactions were removed whereas the number was 34 for GSE13276 dataset. Finally, five different GBM-specific metabolic models were reconstructed by GIMME. Then, flux balance analysis (FBA) (Orth et al., 2010) for maximizing biomass growth rate as primary objective function with subsequent minimization of Euclidean norm of internal fluxes was applied to five different GBM models obtained by GIMME algorithm. This dual objective function framework was shown to give better results (Cakir et al., 2007; Tarlak et al., 2014) since it ensures minimal use of enzyme resources to achieve the primary objective. GBM models corresponding to the five transcriptome datasets and the growth-implemented model are available in SBML format in Supplementary File 2.

### MADE

In addition to GIMME, we used MADE algorithm to see if there is any difference between the methods that maps transcriptome data on metabolic models to obtain conditionspecific models. MADE (Metabolic Adjustment by Differential Expression) algorithm uses the expression levels of significantly changed genes or proteins to generate a functional metabolic model that most accurately recapitulates the expression dynamics (Jensen and Papin, 2011). MADE eliminates the probable problems of arbitrary user-specified threshold, as employed by GIMME, by using statistically significant changes in gene expression measurements between two conditions to determine highly and lowly expressed genes (Blazier and Papin, 2012). Since MADE requires a control group for the analysis, only transcriptome data GSE13276 was used in the MADE-based analysis. MADE algorithm requires three inputs to generate

a context-specific flux distribution, which are a genome-scale metabolic model, fold changes and p-values of gene expression levels between the compared conditions (Jensen and Papin, 2011). Fold changes and p-values calculated by student's t-test were based on the white matter data as a control group. The objective function used by MADE algorithm for the iMS570<sup>g</sup> was biomass growth reaction, as used in GIMME. MADE algorithm was used via TIGER (Toolbox for Integrating Genome-scale Metabolism, Expression, and Regulation; Jensen et al., 2011) toolbox under MATLAB environment. MADE uses its own flux calculation algorithm which is based on mixed integer linear programming. Both GIMME and MADE were run in default settings, and GUROBI optimizer (http://www.gurobi.com) was used as a solver in both tools.

# Constraints Reflecting Physiology of Glioblastoma Multiforme (GBM)

GIMME was run with the constraints which reflect basic characteristics of GBM. The reactions defining glutamine exchange from astrocyte to neuron (r95) and glutamine release (r580) were constrained as zero due to the fact that glutamine exchange between astrocytes and neurons in healthy brain is perturbed within GBM (Marin-Valencia et al., 2012). Glycogen uptake (r575) and ketone body metabolism (r608, r609), which are used as alternative pathways in case of low activity of the glucose metabolism (Cakir et al., 2007), were also constrained to zero in our study. NH<sup>3</sup> exchange reaction (r607), which is defined only as uptake reaction in iMS570, were changed to a reversible exchange reaction to allow NH<sup>3</sup> release as observed in the GBM (DeBerardinis et al., 2007). Also a ratio, not an absolute value, of 94/6 was defined as a constraint for relative glucose, glutamine, and oxygen uptake rates of astrocytes and neurons considering the relative amounts of the two cell types in brain white matter (see below). GIMME generated GBM specific metabolic networks based on these constraints. Afterwards, corresponding flux distributions were calculated by employing more specific additional constraints as follows: Glucose, glutamine and oxygen uptake rates were fixed to the experimental flux values indicated in the study of DeBerardinis et al. (2007), which are 0.852 mmol/gDW/h, 0.080 mmol/gDW/h, and 0.272 mmol/gDW/h respectively. These uptake rates were distributed among astrocyte and neuron as 94 and 6% respectively, according to the relative amount of the cell types in white matter (Azevedo et al., 2009). Furthermore, the upper bound of the uptake rates of the amino acids other than glutamine were constrained to one tenth of the glutamine uptake rate because Yang et al. (2009) observed that GBM cells consume glutamine at a rate at least 10-fold higher than any other amino acids. Complete list of the constraints used in the models are also given in Supplementary File 1. All the above mentioned constraints were directly used as input to MADE algorithm since, in contrast to GIMME, it calculates flux distribution as a direct output. Although the contribution of neuronal reaction rates is usually considerably lower, this contribution was accounted by summing up astrocytic fluxes with neuronal counterparts to evaluate the phenotype of GBM metabolic models in the results part.

# RESULTS

The GBM metabolic models used in this study were reconstructed by integrating GBM gene expression data with the generic genome-scale brain metabolic model, iMS570<sup>g</sup> (See Materials and Methods). Five GBM metabolic models were created by GIMME, three for the comparison of Mes, PN, and Pro subtypes of GBM, one to compare the effect of different microarray platforms (The GBM metabolic model obtained using a different microarray platform but from the same dataset as GBM subtypes), one to compare the effect of different datasets (The GBM metabolic model obtained using the same microarray platform as GBM subtypes but from a different dataset). In-silico GBM phenotypes obtained via metabolic modeling were compared with the literature based experimental results. Metabolic fluxes calculated for healthy brain in resting state (Sertbas et al., 2014) were also used for comparison.

# GBM Subtypes

The GBM subtypes, Mesenchymal (Mes), ProNeural (PN), and Proliferative (Pro) GBM, are classified based on the clustering of gene expression profiling (Phillips et al., 2006; Lee et al., 2008; Verhaak et al., 2010; Huse et al., 2013). In-silico GBM subtype metabolic models for Mes, PN and Pro types were reconstructed by GIMME algorithm (Becker and Palsson, 2008; see Materials and Methods) and used in the calculation of metabolic fluxes. Key fluxes and flux ratios are presented in **Table 1**. The results in **Table 1** show high qualitative and quantitative agreements between the flux predictions and the literature results. Results reveal subtle differences between the GBM subtypes in terms of simulated metabolic flux phenotypes. However, all GBM subtype metabolic models exhibit the same behavior in terms of active flux routes. This flux routing, as shared by the three subtypes, are summarized in **Figure 2**. The simulation results depicted in the figure are in agreement with the major properties of GBM metabolic phenotypes reported in literature. The results confirm the study which reports similar metabolic characteristics for different GBM types derived from independent human tumors with different driver mutations (Marin-Valencia et al., 2012). Detailed metabolic remodeling observed in the GBM subtype metabolic models is discussed below.

# Aerobic Glycolysis and Pyruvate Branch Point

One of the most known metabolic alterations observed in cancer metabolism is related to aerobic glycolysis, which is also called Warburg effect. This phenomenon is characterized by a high rate of glucose consumption, which is mostly


TABLE 1 | GIMME-derived key fluxes and flux ratios for GBM subtype metabolic models, Mes, PN, and Pro.

Rate units of metabolic reaction fluxes and growth rates are in mmol/gDW/h and 1/h respectively. In-silico flux values and ratios are compared. Corresponding reactions for reaction IDs can be found in Supplementary File 1.

metabolized in glycolysis rather than in mitochondrial oxidative phosphorylation even in aerobic conditions, resulting in a high rate of lactate production (Warburg, 1956; Wolf et al., 2010). All in-silico GBM subtype metabolic phenotypes computed in this study exhibited the Warburg effect with a high rate of lactate production and comparatively low tricarboxylic acid (TCA) cycle activity. Flux values of the lactate production (r11+ r56) by GBM metabolic models were around 1.68 mmol/gDW/h for all GBM subtypes (**Table 1**).

Unlike initial cancer studies which reports that the glycolytic phenotype in cancer is due to a permanent impairment of mitochondrial oxidative phosphorylation (Zheng, 2012), both recent in-vitro and in-vivo studies demonstrate that oxidative metabolism in GBM is more active than thought (Maher et al., 2012; Marin-Valencia et al., 2012). <sup>13</sup>C-labeled nutrient experiments show that glucose is metabolized through pyruvate dehydrogenase rather than pyruvate carboxylase in GBM cells. Acetyl-CoA produced from pyruvate dehydrogenase reaction then enters the TCA cycle (Maher et al., 2012; Marin-Valencia et al., 2012). Our results, in agreement with the literature, show an active flux for pyruvate dehydrogenase reaction (r13+ r57), with flux values 0.073, 0.055, and 0.073 mmol/gDW/h for Mes, PN and Pro respectively. GBM subtype metabolic models also exhibit a much lower flux in pyruvate carboxylase reaction (r12) (**Table 1**) with respect to the pyruvate dehydrogenase reaction flux, which confirms the fact that glucose metabolism does not significantly contribute to anaplerosis in GBM cells (DeBerardinis et al., 2007). Acetyl-CoA produced from pyruvate dehydrogenase and oxaloacetic acid (OAA) generate citrate via citrate synthase reaction, which is the first reaction of the TCA cycle (**Figure 2**). In all GBM subtype models, citrate synthase reaction (r25+ r69) is active but possesses a very low flux value compared to the healthy brain metabolic model (**Table 1**). Furthermore, following reactions of TCA cycle, which are conversion of citrate to isocitrate (r26−27+ r70−71) and then alpha-ketoglutarate (α-KG) (r30−32+ r73−75), carried much lower flux than the citrate synthase reaction (see Supplementary File 1 for complete flux distributions of the GBM models). One of the most known features in malignant gliomas including GBM is the mutation in isocitrate dehydrogenase (IDH) gene. While wild type IDH1 and IDH2 convert isocitrate to α-KG resulting in NADPH production, mutant IDH2 and especially IDH1 convert isocitrate to 2-hydroxyglutarate known as an oncometabolite, without producing NADPH (Dunn et al., 2012). Although the GBM metabolic models do not include mutant IDH genes and the related reaction, low flux values for conversion of isocitrate to α-KG (r31−32+ r74−75) can be explained by insufficiency in the wild type isocitrate dehydrogenase genes.

### Glutaminolysis

Glutaminolysis is one of the key pathways in GBM since it provides glutamine as an alternative carbon source for TCA cycle (Wolf et al., 2010; Ru et al., 2013). Several experimental studies were performed to reveal the role of glutaminolysis in brain tumors (Wise et al., 2008; Yang et al., 2009; Chinnaiyan et al., 2012). Glutamine is not only the nitrogen source for nucleotide synthesis or maintenance of non-essential amino acid pools, but also the carbon and energy source which can replenish the TCA cycle intermediates in GBM cells (DeBerardinis et al., 2007; Maher et al., 2012; Ru et al., 2013). Glutamine is converted to glutamate by glutaminase reaction (r96+ r97) and then to α-KG by the anaplerotic glutamate dehydrogenase reaction (r89+ r90+ r92+ r93) in order to replenish the TCA cycle intermediates (see **Table 1** for flux values).

Glutamate derived α-KG is also produced by transamination reactions (DeBerardinis et al., 2007; Yang et al., 2009) while nonessential amino acids such as aspartate and alanine are produced (see Supplementary File 1 for flux values of aspartate and alanine metabolisms). After replenishing α-KG by glutaminolysis, turnover of the TCA cycle is completed by converting α-KG to succinate, fumarate, malate and OAA respectively. All GBM subtype metabolic models exhibit active and similar flux values for the conversion of α-KG to malate (r33−<sup>36</sup> and r76−79, see Supplementary File 1 for detailed flux values). Malate is both converted to OAA (r37+ r38+ r80) to complete turnover of the TCA cycle, and converted to pyruvate by malic enzyme reaction (r39+ r82), resulting in NADPH production. All GBM subtype metabolic models exhibit similar flux value (around 0.066 mmol/gDW/h) for malic enzyme reaction. In agreement with this finding, labeled glutamine experiments showed that labeled carbon was observed in lactate derived from glutamine through malic enzyme reaction (DeBerardinis et al., 2007). Ratio of the contribution of glutaminolysis and glycolysis to pyruvate pool, (r<sup>39</sup> + r82)/(r<sup>10</sup> + r55), was around 1/25 for GBM subtypes. Flux values for the conversion of the malate to OAA by malate dehydrogenase reaction (r<sup>37</sup> + r<sup>38</sup> + r80) were around 0.012 for GBM subtypes. As a result, major sources of the acetyl-CoA and OAA pool used in TCA cycle were respectively the pyruvate dehydrogenase from glycolysis and anaplerotic flux from glutaminolysis, which was found to be consistent with both in-vitro and in-vivo experiments (DeBerardinis et al., 2007; Yang et al., 2009; Maher et al., 2012; Marin-Valencia et al., 2012). The other source for OAA pool, pyruvate carboxylase, was found to have very low flux in our results, in accordance with the literature. Although, the turnover of the TCA cycle can be completed in GBM subtype metabolic models, ATP production fluxes from glycolysis were considerably higher than ATP production fluxes from oxidative phosphorylation pathway. ATP production fluxes from glycolysis (r<sup>7</sup> + r<sup>10</sup> + r<sup>52</sup> + r55) and oxidative phosphorylation (r<sup>45</sup> + r88) were around 2.7 and 0.4 mmol/gDW/h respectively, which shows the phenomena that although the turnover of TCA cycle can be completed, major energy source is aerobic glycolysis in GBM cells. The flux value of ATP production from oxidative phosphorylation for GBM subtypes were calculated to be 3-fold lower than the value calculated by the healthy brain metabolic model (1.236 mmol/gDW/h; Sertbas et al., 2014).

# Precursors for Tumor Proliferation

In addition to energy metabolism, increased aerobic glycolysis and glutaminolysis also provide macromolecule precursors required for cell proliferation in tumors (Wolf et al., 2010; Chinnaiyan et al., 2012). Ribose-5-phosphate (R5P), produced through pentose phosphate pathway (PPP), is used as a nucleotide precursor. R5P is also used in cancer diagnosis as a tumor biomarker for its excess molecular level (Iqbal and Bamezai, 2012). Flux values of the R5P isomerase reaction (r<sup>21</sup> + r65) producing R5P were higher for all GBM subtype metabolic models than healthy brain metabolic model (**Table 1**). Fatty acid synthesis relies on citrate exported from the mitochondria to cytoplasm (DeBerardinis et al., 2007; Wolf et al., 2010). The exported citrate from TCA cycle is converted to acetyl-CoA by ATP citrate lyase reaction (r<sup>28</sup> + r72), which is the precursor for fatty acid, thereby lipid synthesis (**Table 1**). All GBM subtype models had a higher flux value than healthy brain metabolic model for the acetyl-CoA carboxylase reaction (r289), which is the first committed step for the fatty acid synthesis (**Table 1**). In addition to fatty acid precursors, lipid synthesis in proliferative GBM cells requires a large amount of NADPH since it is the electron donor for fatty acid synthesis (Wolf et al., 2010; Ru et al., 2013). The sources of the NADPH in the GBM models are oxidative arm of the PPP (r<sup>17</sup> + r61), anaplerotic reaction through glutaminolysis (r<sup>89</sup> + r92) and malic enzyme reaction (r<sup>39</sup> + r82), which are all active for GBM subtype metabolic models (see Supplementary File 1 for flux values). A high enough NADPH supply by malic enzyme flux was reported for fatty acid synthesis, together with a glutaminolytic flux higher than PPP flux for NADPH generation (DeBerardinis et al., 2007). Our results report about 50% contribution by the glutaminolysis and 25% contribution by malic enzyme and PPP, in agreement with literature.

We found around 20% less growth rate in PN compared to the other subtypes (**Table 1**), which is in perfect agreement with the clinical observation that patients with PN subtype GBMs have longer survival (Lee et al., 2008; Verhaak et al., 2010). Experimental growth rates derived from doubling time (td) using the formula ln2/t<sup>d</sup> (Stensjoen et al., 2015) for GBM cells were also used to compare with in-silico derived growth rates. Growth rates of GBM subtype metabolic models are in the range obtained by both in-vitro and in-vivo experimental studies (Perego et al., 1994; Pennington et al., 2006; Wang et al., 2009; Stensjoen et al., 2015).

# Effect of Platform Difference and Dataset Difference on Simulation Results

In order to demonstrate the robustness of the results, transcriptome data from a different microarray platform (GPL570) but from the same dataset (GSE13041) and from the same platform (GPL96) but from a different dataset (GSE13276) were additionally used to derive GBM-specific metabolic models and calculate corresponding flux distributions (see Materials and Methods for details). This is an important issue to be considered to validate our results since platform or laboratory differences

may cause serious reproducibility problems in microarray experiments (Draghici et al., 2006). No sub-type differences were accounted in these calculations. The calculated fluxes are depicted in **Figure 3**.

When the flux ratios reported in **Table 1** were calculated for the new dataset from the same platform, a very similar profile to PN subtype results was observed, with the same insilico growth rate. A classification analysis of the transcriptome data of this dataset with the PN dataset via Fisher discriminant analysis method revealed that the data of the new dataset had an acceptable degree of similarity to the PN data at transcriptome level. Based on these results, it is shown that different GBM datasets give consistent results in terms of the calculated flux phenotypes (**Figure 3**).

On the other hand, the use of data from the different microarray platform resulted in slightly different quantitative results, albeit not deviating from the flux rerouting behavior depicted in **Figure 2**. The differences include a higher growth rate (0.00118 1/h) and a higher flux to lipid metabolism through Acetyl-CoA. The lactate production rate was lower (1.51 mmol/gDW/h) then the GBM subtypes, whereas pyruvate dehydrogenase reaction and citrate synthase reaction fluxes were higher than the GBM subtypes, which were 0.152 and 0.070 mmol/gDW/h respectively (**Figure 3**). Furthermore, the conversion rates of citrate to α-KG (r26−<sup>27</sup> + r70−<sup>71</sup> and r30−<sup>32</sup> + r73−75) and ATP production flux via oxidative phosphorylation (r<sup>45</sup> + r88) were higher than GBM subtypes (see Supplementary File 1 for complete flux distributions of the GBM models). This shows that the metabolic model obtained using the expression data from the different platform gives TCA cycle flux more active than other in-silico models (**Figure 3**).

# Effect of Transcriptome-Based Model-Generation Algorithms on Simulation Results

The data of GSE13276 includes also data for a reference state, which is required for MADE simulations. Therefore, this dataset was also used by MADE algorithm and compared with the results obtained by GIMME from the same dataset to enable a validation of GIMME-based flux results (see Materials and Methods). The comparison of resulting flux distributions allowed to check if there is any difference between the algorithms mapping gene expression data to the brain metabolic network, iMS570<sup>g</sup> . Growth rate, the objective function, was calculated to be the same (0.057 1/h) by both GIMME and MADE. Although there are some significantly different flux values for same reactions generated by GIMME and MADE, MADE-based metabolic model obeys the GBM metabolic remodeling depicted in **Figure 2**. For instance, citrate synthase reaction, the first reaction of the TCA cycle, is less active in MADE-based metabolic model. Citrate lyase reaction (r<sup>28</sup> + r72) producing acetyl-CoA as the precursor for fatty acid is active for all GIMME-based metabolic models; whereas flux value of this reaction is zero for MADE. TCA cycle behavior differs quantitatively for GIMME and MADE-based models. After maximizing biomass growth rate as objective function, the minimization of Euclidean norm of internal fluxes was applied in GIMME to narrow flux ranges of reactions due to alternate optima (Cakir et al., 2007; Tarlak et al., 2014), which shows more realistic results. MADE algorithm automatically generates context-specific flux distribution and we could not apply a second optimization step to minimize the Euclidean norm of internal fluxes. This can be the reason of the significantly different flux values for same reactions generated by GIMME and MADE (for detailed flux distributions see Supplementary File 1).

# DISCUSSION

This study provides GBM-specific genome scale metabolic models, derived from a brain-specific metabolic network. A growth reaction for tumor proliferation is implemented based on the lipid and protein content of white matter, where the GBM arises. By incorporating appropriate constraints from the literature, different GBM datasets were shown to predict similar metabolic flux reroutings, both validating our results and providing a proof of data consistency over transcriptome datasets. Moreover, the effect of different computational methods to incorporate transcriptome data with genome-scale models was investigated, and the two different methods, GIMME and MADE, which differ considerably in terms of dealing with the gene expression data, give similar qualitative results. Although there are several studies to apply genome-scale constraint-based modeling to cancer metabolism, an analysis of brain tumors with this approach is scarce. Only a recent study analyzed GBM with the FBA approach, by using a metabolic model with 147 genes and 12 pathways, without the incorporation of gene expression data (Bhowmick et al., 2015). Our work provides a much more comprehensive coverage of brain tumor metabolism. Correct predictions of flux distributions in glycolysis, glutaminolysis, TCA cycle and lipid metabolism discussed in the paper validate the reconstructed GBM specific models for further use of these models in future to simulate more specific metabolic patterns for GBM, or to predict drug targets.

When the constraints are directly applied to ˙IMS570<sup>g</sup> without any incorporation of transcriptome data constraint, the calculated fluxes show some basic characteristics of GBM such as high glycolysis rate and increased lactate production

# REFERENCES


and reduced TCA cycle activity, but missed to capture other basic GBM remodeling patterns such as the contribution of glutaminolysis to TCA cycle and a lower activity of oxidative phosphorylation. This meant higher contribution of glutamate and glutamine to growth for the purely constrained model, leading to an almost doubled growth rate compared to GBM subtype metabolic models. Therefore, some GBM patterns are observed solely because of the change in experimental constraints to GBM-specific values and change in the objective function (increased flux toward biomass precursors), but it is the use of measurement and gene expression constraints together that leads to a better prediction of GBM flux remodeling. This also shows the importance of incorporating gene expression data for flux calculations.

Another commonly used method for constraint-based genome-scale metabolic modeling is the sampling of solution space. The solution space is randomly uniformly sampled for a high number of flux vectors rather than searching the space for an optimum flux vector (Thiele et al., 2005; Megchelenbrink et al., 2014). The approach is especially preferred for the analysis of mammalian metabolism. Here, we re-analyzed all GIMME-derived GBM-specific metabolic models with the sampling approach by constraining the growth rate between the optimum value and 80% of the optimum. Resulting flux values and flux reroutings were similar to the values reported in **Table 1**, **Figure 2** (results not shown). Interestingly, the flux values obtained for the effect of a different microarray platform (GPL570) were more similar to the values obtained for the GBM subtypes (**Figure 3**). For example, a lower flux pentose phosphate pathway was obtained for this platform with sampling approach.

# AUTHOR CONTRIBUTIONS

TÇ conceived and designed the study, EO performed all simulations, TÇ, EO analyzed the results and wrote the manuscript.

# ACKNOWLEDGMENTS

TÇ acknowledges The Turkish Academy of Sciences - Distinguished Young Scientists Award Program (TUBA-GEBIP) for financial support.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnins. 2016.00156

cells make the human brain an isometrically scaled-up primate brain. J. Comp. Neurol. 513, 532–541. doi: 10.1002/cne.21974


adaptation in Escherichia coli. BMC Syst. Biol. 6:148. doi: 10.1186/1752-0509- 6-148


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Özcan and Çakır. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Forward Individualized Medicine from Personal Genomes to Interactomes

#### Xiang Zhang<sup>1</sup> \*, Jan A. Kuivenhoven<sup>2</sup> and Albert K. Groen1, <sup>3</sup>

*<sup>1</sup> Department of Pediatrics, Center for Liver Digestive and Metabolic Diseases, University of Groningen, University Medical Center Groningen, Groningen, Netherlands, <sup>2</sup> Section Molecular Genetics, Department of Pediatrics, University of Groningen, University Medical Center Groningen, Groningen, Netherlands, <sup>3</sup> Department of Laboratory Medicine, Center for Liver Digestive and Metabolic Diseases, University of Groningen, University Medical Center Groningen, Groningen, Netherlands*

When considering the variation in the genome, transcriptome, proteome and metabolome, and their interaction with the environment, every individual can be rightfully considered as a unique biological entity. Individualized medicine promises to take this uniqueness into account to optimize disease treatment and thereby improve health benefits for every patient. The success of individualized medicine relies on a precise understanding of the genotype-phenotype relationship. Although omics technologies advance rapidly, there are several challenges that need to be overcome: Next generation sequencing can efficiently decipher genomic sequences, epigenetic changes, and transcriptomic variation in patients, but it does not automatically indicate how or whether the identified variation will cause pathological changes. This is likely due to the inability to account for (1) the consequences of gene-gene and gene-environment interactions, and (2) (post)transcriptional as well as (post)translational processes that eventually determine the concentration of key metabolites. The technologies to accurately measure changes in these latter layers are still under development, and such measurements in humans are also mainly restricted to blood and circulating cells. Despite these challenges, it is already possible to track dynamic changes in the human interactome in healthy and diseased states by using the integration of multi-omics data. In this review, we evaluate the potential value of current major bioinformatics and systems biology-based approaches, including genome wide association studies, epigenetics, gene regulatory and protein-protein interaction networks, and genome-scale metabolic modeling. Moreover, we address the question whether integrative analysis of personal multi-omics data will help understanding of personal genotype-phenotype relationships.

Keywords: personalized medicine, interactome, gene regulatory networks (GRN), protein-protein interaction (PPI), genome-scale metabolic models, integrative genomics, network medicine

# 1. INTRODUCTION

Humans share the same genes but do not have identical DNA sequences. The latest 1000 Genomes Project reported over 84,000,000 single nucleotide polymorphisms (SNPs), 3,000,000 short insertions/deletions, and 60,000 structural variants in 2504 subjects from 26 populations, by applying whole genome sequencing as well as exome sequencing and microarray genotyping technologies (1000 Genomes Project Consortium et al., 2015). While there are large differences in the presence of both rare and common variants, it has been reported that every subject

### Edited by:

*Adil Mardinoglu, Chalmers University of Technology, Sweden*

### Reviewed by:

*Anshu Bhardwaj, Council of Scientific and Industrial Research, India Hyun Uk Kim, Korea Advanced Institute of Science and Technology, South Korea*

> \*Correspondence: *Xiang Zhang x.zhang01@umcg.nl*

#### Specialty section:

*This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology*

Received: *12 September 2015* Accepted: *16 November 2015* Published: *09 December 2015*

#### Citation:

*Zhang X, Kuivenhoven JA and Groen AK (2015) Forward Individualized Medicine from Personal Genomes to Interactomes. Front. Physiol. 6:364. doi: 10.3389/fphys.2015.00364* carries around 250–300 loss-of-function variants that lead gene products to having less or no function (1000 Genomes Project Consortium et al., 2010, 2012; UK10K Consortium et al., 2015). Nowadays, whole genome sequencing allows the determination of the entire DNA sequence of an individual, and the resulting genomic information is believed to enable prediction of disease risk and optimization of treatment outcome (Sadee, 2011). In practice, predicting disease phenotypes from genetic sequences is extremely challenging because the genotypephenotype relationship is far more complex than anticipated. A single gene can be associated with multiple disease phenotypes while a single disease phenotype can be caused by mutations in multiple genes (Barabási et al., 2011). Importantly, mutations do not have identical effects on individuals due to the individual variation in interaction between genes, proteins, metabolites and environmental factors (Barabási et al., 2011; Kathiresan and Srivastava, 2012).

The complete set of (physical) interactions between molecules, such as genes, proteins and metabolites is known as the interactome (Cusick et al., 2005). In this review, we focus on the interactome in human cells. If we consider genome sequences as stills and phenotypes as a movie, then there must be a biological system which serves as a projector. It is indeed proposed that the interactome that acts as the projector and eventually translates the phenotypic effects determined by both genotypes and environmental factors (**Figure 1**). Vidal et al. (2011), Emmert-Streib et al. (2014) proposed that most disease phenotypes may be caused by the perturbation of the interactome, in which the products of disease genes were found to interact with each other and cluster as modules (Ghiassian et al., 2015; Menche et al., 2015). These disease modules may overlap each other, explaining the shared associated genes and clinical symptoms of different diseases (Ghiassian et al., 2015; Menche et al., 2015).

To understand the projector function of the interactome, one must capture all molecular components involved in cellular functions. With the rapid development of omics technologies, it is now possible to readily profile up to 19,797 protein-coding genes, 79,795 protein-coding transcripts, 30,057 proteins, and 4229 metabolites (Psychogios et al., 2011; Harrow et al., 2012; Kim et al., 2014). Since individuality is present in the genomes, epigenomes, transcriptomes, proteomes, and metabolomes, each cell type in every human subject will have a different interactome (Feinberg et al., 2010; Montgomery and Dermitzakis, 2011; Suhre et al., 2011; Forler et al., 2014). In contrast to nonindividualized medicine, personalized medicine attempts to address such subject-specific differences with respect to diagnosis and treatment (Topol, 2014). This review aims to give an overview of bioinformatic and network modeling approaches that can be used to develop individualized medicine.

# 2. GENOME-WIDE ASSOCIATION STUDIES, EPIGENETICS AND INDIVIDUALIZED MEDICINE

Genome-wide association studies (GWAS) have identified a great number of common single nucleotide polymorphisms (SNPs) that are statistically associated with complex disease phenotypes. The National Human Genome Research Institute (NHGRI) GWAS catalog (www.genome.gov/gwastudies/) includes 1751 curated publications of 11,912 SNPs (Welter et al., 2014). Besides disease-associated SNPs, GWAS also identified SNPs associated with drug efficacy and toxicity, fueling the development of pharmacogenomics and guiding individualized therapies (Sadee, 2011; Crews et al., 2012; Low et al., 2014). The Pharmacogenomics Knowledgebase

is a complex network constituted by gene regulatory network, protein interaction network, and metabolism.

(PharmGKB, http://www.pharmgkb.org/) (Hewett et al., 2002; Altman, 2007) is a literature-based database which provides useful annotations on genes involved in pharmacokinetics (how the drug is absorbed, distributed, metabolized and eliminated) and pharmacodynamics (how the drug acts on its target and its mechanism of action). In the current release of PharmGKB, curated evidence for 1073 human genes involved in drug response is present.

Epigenetics has been shown to play a key role in the crosstalk between environment and genome, pointing toward the notion that epigenetic marks might explain in part the role of the environment in disease development (Bjornsson et al., 2004; Rivera and Ren, 2013). Major epigenetic alterations include DNA methylation, histone modification, and chromatin remodeling (Rasool et al., 2015). A total number of 127 reference human epigenomes are available on the website of the Roadmap Epigenomics Project (http://www.roadmapepigenomics.org/), including epigenetic landscapes of 111 primary cell and tissue types as well as 16 cell lines (Roadmap Epigenomics Consortium et al., 2015). Due to epigenetic modifications, cells can exhibit different phenotypes in response to various environmental factors, such as nutritional changes and oxidative stress. Feinberg (2007) defined this ability as phenotypic plasticity, whose abnormality is linked to diseases, such as cancers, neurodegenerative and autoimmune disorders (Howell et al., 2009). By integrating GWAS SNPs with epigenetic annotations, Farh et al. (2015) identified that 90% of potentially causal variants of autoimmune diseases are non-coding and 60% map to enhancers of immune cells.

In general, information deriving from GWAS (**Table 1**) and epigenetics provide possible etiological pathways rather than the exact molecular mechanisms underlying diseases. Burke and Korngiebel (2015) pointed out that although dramatic progress has been made in genomics research, there is still a gap between genomic knowledge and clinical application. To fill such gap, an accurate understanding of the genotypephenotype relationship, which is hierarchically bridged by DNA, RNA, protein, metabolite and flux, must be developed (**Figure 2**). The integrative personal omics profile (iPOP) study (Chen et al., 2012) was the first example of individualized medicine attempting to overcome the gap by combining omics data sets. Over a 14-month period which also included two viral infections (HRV: human rhinovirus and RSV: respiratory syncytial virus), dr. Michael Snyder not only profiled his whole genome, but also

TABLE 1 | Major SNP-trait association databases.


Frontiers in Physiology | www.frontiersin.org December 2015 | Volume 6 | Article 364 |

the transcriptomes of his PBMCs (Peripheral Blood Mononuclear Cells) at 20 different time points, proteomes from PBMCs and serum across 14 time points, and metabolomes of his serum sampled 17 time points, respectively. Integration of the data sets revealed the great potential of the individualized approach. In particular, the genetic variant information of dr. Snyder indicated that he is at risk for developing coronary artery disease, basal cell carcinoma, hypertryglyceridemia, and type 2 diabetes, At the same time, he was found carrying variants that are associated with response to glucose lowering drugs, rosigitazone and metformin. Interestingly, his time series measurements of transcriptome, proteome, and metabolome across healthy states, response to RSV infection, and recovery, enabled the authors to identify an alteration of the insulin signaling response following the RSV infection (Chen et al., 2012).

The iPOP study also provided us with some important insights on omics-based individualized medicine. First of all, as sequencing technologies vary considerably from each other due to sensitivity, accuracy, coverage and resolution, the measurements may contain systematic errors. Fortunately, since the human genome is constant over time, profiling with multiple DNA sequencing technologies is a way to improve the accuracy of genetic variant detection in an individual genome. As shown in the iPOP study (Chen et al., 2012), a genetic variant in the protein-coding genes can be trusted, if it is captured by the whole genome sequencing as well as whole exome sequencing. Same as above, we can also trust a genetic variant in the non protein-coding genes, if it is identified by different whole genome sequencing platforms. In contrast to the genome which is static, transcriptome, proteome, and metabolome are more dynamic and changes in their patterns represent the most valuable information for individualized medicine. To minimize systematic errors, the personal transcriptomes, proteomes, and metabolomes should be measured with standardized highthroughput methods at different time points and compared longitudinally. The longitudinal design also allows to perform statistical analysis with a single sample through applying wellestablished time-series data analysis techniques, such as Fourier spectral analysis and autocorrelation calculations (Chen et al., 2012). However, we have to admit that although the cost of sequencing technologies has dramatically decreased, sequencing with different platforms or multiple time points is unlikely to be performed for more than economic reasons only. In addition, the large volume of omics data sets will require substantial investments in data storage and management.

Topol (2014) rightfully indicated that individualized medicine needs translating large-scale omics data sets into useful knowledge. The approaches of omics data analysis can be roughly categorized as bioinformatics and network-based. Bioinformatics-based approaches often use statistical techniques to assess significant difference or association in the omics data. Their biological interpretation mainly relies on annotations in the community databases. Due to the chosen scope of this review, we are not going into details of these approaches. Network-based approaches, on the other hand, are mainly used to integrate multi-omics simultaneously and the network itself is subsequently used to explore biological insights. In general, network-based approaches first reconstruct biased or unbiased

networks in silico, and then use the reconstructed network to interpret the omics data. A biased network indicates that prior biological knowledge is incorporated, whereas an unbiased network is purely data-driven.

Network-based approaches enable us to link genotype to phenotype, and vice versa. The constructed networks can be viewed as maps, in which we can locate GWAS results and improve our understanding the roles of genetic/epigenetic alterations in disease predisposition (Califano et al., 2012; Ghiassian et al., 2015). At the same time, these maps can also help us tracking back molecular mechanisms of given clinical phenotypes. Like what has been shown by Bartel et al. (2015), the "human blood metabolome-transcriptome interface," a network constructed based on the correlation between serum metabolomes and whole blood transcriptomes of 712 subjects, can identify active pathways/modules with concentrations of blood cholesterol and triglycerides. In the next sections, we focus on three types of network-based approaches, namely gene regulatory network, protein-protein interaction networks, and genome-scale metabolic modeling and discuss them in a schematic manner: i.e., (1) definition and generation; (2) usage and results; (3) strength and weakness. We also discuss their applicability for individualized medicine.

# 3. GENE REGULATORY NETWORKS

# 3.1. What are Gene Regulatory Networks?

Thousands of gene products are produced from the human genome to support cell function and survival. The protein-coding genes can induce protein synthesis, whereas the non proteincoding genes encode noncoding RNAs (ncRNAs) as their gene products. Gene regulatory networks (GRNs) ensure proper levels of gene products present at the right time in the cell (Karlebach and Shamir, 2008). In the GRN, nodes represent the genes and edges indicate the interactions between gene products.

# 3.2. How are GRNs Generated?

Similar to gene coexpression networks, GRNs are statistically inferred from a large number of gene expression data sets. However, gene coexpression networks and GRNs are fundamentally different from each other. Pearson's correlation coefficient is used to infer coexpression networks, meaning that there is always a direct interaction for any pair of genes when their expressions are statistically correlated (Stuart et al., 2003). In contrast, GRNs are inferred mainly based on mutual information, which explicitly specifies direct or indirect interaction for each pair of genes. Mutual information defines how much information one random variable X provides about another random variable Y (Cover and Thomas, 2006). For GRNs, the random variables refer to the gene expression levels. Almost all major algorithms developed for GRN inference are mutual information-based and include ARACNe (Algorithm for the Reconstruction of accurate Cellular Networks) (Basso et al., 2005; Margolin et al., 2006), CLR (Context Likelihood of Relatdeness) (Faith et al., 2007), MRNET (Meyer et al., 2007), RN (Relevance Network) (Butte and Kohane, 2000), C3Net (Altay and Emmert-Streib, 2010a), and BC3Net (de Matos Simoes and Emmert-Streib, 2012). Different inference algorithms above were used to reconstruct human B cell GRNs and found the networks contained consistent biological information (Altay and Emmert-Streib, 2010b; de Matos Simoes et al., 2013). We refer readers to a recent review (Emmert-Streib et al., 2014) for more general concepts of GRN inference and applications. In this review, we focus on ARACNe since it is the most widely used method. ARACNe makes use of two steps to infer a genome-wide GRN (Basso et al., 2005). First, ARACNe assesses all the pair of genes by calculating their mutual information. Then, ARACNe discriminates whether the pair of genes are directly linked or separated by any other genes through applying a well-known property of mutual information called the data processing inequality (Basso et al., 2005; Cover and Thomas, 2006).

# 3.3. What are GRNs Used for?

The rationale of the GRN lies in the idea that genetic/epigenetic alterations contribute to disease phenotypes by inducing changes in a finite number of regulatory bottlenecks, i.e., transcription factors (TFs; Lefebvre et al., 2010; Califano et al., 2012). ARACNe-inferred GRNs are used for identification of the crucial TFs (also called master regulators) that affect the transition from healthy to diseased states and vice versa. The identified master regulators then serve as starting points to search for the driver genetic/epigenetic alterations upstream.

# 3.4. What Has Come Out?

Lefebvre et al. (2010) applied ARACNe to infer a human Bcell specific GRN from 254 B-cell microarray expression profiles representing 24 distinct phenotypes. The ARACNe-inferred Bcell GRN was subsequently used to identify MYB and FOXM1 as the master regulators of B-cell proliferation. Similarly, an ARACNe-inferred glioblastoma GRN was created and used by Chen et al. (2014) to identify two master regulators, C/EBPβ and C/EBPδ that are known to be involved in mesenchymal subtype of glioblastoma patients (Carro et al., 2010). Furthermore, by combining the genetic variants from the same glioblastoma patients, the authors identified that KLHL9 deletions are upstream of the two identified master regulators and act as driver mutations (Chen et al., 2014).

# 3.5. Strengths and Weaknesses

One of the major advantages of ARACNe-inferred GRNs is that with whole genome microarray or total RNA sequencing, the entire genome can actually be included in the ARACNe-inferred GRNs. Moreover, since it has been shown that the interactions inferred by the ARACNe algorithm are very likely to represent real biophysical and biochemical interactions (Basso et al., 2005; Lefebvre et al., 2010), ARACNe-inferred GRNs are suitable to explore all the possible interactions related to ncRNAs. This represents an important feature of ARACNe-inferred GRNs, as more or less 90% of the human genome is being transcribed, but only about 3% encodes protein. It is known that long noncoding RNAs (lncRNAs) can interact with DNA and proteins (Quinodoz and Guttman, 2014), and some lncRNA interactions are related to human diseases. For example, Hirata et al. (2015) reported that interaction between lncRNA MALAT1 and histone-lysine N-methyltransferase EZH2 is involved in renal cell carcinoma.

The major drawback of ARACNe is that a large number (≥100) of gene expression profile data covering a broad range of phenotypes is required to infer the target GRNs (Basso et al., 2005; Margolin et al., 2006). This is indeed necessary to explore a significant range of gene expression dynamics in order to obtain adequate mutual information for inferring GRNs (Margolin et al., 2006). Obviously, in practice it is costly and time-consuming.

# 4. PROTEIN-PROTEIN INTERACTION NETWORKS

# 4.1. What are Protein-protein Interaction Networks?

Proteins exert their function through interactions with other molecules (e.g., DNA, RNA, proteins, and metabolites). For instance, signal transduction is mediated through protein-protein interactions (PPIs), whereas gene expression (transcription factor-DNA) and metabolism (enzyme-substrate interaction) are mediated by protein-DNA and proteinmetabolite interactions, respectively (Sevimoglu and Arga, 2014). PPIs can also refer to formation of dimers, multi-protein complexes or supramolecular assemblies (e.g., actin filaments). Since some proteins are shared by different PPIs, individual PPIs are interconnected. In the PPI network, nodes represent genes whereas edges refer to physical interactions of the respective proteins.

# 4.2. How are PPI Networks Generated?

There are three main resources of generic human PPI networks. The first resource is from the literature mining. We listed six primary databases (**Table 2**) that store and combine experimentally supported PPIs from small-scale studies. The second resource is derived from large-scale yeast-two-hybrid (Y2H) screening. In 2005, the first generation of Y2H-based human PPI network, HI-I-05, was introduced and included 2700 high-quality binary PPIs among 1705 proteins (Rual et al., 2005; Stelzl et al., 2005). In 2014, the second generation of Y2H-based human PPI network, HI-II-14, was released (Rolland et al., 2014). This time 13,944 PPIs were identified among 4303 proteins. Both HI-I-05 and HI-II-14 can be downloaded (http://interactome.dfci.harvard.edu/H\_sapiens/). In addition to the Y2H system, affinity-purification mass spectrometry (AP-MS) is also developed to profile PPIs in human cells (e.g., human HEK293T, Huttlin et al., 2015). Compared to Y2H which is


mainly used to identify binary interactions between two proteins, AP-MS is more focusing on deciphering the composition of protein complexes. The third resource of the human PPI network is the computational prediction, in which machine learning algorithms are applied to calculate the likelihood of interactions between two proteins based on the known interactions in the databases (**Table 2**). STRING (Search Tool for the Retrieval of Interacting Genes, http://string-db.org/) (Snel et al., 2000) is such a web-server including known and predicted protein interactions of over 2000 organisms. In addition to STRING, databases, such as PIPs (http://www.compbio.dundee. ac.uk/www-pips/) (McDowall et al., 2009) and hPRINT (human Predicted Protein Interactome) (Elefsinioti et al., 2011) also predict PPIs without priori experimental evidence. The hPRINT results can be retrieved in STRING as well (Franceschini et al., 2013).

Human proteome studies have shown distinct proteome profiles in different cell and tissue types (Kim et al., 2014; Uhlén et al., 2015). This makes it necessary to specify PPI networks in the target cell and tissue (Schaefer et al., 2013). TissueNet database (http://netbio.bgu.ac.il/tissuenet/) provides such context-specific PPI networks for 16 human tissues (Barshir et al., 2013). A generalized way to construct such contextspecific PPI networks is introduced by Magger et al. (2012), who developed a specific algorithm integrating context-specific gene expression data (proteomics or transcriptomics). Gene expression data are used to assess the probability of PPIs in the generic PPI network. If a gene is not expressed, the algorithm can either remove the gene from the generic PPI network or reduce the weight of the interactions associated with the gene.

# 4.3. What are PPI Networks Used for?

Human PPI networks are used to identify genes, proteins and subnetworks associated with diseases (Sevimoglu and Arga, 2014). They are also used to systematically characterize PPI network perturbations associated with disease mutations. The PPI network perturbations include complete loss of gene products or alteration of PPI arrangement (Zhong et al., 2009; Sahni et al., 2013).

# 4.4. What Has Come Out?

Goehler et al. (2004) generated a PPI network for Huntington's disease by using the Y2H. From there, they identified GIT1, a G protein-coupled receptor kinase-interacting protein, which directly interacts with huntingtin and turns out to enhance huntingtin aggregation. Based on the generic human PPI network derived from HPRD (Human Protein Reference Database; Keshava Prasad et al., 2009), Jia and Zhao (2014) focused on PPI subnetworks that contain multiple genes frequently mutated in lung adenocarcinoma and melanoma patients. The results showed that the driver mutations interrupted the PPIs that are involved in signaling pathways (e.g., EGF receptor signaling pathway) and biological processes (e.g., DNA repair systems; Jia and Zhao, 2014). Based on the Y2H protein interaction assays, Sahni et al. (2015) reported that common SNPs from healthy subjects rarely affected PPIs, but around 60% of human diseaseassociated missense mutations perturbed PPIs. Furthermore, they also noticed that different mutations in the same gene can influence different PPIs.

# 4.5. Strengths and Weaknesses

Unlike the ARACNe-inferred GRNs, in which the interactions are statistically inferred from the gene expression levels, PPI networks derived from the literature or Y2H screening are experimentally supported. Therefore, perturbations in PPI networks can be used with confidence to elucidate the molecular basis of diseases as described in the examples given above.

A weakness of the PPI networks is incomplete coverage. According to the up-to-date GENCODE release 23 (), there are 19,797 protein-coding genes in the human genome. The number of genes covered by the most comprehensive human PPI network, HI-II-14 (Rolland et al., 2014), is only 3146 which suggests that there is still a long way to go. In addition, PPIs are often evaluated under unphysiological conditions, leading to false positive and negative PPIs included in generic PPI networks (Schaefer et al., 2013). Kuchaiev et al. (2009) reported that the false positive and negative rate of Y2H could be as high as 64 and 71%, respectively.

# 5. GENOME-SCALE METABOLIC MODELS

# 5.1. What are Genome-scale Metabolic Models?

Metabolites are implicated in maintenance of cellular functions and production of building blocks (e.g., purines and pyrimidines) for macromolecular biosynthesis. Computational biologists have reconstructed all metabolic reactions into one large network and name it "genome-scale metabolic model." GEMs and GSMMs are typically used as abbreviations in the literature.

# 5.2. How are GEMs Generated?

In general, GEMs are constructed by using enzyme-mediated reactions, transporters and intermediary metabolites (Bordbar et al., 2014). The first landmark studies in this field emerged in 2007 when Recon1 (Duarte et al., 2007) and EHMN (Edinburgh Human Metabolic Network) (Ma et al., 2007) were manually reconstructed based on genomic and experimental data in the literature. These two human metabolic networks were merged into the HMR (Human Metabolic Reaction) database (Agren et al., 2012). In 2010, a human hepatocyte-specific metabolic network, HepatoNet1, was reconstructed based on experimental evidence for presence of metabolic reactions in human hepatocytes (Gille et al., 2010). The experimental evidence was manually curated based on information from over 1500 scientific articles. In 2013, the continuing development of Recon1, EHMN, and HepatoNet1 leads to the release of Recon2 (Thiele et al., 2013). A year later, another reconstruction of human hepatocytespecific metabolic network, iHepatocytes2322, together with a new release of the Human Metabolic Reaction database, HMR2, were published (Mardinoglu et al., 2014).

Recon2 (Thiele et al., 2013) and HMR2 (Mardinoglu et al., 2014) represents all current knowledge of global human metabolism. Since different cell/tissue types may harbor synonymous enzymes to catalyze the same reaction and different metabolic pathways may result in the same product (Uhlén et al., 2015), it is important to reconstruct cell/tissue type specific GEMs to characterize the metabolism of target cells and tissues. For this purpose, algorithms, such as tINIT (task-driven Integrative Network Inference for Tissues) (Agren et al., 2014), GIMME (Gene Inactivity Moderated by Metabolism and Expression) (Becker and Palsson, 2008), and mCADRE (metabolic Contextspecificity Assessed by Deterministic Reaction Evaluation) (Wang et al., 2012) are used to generate cell/tissue type specific GEMs from the generic GEMs (e.g., Recon2 or HMR2). These algorithms use abundances of transcripts and proteins to estimate the probability of presence of enzymes in the generic GEMs. We refer readers to an excellent review (Machado and Herrgård, 2014) for more details on the differences between the various algorithms.

# 5.3. What are GEMs Used for?

Human GEMs, especially cell/tissue type specific GEMs, are mainly used as scaffolds to analyze transcriptomics data obtained from patient samples, in order to identify the metabolic pathways and metabolite biomarkers that are related to disease pathogenesis.

# 5.4. What Has Come Out?

Using the tINIT algorithm with proteomics and transcriptomics data of human myocytes, Väremo et al. (2015) reconstructed a myocyte-specific GEM, iMyocytes2419, which made it possible to reveal that type 2 diabetes patients show extensive transcriptional changes in reactions involved in pyruvate oxidation, branched-chain amino acid catabolism, and tetrahydroflate metabolism. Mardinoglu et al. (2014) applied iHepatocytes2322 and their previously developed Reporter Metabolite algorithm (Patil and Nielsen, 2005) to analyze transcriptomics data of patients with non-alcoholic fatty liver disease, and identified that concentrations of chondroitin and heparan sulfates may represent novel biomarkers for diagnosing non-alcoholic steatohepatitis. Similar GEM-based analyses have been performed to study diseases such as, Alzheimer's disease (Lewis et al., 2010), obesity (Mardinoglu et al., 2013), and cancer (Agren et al., 2014; Yizhak et al., 2014).

# 5.5. Strengths and Weaknesses

In our opinion, the major advantage of GEMs is that it allows to study global metabolic flux distributions. The rate of the metabolic reactions in a pathway (metabolic flux) is determined by many aspects, such as protein concentration, protein interaction (signal transduction), enzyme kinetics and metabolite concentrations (Winter and Krömer, 2013). Therefore, metabolic fluxes can be considered as the ultimate outcome of cellular regulation at different levels (Nielsen, 2003). When listing all the reactions as well as their corresponding flux values under a particular condition, one can construct a metabolic flux distribution that represents a particular cellular phenotype in detail.

Currently, <sup>13</sup>C stable isotope labeling is the most popular experimental method to measure in vivo fluxes (Blank and Ebert, 2013). By performing <sup>13</sup>C fluxomic experiments, Murphy et al. (2013) noticed that different levels of oncoprotein MYC can induce distinct metabolic flux distributions in P493-6 B cells. They showed that high MYC cells as rely more heavily on amino acids and mitochondrial oxidative metabolism than low MYC cells. <sup>13</sup>C fluxomics also revealed distinct metabolic flux distributions in different cell lines. Niklas et al. (2011) reported that human neuronal AGE1.HN cells had lower flux rates (around 2.3% of the glucose uptake flux) in the pentose phosphate pathway than other cell lines, such as HEK-293 cells (15%) and hybridoma cells (20%). These <sup>13</sup>C fluxomic studies illustrate that various biological conditions can induce distinct metabolic flux distributions.

However, <sup>13</sup>C fluxomics cannot deliver us a complete picture of flux distributions in the metabolic network, since only a small number of reactions can be measured. Here, GEMs provide a means to estimate metabolic flux distributions under different conditions relying on a limited number of exchange fluxes, i.e., fluxes of substrates entering the cells and the fluxes of metabolites that are secreted from the cells. It is beyond the scope of this review to explain the related mathematical theory, but we recommend the article by Rossell et al. (2011), in which they formulated how to compute complete set of fluxes from the exchange fluxes.

Bordel et al. (2010) introduced a random sampling method which can calculate means and standard deviations for each flux in the GEM under different experimental conditions, when a limited number of measurements of exchange fluxes are given. By integrating changes in gene expression between different conditions, metabolic reactions can be classified as either transcriptionally regulated (significant changes in both flux and gene expression levels), post-transcriptionally regulated (significant changes in gene expression levels but not flux), or metabolically regulated (significant changes in flux but not gene expression levels). This random sampling method was applied together with the adipocyte-specific GEM, iAdipocytes1809, and helped identifying the fluxes of glucose uptake, fatty acids uptake, oxidative phosphorylation, mitochondrial and peroxisomal βoxidation, fatty acid metabolism and tricarboxylic acid cycle as being differentially down regulated in obese subjects (Mardinoglu et al., 2013). Gavai et al. (2015) developed a novel algorithm called Lsei-FBA (Lesat-squares with equalities and inequalities Flux Balance Analysis), and identified the fluxes of glycolysis and oxygen uptake as being decreased in brains of Alzheimer's disease patients (29 and 46%, respectively) compared to healthy subjects. Similar to the random sampling method, Lsei-FBA also requires tissue-specific GEMs, and measurements of gene expression as well as exchange fluxes.

The second biggest advantage of GEMs is that up to now it is currently the only platform that can integrate genomics, transcriptomics, proteomics, metabolomics, and fluxomics data. Yizhak et al. (2010) integrated quantitative proteomics and metabolomics with a GEM of the human erythrocyte, and predicted metabolic flux distributions in red blood cells. The flux distribution predictions were found to be consistent with the simulations made by a detailed kinetic model of human red blood cells. Bordbar et al. (2012) analyzed transcriptomics, proteomics, and metabolomics data sets of LPS-stimulated RAW 264.7 cells with a GEM of the RAW 264.7 cell line, and identified a suppressive role for de novo nucleotide synthesis in macrophage activation.

Last but not the least, it has been shown by Uhlén et al. (2015) that the minimum requirement of generating a cell/tissue type specific GEM is a single RNA sequencing profile.

Naturally, GEMs also have their limitations. First of all, although novel metabolite biomarkers for various diseases have been predicted by using cell/tissue type specific GEMs, few of them have been validated in humans, because of either technical limitation of measuring the metabolites in question or difficulty of accessing the patient materials. Secondly, since GEMs focus on metabolic enzyme-coding genes, reactions and pathways, GEMs cannot be used to study signal transduction pathways. Lastly, GEMs do not contain detailed kinetics of enzymes and produce metabolic flux distributions only under steady state conditions.

# 6. THE FUTURE OF INDIVIDUALIZED MEDICINE

# 6.1. Role for GRNs

Regarding individualized medicine, longitudinal transcriptomics derived from cells/tissues of an individual including healthy and diseased states are the ideal resources to assemble an individualized GRN. Zoppoli et al. (2010) introduced TimeDelay-ARACNe to infer GRNs specifically from timecourse data. Such ARACNe-inferred GRN provides a personalized map, with which one can locate the genetic mutations identified in the one-dimensional genome sequences in a multi-dimensional network. By integrating gene differential expression information between healthy and diseased states, one can also identify the crucial transcription factors controlling the phenotype transition. Taken together with the network location information, one can make the most of the personal genomic information and further prioritize the damaging effect of genetic mutations.

# 6.2. Role for PPI Networks

PPI networks are proposed playing a role in buffering the impact of genetic mutations and environmental challenges (Forler et al., 2014; Garcia-Alonso et al., 2014). This opinion has been investigated by Garcia-Alonso et al. (2014), who built up a human PPI network by merging generic PPI networks derived from three public databases (BioGRID, Stark et al., 2006, IntAct, Orchard et al., 2014, and MINT, Licata et al., 2012). They used the reconstructed PPI network to study the effect of genetic variants predicted to be deleterious in the subjects participating in the 1000 Genomes Project, 252 healthy Spanish individuals, and 41 chronic lymphocytic leukemia patients. Interestingly, most of the potentially damaging genetic variants in healthy individuals were located in peripheral regions of the PPI network and did not really perturb the structure of the PPI network. However, when investigating the somatic variants that were predicted to be deleterious in chronic lymphocytic leukemia patients, they noticed that these mutations tended to be in internal regions of the PPI network. The above study indicates that PPI networks can help to identify whether genetic variants may be disrupting PPIs and hence may be important in explaining diseases.

# 6.3. Role for GEMs

GEMs have already been used successfully for individualized medicine. Agren et al. (2014) reconstructed personalized GEMs for 6 hepatocellular carcinoma patients based on proteomics data, and used these models to identify potential anticancer drug targets for the individual patients. Yizhak et al. (2014) reconstructed personalized GEMs for breast and lung cancer patients based on gene expression measurements obtained from biopsy samples. These personalized GEMs were used to predict the cancer cell growth rate, which was used to infer patient survival.

For successful individualized medicine, it should be realized that it is important to integrate information of cell/tissue type specific GEMs, in an attempt to capture whole-body metabolism. Urine, plasma, and serum are the most common samples from human subjects for diagnostic purpose (Nicholson et al., 2012). Metabolic measurements based on these samples are the results of the crosstalk of many organs and can be regarded as serving the readouts of whole-body metabolism.

Bordbar et al. (2011) build a multi-tissue GEM by integrating adipocyte, hepatocyte and myocyte-specific GEMs via a blood compartment. The assembled multi-tissue GEM was used to study the metabolic differences between non-type 2 diabetes obese and type 2 diabetes obese individuals. They reported that type 2 diabetes obese individuals lack activity in reactions catalyzed by lactate dehydrogenase, catalase and cysteine dioxygenase, comparing to the non-type 2 diabetes obese subjects. Besides integrating metabolism of different tissues and cells, the human gut microbiome is also considered important for whole-body metabolism (Mardinoglu and Nielsen, 2015). Shoaie et al. (2015) reconstructed five GEMs for five representative bacteria in the human gut, including Bacteroides thetaiotanmicron, Eubacterium rectale, Bifidobacterium adolescentis, Faecalibacterium prausnitzii, and Ruminococcus bromii. These GEMs were used to study 45 overweight and obese individuals who were subjected to an energy-restricted, high-protein diet intervention for 6 weeks. The authors reported that the diet intervention decreased the gut microbiota production of short chain fatty acids (acetate, butyrate, and propionate) and amino acids (e.g., alanine, proline and glycine etc.).

# 6.4. Concluding Remarks

Due to the central role of the interactome in cellular functions, we think that the roadmap of individualized medicine is moving from human genomes to interactomes. However, construction of a complete human interactome is extremely complex and it might take at least another decade (Menche et al., 2015). This review shows that GRNs, PPI networks, GEMs can characterize part of the interactome in cells. Integrating different type of networks may contribute to better understanding of the interactome, and ultimately realizing true individualized medicine.

# AUTHOR CONTRIBUTIONS

XZ wrote the manuscript. JK edited the manuscript. AG edited the manuscript.

# REFERENCES


# FUNDING

This work was supported by grants CVON-Genius (CVON2011- 19) and RESOLVE (FP7 305707).

whole-body systems physiology. BMC Syst. Biol. 5:180. doi: 10.1186/1752-0509- 5-180


regulation from a compendium of expression profiles. PLoS Biol. 5:e8. doi: 10.1371/journal.pbio.0050008


experiments and 13c metabolic flux analysis. J. Biosci. Bioeng. 112, 616–623. doi: 10.1016/j.jbiosc.2011.07.021


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Zhang, Kuivenhoven and Groen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Building SuperModels: emerging patient avatars for use in precision and systems medicine

Sherry-Ann Brown\*

*Division of Cardiovascular Diseases, Department of Medicine, Mayo Clinic, Rochester, MN, USA*

Keywords: computational, avatars, virtual patient, precision medicine, systems medicine, supermodels, discipulus, electronic health record

# INTRODUCTION

#### Edited by:

*Adil Mardinoglu, Chalmers University of Technology, Sweden*

#### Reviewed by:

*Xiaogang Wu, Institute for Systems Biology, USA Tunahan Cakir, Gebze Technical University, Turkey Enrico Capobianco, University of Miami, USA Gustavo Glusman, Institute for Systems Biology, USA*

> \*Correspondence: *Sherry-Ann Brown brown.sherryann@mayo.edu*

#### Specialty section:

*This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology*

Received: *05 August 2015* Accepted: *20 October 2015* Published: *06 November 2015*

#### Citation:

*Brown S-A (2015) Building SuperModels: emerging patient avatars for use in precision and systems medicine. Front. Physiol. 6:318. doi: 10.3389/fphys.2015.00318* The era is approaching when each individual can be mapped to a patient avatar—not a life-sized 3D blue form of the patient filled with physical substance (as in the movie "Avatar"), but a hologram of the patient that simulates key medical components. Patient avatars will be composed of computational models combined with various data types and analytics to form what might be called SuperModels (**Figure 1**). These SuperModels (comprehensive virtual representations of the patient, not fashion models) will be important to help realize visionary precision medicine initiatives that have recently been announced (Collins and Varmus, 2015; Nature Biotechnology, 2015).

Precision medicine tailors prevention, diagnosis, therapeutics, and prognosis for each patient (Garay and Gray, 2012; Highnam and Mittelman, 2012; Mirnezami et al., 2012). Related to precision medicine is systems medicine (Auffray et al., 2009; Capobianco, 2012; Emmert-Streib and Dehmer, 2013; Wolkenhauer, 2013), which leverages systems biology (Noble, 2008) for clinical application, with resulting data termed "systems medicine data" (Brown et al., 2015b). Systems biology studies the behavior of organisms or cells as whole systems, and uses various advances in biotechnology, including genomics, transcriptomics, proteomics, metabolomics, methylomics, microbiomics, and elucidation of cellular interaction networks by network biology (**Figure 1**, top left panel labeled X). Often, systems medicine data from these various advances can be modeled and simulated with complementary computer science, mathematics, chemistry, physics, and engineering concepts in computational biology.

A variety of fields have used computational models as virtual surrogates of specific portions of patient physiology. These individual models can be considered computational avatars of a subset of the patient's organic identity. This is akin to cancer avatars in mice, which involve mouse models mapped to individual patients, for example, by injection of tumor cells from a particular individual. These cancer avatars facilitate personalized study of the pathophysiology and response to drugs of a particular patient's cancer cells. Similarly, biomathematical or computational cancer avatars simulate the micro-environment of individualized cancer cells.

Beyond such in silico exemplars, a computational avatar may also be thought of as any finite representation of a specific portion of the patient, that harnesses computing power. This includes electronic health records (EHR), patient portals, and a variety of other precision medicine tools. However, this paper focuses on individual biomathematical models as computational avatars that can be incorporated into comprehensive patient avatars for use in precision and systems medicine. The following section describes a non-exhaustive sample of biomathematical models, including genome-scale metabolic models (GEMs) that use computational approaches to integrate omics data (Yizhak et al., 2015). These computational avatars can serve as ingredients for SuperModels, forming the first portion (X→Z) of a positive feedback loop in **Figure 1**.

FIGURE 1 | Building SuperModels for precision and systems medicine, with incorporation of computational avatars. *Input*: A feedforward loop (top left panel X connects to patient avatar Z both indirectly via bottom left panel Y and also directly, by large blue arrows) is fueled by two disparate paths. In the top left panel labeled X, various systems medicine data types, along with network biology and emerging computational models (including genome-scale metabolic models, Agren et al., 2014; Yizhak et al., 2015), provide input for patient representations such as the digital human construct being developed by The Discipulus Project. In the bottom left panel labeled Y, systems medicine data are integrated into patients' mobile health (mHealth) technologies and electronic health records (EHR) (Brown et al., 2015b). mHealth and EHR data are coupled with external knowledge [e.g., from medical societies guidelines and Food and Drug Administration (FDA) approvals] by cognitive machines such as Watson, and analytics are employed to process multi-omics (integrated personal omics profile; iPOP) and patient similarity algorithms. The path labeled Y is already in progress with EHR data, independent of digital human constructs described in the path labeled X. Paths X and Y can be bridged by locally supervised metric learning (LSML) similarity measures and similarity network fusions (SNF), for synergistic creation of SuperModels to produce results that cannot be obtained from path X or path Y alone. *Output*: High-yield predictive, preventive, and personalized data indicate patients at low/high risk for disease/adverse effect development. Individualized therapeutic plans can therefore be devised, also guided by the patient's likelihood of being a responder or non-responder to specific medications. Provision of personalized data should be in the context of systems medicine counseling, integrating genetic counseling with information about various forms of systems medicine data (Brown et al., 2015b). *Iteration*: The curved gray arrow linking output to input represents using outcome observations to iterate and refine SuperModels at all stages of development, to guide precision medicine. <sup>1</sup>Hood and Flores (2012), <sup>2</sup>Brown et al. (2015b), <sup>3</sup>Barabási and Oltvai (2004), <sup>4</sup>Schuyler et al. (2011), <sup>5</sup>Plotkin et al. (2013), <sup>6</sup>Bikson et al. (2012a), <sup>7</sup>Brown and Loew (2014), <sup>8</sup>Brown et al. (2015a), <sup>9</sup>Henderson et al. (2014), <sup>10</sup>Tortolina et al. (2012), <sup>11</sup>Stamatakos et al. (2010), <sup>12</sup>El-Kareh and Secomb (2000), <sup>13</sup>Utsler (2015), <sup>14</sup>Agren et al. (2014), <sup>15</sup>Yizhak et al. (2015), <sup>16</sup>The Discipulus Project (2013), <sup>17</sup>Kullo et al. (2013), <sup>18</sup>Steinhubl et al. (2015), <sup>19</sup>Zhang et al. (2014), <sup>20</sup>Chen et al. (2012), <sup>21</sup>Savage (2014).

# COMPUTATIONAL AVATAR EXEMPLARS

Several biomathematical models focus on understanding mechanism and prediction of pathophysiology progression, as well as delivery, efficacy, and adverse effects of therapeutics, such as deep brain stimulation or chemotherapy. Many of these computational models can replicate biomedical and/or electrophysiological properties of brain, cancer, and heart cells, personalized for each patient. Strategies and predictions for survival or for safer and more efficacious and well-timed therapy are studied and influence care of neurological and cardiovascular disorders and cancer, among others.

# Brain

Computational models of the brain have been developed, e.g., for amyotrophic lateral sclerosis (ALS). These models use known familial ALS mutations to predict functional implications or patient survival, based on mechanical properties of the mutant proteins that would be nearly impossible to produce experimentally (Schuyler et al., 2011; Plotkin et al., 2013). Thus, computational power is harnessed to predict survival and function in ALS.

Customized computational models of transcranial direct current stimulation (tDCS) have also been created (Bikson et al., 2012a,b). This non-invasive electrotherapy limits anatomic and temporal exposure to electricity in specific brain regions, minimizing side effects that could otherwise be experienced from pharmacotherapy (Bikson et al., 2012b). The tDCS computational avatars provide opportunities to individualize therapy for stroke, Parkinson's disease, and treatment-resistant depression, among others (Fregni et al., 2006; Truong et al., 2013; Tortella et al., 2015). Using these computational avatars to customize tDCS for patients at extremes of age or those with skull defects or brain damage could become standard tools to guide trials and therapy (Bikson et al., 2012a).

Another example of integrating modeling with experimental observations and clinical findings lies in spinocerebellar ataxia (SCA). The SCA modeling suite explains, interprets, or predicts experimental results in mouse models, post-mortem human brains, and peripheral blood samples from living patients (Brown and Loew, 2012, 2014). The suite is an example of the utility of computational systems biology in translational medicine (Brown et al., 2015a).

# Cancer

A number of computational avatars have also been designed for precision cancer care. These include models for colon cancer that synthesizesmathematical modeling, omics, other systems biology approaches, and pharmacologics to produce personalized molecular imprints aimed at predicting the right diagnostics and prescriptions (Tortolina et al., 2012; Henderson et al., 2014; American Cancer Society, 2015). With further study, these could be used to tailor therapy. The ContraCancrum project has moved in this direction with clinical trials illustrating the utility of computational models for lung carcinoma and other cancers (Stamatakos et al., 2010). Lung cancers account annually for ∼15 and ∼30% of all new cancer cases and deaths, respectively; colon cancer accounts for ∼10% of all new cancer cases and deaths annually (American Cancer Society, 2015). Thus, computational avatars have tremendous potential for individualizing care of cancers that account for a great proportion of morbidity and mortality in adult patients.

Adverse drug effects (ADE) on the heart or other organs limit the administration of optimal pharmacologics for cancer care (Vejpongsa and Yeh, 2014). As an example to counteract this, mathematical models predict the optimal modes of doxorubicin delivery (El-Kareh and Secomb, 2000) for breast cancer, which annually accounts for 30 and 15% of all new cancer cases and deaths, respectively (American Cancer Society, 2015). Consistent with model predictions, liposomal delivery has subsequently been studied in a number of clinical trials, which have shown superior toxicity profiles compared to standard non-liposomal delivery for breast cancer (Lao et al., 2013). Computational avatars can therefore be used to predict and hopefully prevent ADE in cancer care.

Additional avatars may focus on pancreatic and hepatocellular cancer. Recent proteomic results implicated Glypican-1 as an unparalleled near perfect non-invasive diagnostic and screening tool detection of early pancreatic cancer (Melo et al., 2015). Addition of this and other biomarkers and systems medicine data to computational avatars will help guide safe, early, and effective cancer therapy. For example, personalized computational models based on proteomics data have predicted potential drugs to treat hepatocellular cancer, one of which has already been validated experimentally (Agren et al., 2014).

# Heart

Computational avatars have also been composed for the heart. The cardioid project from the International Business Machines (IBM) Corporation uses advanced computing to compose individualized 3D models of the heart (Utsler, 2015). The system is devised to predict the risk of sudden cardiac death due to Torsade de pointes or similar arrhythmic complications. These arrhythmias are another form of cardiotoxicity, in this case related to prolongation of the QT interval (distance between the start of the wave labeled "Q" and the end of the wave labeled "T" on an electrocardiogram), induced by drugs (e.g., antibiotics). Cardioid is thought to be the world's most detailed real-time human heart simulation (Lawrence Livermore National Security, 2015). Cardioid complements prior human heart models, and expands the capabilities of avatars to geometric point-of-care. Such an achievement resulted from computational, natural, and life sciences teamwork among computational biologists, physicists, and mathematicians.

These examples of computational avatars for the brain, heart, and cancer provide evidence for established units, which can 1 day be merged (e.g., with immersive virtual environment technology for 3D animated photorealistic virtual representations of the self, Fox et al., 2009) to form SuperModels for precision and systems medicine.

# BUILDING SUPERMODELS

Computational avatars can be integrated with EHR or patient portals to build SuperModels, merging with clinical information about past medical history and diet and lifestyle habits, as well as measurements from wearable sensors, mobile health (mHealth) technologies (Steinhubl et al., 2015), and telemedicine (left section of X→Y in **Figure 1**). Some have termed a similar concept proposing patient mapping by integration of computational models with EHR information, and inviting incorporation of other biotechnological tools, as the "digital patient," "virtual patient," "medical avatar," or "patient avatar" (The Discipulus Project, 2013). Digital patient platforms, similar to Discipulus (The Discipulus Project, 2013), will use 3D scanning to produce a virtual geometric and physiologic view of the patient. MRI and CT scan results will guide reproduction of individualized anatomy, organ structure, and temporal blood flow. This paper proposes that all of this information and all of these technologies can ultimately be amalgamated with knowledge sources (such as medical society guidelines documents) and analytics to create SuperModels as the most advanced patient avatars. This forms the second portion (X→Y→Z) of a positive feedforward loop in **Figure 1**.

Systems medicine EHR data can provide input for natural-language processors, such as Watson (Savage, 2014; **Figure 1**, bottom panel labeled Y). Cognitive machines like Watson assimilate patient information to tailor medical recommendations, guidelines, and treatment options to the individual (Savage, 2014). Watson is outfitted with virtual advisors trained by medical experts, to assist with personalized risk factor identification and associated recommendations. Cognitive machines and analytics are also employed to use machine learning and natural language processing in patient similarity algorithms that yield a cohort of patients similar to a target patient, stratified by medical conditions of most concern to the engaged target patient in participatory medicine. Integrative personal omics profile (iPOP) longitudinal analysis can also combine the various omics data integrated in the EHR to uncover extensive, dynamic changes over time across healthy and diseased conditions for the target patient (Chen and Snyder, 2013), and for similar patients.

Bridging the two paths (labeled X and Y in **Figure 1**) to create SuperModels can be achieved with implementation of methodologies such as locally supervised metric learning (LSML) similarity measures and similarity network fusions (SNF). LSML and SNF facilitate personalization and prediction for risk factor profiles and computational avatars by constructing networks of patient samples for a variety of available data types, and efficiently fusing data types into one representative network that captures the full pathophysiological spectrum, respectively, while harnessing the power of complementarity in the data (Wang et al., 2014; Ng et al., 2015). Both LSML and SNF substantially outperform single data type analysis, and models created from global datasets that do not address patient similarity, respectively, while establishing integrative pathways (Wang et al., 2014; Ng et al., 2015). Synergistically not additively combining EHR integration, knowledge sources, and analytics with systems medicine data, network biology, computational models, and digital human constructs in this way produce a novel modeling perspective that can be considered the advent of SuperModels. Emergent properties of such a powerful combination are the epitome of systems medicine.

SuperModels can be interrogated to determine whether an individual might be at low or high risk for developing serious side effects to certain medications, or whether a patient is likely to respond—or not respond—to chemotherapy, for example. SuperModels will therefore in part serve as a clinical decision support tool for shared decisionmaking, supporting patient engagement in participatory medicine. Participatory medicine, which advocates for patient input and education in all phases of their individualized

# REFERENCES

Agren, R., Mardinoglu, A., Asplund, A., Kampf, C., Uhlen, M., and Nielsen, J. (2014). Identification of anticancer drugs for hepatocellular carcinoma through care, is a component of P4 (predictive, preventive, personalized, and participatory) medicine, which has been proposed as the clinical face of Systems Medicine (Hood and Flores, 2012).

# CHARTING A COURSE FORWARD

Development of SuperModels will require worldwide partnerships in academia and industry, for creation, education (e.g., https://sems.uni-rostock.de/reproducible-and-citable-dataand-models/, implementation, and troubleshooting challenges. Accordingly, large interdisciplinary consultation meetings and online fora like those of the Discipulus project initiative will become the norm, bringing together clinicians, scientists, mathematicians, bioengineers, technologists, and patients (The Discipulus Project, 2013), and may ultimately engage crowd sourcing. Difficulties, such as assuring accuracy postdata-processing (Capobianco, 2012) involving (1) merger of multi-scale noisy biased data sets with small sample sizes and large amounts of measured data (Wang et al., 2014), (2) harmonization of whole-body pharmacokinetics and pharmacodynamics with cellular network and tissue-level models (Agren et al., 2014) and diverse systems medicine data types to form digital human constructs (**Figure 1**, top left panel labeled X), (3) robust cross-validation of highly complex model findings including temporal features of more diversified disease targets (Ng et al., 2015), and (4) intercalation with analytics (**Figure 1**, bottom left panel labeled Y), along with other systems medicine challenges (Capobianco, 2012) that may be encountered when building SuperModels (**Figure 1**, X→Z and Y→Z), will most effectively be addressed through collaborative efforts. These and other principles, including ones for efficiency and cost-effectiveness, will be needed to guide the use of SuperModels in systems medicine (see companion paper in Frontiers in Genetics, Brown, in review), along with ethical and other considerations for EHR integration (Kullo et al., 2013).

# AUTHOR CONTRIBUTIONS

SB conceived of, analyzed, designed, drafted, critically revised, approved, and agreed to be accountable for this submitted work.

# ACKNOWLEDGMENTS

The author is grateful to Dr. Iftikhar Kullo of Mayo Clinic in Rochester, Minnesota for helpful discussion, to Dr. Joerg Herrmann of Mayo Clinic in Rochester, Minnesota for figure assistance, and to both for reading the manuscript.

American Cancer Society (2015). Cancer Facts and Figures 2015. Atlanta, GA: American Cancer Society.

personalized genome-scale metabolic modeling. Mol. Syst. Biol. 10, 721. doi: 10.1002/msb.145122


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Brown. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Kinetic Studies to Elucidate Impaired Metabolism of Triglyceride-rich Lipoproteins in Humans

Martin Adiels 1, 2, Adil Mardinoglu3, 4, Marja-Riitta Taskinen<sup>5</sup> and Jan Borén<sup>1</sup> \*

<sup>1</sup> Department of Molecular and Clinical Medicine/Wallenberg Laboratory, University of Gothenburg, Gothenburg, Sweden, <sup>2</sup> Health Metrics Unit, Sahlgrenska Academy, University of Gothenburg, Gothenburg, Sweden, <sup>3</sup> Department of Biology and Biological Engineering, Chalmers University of Technology, Gothenburg, Sweden, <sup>4</sup> Science for Life Laboratory, KTH - Royal Institute of Technology, Stockholm, Sweden, <sup>5</sup> Heart and Lung Centre, Helsinki University Hospital and Research Programs' Unit, Diabetes & Obesity, University of Helsinki, Helsinki, Finland

To develop novel strategies for prevention and treatment of dyslipidemia, it is essential to understand the pathophysiology of dyslipoproteinemia in humans. Lipoprotein metabolism is a complex system in which abnormal concentrations of various lipoprotein particles can result from alterations in their rates of production, conversion, and/or catabolism. Traditional methods that measure plasma lipoprotein concentrations only provide static estimates of lipoprotein metabolism and hence limited mechanistic information. By contrast, the use of tracers labeled with stable isotopes and mathematical modeling, provides us with a powerful tool for probing lipid and lipoprotein kinetics in vivo and furthering our understanding of the pathogenesis of dyslipoproteinemia.

#### Edited by:

Ioannis P. Androulakis, Rutgers, The State University of New Jersey, USA

### Reviewed by:

Litao Sun, The Scripps Research Institute, USA Hisham Bahmad, American University of Beirut, Lebanon

> \*Correspondence: Jan Borén jan.boren@wlab.gu.se

#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology

Received: 31 August 2015 Accepted: 03 November 2015 Published: 20 November 2015

#### Citation:

Adiels M, Mardinoglu A, Taskinen M-R and Borén J (2015) Kinetic Studies to Elucidate Impaired Metabolism of Triglyceride-rich Lipoproteins in Humans. Front. Physiol. 6:342. doi: 10.3389/fphys.2015.00342 Keywords: very low density lipoproteins, apoB, multocomnpartmental modeling, kinetics, stable isotopes

# INTRODUCTION

The most abundant lipids in plasma are: triglycerides, cholesterol, cholesterol esters, and phospholipids. Since lipids are water-insoluble they have to be transported in lipoprotein particles. They consist of a hydrophobic core of triglycerides and cholesterol esters, shielded from the water by a surface monolayer of phospholipids, unesterified cholesterol, and specific proteins (Mahley et al., 1984). The protein components of the lipoprotein are known as apolipoproteins (apo) (Mahley et al., 1984). The amount of lipids and proteins in the lipoprotein particles affect their density—the lower the density of a lipoprotein, the more lipids it contains relative to protein. Depending on function and hydrated density, the lipoproteins are traditionally divided into four major classes. These are chylomicrons, very low-density lipoproteins (VLDL), low-density lipoproteins (LDL), and high-density lipoprotein (HDL).

Chylomicrons and VLDL particles are the major carriers of triglycerides in the circulation. Chylomicrons are synthesized in the intestine and carry dietary lipids absorbed by the intestine. VLDL particles are synthesized by the liver. The function of these lipoprotein particles is to transport and deliver triglycerides to adipose tissue and muscles. Elevated triglycerides in plasma are associated with increased risk for cardiovascular disease (CVD).

In order to prevent and treat disturbances in metabolism of triglyceride-rich lipoproteins, it is necessary to clarify the underlying mechanism(s). The hypertriglyceridemia can either be caused by increased secretion, conversion, or catabolism of lipoprotein particles of triglyceride-rich lipoproteins. Although static s measurements of plasma lipids and functional assays may give some information, in the end, it is necessary to study the true unit of function (the integrated metabolic pathway) to understand the complexity of lipoprotein metabolism (Chan et al., 2004a,b; Adiels et al., 2008). Therefore, kinetic studies with stable isotopes are critical to explore and clarify the pathophysiology of lipid disorders in humans (Chan et al., 2004a,b; Adiels et al., 2008). The aim of this review is to illustrate how kinetic studies have furthered our understanding of impaired human lipoprotein metabolism.

# METABOLISM OF apoB-CONTAINING LIPOPROTEINS

Triglyceride-rich lipoproteins in the circulation are a mixture of chylomicrons (synthesized in the intestine) and VLDL particles (synthesized in the liver) (**Figure 1**). Each of these lipoproteins contains one molecule of apolipoprotein B (apoB). apoB is a large hydrophobic protein that remains bound to the lipoprotein particles (Segrest et al., 2001). This unique characteristic of apoB makes it possible to use apoB as a tool to trace the intravascular kinetics of the triglyceride-rich lipoproteins.

apoB is present in two different lengths; apoB100 and apoB48. The shorter form apoB48 corresponds to the amino-terminal 48% of apoB100. It is in humans exclusively synthesized in the intestine, and thus present on intestinal-derived chylomicrons and their remnants. The longer form, apoB100 is synthesized in the liver and present on VLDL, IDL, and LDL. Both apoB100 and apoB48 are coded by the same gene, but the shorter apoB48 is generated as a result of a posttranscriptional process called "RNA editing" that converts a cytidine-to-uridine (C-to-U) that generates a stop-codon and thus a truncated form of the fulllength protein. In humans, this posttranscriptional RNA editing occurs in the intestine only but in certain animals such as rodents and dogs, the process occurs also in the liver (Powell et al., 1987).

Why is apoB48 synthesized by this complex mechanism in the Intestine? The explanation is still unclear but it has been shown that apoB48-containing lipoproteins can carry more lipids than apoB100-containing lipoproteins. This is important since the chylomicrons must have capacity to rapidly and efficiently absorb large amounts of dietary lipids and a meal (Hussain, 2014). Once synthesized, chylomicrons are secreted into the lymphatic vessels until they enter into the bloodstream close to the heart. Thus, they are delivered directly to adipose tissue and muscles without first being metabolized by the liver. In the circulation, chylomicrontriglycerides are hydrolyzed by the enzyme lipoprotein lipase, which is present on the endothelial cells in the heart, muscle, and adipose tissue. The released free fatty acids are then taken up these tissues where they are stored or used for energy production. The regulation is LPL is controlled transcriptionally and posttranscriptionally. On the endothelial cells, LPL is activated by apoC-II and inhibited by apoC-III and ANGPTL4.

When the triglycerides are removed from the hydrophobic core of the chylomicrons, they shrink in size and become chylomicron remnants. The smaller chylomicron remnants are cleared from the circulation by the liver (Mahley and Ji, 1999; MacArthur et al., 2007; Williams, 2008; Williams and Chen, 2010). Recent studies have shown that chylomicron remnants are atherogenic and directly involved in atherogenesis. The explanation for this is that chylomicrons in addition to triglycerides also contain some cholesterol esters. When triglycerides are removed from the lipoprotein particle, the cholesterol esters remain. Thus, the chylomicron remnants become enriched in cholesteryl esters.

The liver secretes apoB100-containing VLDL particles. In fact, the liver produces two different forms of VLDL; larger VLDL<sup>1</sup> and smaller VLDL<sup>2</sup> particles. The triglyceride-rich VLDL<sup>1</sup> particles carry most of the plasma triglycerides and have been shown to be the major determinant for the variation of plasma triglycerides in both healthy subjects and individuals with type 2 diabetes (Hiukka et al., 2005; Boren et al., 2014; Taskinen and Boren, 2015). The formation of VLDL<sup>1</sup> is highly associated with fatty liver (Adiels et al., 2008; Boren et al., 2013). Indeed, subjects with non-alcoholic fatty live disease (NAFLD) have increased VLDL<sup>1</sup> production (Adiels et al., 2006b).

VLDL-triglycerides are hydrolyzed by the same LPLdependent mechanisms as chylomicrons, **Figure 1**. As VLDL-triglycerides are hydrolyzed, the density of the lipoprotein increases and they become intermediate density lipoproteins (IDL) and subsequently LDL. Like chylomicrons, also VLDL contains mainly triglycerides but also some cholesterol esters that become enriched in the lipoprotein particle when the triglycerides are removed. Thus, the end-product of the lipolytic cascade, LDL contains mainly cholesterol esters and is the major determinant of cholesterol in plasma. LDL particles are removed from the circulation by the LDL-receptor on the hepatocytes. The importance of the LDL-receptor for clearance of LDL particles is illustrated by the genetic disorder, familial hypercholesterolemia caused by mutation in the LDL-receptor. The disease is associated with hypercholesterolemia and premature CVD.

Patients with obesity and insulin resistance have a characteristic atherogenic dyslipidemia characterized by increased plasma triglycerides, excessive postprandial lipemia (i.e., rise in triglyceride-rich lipoproteins after eating) postprandial hyperlipidemia, and low concentrations of HDL cholesterol (Taskinen, 2003; Adiels et al., 2006a). Interestingly, these lipid disturbances are not isolated abnormalities but metabolically linked to each other (Taskinen, 2005), and they appear years before type 2 diabetes is diagnosed (Taskinen, 2003).

# PRINCIPLES OF TRACER METHODOLOGY

To simplest approach to study the kinetics of a molecule (i.e., the tracee) is to introduce the same molecule (tracer), but labeled into the system (Chan et al., 2004a,b). This process has many advantages but also obvious concerns, especially in human studies. The alternative is to introduce a labeled precursor of the molecule of interest. Ideally, the tracer should be easily detected and quantified, and not affect the system. Usually, kinetic studies are performed in a steady-state, where the rates of input and output for a given unlabeled tracee substance are equal and timeinvariant. Thus, the information provided by the tracer reflects the behavior of the tracee (Barrett et al., 2006). At various times, the amount of tracer is quantified to provide a kinetic curve. Then a mathematical model is constructed to extract all of the information contained in the kinetic curve. By fitting a model to the data, it is possible to calculate the parameters of the model that characterize the flux of molecules between kinetically

homogeneous pools of molecules. Historically, the radioactive isotopes were used as tracers, but today naturally occurring non-radioactive stable isotopes are almost exclusively used in human studies. The technical advances in mass-spectrometry technologies now permit accurate measurement of stable isotopes in smaller samples and in low concentrations (Barrett et al., 2006; Adiels et al., 2010).

A tracer can be introduced into the system either as a bolus injection or as a constant infusion given immediately after a priming dose. The bolus administration of the tracer is suitable to study kinetics of molecules with a relatively slow rate of turnover since the enrichment curves (the tracer/tracee ratios) after a bolus injection corresponds to the impulse response of the system. Furthermore, it also enables determination of newly synthesized

Adiels et al. In vivo Kinetics of VLDL

particles, since the intracellular precursor enrichment is greater at the start of the study. Primed constant infusions require a longer time to achieve a plateau of isotopic enrichment but can be appropriate for kinetic studies of molecules with rapid turnover.

# MODELING OF LIPOPROTEIN METBOLISM

Mathematical modeling enables us to better understand experimental observations. Predictions derived from model predictions can subsequently be experimentally tested. Toward this goal, known pathways are described by a set of differential equations, thereby allowing quantitative estimates to be derived. In multi-compartment models, molecules of interest move among different compartments of a system. Each compartment is assumed to be a homogeneous entity within which the entities being modeled are equivalent. Multicompartment modeling has proven to yield predictions which are as accurate as those made by physiological models, and the data required can be derived from measurements of tracer/tracee ratios of stable isotopes. Multicompartment models differ in how complex they are. Practically, their design is usually a compromise for what is practically feasible. Too simplified models may not adequately describe the kinetic heterogeneity present within the system. On the other hand, it's hard to generate experimental data for overly complex models.

The power of mathematical modeling to describe the metabolic pathways of lipid and lipoprotein metabolism was first demonstrated by Drs. Berman and Zech (Grundy et al., 1979; Zech et al., 1979). Since then, kinetic studies combined with mathematical modeling have been used to clarify the pathogenesis of impaired lipoprotein metabolism in humans linked to accelerated CVD, obesity, and insulin resistance (**Figure 2**). The methodology has also been instrumental in testing how efficiently novel drugs improve the dyslipidemia. However, it's important to emphasize that all models are based on several assumptions and simplifications. Therefore, mathematical modeling does not determine the kinetics of lipids directly; rather, they derive an indirect approximation.

# MODELS TO STUDY VLDL KINETICS

Increased levels of apoB-containing lipoproteins are the most important risk-factor for developing CVD. It's therefore clinically relevant to decipher the pathophysiology of impaired metabolism of apoB-containing lipoproteins. Until today most studies have used steady-state models to elucidate VLDL kinetics in the fasting state. However, it can be argued that these studies are not fully physiologically designed since human are in a postprandial state most of the time awake.

To study postprandial lipid metabolism it is necessary to include also intestinal derived lipoproteins. The majority of studies of intestinal lipid metabolism has been conducted in a "constant feeding" regime where small meals are served frequently to achieve a steady state of the intestinal lipoprotein secretion (Lamon-Fava et al., 2007; Sun et al., 2013; Padilla et al., 2014; Tremblay et al., 2014; Xiao et al., 2014). To date, only a few studies have focused on apoB48 metabolism following a

single meal (Wong et al., 2014a,b). No studies have yet linked hepatic and intestinal lipid metabolism in a combined model, and therefore it is not possible to study the interaction of these lipoproteins as depicted in **Figure 1**. From fat load studies it is obvious that there is a clear interplay between hepatic and intestinal lipoproteins (Adiels et al., 2012). Therefore, such models must include that apoB48- and apoB100-containing particles are cleared from the circulation by common pathways and therefore compete for clearance (Brunzell et al., 1973). During the postprandial phase, other factors need also to be considered such as insulin secretion which may affect hepatic lipoprotein secretion (Lewis et al., 1993; Adiels et al., 2007; Sorensen et al., 2011; Sondergaard et al., 2012).

In steady-state modeling of VLDL kinetics it has been discussed whether VLDL-apoB and VLDL-TG should be modeled in the same integrated model, or if they should be modeled independently. In integrated models the equation for the rate of change of an apoB100 compartment is linked to the rate of change of the corresponding triglyceride compartment size. This procedure of tying together the apoB100 and triglyceride models enhances the precision of the model as a whole. As each particle contains one single copy of apoB100, the model provides an estimate of the lipolytic rates (the loss of triglycerides per time unit), which can then be used as a physiological readout for answering study questions related to dyslipidemia. Drs. Ramakrishnan and Ginsberg recently reported that the VLDL-apoB and VLDL-TG pools in the delipidation cascade have identical rate constants despite different fates and mass distribution (Ramakrishnan and Ginsberg, 2015). These results strongly support integrated steady-state models of VLDL-apoB and VLDL-TG kinetics.

The major predictor of plasma triglycerides are VLDL1-TG. We therefore developed a multicompartment model that allows

FIGURE 3 | Compartment model of VLDL apoB and triglyceride kinetics. The model includes separate modules for leucine and glycerol. Plasma leucine kinetics is modeled using a four-compartment system that drives the synthesis and secretion of apoB into VLDL1 and VLDL2 . Plasma glycerol kinetics is modeled using a two-compartment system connected to fast and slow pathways for triglyceride (TG) synthesis. Plasma apoB and triglyceride kinetics are modeled using a four-compartment hydrolysis chain, in which the kinetics of apoB and triglyceride coupled. For each apoB compartment, there is an equivalent compartment for triglyceride. Triglycerides hydrolyzed from VLDL particles are represented by the dashed arrows, and particles lost from the plasma space are represented by the solid arrows. See Adiels et al. (2005) for additional model details.

the kinetics of triglycerides and apoB100 in VLDL<sup>1</sup> and VLDL<sup>2</sup> to be simultaneously assessed after a bolus injection of glycerol and leucine stable isotopes (**Figure 3**) (Adiels et al., 2005). Analysis of tracer/trace curves of the stable isotopes in VLDL<sup>1</sup> and VLDL<sup>2</sup> was used to derive estimates of kinetic parameters using mathematical modeling. By integrating apoB and triglycerides in the model, the triglyceride: apoB ratio of newly produced VLDL<sup>1</sup> and VLDL<sup>2</sup> particles can be estimated to follow the transfer and removal of lipids (Adiels et al., 2005).

The model can be envisioned as a two-layer model, connected at certain points and is based on the apoB model as described by Packard et al. (1995). The model consists of four parts; plasma leucine, plasma glycerol, the assembly of lipoprotein, and lipoprotein plasma kinetics. The plasma kinetics is modeled by a four-compartment hydrolysis chain, where the apoB and triglyceride kinetics are coupled at the transfer between compartments. Removal from triglyceride compartments consists of both removal of whole particles and removal of triglycerides.

# INDIVIDUAL AND POPULATION KINETICS

Traditionally, parameters for each individual are estimated individually and conclusions are made based on some statistical model applied to the model output. The model complexity ranges from one single compartment models describing VLDL-TG-kinetics (Patterson et al., 2002), to 12 compartments describing the combined apoB and TG kinetics in VLDL1 and VLDL2 (Adiels et al., 2005). The majority of published studies have used the SAAMII software (The Epsilon Group, US) (Barrett et al., 1998).

As models are becoming more complex and includes more unknown parameters, more data is needed to support the estimation of the model parameters. Individual data sets are also sensitive to loss of data and data quality, which is directly reflected in the variability in the estimated parameters.

Modern modeling techniques combine the mechanistic model (describing the system) with a statistical model (describing the populations). This is an extension to non-linear mixed effects models (NLME) where the non-linear model is the set of ordinary differential equations describing the system. Model parameters are described as random variables drawn from a distribution centered round the population mean. Using data from all individuals, population means and variances as well as individual estimates are calculated for all parameters (Beal and Sheiner, 1982).

The major gain using NLME approaches are that generally the estimated variances are smaller compared to the traditional approach, thus statistical power is greatly increased. Furthermore, such methods have shown produce better results also when data are sparse (Denti et al., 2009; Largajolli et al., 2012). Using these techniques we have recently shown that estimation of lipoprotein kinetics parameters can greatly be improved by an NLME approach (Berglund et al., 2012) (leucine subsystem) and (Berglund et al., 2015) (full system).

To estimate the day-to-day variability in VLDL kinetics and other measures, Magkos et al. repeated a kinetic study in 8 obese men on two occasions 2 months apart (Magkos et al., 2011). Using this data they calculated the sample size needed to detect differences (15–35%) between two groups during an intervention. For VLDL-TG secretion, n=15 was needed in each group to detect a difference of 25% using an unpaired study design at a study power of 80%. We compared 15 healthy controls and 15 type 2 diabetic subjects and found that using the traditional approach, n=9 was needed to detect the difference in VLDL<sup>1</sup> secretion. In contrast, using the NLME approach, a sample size of only n=5 was needed to detect the same difference (Berglund et al., 2015).

The limitations for these methods are so far that they are very computational intensive and the methods are not yet implemented in traditional software.

# PATHOPHYSIOLOGY OF DYSLIPIPIDEMIA IN OBESE SUBJECTS

To elucidate the pathophysiology of the dyslipidemia in subjects with abdominal obesity, we recently performed lipoprotein kinetic studies in 46 middle-aged well-phenotyped men and women with abdominal obesity and additional cardiometabolic risk factors to clarify determinants of plasma triglyceride concentration (Borén et al., 2015). The results are summarized in **Figure 4**. The concentration of triglycerides in plasma is determined by the balance between synthesis and removal of VLDL1-TG. Thus, dual metabolic defects are required to

FIGURE 4 | The dysregulation of the VLDL<sup>1</sup> -TG and apoAI metabolism in obese subjects. The plasma triglycerides are determined by the balance between synthesis and removal of VLDL1 -TG, illustrated by a Synthesis pathway and a Clearance pathway. The Synthesis pathway includes liver fat and total fat mass, the two independent predictors of VLDL1–triglyceride secretion rate. In the Clearance pathway, plasma concentration of the LPL inhibitor apoC-III was tightly linked with plasma TG and catabolism of VLDL1 -TG. In addition, apo-CIII was closely related to plasma TG concentration (dotted arrow). This indirect effect of apoC-III on plasma triglycerides is likely explained by effect(s) of apoC-III beyond LPL-independent pathways of triglyceride metabolism (Huard et al., 2005; Sacks et al., 2011). Increased plasma triglycerides are associated with increased HDL catabolism (Pont et al., 2002; Ooi et al., 2005; Chan et al., 2009). Thus, increased plasma triglycerides are closely associated with decreased HDL cholest erol (Verges et al., 2014). Figure modified from Borén et al. (2015). \*\*p < 0.01; \*\*\*p < 0.001.

produce hypertriglyceridemia in obese subjects (Taskinen et al., 2011). To illustrate this, two pathways are shown; a Synthesis pathway and a Clearance pathway. The Synthesis pathway include liver fat and total fat mass since these remained independent predictors of VLDL1–triglyceride secretion rate in a stepwise multivariable regression analysis. In the Clearance pathway, plasma concentration of the LPL inhibitor apoC-III was tightly linked with plasma TG and catabolism of VLDL1-TG. In addition, apo-CIII was closely related to plasma TG concentration. This indirect effect of apoC-III on plasma triglycerides is likely explained by effect(s) of apoC-III beyond LPL-independent pathways of triglyceride metabolism (Huard et al., 2005; Sacks et al., 2011).

HDL cholesterol is closely associated with diabetic dyslipidemia and abdominal obesity (Verges et al., 2014), and low HDL cholesterol is strongly associated with increased cardiovascular risk (Ninomiya et al., 2004). Pharma industry has therefore developed drugs that increase HDL cholesterol. However, clinical studies with these agents have not been successful, and the causative role of HDL is therefore questioned. This indicates that low HDL cholesterol is more a marker of an atherogenic lipoprotein profile. In vivo kinetic studies performed in abdominally obese individuals have shown that low plasma concentration of HDL cholesterol is the consequence of increased HDL catabolism (Pont et al., 2002; Ooi et al., 2005; Chan et al., 2009), and increased plasma triglycerides is closely associated with decreased HDL cholesterol (Verges et al., 2014).

Interestingly, when comparing the Synthesis pathway and the Clearance pathway, indices of catabolism were stronger predictors of plasma triglycerides than parameters of secretion. In a multivariable regression model, VLDL1-TG kinetics explained 76% of the variation in the total plasma triglycerides. Kinetic parameters of VLDL1-TG secretion explained 19–20% only of the variation in plasma triglyceride concentrations in the study subjects. Thus, in subjects with abdominal obesity and dyslipidemia, the VLDL1-triglyceride clearance is a stronger determinant of the plasma triglyceride concentration than increased secretion of VLDL<sup>1</sup> particles. This finding may support combination therapies in subjects with abdominal obesity and dyslipidemia affecting both secretion and catabolism of VLDL1- TG. The results further support apoC-III as a key target for reducing residual cardiovascular risk.

# LIVER FAT ACCUMULATION AND CVD

The observed associations between liver fat plasma triglycerides are in line with earlier studies showing that liver fat content is closely related to triglyceride secretion in different settings and populations (Chan et al., 2010; Taskinen et al., 2011) and a better predictor of triglyceride secretion than intra-abdominal fat (Fabbrini et al., 2009; Magkos et al., 2010). It is also clinically important because of the worrisome increase of non-alcoholic fatty liver disease (NAFLD), defined as hepatic fat accumulation that exceeds 5% of liver weight in individuals who do not consume significant amounts of alcohol (Neuschwander-Tetri and Caldwell, 2003; Vernon et al., 2011). Approximately 25% of adults have NAFLD, and its prevalence increases to 70– 90% among adults with obesity or type 2 diabetes (Ray, 2013). Even though NAFLD may progress to severe liver diseases, the most common cause of death in patients with NAFLD is CVD. Several epidemiological studies indicate that NAFLD is not merely a marker of CVD, but may also be actively involved in its pathogenesis (Targher et al., 2008). Thus, the finding that NAFLD is closely linked to overproduction to of VLDL1-TG that drives an atherogenic dyslipidemia characterized with hypertriglyceridemia, HDL-cholesterol, and postprandial hyperlipidemia provide an important explanation for this and indicates that novel treatments reducing liver fat accumulation might be important in preventing CVD (Boren et al., 2013).

# CONCLUSION

Lipid homeostasis is essential for human health but elevated lipid levels are a risk factor for atherosclerosis and thus can lead to symptomatic CVD. Increased lipid levels can be caused by either increased secretion of atherogenic lipoproteins, and/or impaired clearance of lipoproteins from the circulation. The use of in vivo kinetic studies using stable isotopes and mass spectrometry in combination with the development of mathematical models has been critical in advancing understanding of normal and dysregulated lipid metabolism. However, kinetic studies are timeconsuming, expensive and require a high level of expertise. Thus, they are limited to rather few research groups. Future, development will hopefully enable us to optimize the protocols and increase the statistical power, in particular when data are sparse (Denti et al., 2009; Largajolli et al., 2012).

Also, combining kinetic studies with advanced modeling, as genome-scale metabolic modeling (Mardinoglu and Nielsen, 2015; O'Brien et al., 2015), might provide even deeper understanding of the metabolic changes and clarify the underlying metabolic perturbations in the occurrence of metabolism related disorders. Genome-scale metabolic models (GEMs) are the collection of the biochemical reactions and associated protein-coding genes, and provide a scaffold for integration of the fluxomics as well as other omics data e.g., proteomics, transcriptomics, metabolomics, and lipidomics. Recently, a functional GEM for the hepatocytes in liver tissue has been reconstructed (Mardinoglu et al., 2014) and its use together with the kinetic studies may provide detailed knowledge for

# REFERENCES


understanding the relationship between the genotype-phenotype in different clinical conditions.

# ACKNOWLEDGMENTS

Work by the authors are supported by grants from the Swedish Medical Research Council, the Swedish Foundation for Strategic Research, the Sigrid Juselius Foundation and EUproject RESOLVE (FP7-HEALTH-INNOVATION-2012-1, Nr. 305707).

in abdominal obesity. Multicenter tracer kinetic study. Arterioscler. Thromb. Vasc. Biol. 35, 2218–2224. doi: 10.1161/ATVBAHA.115.305614


suppressed in obese type 2 diabetic men. Diabetologia 55, 2733–2740. doi: 10.1007/s00125-012-2624-z


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Adiels, Mardinoglu, Taskinen and Borén. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Metabolome-Wide Analysis of Stable Isotope Labeling—Is It Worth the Effort?

Daniel Weindl, Andre Wegner and Karsten Hiller\*

Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg

Keywords: stable isotope labeling, non-targeted metabolomics, Metabolic Flux Analysis (MFA), metabolism, 13C

Stable isotope assisted metabolomics techniques have emerged as a valuable tool in systems biology for metabolic flux analysis and pathway discovery (Duckwall et al., 2013; Chokkathukalam et al., 2014; Niedenführ et al., 2015). Traditionally, stable isotope labeling analyses have been highly targeted, meaning that isotopic enrichment of only a small set of metabolites has been analyzed to deduce metabolic fluxes. However, in recent years, tools for the global non-targeted detection, quantification, and computational analysis of isotopic enrichment have become available (Hiller et al., 2010, 2013; Bueschl et al., 2012; Creek et al., 2012; Chokkathukalam et al., 2013; Cho et al., 2014; Huang et al., 2014). In this article, we will discuss whether and how such novel non-targeted stable isotope labeling analyses can be applied for systems-biomedical research.

#### Edited by:

Adil Mardinoglu, Chalmers University of Technology, Sweden

#### Reviewed by:

Alessandro Giuliani, Istituto Superiore di Sanità, Italy Ranjan K. Dash, Medical College of Wisconsin, USA

> \*Correspondence: Karsten Hiller karsten.hiller@uni.lu

#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology

Received: 11 September 2015 Accepted: 06 November 2015 Published: 20 November 2015

#### Citation:

Weindl D, Wegner A and Hiller K (2015) Metabolome-Wide Analysis of Stable Isotope Labeling—Is It Worth the Effort? Front. Physiol. 6:344. doi: 10.3389/fphys.2015.00344

Pathological alterations of cellular processes usually manifest in altered metabolism. Therefore, analysis of metabolism is an ideal entry-point to diagnose or analyze diseases. Metabolic fluxes are the endpoint of cellular regulation and most likely to reflect changes on the genome, transcriptome or proteome level, and hence, are a valuable read-out for biomedical research (Wegner et al., 2015). Since intracellular metabolic fluxes cannot be measured directly, they are probed using stable isotope labeling: An isotopically enriched substrate is applied and the metabolization of this tracer leads to isotopic enrichment in downstream metabolites, depending on the underlying metabolic fluxes (Buescher et al., 2015). These enrichment patterns are analyzed by mass spectrometry (MS) or nuclear magnetic resonance (NMR) and are used to deduce metabolic fluxes (Truong et al., 2014; Young et al., 2014).

Because isotopic labeling patterns are a direct consequence of metabolic fluxes, changes in these patterns indicate metabolic flux changes (Sauer, 2006). Consequently, global analysis of labeling patterns would allow for the global detection of metabolic flux changes. Because isotopic enrichment can be deduced solely from mass spectra (Jennings and Matthews, 2005), differential flux analysis based on differences in these isotopic labeling patterns does not require prior compound identification or a model of the metabolic network. This way, a non-targeted stable isotope labeling analysis can also consider unexpected or unknown compounds and reactions, circumventing current limitations of compound identification. Such a data-driven metabolic flux analysis is perfectly suited to pinpoint disease-specific alterations of cellular metabolism as little information on the experimental outcome is required. However, even in cases where there are specific hypotheses in place, a data-driven analysis may identify previously overlooked features. As such, non-targeted stable isotope labeling analysis can function as a hypothesis generator, aiding the design of subsequent experiments. Other than e.g., <sup>13</sup>C metabolic flux analysis (13C-MFA) or flux balance analysis (FBA) which completely rely on the integrity of the underlying metabolic network and its constraints, such a non-targeted analysis does not require much biochemical a priori information and is, therefore, not biased by the limited knowledge of the metabolic network.

FIGURE 1 | Different approaches for metabolic analyses. (A) Targeted and non-targeted analysis of stable isotope labeling (SIL). Non-targeted detection of isotopic enrichment requires the measurement of an enriched and a non-enriched sample. Targeted approaches usually have higher sensitivity toward target analyses, but non-targeted analyses generally yield information on more compounds. Targeted <sup>13</sup>C-MFA provides absolute flux estimates for a given metabolic network model, but does not provide information on unknown reactions or metabolites. Non-targeted stable isotope labeling analysis can point to novel metabolic reactions and metabolites. Differential mass isotopomer distribution (MID) analysis reveals changes in metabolic fluxes in a data-driven manner. However, the mechanistic interpretation is more complex and requires in-depth biochemical knowledge. (B) Comparison of different metabolic analyses in terms of flux accuracy they provide and to which extent they allow for the discovery of novel reactions or metabolites. Non-targeted methods have a higher discovery potential, but fail to provide accurate absolute flux information.

The experimental effort for a non-targeted stable isotope labeling study is not much higher than for a conventional targeted analysis; the main difference is that a non-enriched sample is required as reference to determine isotopic enrichment in the sample of interest (**Figure 1A**; Weindl et al., 2015). In terms of analytical requirements, the non-targeted analysis is similar to label-free non-targeted metabolomics: The wide range of concentrations and chemical heterogeneity of metabolites pose problems to the metabolome-wide analysis of isotopic enrichment. Therefore, sensitive and robust analytics and powerful ion-chromatographic deconvolution and spectrum matching algorithms are required to accurately quantify also low-intensity isotopic peaks (Wegner et al., 2013). Yet, isotopic labeling patterns are much more robust towards technical variation during sample workup than metabolite levels, because the labeled and unlabeled species are subject to the same biases, thus increasing comparability of different measurements. In comparison to targeted stable isotope labeling analyses, one faces the problem of lower sensitivity because single ion (SIM) or selected reaction monitoring (SRM) cannot be used, as the analytes of interest are not known beforehand. Apart from that, a non-targeted approach would yield information on the same analytes as a targeted analysis, but also include additional compounds.

If non-targeted stable isotope labeling analysis is that powerful, why is not more widely applied? There is a great potential and the experimental effort is manageable, however, the subsequent data interpretation is much more complex. Changes in isotopic enrichment can be detected easily and do not require any biochemical knowledge (Huang et al., 2014). However, interpreting these changes in terms of underlying metabolic flux changes requires a deep understanding of cellular metabolism and is further complicated by the usually high number of compounds, many of which remain unidentified. Furthermore, although algorithms for the non-targeted and quantitative detection of isotopic enrichment in MS data have been around for a couple of years, non-targeted analysis of stable isotope labeling is still hampered by the lack of dedicated and intuitive software tools. For these reasons, only a few studies have successfully applied non-targeted stable isotope labeling analysis to provide novel biological insights, underlining the higher complexity of data interpretation and the need for more advanced tools and algorithms.

In vivo application of non-targeted stable isotope labeling analysis is technically possible, but data analysis is complicated by the metabolic complexity of a whole organism. For that reason, the main field of application would be in vitro setups with e.g., cell lines or patient-derived primary cells which are often applied in biomedical research. In such a context, non-targeted stable isotope labeling has great potential to reveal hitherto overlooked metabolic pathways or metabolic flux changes. Although such non-targeted approaches are not meant to replace targeted methods and cannot provide absolute flux information (**Figure 1B**), they provide means to get closer to the full picture.

In summary, if the aforementioned limitations can be dealt with, non-targeted stable isotope labeling analysis is worth the effort and a valuable addition to the systems biology toolbox. This rather unbiased approach can help to detect flux changes in regions of the metabolic network which were not expected to be affected or are not known yet, thereby deepening the understanding of pathobiochemical mechanisms. This understanding is the basis for the development of novel diagnostic and therapeutic methods that will impact human health. Further improvements in sensitivity and specificity of current algorithms and the development of easy-to-use software tools will help to realize the great potential of non-targeted stable isotope labeling analysis.

# REFERENCES


# AUTHOR CONTRIBUTIONS

AW and DW jointly wrote the manuscript. KH critically revised the manuscript.

# FUNDING

DW and KH are supported by the Fonds National de la Recherche (FNR) Luxembourg (ATTRACT A10/03).


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Weindl, Wegner and Hiller. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Historeceptomic Fingerprints for Drug-Like Compounds

Evgeny Shmelkov 1 †, Arsen Grigoryan1 †, James Swetnam<sup>2</sup> , Junyang Xin<sup>1</sup> , Doreen Tivon<sup>3</sup> , Sergey V. Shmelkov 4, 5 and Timothy Cardozo<sup>1</sup> \*

*<sup>1</sup> Department of Biochemistry and Molecular Pharmacology, New York University School of Medicine, New York, NY, USA, <sup>2</sup> Google Inc., Mountain View, CA, USA, <sup>3</sup> GeneCentrix Inc., New York, NY, USA, <sup>4</sup> Department of Neuroscience and Physiology, New York University School of Medicine, New York, NY, USA, <sup>5</sup> Department of Psychiatry, New York University School of Medicine, New York, NY, USA*

Most drugs exert their beneficial and adverse effects through their combined action on several different molecular targets (polypharmacology). The true molecular fingerprint of the direct action of a drug has two components: the ensemble of all the receptors upon which a drug acts and their level of expression in organs/tissues. Conversely, the fingerprint of the adverse effects of a drug may derive from its action in bystander tissues. The ensemble of targets is almost always only partially known. Here we describe an approach improving upon and integrating both components: *in silico* identification of a more comprehensive ensemble of targets for any drug weighted by the expression of those receptors in relevant tissues. Our system combines more than 300,000 experimentally determined bioactivity values from the ChEMBL database and 4.2 billion molecular docking scores. We integrated these scores with gene expression data for human receptors across a panel of human tissues to produce drug-specific tissue-receptor (historeceptomics) scores. A statistical model was designed to identify significant scores, which define an improved fingerprint representing the unique activity of any drug. These multi-dimensional historeceptomic fingerprints describe, in a novel, intuitive, and easy to interpret style, the holistic, *in vivo* picture of the mechanism of any drug's action. Valuable applications in drug discovery and personalized medicine, including the identification of molecular signatures for drugs with polypharmacologic modes of action, detection of tissue-specific adverse effects of drugs, matching molecular signatures of a disease to drugs, target identification for bioactive compounds with unknown receptors, and hypothesis generation for drug/compound phenotypes may be enabled by this approach. The system has been deployed at drugable.org for access through a user-friendly web site.

Keywords: polypharmacology, molecular docking simulation, gene expression, mechanism of drug action, drug target

# INTRODUCTION

Enormous quantities of "omics" data characterizing both normal and diseased tissues continue to accumulate, leading to the development of increasingly complex molecular biomarkers for diseases. The majority of drugs in current clinical use were discovered by phenotypic screens, leaving their precise mechanism of action unknown. Many if not most of these drugs likely act polypharmacologically (on multiple receptors simultaneously). These two trends result in

#### Edited by:

*Pierre De Meyts, De Meyts R&D Consulting, Belgium*

#### Reviewed by:

*Steven G. Gray, St. James Hospital/Trinity College Dublin, Ireland Irina Kufareva, University of California, San Diego, USA*

\*Correspondence:

*Timothy Cardozo timothy.cardozo@nyumc.org † These authors have contributed equally to this work.*

#### Specialty section:

*This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology*

Received: *18 August 2015* Accepted: *20 November 2015* Published: *18 December 2015*

#### Citation:

*Shmelkov E, Grigoryan A, Swetnam J, Xin J, Tivon D, Shmelkov SV and Cardozo T (2015) Historeceptomic Fingerprints for Drug-Like Compounds. Front. Physiol. 6:371. doi: 10.3389/fphys.2015.00371* a growing knowledge gap between the efforts to mechanistically and genomically characterize diseases on the molecular level and the chemicals used for their treatment (**Figure 1**).

Polypharmacology partly addresses this gap and has gained increasing attention in the field of drug discovery (Peters, 2011). At least some approved drugs exhibit polypharmacological signatures by interacting with multiple targets(Ashburn and Thor, 2004; Keiser et al., 2007, 2009; Mestres et al., 2009; Durrant et al., 2010; Yang et al., 2011). The identification of more of the ensemble of these targets is essential for both understanding the mechanism of drug action and predicting toxicity (Cereto-Massagué et al., 2015). Moreover, the development of compounds that rationally interact with multiple targets is appealing in the case of complex multigenic diseases, such as cancer (Knight et al., 2010) or psychiatric disorders (Metz and Hajduk, 2010; Allen and Roth, 2011; Brown and Okuno, 2012). Improved polypharmacological profiles of a drug can be identified only by a more comprehensive analysis of drug-target interactions on a proteome-wide scale (Xie et al., 2012).

In recent years, growing databases of compound-receptor bioactivities have become available (Wang et al., 2009; Sharman et al., 2011; Gaulton et al., 2012). However, the complete universe of bioactivity scores between putative or actual drugs/compounds and their receptors is still far from approachable. A number of ligand-based and structure-based in silico approaches emerged to address the off-target identification aspect of this issue (Rognan, 2013). Ligand-based approaches are based on an assumption that chemically similar structures are more likely to have similar pharmacological profiles. The idea behind the structure-based off-target identification approaches is based on inverse docking (Chen and Zhi, 2001), where a single compound is docked to multiple targets and the potential biological targets are ranked based on the docking (Chen and Zhi, 2001; Paul et al., 2004; Gao et al., 2008; Yang et al., 2009; Durrant et al., 2010; Li et al., 2010a,b; Grinter et al., 2011).

The combination of in silico target identification methods and growing databases of experimental bioactivity scores improves the feasibility of using these methods to identify a significant subset of the complete ensemble of receptors for known drugs and drug-like compounds by computational approaches. However, a gap would still remain between the polypharmacology of a drug and its pharmacodynamics, i.e., the distribution of its receptor targets in the human body. In order for the affinity of a drug for a given receptor in a given tissue to be a significant factor, the receptor has to be expressed in this tissue. For example, no matter how high the affinity of LSD is for the serotonin 5-HT2a receptor (HTR2A), this drug-target interaction is not physiologically significant in uterine tissue as HTR2A is not expressed there. The true fingerprint of drug action is the totality ("omics") of receptors for which a drug has affinity, weighted by the expression levels of these receptors in the tissues ("histos") across human body. Hence we introduced the term "historeceptomic fingerprint" for the holistic signature of drug action. Thus, here, we aim to develop a novel approach for the identification of historeceptomic fingerprints for any given drug/compound.

# METHODS

# Chemical Library

Chemical structures in Drugable were obtained from three sources: DrugBank, PubChem, and ChEMBL. 1423 approved and 4752 experimental drugs were imported from DrugBank 2.5 via the XML format release. An additional 1,138,288 compounds were imported from the SDF format release of ChEMBL 14. Additionally, PubChem compound identifiers from the SDF release were assigned to 1,006,895 DrugBank or ChEMBL compounds in Drugable on the basis of equal canonical SMILES strings as computed from RDKit (Landrum, 2008). Overall 1,141,434 unique chemical structures are represented in Drugable.

# Compound-Compound Associations

Compound-compound associations were evaluated as a chemical similarity measure between two compounds and derived as Tanimoto distance between their molecular fingerprints as implemented in the RDKit PostgreSQL extension. Briefly, given a molecule, all linear and non-linear fragments of different size were enumerated and hashed into a bit string called a fingerprint. The Tanimoto coefficient, T, for two fingerprints was calculated as the number of bits in which they differ divided by the number of non-zero bits they have in common. The Tanimoto distance was defined as 1—T. Compounds are shown in the "Similar Compounds" section of a compound page if their Tanimoto distance is less than 0.5.

# Protein Library

20,266 Human proteins were imported from the XML release of UniProt into drugable.org.

# Structure Library

3D Structures for the human proteins imported as above were obtained from two sources, the Pocketome (Abagyan and Kufareva, 2009) and ModBase (Pieper et al., 2006). 6857 experimental structures come from the Pocketome and 64,801 homology models are available from ModBase.

Consideration of receptor flexibility is crucial for structurebased drug design and the conformational ensembles of protein receptors derived from Pocketome are a practical alternative to mimic receptor flexibility. However, blindly adding certain conformations to an ensemble may be counterproductive (Rueda et al., 2010). To ensure the high quality of selected conformers, we performed retrospective virtual screening experiments and only structures with high separation power of known ligand binders from decoys were selected. Initially, for a benchmark screen, pockets on Pocketome human proteins (**Table 1**) were screened against a custom chemical library consisting of compounds solved crystallographically with several proteins and 100 random chemical decoys in order to measure the docking quality of the pockets. Having established that only the highest quality pockets could produce accurate docking scores, a subset of 6857 high quality X-ray conformations of 570 human protein targets from Pocketome was imported into the data warehouse. The 4.2 billion scores generated for Pubchem Bioassay, ChEMBL, and DrugBank compounds against these 6857 high quality pockets on 570 protein targets from the Pocketome have been



*TP and TN are the numbers of positive and negative bioactivity values available for a given pocket on a protein. Since estimation of AUC for pockets with a very small number of bioactivity values may not be fair, we also provide estimates obtained on pockets with at least 5 positive and 5 negative bioactivity values.*

integrated into the drugable.org historeceptomics system. Where there are multiple conformations for a pocket, the best score was retained. An additional complete matrix of docking scores of 4313 unique chemotypes from drugbank against ModBase homology model database is available in raw form from the authors. As a complete matrix, this data can be used for routine mathematical transformations to study symmetries and trends in the data that relate to polypharmacology. In all, docking to the largest possible set of pockets representing the druggable human genome was evaluated in this study.

# Pharmareceptomic (Bioactivity or Docking) Scores

In order to score the probability of interaction of compounds to a comprehensive set of protein targets, we used the largest available set of experimentally obtained bioactivities and in silico predicted compound-protein docking associations.

# Source of In vitro Binding Data

1,062,908 experimental compound-protein binding affinity measurements were downloaded from ChEMBL 14 PostgreSQL release. We used only binding measurements annotated with a confidence score ≥7, "assay type" field of "B," or direct protein-ligand binding, and "standard\_type" field of "Kd," "Ki," or "Potency." All compound-protein associations obtained from ChEMBL are linked to their original scientific publications in PubMed where data was available from ChEMBL.

# Source of In silico Docking Data

More than four billion compound-protein associations were derived from in silico docking experiments. The AutoDock docking program was used for the docking calculations and all the parameters were set to default values. AutoDock addresses the docking issue as a global optimization problem of an energy function, implementing an iterated local search global optimizer, using the Broyden-Fletcher-Goldfarb-Shanno criterion for local search (Trott and Olson, 2010).

Target Structure Preparation: The approach is intended to be proteome wide. Therefore, many targets with unknown biological function are expected to be available from structural genomics efforts for this approach. In order to simulate the realistic situation wherein the specific functional site on a new crystallographically resolved target receptor with unknown biological function is unknown, we rendered pockets on all receptors blindly based only on the structure coordinates and randomly selected one pocket per receptor. This pocket was then defined as the binding site for docking. Receptors were then setup by deleting the chains, heteroatoms, and prosthetic groups not involved in the binding site definition using ICM Browser (Molsoft LLC, La Jolla CA). Protein atom types were assigned, and hydrogen atoms and missing heavy atoms were added. The added or zero occupancy side chains and polar hydrogen atoms were optimized and assigned the lowest energy. Tautomeric states of histidines and the rotations of asparagine and glutamine side chain amidic groups were optimized to improve the hydrogenbonding patterns. The cognate ligands were deleted from the complexes only after hydrogen optimization. Following this receptor preparation, we used the prepare\_receptor4.py script (a part of the AutoDock Tools distribution) with default settings to convert the PDB models produced by ICM to the native PDBQT format of AutoDock.

Ligand Structure Preparation: For each compound, bond orders, tautomeric forms, stereochemistry, hydrogen atoms, and protonation states were assigned automatically by the AutoDock chemical conversion procedure. Each ligand was assigned the modified X-Score force field atom types and charges implemented in Arg. Canonical SMILES of each ligand to be screened were matched to the appropriate PubChem 3D structure (Bolton et al., 2011) to be used as a starting conformation for AutoDock docking.

After each docking simulation a stack of diverse binding poses was generated, and their respective docking scores were evaluated using the AutoDock scoring function (Trott and Olson, 2010). Three docking runs were performed for each compound-pocket pair; all binding poses accumulated after each run were merged in a single conformational stack and ranked based on their binding scores; finally, the conformation with the best docking score was retained.

# Predicted Pharmareceptomics Score (Probability) of Compound-Target Interaction

In our approach, the pharmareceptomics score is equivalent to the estimated probability that the compound will interact with the target at a physiologically significant level. For experimental bioactivities, the pharmareceptomics score is set equal to experimental affinity. For docking scores, we used the relationship between binding affinity and docking score published in Husby et al. (2015) to estimate a pharmareceptomics score from a docking score.

# Protein Target–Gene Expression Associations

Gene expression patterns of protein targets from a diverse set of tissues and cell types were derived from the "GeneAtlas U133A, gcrma" dataset (Su et al., 2004) via the BioGPS web-tool (http:// biogps.org/, accessed on 5/7/2013; Wu et al., 2009, 2013). If for a given gene, data from multiple probes/experiments were available, the mean of those values was used. For each target protein, the level of expression in each tissue was normalized with regard to its level of expression in all tissues of the dataset and projected into the Z-score.

# Data Access

The system ("Drugable") is accessible via user-friendly interface at http://drugable.org/. A flexible free-text search index is available for common names of compounds and targets, medical conditions, etc. Chemical drawer allows user to search by chemical similarity or substructure.

For example when searching Drugable by compound common name, the user is presented with compound chemical structure, compound information (Number of Hydrogen Bond Donors and Hydrogen Bond Acceptors, Number of Rotatable Bonds, Number of Rings, Walden-Crippen LogP, Indication, Pharmacology, Mechanism of Action etc.), and a table of compound-protein associations (experimentally derived and/or predicted by in silico docking experiments) available for this specific compound. The resulting table gives a list of protein targets of the compound of interest with reported or predicted affinity, including protein target UniProt accession ID, the measured activity value and type or docking score. Note that all the experimentally obtained activities are displayed in nM. In addition, a list of compounds that are chemically similar to the compound of interest is also presented. Furthermore, tissue-specific levels of expression for all genes, correspond to the protein targets of the compound of interest, are presented as a heat map.

Alternatively, a user may want to search for a particular protein of interest. In this case, the user is presented with details of the protein target, such as X-ray structure (if available), protein name synonyms, gene names, organism this protein belongs to, and UniProt accession ID.

Furthermore, users may search for a medical condition of interest. In this case user is presented with a list of drugs/druglike compounds as well as protein targets associated with this medical condition.

# RESULTS

# Generation of Bioactivity Scores

First, we generated bioactivity probability scores for the compound-receptor pairs by executing the largest computational molecular docking reported to date (see Section Methods). A benchmark docking screen was performed against 3D structural models of human proteins (**Table 1**). The mean area under the receiver operating curve (AUC) for benchmark docking was 0.59 (with about 23% of structures having separation power above 0.7) when performed on 3D structural models from Pocketome (Kufareva et al., 2012), but only 0.53 (with 8.5% of structures above AUC of 0.7) on ModBase (Pieper et al., 2006) homology models proteins (**Table 2** and Supplementary Table 1).

TABLE 2 | Number of receptors from the benchmark study with AUC above a certain threshold.


This result suggests that only the docking scores achieved with the highest quality Pocketome pockets should be included in our "omics" set of compound-receptor scores, which are used to predict mechanistic signatures solely from chemotype. The Pocketome currently includes 6857 pockets derived from high quality crystallographic structures of 570 target human proteins. Therefore, for our "omics" set we docked over 600,000 unique non-overlapping chemical structures from PubChem Bioassay, ChEMBL, and DrugBank against these 6857 pockets for a total of 4.2 billion pairwise docking scores between compounds and targets. These "omics" in silico docking scores together with the compound–receptor affinities obtained experimentally, constitute the bioactivity scores data set (**Figure 2**), which comprise a significant fraction of the druggable targets encoded in the human genome, by one estimate to be around 4000 targets (Reardon, 2013).

## Generation of Historeceptomic Scores

To address the issue of physiological significance of drug targets detected in the first step, we endeavored to calculate a tissue-specific (historeceptomic) compound-receptor score (**Figure 3A**). Tissue-specific gene expression data on protein targets were obtained from the BioGPS database. The level of expression of each receptor in every tissue was normalized with regard to its expression level in all tissues of the dataset by calculating its standard score (Z-score, see Section Methods). Each compound-receptor association in each tissue was scored by integrating their bioactivity with the receptor expression in a given tissue as follows:

$$Hs = -\log\_{10} Ps \times Z,$$

where Hs is a historeceptomic score, Ps is a bioactivity score, and Z is the gene-expression Z-score.

By this method, for any given drug/compound, thousands of historeceptomic scores can be generated, but only a tiny fraction of these, which measure the probability that the compound will affect the receptor in a physiologically significant way, are important. The average drug may have hundreds of low affinity receptors, resulting in a set of scores numbering in the tens of thousands across all tissues in the human body. To identify the physiologically significant compound-receptor interactions out of the large number of all on-/off-target interactions of a given compound, we used the generalized extreme Studentized deviate test as a statistical novelty detection approach using the α = 0.0001 level of significance (**Figure 3B**). Statistically significant historeceptomic scores of a given drug/compound form its historeceptomic fingerprint.

Fingerprints were pre-calculated for all known drugs into an integrated system suitable for searching with any chemical structure to find its historeceptomic fingerprint. The system includes the 4.2 billion docking scores with experimental affinity scores in a graph linking drugs/compounds to protein targets in order to maximize the sensitivity of target detection for any drug.

## Illustrative Use Case

Historeceptomics fingerprints may specifically localize in vivo significant mechanisms of action of a polypharmacologic drug, translating purely molecular data into a clinically interpretable profile. An example is shown in **Figure 3A**. Lysergic acid diethylamide (LSD) is a hallucinogenic drug in humans, which makes it difficult to study in animal models, as many

hallucinations are only represented internally and can only be communicated verbally. We calculated the historeceptomics profile for LSD. In this case, the inputs into our system were only molecular in nature: the affinity scores and the expression data. We did not use docking in this example. Our historeceptomics approach identified the 5HT2A receptor in the prefrontal cortex (PFC) as the most significant of tissue-target pair associated with the phenotype induced by LSD. Independently, we analyzed the preclinical and clinical literature on LSD targets, which is exclusively non-molecular data. The textbook and literature consensus from animal neuroperturbation studies, pharmacologic studies and clinical neuroimaging is that 5HT2A is the primary molecular target of LSD, and that, specifically, its activity in the PFC is responsible for its effects. Thus, there are many non-molecular clinical and translational papers in the literature, none of which were input to our system, that clearly establish 5HT2A specifically in the PFC not only as a key pathway for LSD psychosis, but also as the epicenter of the very similar psychoses seen in human schizophrenia (Arvanov et al., 1999; Vollenweider and Geyer, 2001; Muschamp et al., 2004; Nichols, 2004). The historeceptomics approach predicts this finding independently of animal or clinical studies.

# DISCUSSION

This report takes on the two major challenges of precisely describing the holistic pharmacodynamics of drugs. First, we expanded the graph of experimental scores linking drugs/compounds to protein targets, which has been used in prior methods such as SEA (Keiser et al., 2007), to include the data from the largest computational molecular docking of compounds to protein pockets yet reported. This should increase the sensitivity of target detection. Second, we addressed, for the first time, the systematic integration of bioactivity/docking scores between drugs/compounds and proteins with the expression patterns of those proteins in human tissues, thus mapping the pharmacology of drugs into human physiologic space.

The integration of bioactivity/docking scores of compoundreceptors with the expression patterns of those receptors in human tissues increases the specificity of the results by eliminating noise and selecting only physiologically significant drug-target interactions. Thus, although for many models/pockets the docking scores correlate only moderately with affinity due to the limited ability to take induced fit into account, this lack of specificity is abrogated by our integration of the gene expression such that many false positives are likely to be culled. While sensitivity is low, it can be steadily improved from our pioneering prototype by (1) improved binding site (pocket) selection methods and (2) natural growth and improved curation of the crystallographic and bioactivity databases.

There are 20,198 reviewed human proteins in UniProt, of which 4300 have human crystal structures in the PDB (21.3% of total). An additional 20–30% of these can likely be modeled reliably by homology. Thus, up to 50% of the "proteome" might already be surveyed by docking. Estimates of the druggable genome range from 8 to 12 thousand targets. The existing structures are probably highly enriched in these targets so, one can speculate that 40–50% of the druggable genome is already accessible by docking. These are highly speculative estimates, but since the number of crystal structures and the power of computation is growing rapidly, it is not unreasonable to speculate that a low resolution representation of the majority of the druggable genome could be available for docking soon.

The system has been deployed for access through a userfriendly web site: drugable.org. For compounds resulting from phenotype screens, where their mechanism of action is not known, searching the site can identify possible mechanisms of action. Similarly, where the tissue pattern of a disease is known, drug activity detected by our approach in tissues not included in the pattern could be suggestive of the mechanism of the adverse effects of a drug. Since the historeceptomic fingerprints contain both a specific pattern of targets and a specific pattern of tissues, they could potentially be matched to complex biomarkers of disease derived from exhaustive molecular profiling, which can have a similar gene-tissue signature. Our novel approach thus potentially fills a currently

# REFERENCES


existing gap between burgeoning "omics" data and drugs/druglike compounds (**Figure 1**).

# AUTHOR CONTRIBUTIONS

ES, AG, JS designed the study, performed experiments, analyzed data, and assisted in writing the manuscript; JX and DT performed experiments and analyzed data; SS. designed the study, analyzed data, and wrote the manuscript; TC conceived and designed the study, analyzed data, and wrote the manuscript.

# FUNDING

This work was supported by an American Recovery and Reinvestment Act grant from the National Library of Medicine (RC LM010994 to TC). Additional support was provided by Google, Inc. via their Exacycle for Visiting Faculty award program (to JS and TC).

# ACKNOWLEDGMENTS

We thank Ruben Abagyan, Maxim Totrov, Yuval Kluger, Fabio Parisi and Francesco Strini for helpful discussions. We thank the Visiting Faculty for Exacycle program for support of the work of JS and Daniel Belov, Chris van Arsdale, David Konerding and Daniel Meredith of Google Inc. for technical assistance.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys. 2015.00371

Proteins 43, 217–226. doi: 10.1002/1097-0134(20010501)43:2<217::AID-PROT1032>3.0.CO;2-G


visualization of pharmacological data. Nucleic Acids Res. 39, D534–D538. doi: 10.1093/nar/gkq1062


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Timothy Cardozo and Sergey V. Shmelkov are shareholders in Genecentrix Inc., a company founded after this study was completed, but based, in part, on this research. No funding was received from Genecentrix Inc. for the present study.

Copyright © 2015 Shmelkov, Grigoryan, Swetnam, Xin, Tivon, Shmelkov and Cardozo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.