# SYSTEMS BIOLOGY AND BIOINFORMATICS IN GASTROENTEROLOGY AND HEPATOLOGY

EDITED BY : Steven Dooley, Kai Breuhahn and Andreas Teufel PUBLISHED IN : Frontiers in Physiology

### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-349-4 DOI 10.3389/978-2-88963-349-4

### 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<EMAIL_ADDRESS>
# SYSTEMS BIOLOGY AND BIOINFORMATICS IN GASTROENTEROLOGY AND HEPATOLOGY

Topic Editors: Steven Dooley, University of Heidelberg, Germany Kai Breuhahn, Heidelberg University, Germany Andreas Teufel, Medical Faculty Mannheim, University of Heidelberg, Germany

Citation: Dooley, S., Breuhahn, K., Teufel, A., eds. (2020). Systems Biology and Bioinformatics in Gastroenterology and Hepatology. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-349-4

# Table of Contents

*05 Editorial: Systems Biology and Bioinformatics in Gastroenterology and Hepatology*

Peter L. M. Jansen, Kai Breuhahn, Andreas Teufel and Steven Dooley

*09 The Potential for Circulating Tumor Cells in Pancreatic Cancer Management*

Michael Pimienta, Mouad Edderkaoui, Ruoxiang Wang and Stephen Pandol

Federico Pinna, Michaela Bissinger, Katharina Beuke, Nicolas Huber, Thomas Longerich, Ursula Kummer, Peter Schirmacher, Sven Sahle and Kai Breuhahn

*60 Model Based Targeting of IL-6-Induced Inflammatory Responses in Cultured Primary Hepatocytes to Improve Application of the JAK Inhibitor Ruxolitinib*

Svantje Sobotta, Andreas Raue, Xiaoyun Huang, Joep Vanlier, Anja Jünger, Sebastian Bohl, Ute Albrecht, Maximilian J. Hahnel, Stephanie Wolf, Nikola S. Mueller, Lorenza A. D'Alessandro, Stephanie Mueller-Bohl, Martin E. Boehm, Philippe Lucarelli, Sandra Bonefas, Georg Damm, Daniel Seehofer, Wolf D. Lehmann, Stefan Rose-John, Frank van der Hoeven, Norbert Gretz, Fabian J. Theis, Christian Ehlting, Johannes G. Bode, Jens Timmer, Marcel Schilling and Ursula Klingmüller

Juan Tang, Ping Hu, Yansen Li, Tin-Tin Win-Shwe and Chunmei Li

*150* LiverSex *Computational Model: Sexual Aspects in Hepatic Metabolism and Abnormalities*

Tanja Cvitanović Tomaš, Žiga Urlep, Miha Moškon, Miha Mraz and Damjana Rozman

*162 Understanding the Representative Gut Microbiota Dysbiosis in Metformin-Treated Type 2 Diabetes Patients Using Genome-Scale Metabolic Modeling*

Dorines Rosario, Rui Benfeitas, Gholamreza Bidkhori, Cheng Zhang, Mathias Uhlen, Saeed Shoaie and Adil Mardinoglu

*176 The Diurnal Timing of Starvation Differently Impacts Murine Hepatic Gene Expression and Lipid Metabolism – A Systems Biology Analysis Using Self-Organizing Maps*

Christiane Rennert, Sebastian Vlaic, Eugenia Marbach-Breitrück, Carlo Thiel, Susanne Sales, Andrej Shevchenko, Rolf Gebhardt and Madlen Matz-Soja

*195 Causality Analysis and Cell Network Modeling of Spatial Calcium Signaling Patterns in Liver Lobules* Aalap Verma, Anil Noronha Antony, Babatunde A. Ogunnaike, Jan B. Hoek

and Rajanikanth Vadigepalli

Changzheng Guo, Daming Sun, Xinfeng Wang and Shengyong Mao

*235 IL-1*b *and TNF*a *Differentially Influence NF-*κ*B Activity and FasL-Induced Apoptosis in Primary Murine Hepatocytes During LPS-Induced Inflammation*

Julia Rex, Anna Lutz, Laura E. Faletti, Ute Albrecht, Maria Thomas, Johannes G. Bode, Christoph Borner, Oliver Sawodny and Irmgard Merfort

Sara Zafarnia, Anna Mrugalla, Anne Rix, Dennis Doleschel, Felix Gremse, Stephanie D. Wolf, Johannes F. Buyel, Ute Albrecht, Johannes G. Bode, Fabian Kiessling and Wiltrud Lederle

# Editorial: Systems Biology and Bioinformatics in Gastroenterology and Hepatology

### Peter L. M. Jansen1 , Kai Breuhahn2 , Andreas Teufel 3 and Steven Dooley 4 \*

*1 Emeritus Professor of Hepatology, Amsterdam University Medical Center, Amsterdam, Netherlands, 2 Institute of Pathology, University Hospital Heidelberg, Heidelberg, Germany, 3 Division of Hepatology, Division of Clinical Bioinformatics, Department of Medicine II, Medical Faculty Mannheim, Heidelberg University, Heidelberg, Germany, 4 Division of Molecular Hepatology, Department of Medicine II, Medical Faculty Mannheim, Heidelberg University, Heidelberg, Germany*

Keywords: gastroenterology, hepatology, systems biology, systems medicine, mathematical modeling

### **Editorial on the Research Topic**

**Systems Biology and Bioinformatics in Gastroenterology and Hepatology**

### HOW SYSTEMS MEDICINE IMPROVES OUR UNDERSTANDING OF COMPLEX GASTROENTEROLOGICAL DISEASES

Traditional medical research gained tremendous improvements in detection and treatment of acute and chronic metabolic and inflammatory diseases as well as cancer. Especially the field of hepatology and gastroenterology has significantly benefitted from these advances. Indeed, the discovery of basic molecular and cellular disease mechanisms in the last 60 years led to the development of reliable diagnostic tests and effective therapies. For instance, the discovery of the hepatitis B and C viruses (HBV, HCV) led to powerful diagnostic tools, antiviral drugs, and an HBV vaccine (Szmuness et al., 1980, 1981; André, 1990; Lau and Wright, 1993). Indeed, mass vaccination in Taiwan led to a significant reduction of HBV prevalence and hepatocellular carcinoma incidence (Chang et al., 1997, 2016). Furthermore, the development of direct-acting antiviral drugs allows the eradication of HCV (Das and Pandya, 2018). These achievements occurred in a relatively short period of time. For example, HCV was discovered in 1989 and the first effective antiviral therapy for one genotype was developed only 20 years later (Boettler et al., 2019; Viganò et al., 2019; Zajac et al., 2019). Equally, for many monogenetic liver diseases reliable diagnostic tests exist (Lammert, 2016; Weber and Lammert, 2017). Although effective pharmacotherapy for some of these diseases exists (Wilson's disease), comparable treatments are not available for others (e.g., progressive familial intrahepatic cholestasis, PFIC). In this context, translation of gene knock-out mouse models emerged as powerful tool to understand the underlying disease processes (Liu, 2013) and may eventually lead to successful gene therapy.

In contrast to mono-factorial diseases caused by viruses or individual gene alterations, the situation is quite different in more complex multi-factorial diseases. For example non-alcoholic fatty liver disease (NAFLD) was first described by the pathologist Jurgen Ludwig in 1980 (Ludwig et al., 1980) but no effective therapy is currently available, 40 years later (Altinbas et al., 2015; Gottlieb et al., 2019). Despite all recent success in the development of treatments of HBV and HCV, targeting multi-factorial metabolic liver disease like NAFLD or the more serious non-alcoholic steatohepatitis (NASH) is still difficult. It remains unclear if targeting liver disease alone will pave the way to success or if broader approaches will ultimately lead to a decrease in patient's mortality and morbidity. Importantly, a plethora of reasons significantly complicates the systematic analysis of complex diseases:

### Edited and reviewed by:

*Stephen J. Pandol, Cedars-Sinai Medical Center, United States*

> \*Correspondence: *Steven Dooley*

*steven.dooley@ medma.uni-heidelberg.de*

### Specialty section:

*This article was submitted to Gastrointestinal Sciences, a section of the journal Frontiers in Physiology*

Received: *04 October 2019* Accepted: *07 November 2019* Published: *22 November 2019*

### Citation:

*Jansen PLM, Breuhahn K, Teufel A and Dooley S (2019) Editorial: Systems Biology and Bioinformatics in Gastroenterology and Hepatology. Front. Physiol. 10:1438. doi: 10.3389/fphys.2019.01438*

The multi-scale and multi-stage complexity of NAFLD and NASH, and the necessity to perform one or more liver biopsies to stage and monitor the disease, form a considerable obstacle in the development and application of targeted precision medicine. For instance, the anti-oxidant vitamin E is a recommended therapy in non-diabetic patients with biopsy-proven NASH (Sanyal et al., 2010; Chalasani et al., 2018). In real life only a minority of patients receives this drug (Ratziu et al., 2012). The FXR-agonist obeticholic acid also showed positive results in biopsy-proven non-cirrhotic patients with NASH (Neuschwander-Tetri et al., 2015). The mechanism of action of this drug in NAFLD is unclear and long-term effects and safety need to be assessed. Treatment of cirrhotic patients with this drug cannot be recommended. For successful wide-scale pharmacotherapy programs, non-invasive disease biomarkers are clearly needed.

NAFLD represents a multi-scale disease, in which the entire metabolic program of liver hepatocytes (incl. structures at the molecular level and subcellular organelles), non-parenchymal cells (incl. cholangiocytes, endothelial cells, Kupffer cells, and hepatic stellate cells), and the blood stream (incl. the presence of immune cells and sub-cellular blood components) are critically involved in disease development and progression. In addition, there is dysfunctional temporal communication between organs such as liver, gut, brain, pancreas and fatty tissues. Indeed, NAFLD develops in 20–30 years, starting from "simple" steatosis and progresses to pronounced liver cirrhosis. These dynamic spatial and temporal changes significantly increase the level of complexity and further complicate biomarker and drug development. For instance, pharmacotherapy targeting steatosis may be more effective in early disease stages while drugs that act on inflammation and fibrosis are more suitable at later stages. Once advanced cirrhosis with profound architectural changes of the liver and portal hypertension is established, effective pharmacotherapy becomes even more difficult.

Are computational approaches and systems medicine the solution for complex diseases? Dynamic processes can be described mathematically with a set of differential equations. With a number of these equations, scientists can generate computational models, which can describe the time-resolved behavior of molecular reactions and cellular processes (Schliess et al., 2014; Meyer et al., 2017; Berndt et al., 2018; Hoehme et al., 2018; Lucarelli et al., 2018; Poloznikov et al., 2018; Kockerling et al., 2019). In this process, experimentalists provide quantitative and semi-quantitative data derived from in vitro and in vivo models to feed these mathematical constructs. Once a reliable and robust computational model is established, the model can be used for in silico research, an approach that has the potential to save laboratory animals and to protect people and patients before a drug is used or tested in a clinical setup.

This model-based gain of knowledge leads to a process of iteration and re-iteration between theoretical and experimental scientists until the mathematical model is a reliable proxy of the in vivo situation. Sometimes predictions cannot perfectly reflect the processes observed in living cells or organisms; however, these complications can also lead to new scientific knowledge. One example within one of the biggest systems biology consortia (LiSyM, see below) was the finding that ammonia detoxification was less affected by damaging the centrizonal glutamine synthase-containing hepatocytes than predicted. These unexpected findings led to the discovery of a novel ammonia detoxification pathway (Schliess et al., 2014).

Since 15 years the German Ministry of Education and Research (BMBF) fosters systems biology and systems medicine by supporting the collaboration of multidisciplinary research groups, including biologists, clinical researchers, and mathematicians, working on liver physiology and liver diseases, including NAFLD. The research network HepatoSys was launched in 2004 to study the processes in liver cells with a systems biology approach. It was Europe's first funding measure in this field. The follow-up project the "Virtual Liver Network (VLN)" took the systems biology liver research to the next biological level. Drawing on the findings at cellular level, the network examined the processes for the whole organ. The initiative was the first systems biology network that focused on an entire organ. The current funding activity "Research Network Systems Medicine of the Liver—LiSyM" builds again on the results produced by HepatoSys and VLN. LiSyM aims to transfer the computational models into clinical application for use as diagnostic tools to assist doctors in choosing the most appropriate therapy. LiSyM and the proceeding initiatives have been successful over the years in developing computational models that help theoreticians and experimentalists to discover new aspects of signaling pathways and mechanisms to test the therapeutic potential of new molecules or biological agents in vivo as in silico (www.lisym.org).

However, opportunities to discuss the diverse facets of systems biology from data generation, utilization of mathematical models, and data integration among experts in the field remain rare. In this regard, we were happy that Frontiers in Physiology provided a platform for such urgently needed discussions and visibility beyond the German networks. The collection of articles in this special issue of Frontiers in Physiology provides examples of the current status of research in gastrointestinal diseases, including NAFLD, alcoholic hepatitis, viral hepatitis, liver fibrosis, and liver cancer, applying systems biology at the level of cells, zones, tissues, networks, and with regard to systemic consequences.

The issue comprises 20 articles from more than 170 authors. From June 2017, where the first article was accepted, until July 30 2019, the manuscripts of the issue have nearly 46,000 views. In more detail, 1 review and 19 original articles are included with nearly half of it (9 in total) from participants of the BMBF-funded networks VLN and LiSyM.

Twelve contributions are related to liver, 3 to colon and gastrointestinal tract and 1 to pancreas, again highlighting

### REFERENCES

the predominant role of the German network in this new scientific field. Experimental data of the contributions include gene expression arrays (4), metabolomics data (6), proteomics data (3), imaging (2), and signal transduction pathways (6). The modeling type of the manuscripts include high throughput data and bioinformatics in 10 cases and mathematical modeling in 9 contributions. We believe that this initiative successfully provided a platform for researchers and clinicians who are interested in systems medicine with focus on gastroenterology and hepatology.

### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

the spatiotemporal phenotype of early hepatocellular carcinoma. Bull. Math. Biol. 80, 1134–1171. doi: 10.1007/s11538-017-0375-1

**Conflict of Interest:** 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 © 2019 Jansen, Breuhahn, Teufel and Dooley. 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) and the copyright owner(s) 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.

# The Potential for Circulating Tumor Cells in Pancreatic Cancer Management

### Michael Pimienta1, 2, Mouad Edderkaoui 2 , Ruoxiang Wang2 and Stephen Pandol 2 \*

1 University of California, San Diego School of Medicine, La Jolla, CA, United States, 2 Cedars-Sinai Medical Center, Basic and Translational Pancreas Research, Los Angeles, CA, United States

Pancreatic cancer is one the most lethal malignancies. Only a small proportion of patients with this disease benefit from surgery. Chemotherapy provides only a transient benefit. Though much effort has gone into finding new ways for early diagnosis and treatment, average patient survival has only been improved in the order of months. Circulating tumor cells (CTCs) are shed from primary tumors, including pre-malignant phases. These cells possess information about the genomic characteristics of their tumor source in situ, and their detection and characterization holds potential in early cancer diagnosis, prognosis, and treatment. Liquid Biopsies present an alternative to tumor biopsy that are hard to sample. Below we summarize current methods of CTC detection, the current literature on CTCs in pancreatic cancer, and future perspectives.

### Edited by:

Steven Dooley, Heidelberg University, Germany

### Reviewed by:

Matthias Löhr, Karolinska Institutet, Sweden Ralf Jesenofsky, Heidelberg University, Germany

> \*Correspondence: Stephen Pandol<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Gastrointestinal Sciences, a section of the journal Frontiers in Physiology

Received: 30 January 2017 Accepted: 22 May 2017 Published: 02 June 2017

### Citation:

Pimienta M, Edderkaoui M, Wang R and Pandol S (2017) The Potential for Circulating Tumor Cells in Pancreatic Cancer Management. Front. Physiol. 8:381. doi: 10.3389/fphys.2017.00381 Keywords: circulating tumor cells, pancreatic cancer, pancreatic adenocarcinoma, metastasis, epithelial to mesenchymal transition, mesenchymal to epithelial transitions

# INTRODUCTION

Virtually all cancers have the potential to metastasize, and metastatic disease comes about from a series of events involving the interplay between primary tumor cells and their microenvironment. The end result is the dissemination and growth of tumor cells in new tissue environments. First described in the literature in the nineteenth century by British surgeon James Paget, metastatic disease largely remains an unsolved worldwide public health concern today—metastasis accounts for more than 90% of cancer-related deaths (Ashworth, 1869; Spano et al., 2012).

We now know that metastasis is an "extremely complex" multistep process. Tumor cells must advance through an invasion-metastasis cascade. In order to produce clinically detectable lesions, primary tumor cells need to progressively intravasate through the basal membrane into the systemic or lymphatic circulation, survive in the circulatory environment, adhere to vessel walls, extravasate into a foreign tissue site, and adapt, survive, and proliferate in their new microenvironment (Fidler, 2003).

The potential for mobilization and invasion are critical to the process of intravasation. Cells should be able to degrade the extracellular matrix (ECM) and secrete proteolytic enzymes to facilitate migration and intravasation into the circulatory system. These epithelium-derived cancer cells are thought to undergo a morphological change of epithelial-to-mesenchymal transition (EMT; Leber and Efferth, 2009; Dhamija and Diederichs, 2016). At metastatic sites, neoplastic cells should also infiltrate the endothelium in order to colonize new tissues and be able to induce neo-angiogenesis to ensure sufficient blood supply to the newly formed tumor in order to maintain metabolic needs (Leber and Efferth, 2009; **Figure 1**). The preferred growth and survival of cancer cells at certain metastatic sites is less understood, although a reverse behavioral change, mesenchymal-to-epithelium transition (MET), is hypothesized to take place for this to occur (Dhamija and Diederichs, 2016).

Metastasis itself is a highly inefficient process, as each step in the metastasis cascade may play a limiting role in disease progression, that is—if one fails then all fails. Only few cancer cells are able to go on to form malignant secondary tumors. Animal studies on the kinetics of each step have shown that postextravasation steps create the largest barrier to metastasis. An early report, for example, found that 80% of injected cancer cells survive in circulation and extravasate into distant tissues, but only 1 in 40 cells formed micrometastasis and 1 in 100 micrometastasis actually progressed to macroscopic tumors (Luzzi et al., 1998). Further, studies found similar results, showing the high efficiency of extravasation and survival in the circulation to be independent of cell's malignant potential (Koop et al., 1995, 1996; Cameron et al., 2000). These findings suggest that the growth of circulating tumor cells (CTCs) in a new microenvironment is a key step in metastatic tumor formation.

Materials shed from tumors are being investigated for their potential use in diagnosis, prognosis, and management of cancer. CTCs, cell-free circulating tumor DNA (ctDNA), and tumor cell produced exosomes (oncosomes) all hold promise in the current cancer research for future clinical use in diagnosis and management. Oncosomes, nanovesicles actively shed from most types of cancers in large numbers, like CTCs, are a dynamic source of information regarding the genomic characteristics from the parent tumor of which they are release (Deneve et al., 2013; Pantel and Alix-Panabieres, 2013; Speicher and Pantel, 2014). Though each has their own technical challenges, these materials are now being isolated and detected in peripheral blood of patients with many types metastatic cancers and are leading the way toward use of "liquid biopsies."

# CTCS IN MODERN CANCER RESEARCH

While, it is estimated that only 0.01% of CTCs have metastatic potential, the clinical importance of CTCs in "modern cancer research" over the past two decades has become increasingly apparent (Zhe et al., 2011). Although the presence of CTCs in blood could not exclusively indicate clinically significant macro-metastases, partly due to metastatic inefficiency, it surely indicates the presence of malignant tumors in situ. More and more data suggests that solid tumor shedding occurs early in disease. A recent study detected disseminated tumor cells (DTCs) in bone marrow of a mouse breast cancer model during pre-malignant stages, reinforcing the idea of early spread of tumor cells to distant organs. Fast, specific, and sensitive detection of CTCs may have potential to enhance diagnosis, treatment, and cancer monitoring. CTCs may additionally be exploited for genotypic and phenotypic abnormalities representative of the tumor in situ.

CTC detection in peripheral blood has been reported in a number of cancer types, such as lung (Zhang et al., 2014), metastatic breast (Riethdorf et al., 2007), prostate (Hu et al., 2013), colorectal cancer (Kuboki et al., 2013), and gastrointestinal and biliary cancers (Al Ustwani et al., 2012; Tsujiura et al., 2014). Many of these cancers are diagnosed at late stages, resulting in high rates of mortality. Advancements in CTC collection, enrichment, and characterization have led to increased interest in the clinical use of CTCs. Studies in various organ systems have consistently shown that CTCs rarely exist in the blood of healthy subjects, consolidating their utility in the clinical laboratory (Sastre et al., 2008; Hou et al., 2013; Tsai et al., 2016), while supporting the premise of CTCs with the potential of a powerful biomarker.

Various available CTC detection technologies have expanded their use from simply diagnostic markers to tools to evaluate overall survival, risk of metastasis, and response to therapy. Recently, using the Veridex Cell Search System, the only FDA approved method for CTC enumeration in whole blood, Weissenstein et al., found a strong correlation between median overall survival in metastatic breast cancer patients with <5 CTC/7.5 ml vs. those with ≥5 CTC (p = 0.00006; Weissenstein et al., 2012). In a review of multicenter study of 1,358 individuals, Miller et al. found a highly significant median overall survival in favorable CTC counts vs. patients with unfavorable CTC counts with metastatic breast cancer, metastatic colorectal cancer, and prostate cancer (p < 0.0001). Additionally, patients enrolled in therapies that decreased CTC counts displayed improvements in overall survival, pointing to the utility that CTC analysis in the response to anti-cancer treatments holds (Miller et al., 2010). Studies have also taken advantage of the current available methods to correlate tumor dissemination and stage to CTC count. Hiraiwa et al. showed that CTC counts were higher in metastatic patients than in non-metastatic esophageal, gastric, and colorectal cancers and were significantly correlated to advanced tumor stages. High CTC count, defined as 2 or more CTCs per 7.5 ml in this study, was linked to pleural and peritoneal dissemination (Hiraiwa et al., 2008). Further, validation in clinical settings will establish CTC detection as a marker for sensitive and non-invasive cancer diagnosis, treatment evaluation and prognosis.

# CTC ENRICHMENT AND DETECTION

CTCs are rare cells, detected in numbers ranging from 1 to 10 s per ml whole blood, among billions of red blood cells and millions of leukocytes. CTC detection and isolation remain being technologically challenging (Joosse and Pantel, 2013). In the infancy of this rapidly growing field, current CTC detection and analysis rely mainly on various methods of enrichment.

CTC enrichment approaches exploit the unique biological and/or physical properties of this specified tumor cell type among vast numbers of peripheral blood cells, in order to increase CTC recovery by many orders of magnitude. Immunoaffinity platforms select for CTCs based on the expression of specific surface antigens, through either positive or negative selection. death.

Since most carcinomas express epithelial markers, epithelial cell adhesion molecule (EpCAM) is most commonly used in antibody-based positive selection with commercial technologies. Negative selection depletes other mononucleated cells through anti-CD45 antibody use. Established collection methods have employed these proteins to attract and adhere CTCs to columns, microposts, or magnetic apparatus (Alix-Panabières and Pantel, 2014).

Two early platforms, magnetic-activate cell sorting system (MACS), and Dynabeads use magnetic fields for attract CTCs to anti-EpCAM antibody coated magnetic microbeads (Nagrath et al., 2016). Similarly, CellSearch (Janssen Diagnostics) uses anti-EpCAM conjugated ferrofluid nano-particles to immunomagnetically capture CTCs, which may then be differentiated from contaminating leukocytes based on positive cytokeratin or EpCAM staining and negative CD45 staining (Hayes et al., 2006). Despite much progress in platform development, CellSearch remains the only FDA approved method of whole blood CTC enrichment and enumeration. Another platform, MagSweeper (Stanford University), uses magnetic rods, stirred through diluted blood samples, to attract CTCs pre-labeled with EpCAM-magnetic beads (Talasaz et al., 2009). This platform was one of the first to enrich CTCs with a notably higher purity than its predecessors. Importantly, it has the ability to isolate live CTCs without perturbing gene expression throughout the enrichment process, providing viable CTCs for analysis (Krebs et al., 2014).

Microfluidic devices, which allow separation of CTCs, from small fluid volumes under laminar flow, are promising technologies. Nearly 10 years ago, Nagrath et al. was able to selectively and efficiently isolate CTCs from whole blood of 115/116 (99%) cancer patients using anti-EpCAM-coated posts with this "CTC-Chip" platform, eliminating the need for pre-labeling or sample processing (Nagrath et al., 2007). New microfluidic approaches have appeared since then, including methods taking advantage of physical properties as well, making the platform applicable to isolation of CTCs that lack or have down regulated EpCAM expression (Stott et al., 2010; Ozkumur et al., 2013).

Positive and negative CTC enrichment is also tested based on physical properties alone. Tumor cells and CTCs are generally thought to be larger (>8µm) than hematologic cells (Vona et al., 2000, 2004; Hosokawa et al., 2013). Size-based filtration methods, such as isolation by size of epithelial tumor cells (ISET) through membrane filters with size exclusive pores, have been previously used to isolate individual CTCs and achieved higher sensitivity than CellSearch (Hofman et al., 2011; Hou et al., 2011).

As CTCs extravasate and intravastate the circulation, they undergo massive deformations in their structure due to mechanical forces which they endure. Cancer cells are known to be more deformable than normal cells, a quality that is correlated to their metastatic potential and is exploited by some enrichment platforms (Byun et al., 2013; Park et al., 2016). Additionally, platforms utilizing tumor cell property differences in electrical charge and density have been reported (Müller et al., 2005; Fabbri et al., 2013; Yoo et al., 2016).

While many of the platforms available have been able to detect CTCs in blood samples, these small peripheral blood collections, of several milliliters, may not be representative of the of the entire patient blood volume. Recently, an antibody coated medical wire (CellCollector; Gilupi GmbH) capable of detecting EpCAM and cytokeratin positive CTCs in vivo was introduced. In a recent study, a 30 min incubation period, in which the wire was exposed to circulating blood in the arm vein of lung cancer patients, showed over a 2-fold increase in CTC detection in comparison to CellSearch (Gorges et al., 2016b). Further, efforts to enhance the biocompatibility of these wire coatings, have been employed to maximize functionality for downstream processes such as sequencing analyses of captured cells (Scherag et al., 2017).

Once enriched, CTCs are detected/confirmed through various techniques. Immunocytological and molecular approaches are the most commonly employed. Immunocytochemistry (ICC) may differentiate CTCs from contaminating cells through biomarker detection. Such biomarkers can be specific for nuclear content, epithelial proteins (i.e., cytokeratins), and hematopoietic markers (i.e., CD45). A common immunocytological CTC definition, currently used by CellSearch and other platforms, is a Nucleus+/CK+/CD45− cell. However, it should be noted that CTC designation depending primarily on epithelial marker expression may lead to false negatives by failing to detect CTCs that have undergone EMT (Lustberg et al., 2012). As cell phenotypes can vary in different malignancies, the heterogeneity of CTCs pose barriers to efficient and thorough detection by liquid biopsies.

Reliance on epithelial markers, which most epithelial carcinomas express, for enrichment and identification, fails to capture subpopulations of CTCs, such as mesenchymal cells, that may harbor clinically important information. Currently, cancer studies in breast and prostate have already demonstrated that mesenchymal marker expression by CTCs is associated with poorer survival (Aktas et al., 2009; Yokobori et al., 2013). Recently, negative depletion strategies that enrich CTCs in phenotype-independent ways have been introduced in an effort to solve this problem and enhance detection. Immunostaining of CD45 and cell sorting with flow cytometry was used to enrich the breast cancer CTC population (Lara et al., 2004). Multimarker Immunomagnetic Negative Depletion Enrichment of CTCs (MINDEC), relies on depletion of non-CTCs as opposed to targeting specific properties of CTCs (Lapin et al., 2016). This technique is based on a multi-marker antibody cocktail (CD45, CD16, CD19, CD163, and CD235a/GYPA) to target various contaminating blood classes. This technique has shown high enrichment efficiency of both epithelial and mesenchymal CTCs, with better hematopoietic depletion than CD45 alone. Additional, novel cell surface marker-independent techniques have been shown to effectively detect CTCs in epithelial and non-epithelial malignancies in the absence of cell surface tumor markers. For example, a novel method introduced by Zhang et al., selectively labeled CTCs through GFP expression in human samples and cancer cell lines transfected with tumor selective replicating HSV-1 with a high detection efficiency (Zhang et al., 2016).

**Table 1** shows various CTC isolation methods used in the last few years.

Newer methods have employed combinations of epithelial, mesenchymal, tumor-specific, and tissue-specific marker expression (Pantel and Alix-Panabieres, 2013).

Additionally, nucleic acid-based technologies have provided an alternate avenue (Yu et al., 2011), as improvements in nonfixating enrichment procedures have allowed for the use of RT-PCR and qRT-PCR to amplify single or multiple gene transcripts for CTC detection.

Most recently, emerging single-cell sequencing techniques have shifted the field toward individual CTC analysis of genetic alterations associated with tumor mechanisms, clinical outcomes, therapy response, and drug targets and resistance. The usefulness of genomic analyses, however, is limited by heterogeneity between cancer subtypes, presenting barriers toward finding universal markers. Similarly, all the current enrichment, detection, and analysis techniques available harbor their own technical challenges and limitations. Most of these are outside the scope of this discussion.

## CTCS IN PANCREATIC CANCER DIAGNOSIS

Pancreatic cancer is one of the deadliest malignancies. Pancreatic ductal carcinoma (PDAC) makes up the majority of pancreatic cancers. While advancements in the treatment of other cancer types may have led to significant improvement in patient survival, advancements in pancreatic cancer research have not been met with the same success. PDAC incidence has remained stable over the last 30 years and the lack of fruitful therapies and new/useful diagnostic methods have yet to be changed in pancreatic cancer (Ryan et al., 2014). With current therapies improving survival outcomes by only a few months, pancreatic cancer patients face a 5 year survival rate of only 7%. Gemcitabine, the first-line of PDAC therapy, only modestly improves survival in advanced pancreatic cancer, while the clinical benefit of combinational-targeted therapies (Erlotinib+Gemcitabine) has proven to have only slight benefit, increasing overall survival by less than a month (Burris et al., 1997; Moore et al., 2007). Recent work in combinational chemotherapeutics has led to a promising approach, FOLFIRINOX (Oxaliplatin, Irinotecan, Leucovorin, and 5-fluorouracil), which has almost doubled survival in metastatic pancreatic adenocarcinoma patients to 11.1 months compared to 6.8 months with single-agent gemcitabine (Conroy et al., 2011). In another promising approach, albumin bound paclitaxel (Abraxane) plus standard gemcitabine therapy increased overall survival to 8.5 months compared to standard gemcitabine therapy with an overall survival of 6.7 months (Von Hoff et al., 2013). The toxicity associated with these regimens is unfavorable and should be used in patients with good performance status. Unfortunately, 5 year survival rates remain

### TABLE 1 | CTC isolation techniques.

relatively unchanged. It is estimated that by 2030 pancreatic cancer will be the second leading cause of cancer-related deaths (Rahib et al., 2014; Dawson and Fernandez-Zapico, 2016).

The idea of a liquid biopsy, which could reveal diagnostic and prognostic information about a patient's state, has been gaining much traction in the past 10 years. In one of the earliest studies, with 12 types of metastatic carcinomas in 964 patients, CTCs were successfully detected in patients with pancreatic cancer using CellSearch, albeit in lower numbers than the other cancers (Allard et al., 2004). Below, we summarize the studies of CTCs in the diagnosis, staging, and prognosis for pancreatic cancer patients.

Pancreatic cancer is a fast progressive disease and its early diagnosis is challenging. Initial pancreatic cancer diagnosis depends largely upon symptoms, which would only appear late when tumor have fully progressed and are not specific to be recognized at early stages. Due to the pathobiology and aggressiveness of PDAC, by the time anorexia, early satiety, pain, and weight loss start present, the disease has already progressed, leaving little room for a favorable prognosis. Additionally, of the 15% of patients seeking medical care 6 months prior to diagnosis, 25% have symptoms resembling upper abdominal disease that may lead to misdiagnosis (DiMagno, 1999). Affirmative diagnosis is made by tissue biopsies obtained by surgery, image guided CT biopsy, or fine needle aspiration through endoscopic ultrasound (EUS-FNA). Despite its widespread use, EUS-FNA does have diagnostic drawbacks, specifically a sensitivity range from 75 to 94% and a specificity of 78 to 95%, with low but lethal complications such as pancreatitis and bowel perforation (Court et al., 2015; Bournet et al., 2016). For CTC detection to be adopted to pancreatic cancer diagnosis, it must be useful for early diagnosis and/or monitoring treatment responses. A key performance milestone necessary for the implementation of CTC technologies is an understanding of disease stage at which CTCs can be detected. At the same time, it should be kept in mind that CTC research in pancreatic cancer is at nascent stage, while CTC detection methods and criteria vary largely between studies. It is important to critically analyze the markers available and to characterize CTCs at different stages of cancer progression.

CTC detection has been explored in early diagnosis of various cancers and CTCs have been detected prior to tumor detection by traditional methods. A recent study, for instance, found that CTCs could be detected 1–4 years before lung cancer became detectable through CT-scan screening in the same patients (Ilie et al., 2014). For breast cancer, the American Society of Clinical Oncology (ASCO) has already approved the use of CTCs as a tumor marker, creating new directions for early breast cancer diagnosis (Harris et al., 2007).

Similar results have been obtained in pancreatic cancer. In a mouse model of PDAC, Rhim et al. found that inflicted pancreatic cells underwent EMT early during cancer development. These cells with EMT were predicted to represent early cancer cells, as the extent of EMT correlated well to invasive properties in tumor cells, facilitated their intravasation to circulating the blood and hepatic seeding, prior to the manifestation of primary tumors (Rhim et al., 2012). Additionally, blood samples from patients with pancreatic cystic lesions were detected to contain pancreas epithelial cells, at a time prior to cancer diagnosis. These findings suggested that pancreatic cell appearance in circulating blood precedes in situ tumor formation, and detection of CTCs could be an early biomarker for PDAC early diagnosis.

CTC-specific gene expression has been explored as surrogate markers for early cancer detection. Such studies mostly detect CTCs via detecting their expression of epithelial proteins. Reverse transcription-coupled polymerase chain reaction (RT-PCR) could be used to examine potential epithelial markers in CTCs derived from tumors of the epithelia. Soeth et al., evaluated cytokeratin 20 (CK20) detection, through RT-PCR detection from the marrow and venous blood in pancreatic ductal carcinoma patients. CK-20 positivity was detected in cells of 52 of 154 patients in venous blood, where higher CK20 level was correlated to UICC-tumor stage (Soeth et al., 2005). In a study of 34 pancreatic cancer patients prior to treatments, de Albuquerque et al used immunomagnetic enrichment for CTCs in peripheral blood based on mucin-1 and EpCAM expression. Subsequently, multi-marker RT-PCR analysis was used to detect tumor-associated transcripts, including KRT19, MUC1, EpCAM, CEACAM5, and BIRC5. CTCs with at least one marker in peripheral blood were detected in 47.1% patients prior to undergoing treatment (de Albuquerque et al., 2012). Detection efficacy was increased with the use of multiple markers as opposed to single marker use, indicating differential gene expression among CTCs from the same cancer patient.

Zhang et al., used an alternative strategy by enriching and identifying CTCs through a combination of CD45 and CK with a FISH-CEP8 probe in 22 pancreatic cancer patients. CTCs were detected from 15 of the patients, with CTCs ranging from 0 to 60 cells/3.75 ml of blood. In comparison, healthy controls, and patients with benign pancreatic tumors were negative for detection of CTCs, and sensitivity and specificity of CTC detection in pancreatic diagnosis were determined to be 68.18 and 94.87%, respectively, when using 2 cells/3.75 ml as cutoff. CTC-positive patients exhibited metastasis and poorer survival rates upon a 1.5 year follow-up. CTC positivity did not correlate significantly to CA19-9 levels of the in situ tumor. It is well-known that, though CA19-9 has a high sensitivity and specificity in advanced pancreatic cancers, its diagnostic usefulness is questionable to diseases at early asymptomatic stages. The combination of CA19-9 and CTC positivity in the study above increased detection rates from 68.18 to 77.3% (Zhang et al., 2015). The results from this study suggest CTCs as biomarkers for the diagnosis of early stage pancreatic cancers, in asymptomatic patients, and from patients with normal CA19-9 plasma levels. A later study by Xu et al., used a similar approach in 40 patients and dramatically high detection rates in PC patients (90%). Diagnostic rates increased to 97% when combining CTC ≥ 2 and CA19-9 > 37 µmol/L as a cutoff. Identification of chromosomal instability in CTCs, characterized as chromosome 8 triploids, showed a significantly statistic prognostic correlation. Patients with triploid CTCs < 3 displayed both higher 1 year and overall survival compared to those with ≥ 3 (Xu et al., 2017).

Doublecortin-like kinase 1 (DCLK1) may be another marker for CTC detection in early PDAC stages. Prior work suggests that Dclk1 marks stem cells, being able to differentiate cancer from normal stem cells. Interestingly, this marker protein is overexpressed in pancreatic and colorectal cancers (Nakanishi et al., 2012; Bailey et al., 2014). Qu et al., found elevated serum DCLK1 levels in stage I and II PDAC patients relative to controls and a decline in stage III and IV patients to levels similar to those seen in control patients. Diagnostic utility analysis showed both DLCK1 sensitivity and specificity to be significant in stage I and II patients. Furthermore, the investigators evaluated DCLK1 in the KPC mouse model, finding the serum DCLK1 levels to be significantly elevated as early as 5 weeks in KPC mice, compared to control mice. Over 50% of CTCs isolated from KPC mice whole blood were DLCK1+, suggesting a possible biomarker to be used in conjunction with CTCs for detection of early stage pancreatic cancer (Qu et al., 2015).

Kulemann et al., investigated the usefulness of CTC detection in pancreatic cancers from both localized and advanced stages. Peripheral blood from PDAC patients was used to capture CTCs for cytological and KRAS mutational analysis using the ScreenCell isolation method (Kulemann et al., 2015). High CTC detection efficiency as low as 2 cells/ml was calibrated by spiking experiments with healthy donor blood. CTC KRAS mutations were identified in 8 of 11 PDAC patients (73%). This is in sharp contrast to conventional biopsy-based diagnosis for the same patients, by which 2 of 11 samples (18%) were cytologically categorized as negative/non-diagnostic, while the rest exhibited abnormal morphology (18%) or were categorized as suspicious (64%). Moreover, the authors found no difference in the detection rate between early and advanced diseases, suggesting that CTCs are disseminated from primary tumors early in disease development and can be used to diagnose pancreatic cancer at initial stages where curative surgery may be available. It should be pointed out that this finding is contrary to those shown by Soeth et al. (2005), in which significant stage dependent differences were observed in CTC detection. Further studies are needed to clarify whether the difference is due to the differences in CTC detection strategies.

### CTCS IN METASTASIS AND EARLY CANCER

Early local invasion and metastasis are prominent factors in the poor prognosis of pancreatic cancer, as most patients are found with metastatic disease at diagnosis (DiMagno et al., 1999; Pandol et al., 2009). The impact of CTCs on pancreatic cancer metastasis, recurrence, and prognosis has been investigated.

A recent 9-cohort meta-analysis of separate studies using CellSearch and RT-PCR detection methods, involving 623 pancreatic cancer patients altogether, revealed associations between CTC detection and poor prognosis. Out of 623 patients, 268 (43%) were classified as CTC positive and displayed poor progression-free survival and worse overall survival than those in the non-CTC group (Han et al., 2014).

Using the CellSearch enrichment method, Kurihara et al., investigated the utility of CTCs in peripheral blood as a marker of clinical outcomes in 26 patients with pancreatic cancer. CTC positivity was found in 11 of 26 pancreatic cancer patients (42%). The authors demonstrated a significant difference in median survival times between CTC positive and negative patients, for 110.5 and 375.8 days, respectively (P < 0.001; Kurihara et al., 2008). Given that the detection method was shown to have 100% specificity, with no CTCs detected in the non-cancer groups, CellSearch detection strategy may not be sensitive enough, as no CTCs were detected from the other 58% pancreatic cancer patients. These results suggest the need of developing more sensitive methods to detect positive CTCs in all pancreatic cancer cases.

Similarly, de Albuquerque et al found a correlation between CTC positivity (47% of patients) and median progressionfree survival (PFS). Patients with at least one tumor-associated transcript found in their CTCs, enriched from peripheral blood using immunomagnetic EpCAM and mucin 1 detection, had a PFS of 66.0 vs. 138.0 days in those who did not. Intriguingly, CTC enumeration was found to have no correlation to clinicopathological features of the disease, including metastasis status and tumor stages (de Albuquerque et al., 2012).

Bidard et al. studied CTC detection rates in a subset of 79 patients with locally advanced pancreatic carcinoma enrolled in the LAP 07 trial. The primary study (LAP 07) assessed the effect of subsequent chemotherapy vs. chemo-radiotherapy continuation on overall survival in patients whose disease was controlled after 4 months of chemotherapy alone. The patient subgroup was screened for CTCs using CellSearch technology prior to chemotherapy administration and 2 months after treatment. While the investigators found that CTC positivity was not prognostic of PFS, they found it to be an independent prognostic factor associated with poor tumor differentiation and shorter overall survival (Bidard et al., 2013).

Bissolati et al., used the same CellSearch technique with systemic and portal vein blood of 20 patients undergoing pancreatic resection. No significant differences in both overall survival and disease-free survival between CTC-positive and -negative groups. The authors did, however, find a higher incidence of liver metastasis upon a 2 and 3 year follow up in the CTC-positive portal vein group (Bissolati et al., 2015).

Similarly, an early study investigated CTC positivity in 67 intraoperative patients with biliary-cancer. Molecular detection of CEA mRNA-positive CTCs from peripheral, central, and portal veins via RT-PCR was associated with a significant incidence of hematogenous metastases compared to CTCnegative patients (37.5 vs. 11.4%; Uchikura et al., 2002). Nonetheless, whether surgical resection of the pancreas itself may contribute to tumor cell shedding remains to be addressed. Pancreaticoduodenectomy (PD), involving the pancreatic head, and distal pancreatosplenectomy (DPS), involving the pancreatic body and tail, are standard surgical procedures. Both require necessary mobilization of the pancreas, and may lead to CTC dissemination via the portal vein to increase the risk of liver metastasis (Kuroki and Eguchi, 2017).

**Table 2** represents data showing association between CTC presence and pancreatic cancer stage and outcome.

Chausovsky et al., used RT-PCR to examine the usefulness of CK20 expression in CTCs in the diagnosis of metastatic lung, stomach, colon, and pancreatic cancers (Chausovsky et al., 1999), since CK20 has been shown to not be transcribed in cells of hematopoietic lineage (Burchill et al., 1995). Chausovsky concluded that CK20 is a potential biomarker for detecting metastasis, with a sensitivity of 22/28 (78.6%) in patients with metastatic pancreatic cancer. Cytokeratins are used conventionally to characterize cancer cells of epithelial origin (Cooper et al., 1985; Lane and Alexander, 1990). Combination of cytokeratin and additional gene expression may improve the efficacy of CTC detection.

Poruk et al., assessed the potential of CTCs as biomarkers in 50 patients prior to surgical resection, based on EMT-related epithelial and mesenchymal marker expression. CTCs were acquired from blood samples through the method of ISET. CTCs were further identified by immunofluorescence staining with antibodies against pan-cytokeratins and the mesenchymal cell protein, vimentin. This analysis found that 78% of patients had CTCs expressing cytokeratin and 67% co-expressed vimentin, while no CTCs were found to express vimentin only. The authors found a significant association between cytokeratin only positive CTCs and worse survival. Interestingly, co-expression of vimentin was predictive of recurrence (p = 0.01). Of the patients diagnosed with metastatic cancer at the time of surgery, all the CTCs were positive for dual staining (Poruk et al., 2016). These findings indicate the involvement of EMT mechanism in metastatic progression. EMT would render CTCs heterogeneous and multi-marker analysis would have to be employed in order to ensure a comprehensive detection of all CTCs in a patient blood sample.

A recent study enumerated CTCs independently of surface marker status using a GFP expressed tumor selective Herpes Simplex Virus replicated based on telomerase activity. Transfected cells of 290 samples of patients with different solid tumors were examined and CTCs were detected in patients with epithelial and non-epithelial tumors from as little as 4 ml of blood. PC patients had a positive CTC detection rate of 88.2% across various stages and had the highest average number of cells identified per samples (43.1). Additionally, CTC detection rates increased to 100% in PC patients with regional lymph node metastasis but no distant metastasis (N+M0), further supporting the use of CTCs as a biomarker in disease progression (Zhang et al., 2016). Other recent phenotypic-independent enrichment platforms have shown some success in CTC enumeration regardless of epithelial or mesenchymal surface proteins. Negative selection of hematopoietic cells in blood samples of

PC patients using MINDEC showed a CTC detection rate of 71%. Further, characterization of the enriched cells showed the presence of both epithelial and mesenchymal CTC populations. While the high rate of positivity in this proof-of-principle study, in comparison to previous phenotype specific platforms, could be due to the authors use of patients with metastatic disease only, its ability to detect both epithelial and mesenchymal cells, in addition to CTC clusters, marks a progressive trend toward comprehensive detection of both epithelial and mesenchymal CTCs with one technique (Lapin et al., 2016). More importantly, both CTC surface marker-independent enrichment techniques allow for the viability of collected cells to be subsequently used for downstream genetic analysis without compromise from high background leukocyte levels.

Recent works characterizing EMT found CTCs positive for both epithelial and mesenchymal markers in peripheral blood of breast cancer patients (Yu et al., 2013). Studies in mouse models have provided insight into the composition of CTCs in pancreatic cancer. Single-cell RNA sequencing revealed the expression of both epithelial and mesenchymal markers in KPC LSL-KrasG12D, Trp53flox/flox or +, Pdx1-Cre (KPC) mouse pancreatic tumors. Moreover, the authors observed substantial loss of the classical E-cadherin expression, suggesting that some CTCs of epithelial lineage could indeed adopt a partial mesenchymal stromal phenotype through EMT, while retaining other epithelial features such as cytokeratin expression (Ting et al., 2014).

Different from epithelial cells, most mesenchymal stromal cells harbor certain stem cell properties, being able to be induced to differentiate into more mature cells (Zhau et al., 2011). The EMT phenotype is usually associated with expression of cancer stemness markers (Kong et al., 2011). Compared to other cancers, however, very little is known about stemness in pancreatic cancer CTCs.

In breast cancer, expression of cancer stem cell markers in CTCs is a sign of increased the metastatic ability (Papadaki et al., 2014). The expression of cancer stemness marker ALDH1 on CTCs was found to correlate to the stage of the disease and to the expression of EMT markers vimentin and fibronectin in prostate cancer patients (Raimondi et al., 2011). A study by Barrière et al. (2012) aimed at the detection of CTCs endowed with mesenchymal and/or stem cell characteristics, at the time of initial diagnosis with breast cancer, found that EMT and cancer stemness occur in the primary tumors and are associated with an enhanced ability for tumor cells to intravasate in the early phase of cancer development.

Multiplex transcriptome profiling of single CTCs revealed presence of sub-populations of CTCs expressing multiple procancer transcripts including cancer stem cell markers such as CD44 and CD24 (Gorges et al., 2016a). So far multiple markers have been used to detect CTC stem cell properties in CTCs, including CD44, CD133, CXCR4, ABCG2, and ALDH1. Other markers used uniquely for pancreatic cancer CTCs include CD24 and c-Met (Yang et al., 2015). Whether CTCs with mesenchymal or stem cell characteristics may be used as a marker for aggressiveness of the disease remains to be evaluated in future studies.

### CHALLENGES

There has been much progresses over the last decade to overcome the initial barriers of CTC research in the laboratory and clinic. Significant technological development has been made for CTC detection, enrichment, and molecular characterization. On the other hand, CTC research in gastrointestinal cancers had a late start relative to other human malignancies (Allard et al., 2004). Due to gastrointestinal biology, pancreatic cancer detection via CTCs has its own set of challenges. It is widely proposed that the liver sequesters CTCs as they pass into the systemic circulation via the portal vein (Jiao et al., 2009). The predominant dissemination of pancreatic cancers to the liver have supported this notion. A recent report detected CTCs in 100% (14/14) of portal vein samples of patients with pancreaticobiliary cancers as opposed to under 25% (3/14) when peripheral samples were used for detection (Catenacci et al., 2015), a result consistent with the notion that higher CTC numbers are detected in portal vein as opposed to peripheral blood (Waxman et al., 2014). Furthermore, there was a significant increase in CTC detection in portal blood vs. peripheral blood, with a mean of 125.64 CTCs/7.5 ml as opposed to 0.8 CTCs/7.5 ml (p = 0.01; Catenacci et al., 2015). This study further emphasizes the importance that the collection site plays in CTC detection.

While preliminary studies in various cancers are demonstrating the potential of CTCs in early cancer detection, the continuous data coming out using different platforms (exploiting size, density, charge, surface antigens, etc.) make it challenging to reach a consensus for clinical application. The diversity of methods for CTC enumeration and characterization can confuse the research and the clinical communities. There is a lack of large studies comparing enrichment and detection between CTC detection platforms. A pilot comparative study of 54 pancreatic cancer patients investigated differences between CellSearch and a marker-independent ISET CTC isolation (Khoja et al., 2012). The authors detected significantly more CTCs using ISET in comparison to CellSearch (93 vs. 40%). Similar studies in other cancers have pointed out the discrepancies between platforms (Farace et al., 2011; Hofman et al., 2011).

Discrepancies can partially be attributed to different types of carcinomas and their expression of surface markers. Furthermore, the differences in specificity and sensitivity may lead investigators to adapt different platforms for their specific study. This raises many concerns in regards to the identities of CTCs. As EMT facilitates CTCs with increased capacity for detachment and invasion, the loss of epithelial lineage marker expression makes the identification particularly difficult. We must further explore the limitations that certain platforms create in capture efficiency. For example, EpCAM expression, which CellSearch exploits, is heterogeneous and cleavage has been reported (Maetzel et al., 2009). Limitations of other epithelial cell markers have also been reported, such as the down regulation of CK20 in tumors leading to false-negatives (Vlems et al., 2002; Krebs et al., 2014). It would be ideal for platforms to detect both mesenchymal and epithelial characteristics of CTCs, and the platforms must be carefully validated.

Additionally, not only do detection rates vary by platforms, but also between cancers. There is currently no consensus on the cutoff value for CTC positivity, even within a single platform. Stringent parameters should be set for CTC use not only in detection, but also as a prognostic marker of clinical outcomes in pancreatic cancer. Although it may be difficult due to the differences in enrichment and detection between platforms, standardization across single or multiple platforms is paramount for future incorporation into the field.

Due to the complex nature of the metastatic process, disseminated cells may be clinically silent for long durations. In breast cancer, cytological assessment suggests that CTCs actively undergoing mitosis are most common in late-stage disease and have prognostic value (Adams et al., 2016). One study found that aberrant VCAM1 expression, a common complication of breast cancer, was crucial for the transition from dormancy to overt metastasis (Lu et al., 2011). We must continue to explore ways to stratify CTCs in ways that will allow us to distinguish indolent micrometastasis from aggressive CTCs prior to clinically significant metastasis in pancreatic cancer patients.

# CONCLUSION AND FUTURE DIRECTIONS

Literature evaluating the diagnostic and prognostic role of CTCs in cancer is continuously being reported. Many studies in different malignancies have shown clear associations of CTCs with clinical cancer progression. Much of the current research is now shifting to CTC characterization in order to select appropriate therapies for individuals based on the gene signatures of the CTCs and to measure response to therapies. For example, CTC count now outperforms traditional response evaluation methods in patients with metastatic castrationresistant prostate cancer (Onstenk et al., 2016). With reports estimating the half-life of CTCs to be on the order of hours, their detection can provide a current representation of the malignancy (Meng et al., 2004; Stott et al., 2010).

Future investigations should thoroughly explore CTC response to pancreatic cancer treatments. Furthermore, ex vivo CTC culture and expansion experiments can improve our understanding of the mechanisms of dissemination and escape from dormancy. Single-cell sequencing with next-generation sequencing platforms is paving the way toward understanding the genetic makeup of CTCs and the clinical significance of their genomic alterations (Alix-Panabières and Pantel, 2014). A recent study in a pancreatic cancer mouse model used single-molecule RNA sequencing of CTCs to identify Wnt2 as an up-regulated gene in pancreatic cancer CTCs, which is implicated in celldeath suppression and cancer dissemination. In addition, the authors observed the same Wnt2 signaling aberrations in CTCs of 5/11 patients with metastatic pancreatic cancer (Yu et al., 2013). Such studies have the potential to improve our current clinical management, especially ones exploring new drug targets involved in cancer spread. CTC use as a biomarker is currently being investigated in over 360 open clinical trials registered on ClinicalTrials.gov (Alix-Panabières and Pantel, 2014).

Considering that methods have been developed that have the possibility of being used in the diagnosis, stratification of patients and monitoring of therapy, next efforts require a focus on validation of leading methods for aiding clinical care. For the methods chosen, the validation of the method for certification in clinical use followed by well-designed studies to show utility of the method in the clinical setting are necessary for approval of test for clinical application.

There are challenges in the pancreatic cancer field for development of a test that has utility in early diagnosis or choice of chemotherapy. For example, in the area of early diagnosis, a population at increased risk is needed to show performance of the method in detecting pancreatic cancer at an earlier stage than that achieved with current approaches. Currently, the Consortium on Chronic Pancreatitis, Diabetes, and Pancreatic Cancer (CPDPC) in the United States is developing a protocol for this purpose choosing patients with diabetes after age 50 as the best high risk group to use to identify early diagnostic biomarkers (http://cpdpc.mdanderson.org). Because about 1% of these patients will be found to have pancreatic cancer over 3 years after the diagnosis of diabetes, the CPDPC has determined that the study will require enrollment of 10,000 subjects. At present, the inclusion of measurement of CTCs is not being considered because of the technical difficulties involve with CTC measurements. On the other hand, once a proteomic, ctDNA and/or RNA technique is developed to identify patients with early pancreatic cancer, measurements of CTCs can be applied to this group for further characterization including choice of therapy.

Similarly, it is difficult to apply CTC technology to the choice of chemotherapy as the current therapies do not have a significant effect on long term survival. One the other hand, surgery does have significant effects on long term survival in a substantial percentage of patients. Thus, it seems that currently the best situation to develop a validated test for CTC measurements uses patients who are candidates for curative surgery. Hypotheses to be tested should focus on the role of CTC measurements in predicting the outcome of curative surgery

### REFERENCES

and early demonstration of disease recurrence. Certainly, studies that show performance of CTC measurements in determining and monitoring outcome in surgical patients will have important impacts in disease management and are much more feasible than studies designed for early diagnosis.

Another area of significant importance in the field is the determination of the biology of CTCs. As these cells represent the metastatic process which is the key determinant of poor outcome in pancreatic cancer patients, a better understanding of the biology of these cells will be central to advancing our treatments. Are there unique mechanisms in pancreatic cancer that account for its high rate of metastasis? Are there properties of pancreatic cancer CTCs that account for its resistance to therapy? Exploring these questions will require advancing the methods of isolation and propagating these cells so that the biologic experiments including observing their behavior in cell culture and animal models can be performed.

### AUTHOR CONTRIBUTIONS

SP provided oversight and direction for the development of the manuscript. MP did literature searches and drafting of a primary manuscript. ME and RW provided additions and editing for the manuscript.

### FUNDING

Funding was provided in part by NIH grants P01 CA163200 and U01 DK108314.

the ancillary CirCe 07 study to the LAP 07 trial. Ann. Oncol. 24, 2057–2061. doi: 10.1093/annonc/mdt176

colorectal cancer: expression in controls and downregulation in tumour tissue. Mol. Pathol. 55, 156–163. doi: 10.1136/mp.55.3.156

and mesenchymal composition. Science 339, 580–584. doi: 10.1126/science. 1228522

**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.

The reviewer RJ and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.

Copyright © 2017 Pimienta, Edderkaoui, Wang and Pandol. 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.

# Multi-Level Integration of Environmentally Perturbed Internal Phenotypes Reveals Key Points of Connectivity between Them

Nirupama Benis 1 \* † , Soumya K. Kar 1 †, Vitor A. P. Martins dos Santos 2, 3, Mari A. Smits 4, 5 , Dirkjan Schokker 4 † and Maria Suarez-Diez 2 †

1 Host Microbe Interactomics, Wageningen University & Research, Wageningen, Netherlands, 2 Systems and Synthetic Biology, Department of Agrotechnology and Food Sciences, Wageningen University & Research, Wageningen, Netherlands, 3 Lifeglimmer GmbH, Berlin, Germany, 4 Wageningen Livestock Research, Wageningen University & Research, Wageningen, Netherlands, 5 Wageningen Bioveterinary Research, Wageningen University & Research, Wageningen, Netherlands

### Edited by:

Andreas Teufel, Johannes Gutenberg-Universität Mainz, Germany

### Reviewed by:

Supriyo Bhattacharya, City of Hope Medical Center, United States Hiroshi Ishiguro, Nagoya University, Japan

> \*Correspondence: Nirupama Benis<EMAIL_ADDRESS>
† These authors have contributed equally to this work.

### Specialty section:

This article was submitted to Gastrointestinal Sciences, a section of the journal Frontiers in Physiology

Received: 01 February 2017 Accepted: 23 May 2017 Published: 12 June 2017

### Citation:

Benis N, Kar SK, Martins dos Santos VAP, Smits MA, Schokker D and Suarez-Diez M (2017) Multi-Level Integration of Environmentally Perturbed Internal Phenotypes Reveals Key Points of Connectivity between Them. Front. Physiol. 8:388. doi: 10.3389/fphys.2017.00388 The genotype and external phenotype of organisms are linked by so-called internal phenotypes which are influenced by environmental conditions. In this study, we used five existing -omics datasets representing five different layers of internal phenotypes, which were simultaneously measured in dietarily perturbed mice. We performed 10 pair-wise correlation analyses verified with a null model built from randomized data. Subsequently, the inferred networks were merged and literature mined for co-occurrences of identified linked nodes. Densely connected internal phenotypes emerged. Forty-five nodes have links with all other data-types and we denote them "connectivity hubs." In literature, we found proof of 6% of the 577 connections, suggesting a biological meaning for the observed correlations. The observed connectivities between metabolite and cytokines hubs showed higher numbers of literature hits as compared to the number of literature hits on the connectivities between the microbiota and gene expression internal phenotypes. We conclude that multi-level integrated networks may help to generate hypotheses and to design experiments aiming to further close the gap between genotype and phenotype. We describe and/or hypothesize on the biological relevance of four identified multi-level connectivity hubs.

Keywords: data integration, internal phenotype, transcriptomics, proteomics, metabolomics, microbiota, gastrointestinal tract, systems biology

# INTRODUCTION

The information encoded in the genome (genotype) and the external quantitative traits or characteristics (phenotype) of an organism are linked to each other by several layers of so-called, intermediate (Leuchter et al., 2014; Fontanesi, 2016) or internal (Houle et al., 2010) phenotypes. Several of these internal phenotypic layers are shown in **Figure 1** that visualizes the conceptual relationship between the external phenotype (P), the genotype (G), the environment (E), and the G&E interactions. The epigenome is tightly associated with the genome and represents the programming of gene expression which is not dependent on the DNA code itself. The transcriptome layer represents direct effects of the environment on the gene expression of the

(epi-)genome. Translation of the transcriptome into proteins represents the next internal phenotype. The subsequent layer is represented by complex metabolite profiles. The organismassociated microbiota, especially those in the gut, can be regarded as a separate internal phenotypic layer, because it is not only dependent on the host genome but also heavily influenced by its environment, particularly by nutrition (Schwartz et al., 2012; Montiel-Castro et al., 2013). Although, for several traits the quantitative effects of the environment on the external phenotypes are known (Gentry et al., 2004; Cani et al., 2008; de Wit et al., 2011), the specific effects of the environment on the internal phenotypes are largely unknown. Furthermore, it is obvious to assume that the various layers of internal phenotypes are connected to each other and that their joint profiles ultimately determine the external phenotype (Leuchter et al., 2014; Fontanesi, 2016). Unfortunately, most of these assumptions are not based on solid evidence and at best represent oversimplifications of the dynamic nature of processes involved in determining external phenotypes. It, furthermore, partly explains the knowledge gap that exists between the genotype and the external phenotype.

Therefore, the objective of this study was to develop methodologies to identify components in the internal phenotypic layers that are connected to components in other internal phenotypic layers. To this end, we integrated multi-scale quantitative (–omics) data using a regression approach. The used data sets were derived from a single experiment with inbred mice which were exposed to five different dietary interventions as a means to perturb the different internal phenotypes. With a data-driven approach we were able to identify a large number of potential connections between the various intermediate phenotypes and for several we found proof of causal relationships in literature. We have used networks to represent the identified connections. The molecular components of each internal phenotype (such as genes, metabolites, cytokines, or bacterial groups) are represented as nodes in the network and the identified connections between each data type are represented as links or edges. The results of this study provide a basis to understand how various internal phenotypic layers are connected to each other. The identified connections may be crucial for the identification of causal relationships (Civelek and Lusis, 2014) between various biological scales and to uncover mechanisms involved in determining external phenotypes.

# MATERIALS AND METHODS

### Origin of Data

We used data from an experiment with 6-week old inbred mice that were fed for 4 weeks with six different semi-synthetic diets (Kar et al., submitted). In brief: thirty-six 21-day-old C57BL/6J mice (Harlan Laboratories, Horst, the Netherlands) were divided into 6 groups and housed in pairs with ad libitum access to diet and water. After adaptation for 1 week to a standard diet, the mice were fed semi-synthetic diets containing 300 g/kg (as fed basis) of one of the alternative protein sources for 28 days: soybean meal; casein; partially delactosed whey powder; spray dried plasma protein; wheat gluten meal and yellow meal worm. At the end of the experiment, mice were sacrificed to collect ileal tissue to acquire gene expression data, ileal digesta to study changes in microbiota, blood serum to profile cytokines and chemokines and blood and urine to profile amine metabolites. All procedures were approved by the Animal Experimentation Board at Wageningen University & Research Center (accession number 2012062.c) and carried out according to the guidelines of the European Council Directive 86/609/EEC dated November, 1986. Multi-omics data were obtained with regards to: whole genome gene expression profiles of ileal tissue as measured with Affymetrix GeneChip mouse gene 1.1 ST microarrays (Affymetrix, Santa Clara, CA, USA); community scale microbiota composition of ileal digesta by targeted-amplicon DNA sequencing of the bacterial 16S rDNA V3 region on an Illumina Mi-Seq sequencer; 23 serum cytokine and chemokine concentrations (pg/ml) using a Bio-Rad Mouse 23-plex kit (Bio-Rad, Hercules, CA, USA); and amine metabolic profiles of serum and urine using an ACQUITY UPLC system coupled online with a Xevo Tandem quadrupole mass spectrometer (Waters) operated using QuanLynx data acquisition software (version 4.1; Waters; Kar et al., in preparation). The data from the ileum reflects the local effects of the dietary interventions, the other three data assess the systemic effects.

# Pre-processing and Selection of Data

An overview of the five types of data and their specifics are given in **Table 1**. Supplementary Figure 1 has an overview of all the analytical methods used in this study. Each dataset was preprocessed in a similar way using the R package limma (Smyth, 2005) to find the differentially significant data-points. The data

TABLE 1 | Pre-processing and specificities of each data-type.

Details of the site of sampling and data dimensions before and after pre-processing are indicated. The first number indicates the number of variables in the data and the second number denotes the number of samples.

is first log transformed and then this data is fitted to a linear model using the function lmFit (Phipson et al., 2016) which will give back information on the differences between the data-points in different samples and subsequently different comparisons of and the other, **Y** with ny elements, the response. As an output, the approach produces a matrix Ma(**X,Y**) of size nx × ny representing the relevant correlations between both datasets, so that:

$$ma\_{\overrightarrow{i}} = \begin{cases} \text{0, if } \mathbf{Y}\_{\overrightarrow{j}} \text{ independent of } \mathbf{X}\_{\overrightarrow{i}}\\ \text{or } \begin{pmatrix} \mathbf{X}\_{\overrightarrow{i}}, \mathbf{Y}\_{\overrightarrow{j}} \end{pmatrix}, \text{ if } \mathbf{Y}\_{\overrightarrow{j}} \text{ dependent on } \mathbf{X}\_{\overrightarrow{i}} \end{cases}, \text{ with } \mathbf{i} \in \{1, \ldots, n\_{\overrightarrow{\mathbf{x}}}\} \text{ and } \mathbf{j} \in \{1, \ldots, n\_{\overrightarrow{\mathbf{y}}}\} \tag{1}$$

control vs. treatment. Then we used the function eBayes (Phipson et al., 2016) which applies an empirical Bayes method to compute p-values for a t-statistic under the assumption that only 1% of the data-points are differentially regulated among all the data-points in the samples. This p-value is then subjected to a Benjamini–Hochberg (Benjamini and Hochberg, 1995) multiple testing correction, also known as a False Discovery Rate (FDR).

This analysis was done by comparing the data of each dietary group against the data of the dietary group that received soy bean meal as protein source, which is the most common source of protein in animal diets. The FDR value of the data, is used to gauge significance and data-points that were significant in at least one of the five comparisons of the diets were included in the integration analysis. Except for the Cytokine and Metabolomics Serum (using the amine measurement), all the data-types had some samples thrown out due to quality control. Two types of metabolomics measurements were done on the sampled urine; Amine and Acyl-carnitine, the amine dataset did not have sufficient statistically significant data-points so was discarded. We only work with the Acyl-carnitine measurement in urine.

### Data Integration, Network Generation, and Network Assessment

All significantly different data-points were used in the integration which was initially performed with two datasets at a time, so that from the 5 datasets 10 integrated networks were generated. The integration was performed using the function sPLS (sparse Partial Least Squares) in regression mode with ncomp = 5, from the R package mixOmics (Lê Cao et al., 2009; Dejean et al., 2011; González et al., 2012). The regression mode is used to model causal relationship between variables in both datasets by identifying combinations of variables between both datasets. Weight vectors used in the regression modeling are termed loading vectors. sPLS is used to perform simultaneous variable selection in the two datasets to be integrated and employs LASSO (Least Absolute Shrinkage and Selection Operator) penalization (Tibshirani, 2011) on the loading vectors. This approach requires one data set, **X** with nx elements, to be designated the predictor Where cor Xi , Yj is Pearson's correlation between elements i and j from datasets **X** and **Y**, respectively. The correlation is computed across all available samples (here corresponding to dietary exposures).

Since it is not trivial to determine the predictor and response with biological data, we swapped the two types of data to compute Mb(**Y**,**X**), a matrix of size ny × nx where the roles of **X** and **Y** have been interchanged. Both matrices, Ma and Mb where combined into a final matrix M(**X**,**Y**) size nx × ny using

$$M(\mathbf{X}, \mathbf{Y}) = \
Ma(\mathbf{X}, \mathbf{Y}) \; + \; t(\!\!\!Mb(\mathbf{Y} \!\!X) \!\!\!) \tag{2}$$

where t represents matrix transposition. Thus, non-null elements of the matrix M(**X**,**Y**) represent correlations between data types that have been deemed associated. This matrix can be seen as a weighted adjacency matrix representing a network where two nodes Xi and Yj are connected via an edge if a non-null weight can be assigned to the edge. This weight is represented by the matrix value mij.

To further prune the network of (possibly) spurious interaction two additional thresholds (thl < 0; and thh > 0) were imposed to obtain an unweighted adjacency matrix A(**X**, **Y**)of size nx × ny

$$A\_{ij} = \begin{cases} \begin{array}{c} 1 \text{ if } m\_{ij} \ge \
th\_h \text{ or } m\_{ij} \le \
th\_l\\ 0 \text{ if } \begin{vmatrix} m\_{ij} \end{vmatrix} < \begin{vmatrix} th\_l \end{vmatrix} \text{ and } m\_{ij} < \
th\_h \end{vmatrix} \end{cases} \tag{3}$$

were |x| represents the absolute value. thl and thh where selected for each network so that only top 5% of the highest (positive) and lowest (negative) weights were kept for building the networks.

Networks represented by these adjacency matrix were transformed into the edge-list format, a two column table of the connected nodes in a network were each row represents an edge and visualized in Cytoscape (Shannon et al., 2003; Ono et al., 2015).

For each pair of integrated datasets a null model of the association networks was constructed using a strategy based on random permutations of measured values (Saccenti et al., 2015). Measured data-points were randomly permuted over samples before data integration to obtain randomized datasets that still retained the same value distribution for each variable. The randomized datasets were then used for data integration following the afore mentioned approach thereby generating randomized associations networks. The process was iterated Nit = 1,000 times for each pair of datasets; For each iteration, k, the values of the dynamic cut-offs (thlk and thhk) (5% of the highest and lowest correlation) were recorded. For the 10 pairwise combinations of datasets, the values obtained for thl and thh obtained using the unpermuted dataset, were compared with the distribution of values of thlk and thhk with k = {1,...,Nit} to get networks from the random data to compare to the networks from the biological data.

### Network Merging and Topological Analysis

The 10 networks arising from pair-wise data integration of the 5 data sets were merged in a combined network including all the nodes and edges of the 10 networks. This network is then restricted by only including nodes present in at least two of the separate networks. We used the igraph R package (Csardi and Nepusz, 2006) to further analyze the network, which was treated as non-directed, since no particular directionality was assigned to the edges. We obtained values for the following topological properties of the merged network (Barabasi and Oltvai, 2004; Csardi and Nepusz, 2006; Zhu et al., 2007): Degree: number of neighbors of a given node, that is the number of nodes connected to it. Clustering coefficient of a node is the ratio of the number of connections between the neighbors of a node and the total number of possible connections between said neighbors. Characteristic path length: median of the average distance between a node and all the rest. Network density: ratio between the total number of existing edges and the total number of possible edges (given the number of nodes in the network). Connected components maximal subgraphs in a network such that each node is connected to all the rest by means of network paths. For node level metrics, such as degree or clustering coefficient average values were computed over all nodes. Cytoscape was used for network visualization.

### Literature Mining

To investigate the co-occurrence of the names of the connected nodes in the association network, we used the R package rentrez (Winter, 2016). This package searches for selected keywords in PubMed abstracts while making use of the MeSH (Medical Subject Headings) thesaurus to maximize results via the API from NCBI. The search was not restricted to a specific tissue type or organism. These results were examined, although not exhaustively, to find literature evidence of established relationships between nodes connected through identified edges; these were then considered as true positive search results.

The script used to generate all these results will be made available on request. All the above mentioned operations were performed using existing functions from R packages. The different steps involved are represented in Supplementary Figure 1.

### RESULTS

### Analysis of the Individual Datasets

A dietary intervention was performed on mice where the protein content was changed and multi-omics data were obtained with regard to: whole genome gene expression profiles of ileal tissue (Transcriptomics), community scale microbiota composition of ileal digesta (Microbiota), 24 different cytokine levels in blood serum (Cytokine), and protein-associated metabolic profiles of serum (Metabolomics Serum) and urine (Metabolomics Urine). These data were pre-processed and analyzed separately by fitting a linear model on the data-points and looking for differentially expressed readouts in each treatment vs. the control. Each dataset had its own p-value (corrected for multiple testing with the Benjamini–Hochberg method) threshold, ranging from 0.001 to 0.1 for difference between the tested and reference diets. The highest number of statistically significant entities was found in Transcriptomics. Furthermore, all the measured variables in Metabolomics Urine were found to be significantly different in at least one comparison.

### Pairwise Data Association and Network Generation

We performed the integration by linking two data-types at a time and in such a way that after the pairwise analysis all the observed association data could be combined to build a multi-level interaction network. Therefore, each data-type was integrated with the other four types of data, resulting in 10 association networks. The topological characteristics of all these 10 networks are given in **Table 2** and **Figure 2**, and the network graphs are available in Supplementary Figure 2 as an image. Data Sheet 1 has the networks in a format that can be uploaded into Cytoscape in order to further explore the connectivities of these networks by simply clicking on these nodes. **Table 2** shows the positive and negative thresholds that were used separately for the association network. Connections between pairs of data points with correlation values between the threshold values, i.e., Low Threshold (negative threshold) and High Threshold (positive threshold) as indicated in **Table 2**, were discarded and the corresponding edges removed from the final network. There were two disconnected sub-graphs in five of the networks while the other five have only a single, fully connected graph. Supplementary Figure 3 shows the pattern of changes induced by the diet in three components of the network Microbiota & Transcriptomics.

The largest network, in terms of nodes, is the Microbiota & Transcriptomics network. This seems logical as it represents the most comprehensive datasets and spacial interactions between the two data-types are known to occur. Overall, networks involving Transcriptomics data had higher number of nodes than other networks. The smallest network with 18 nodes and 22 edges was the Metabolomics Urine & Cytokine network.

### Technical Validation of Pairwise Integration Networks by Random Permutation

We performed the same method of integration on the five different data-types after randomly permuting the measured

### TABLE 2 | The 10 individual correlation networks.

Each row represents one of the 10 correlation networks. Low Threshold and High Threshold represent the thresholds used for the correlation values. The 3rd and 4th columns have the number of nodes in the network that belong to the first and second data, respectively. The last column displays the number of connected graphs in the network.

data, this process was iterated a 1,000 times. In this way, the networks obtained from random permutations are considered a null model with no biological information, and used to assess the significance of the results obtained with the non-permuted data. **Figure 3** shows the spread of correlation values for the integration of Metabolomics Serum and Transcriptomics. The thresholds for network reconstruction were selected so that only the 5% highest and lowest correlations were kept. The separation between the values obtained for the integrated data and the randomly permutated datasets indicates the high significance of

the edges in the integration networks. In this way, selection of the 5% highest and lowest correlations and significant limits the number of spurious correlations that could be due to chance alone while retaining maximum information in the networks.

Similar results were obtained for most of the integration networks (Supplementary Figure 3). In three of the networks, there is an overlap between the correlation values from the inferred network and the values arising from the randomly generated networks. The overlaps are in the networks Metabolomics Urine & Microbiota, Metabolomics Urine & Cytokine, and Transcriptomics & Cytokine network. The highest overlap appears in the first two and mostly affects edges with negative correlations.

### Merged Network

All the 10 integration networks (Data Sheet 1) were merged and only nodes linked with nodes of at least two other datatypes were kept (see **Table 3**). The gene expression data has the highest number of nodes in the merged network. However, nodes with the highest degree (number of connecting edges) arise from the microbiota data, with S24-7 having 57 neighbors and Bifidobacterium having 47 neighbors. The merged network encompasses 45 nodes that are connected to all the other types of data. For that reason we denote them "Connectivity hubs" and they are included in **Table 3** and Supplementary Table 1.

# Functional Validation of Merged Network by Text Mining

A PubMed literature search for co-occurrence of linked nodes gave results for 6% of the links corresponding to 37 edges. We further investigated reported causality effects between the nodes in question. Most of the retrieved results are related to metabolites and cytokines measurements whereas a few results confirming causal relationships were found involving gene nodes. We were able to find literature confirmation pertaining to associations for six out of the 10 pair-wise connections between phenotypes, as summarized in **Table 4**. Supplementary Table 2 contains all the PubMed identifiers from the literature mining and Supplementary Table 3 has phrases from a maximum of three PubMed abstracts from the results. Among the nodes with literature results, four are from Microbiota, two from Transcriptomics, 15 from Metabolomics Serum, three from Metabolomics Urine, and six from Cytokines. The node with the highest number of hits in literature is Tnfa which co-occurs 8,563 times with nine metabolites from the Metabolomics Serum data and one bacterial group (Bifidobacterium).

TABLE 3 | Characteristics of the merged network.

Characteristics of the merged correlation network. The number of nodes from each datatype are given in rows three to seven. Between brackets the number of connectivity hubs is indicated.

Of the 30 data-points from all the types of data that have literature results, 15 are connectivity hubs. One such connectivity hub is Glutathione (GSH) which has 21 direct neighbors from four data-types as shown in **Figure 4**. This hub is especially interesting because six of the connected nodes (Carnitine, Tnfa, Il-1b, Il17c, Bifidobacterium, and Dapk2) have textual cooccurrences found by the text mining algorithm. The terms GSH and Tnfa were found 2,231 times in the abstracts of Pubmed indexed articles. Full text inspection shows that some of the connections are causal relationships as one of the connected nodes activates or inhibits the other.

# DISCUSSION

In this study we developed and used a set of computational methods to identify components in internal phenotypic layers that are connected to components in other internal phenotypic layers of an organism. We successfully integrated multi-scale quantitative (-omics) data, derived from a single experiment with inbred mice and which were exposed to five different diets. Here the mice had been exposed to the dietary intervention for 4 weeks. Four weeks is a significant amount of time in the life of mice and previous studies comparing the development of mice and humans (specifically the immune system in Holladay and Smialowicz, 2000) indicate that the development of different systems is much faster in mice than in humans. Hence it is reasonable to assume that the mice have adapted to the new diet in 4 weeks. Since the data originated from an animal experiment that was not designed for the detection of genetically and/or dietarily induced differences in external phenotypes, we only focused on the connectivity between 5 intermediate phenotypic levels. Some studies have reported pairwise data integration of two (Lu et al., 2014; Rajasundaram et al., 2014; Benis et al., 2015) or three data sets (Adourian et al., 2008). But this is, to the best of our knowledge, the first time that an integration of such heterogeneous data-types from different tissues, arising from a single experiment, has been reported. The approach as described here could, in principle, be applied on any number and type of datasets, as long as they are from the same experiment, from samples at the same time-point and have comparable dimensions of differentially regulated data.

The first column shows the types of data that are connected by the edges that were found in the PubMed literature search.

# Internal Phenotypic Data and Pairwise Data Integration

Each used data-type represents a different internal phenotype and a different layer of the system that (co-) drives the manifestation of external phenotypes. We subjected each data-type to a separate analysis in order to correlate only those changes induced by the dietary intervention. Nodes with significantly different values could easily be identified in each of the sampled tissues and fluids (ileum, blood, and urine) thereby representing the local and systemic effects of the interventions and the need of a multi-scale approach.

In order to investigate connections between the five datatypes we used sPLS, an integration method that can be applied to several types of data, two at a time. This method can also handle the dimensionality problem of biological datasets where the number of variables is usually higher than the number of samples. sPLS has been previously used for integration of microbiota with gene expression data (Benis et al., 2015; Steegenga et al., 2016), and measurements on cell wall polysaccharides of fibers with phenotypic characterizations of fibers in cotton balls (Rajasundaram et al., 2014).

We performed pairwise integration of the datasets, resulting in 10 networks with varying spreads of correlation values. Deciding on a threshold to distinguish genuine from spurious correlations is a major bottleneck for the definition of association networks. While a 0.8 threshold (absolute value) has been suggested for gene expression data (Schäfer and Strimmer, 2005), other authors suggested smaller values (0.6) in metabolomics data sets (Camacho et al., 2005). The correlation values greatly depend on the biological dataset under study and its dimensionality. There are several methods to choose a threshold based on the data: use assigned p-values as threshold; use network characteristics of the correlations; or use a percentage of the correlation distribution. When evaluated by Borate et al. (2009) they concluded that threshold selection methods based on network properties such as the clustering coefficient are best for gene co-expression networks. This would not work here because the generated networks always induce connections between data points of different type and as a result they have a zero clustering coefficient for every node. While integrating two types of metabolomics datasets with gene expression of the tissues in which they were measured Adourian et al. (2008) assigned p-values to the correlation values and then set a threshold. Selecting a threshold is further complicated by the possible appearance of spurious correlations due to a common response variable influencing the connecting nodes (A is correlated to B, A is correlated to C, therefore, B and C appear correlated). Regarding gene expression data, multiple methods (reviewed for example in Marbach et al., 2012) have been developed to minimize the number of falsely predicted associations. In this study, we used the top 5% of the correlation values because this dynamic threshold (separate for the positive and negative values) eliminates bias toward the size of the datasets. To further evaluate the impact of the correlation scores we have inspected the correlations between some linked nodes. Supplementary Figure 3 shows an extreme case in which transcript abundance of two genes negatively correlated with the abundance of a bacterial group. This might induce a spurious association between the genes. Spurious associations due to a common response variable influencing the connecting nodes are more likely to appear when both nodes are of the same type. Therefore, to further minimize the number of spurious associations we have focused on associations between different internal phenotypes.

We further validated the observed correlations by comparing them with a null model obtained by randomly permuting the data along the samples (Eguíluz et al., 2005; Saccenti et al., 2015). In the randomly permuted samples we expect all inferred associations to be spurious, as the permutation process destroys any possible correlation between the variables. In that case, even the correlations corresponding to the highest and lowest 5% of the population would be spurious. The values of the correlations deemed significant in the experimental data sets are found to be higher than these false positives. In two of the networks, Metabolomics Urine & Microbiota and Metabolomics Urine & Cytokine (the smallest network), the significance of the negative correlation values could not be established as we observed a considerable overlap between the negative correlation values of this network and the negative thresholds of the random networks. This calls for caution when biologically interpreting these networks. For five of the networks we observed a very clear separation of the random thresholds and the start of the correlation values in the network (Supplementary Figure 4). The other networks showed slight overlaps between the random threshold distribution and the network correlation distribution. This extra validation step reassured us that the observed correlations are rooted in biological phenomena. To our knowledge this technical validation step is not common in current studies of this type.

The edges of the inferred networks, indicate significant computationally-determined correlations between values of connected nodes. Our approach does not require a mechanistic model on how the associations are established and in each network these associations may be caused through entirely different mechanisms. In some cases the associations would be due to causal relationships between the connected nodes, such as increased expression levels of a cytokine gene linked to increased cytokine levels. However, in many cases, the associations could be indirect, mediated by elements that have not been measured in the experimental set up. In a formal mathematical model, they are considered hidden variables. Such would be the case of, for example, the changes in the metabolite levels of urine. These changes might have been caused by the colonic microbiota, in turn affected by the ileal microbiota. Since we only used the ileal microbiota data, we observe correlations between the ileal microbial populations and the urine metabolite levels which could be in reality, indirect relationships mediated by the colonic microbiota.

### Network of Connected Internal Phenotypes

The pair-wise integration method allowed us to merge the 10 individual networks into a single network. Correlations within a dataset were deliberately excluded from this study because we only wanted to focus on connections between different internal phenotypes, where little work has been done. Thus, in the 10 networks, all detected connections are between two different data types and every node has a zero clustering coefficient. However, in the merged network, a non-zero clustering coefficient emerges as a result of nodes connecting to multiple data types (**Table 3**). This emphasizes the biological relevance of this method because the 10 networks were built without any information on crosslinking. Thus, we identified individual nodes that directly or indirectly participate in processes of the other four individual networks. Because they seem to connect different internal phenotypes, we denoted them "Connectivity Hubs." Starting the procedure as developed and applied here with networks with non-zero clustering coefficients (correlating within a dataset) would, however, not alter the connections between internal phenotypes.

# Functional Validations of Phenotype Connections

Results of the text-mining were used to validate some of the identified links. This revealed insights into the mechanistic relationships between the variables predicted to be linked to each other. Thirty-seven of the 577 (6%) computational inferred links have already been described in literature as detected by our text-mining approach, which was not exhaustive because it focused only on text in journal abstracts. This indicates that our method identifies currently known biological interactions. The rest of the predicted links have not been discovered and investigated yet, have not been mentioned in abstracts, or do not exist in the biological system. Furthermore, by inspecting some of the retrieved abstracts and corresponding articles, we were even able to find causal relationships between some of the computational identified nodes where one of the nodes was used as an experimental perturbation and the other node was measured as a response parameter. Some examples are shown in Supplementary Table 3. Several indirect associations were also validated through reports on experiments where nodes, found to be connected in this study, were measured in response to another perturbation. During text-mining, in order to retrieve as many results as possible, search terms were matched against the MeSH thesaurus, irrespective of the organism, and all the synonyms were included in the search. The downside to this approach is the inclusion of several false textual associations. The most striking case is that of the identified association between Glutathione and Il17c. In the literature results, the reported association is between Glutathione and Il17a and not Il17c. Through the thesaurus, Il17c was mapped to Il17 and subsequently to Il17a thereby giving rise to that falsely identified association in literature.

In order to increase the precision and recall of text mining searches, and overcome problems associated to the use of a thesaurus, one needs to move from mining text, to mining the knowledge embedded in the text and the use of data hidden in public databases. Such an approach requires the use of knowledge management tools and representations that can be automatically accessed (Antezana et al., 2009). Semantic web technologies represent a new class of tools that include natural language processing, ontologies, machine learning algorithms and much more to facilitate integration knowledge from heterogeneous sources. The expansion of the use of semantic technologies in the life sciences domain will allow associating concepts such that inferences on causality, regulation, organism, or tissue can be made using high-throughput methods and automated reasoning.

Among the interactions retrieved from the automated literature search, a high prevalence of associations involving cytokines and/or metabolites was observed. In fact, such type of interactions represent 97% of the retrieved results. This probably highlights the extraordinary amount of work that has been done in these types of data in the past. On the opposite extreme, only 8% of the retrieved interactions involved associations between the expression of genes, reflecting the fact that most of the available gene expression data originates from genome-wide techniques. In such type of experiments, papers, especially abstracts, usually report on systems behaviors and pathways and less frequently on the individual behavior or role of individual genes and connected response nodes.

### Validated Connectivity Hubs

Even though we only performed integrations of two datasets at a time, we find data-points (metabolites, cytokines, genes, or microbial groups) that correlate with different types of data. We identified 45 connectivity hubs in the merged network that seem to have associations with all four types of data. More than 30% of them are involved in links that were retrieved in literature. To further support the biological relevance of identified multilevel connectivities we discuss the implications of two of the 15 biologically validated connectivity hubs as examples. The two connectivity hubs were chosen because of the large amount of literature results for these hubs. The first hub, Tnfa has the highest number of literature results among all the nodes in the network and the other hub, Glutathione, has literature validations to the most number of data-types.

Tnfa is a connectivity hub in the merged network, with links to several neighbors belonging to the four other types of data. The position of this cytokine in our merged network shows that it plays a role in processes of the other internal phenotypes. The literature validated links are between Tnfa and two other types of data (Metabolomics Serum, Microbiota). Many of the validated links represent causal relationships. With regards to immune responses and as a drug target, Tnfa has been studied in great detail (Cicha and Urschel, 2015). The un-validated edges show that Tnfa could be a regulator of other internal phenotypes as well, than currently known.

The metabolite Glutathione (GSH) was measured in the serum and in the merged network is a connectivity hub proving that it is vital part of the system that connects several internal phenotypes. Among the 15 connectivity hubs with functionally validated links, GSH is the only one that has validated links to all other data-types based on our literature mining. These results support our claim of GSH being a connectivity hub, a biological component influencing several internal phenotypes. Several PubMed results for GSH are from in-vivo studies where GSH was administered to alleviate symptoms of a disease. Our literature results show that GSH has been studied in relation to all different types of data. Of the six validated links in our merged network, five represent proven causal relationships (see **Figure 4** and discussion of the functional validation). These neighboring nodes in the merged network are mostly related to immune and homeostatic mechanisms. GSH is a tripeptide, ubiquitously distributed in living cells and plays an important role in the intracellular defense mechanism against oxidative stress (Diaz-Vivancos et al., 2015; Couto et al., 2016). It is known that GSH metabolism is very important for the antioxidant and detoxifying action of the intestine. It is also essential for the maintenance of the luminal thiol-disulfide ratio involved in regulation mechanisms of the protein activity of epithelial cells (Iantomasi et al., 1997) which could be important since the intervention is changes in protein. Our results also demonstrate the manifold and central role of GSH when it comes to proteins, peptides and amino acids in nutrition. These observations indicate that the presented merged network represents, at least in part, associations of biological phenomena.

### Potential Relevance of Selected Connectivity Hubs

There are 30 connectivity hubs in the merged network that do not co-occur with their connected nodes in our literature search. However, the prominence of these nodes in our merged network indicates that they could represent potential relevant interactions with components of the other internal phenotypes. In order to demonstrate how the results of this study may be used to hypothesize on functional relationships between different molecular components, we here describe the potential biological relevance of two highly linked connectivity hubs, Tmem72 and S24-7. Both hubs are not yet described in literature abstracts in conjunction with other data-types.

The high number of connectivity hubs in the Transcriptomics layer suggest that the expression of several intestinal genes is involved in many more interactions than currently known. None of the observed Transcriptomics connectivity hubs popped-up in our literature mining results. The most highly connected Transcriptomics node, Tmem72 (Transmembrane Protein 72), has only been studied in the kidney so far (Habuka et al., 2014) and not much information is available on it. But in the merged network this node has 27 links to other data-types (can be visualized in Data Sheet 1), mostly to metabolites from both the metabolomics datasets. Based on this, we hypothesize that Tmem72 is not specific to the kidney and that it has some sort of communication function in intestinal mucosa as well. The fact that Tmem72 is a transmembrane protein is supportive for this. Given its observed links with different microbiota, metabolites, and cytokines, it might be involved in diverse interactions with other internal phenotypes. Based on such an hypothesis, targeted experimental designs may be developed in order to investigate the hypothesized "communication" function of Tmem72 in intestinal mucosal tissue.

The most highly linked node of the merged network is the bacterial family classification, S24-7, suggesting an important role for this species in gut functionality. In some of the inferred individual association networks we already found it to be linked to a high number of nodes. Unfortunately, this node is not represented in literature abstracts together with the here observed neighbors. However, there is compelling literature that shows this microbial classification to be a significant part of the gut microbial community structure (Harris et al., 2014; Jakobsson et al., 2015). This family classification does not have a good functional definition, yet several studies show that it could be an important player in the functionality of the gut (Evans et al., 2014; Harris et al., 2014; Rooks et al., 2014). The latter claims are in line with the high number of neighbors that S24-7 has in our merged network. The current technical inability to cultivate S24-7 is most certainly due to the absence of knowledge on S24-7 interactions. However, a recent in-silico study (Ormerod et al., 2016) shows that S24-7 species have the ability to survive on different types of carbohydrate sources, similar to the genus Bifidobacteria. In the merged network, the connectivity hubs S24-7 and Bifidobacteria, share the highest number of neighbors (directly linked nodes). Among them are 16 genes, and neither S24-7 nor Bifidobacteria have literature results with any of these genes. An enrichment analysis on these shared network gene neighbors shows that they are involved in functions related to linoleic and linolenic acid metabolism (data not shown). It is known that these fatty acids are produced by Bifidobacteria (Teran et al., 2015) and that they are involved in the maintenance of the epidermal barrier function (Muñoz-Garcia et al., 2014). The observation that in our network these genes are shared between S24-7 and Bifidobacteria underscores the here hypothesized importance of S24-7 and indicates that these two bacterial groups are indeed closely related in function as hypothesized before (Ormerod et al., 2016).

From the results described in this paper, we conclude that we successfully developed methodologies to identify components in internal phenotypic layers that are connected to components in other internal phenotypic layers. By integrating multi-scale quantitative (-omics) data using a regression approach, we were able to provide provisional insight into potential ways internal phenotypic layers are connected to each other, including those between local and systemic layers. By a technical and functional validations, we underscored the relevance of our findings. Based on data generated by this type of integrated approaches, hypothesis driven and targeted research may be developed to identify causal relationships between various biological scales in order to diminish our knowledge gap between genotype and external phenotype. In addition, by expanding comparable approaches by incorporating data on genetic diversity and/or variation in external phenotypes, this knowledge gap may be even further closed down. The analysis pipeline that we developed is very general. Here we demonstrated this pipeline with datasets that address only one of the multiple environmental factors that might affect the internal phenotypes, namely the diet. However, the approach is very general and can be adapted to any type or number of data sets describing the impact of other perturbations.

### AVAILABILITY OF DATA AND MATERIALS

Transcriptomics data has been uploaded into GEO with the accession number GSE84442. The microbiota data, the two metabolomics datasets and the cytokine data are available on request. The R scripts using functions from existing R packages are also available on request.

### REFERENCES

### AUTHOR CONTRIBUTIONS

NB performed the data analysis and prepared the manuscript. SK performed the animal experiment and contributed significantly to the biological interpretation of the results. VM contributed to the direction of the analysis and the manuscript. MS was involved in the animal experiment and helped with the direction and critical revision of the manuscript. DS was involved in the data analysis and biological interpretation. MS-D helped with the direction of the data analysis and the manuscript. All authors have read and approved of the final manuscript.

### FUNDING

This work has been financially supported by the Systems Biology Investment Programme of Wageningen University, KB-17-003.02-022.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys. 2017.00388/full#supplementary-material

on time of exposure. Environ. Health Perspect. 108(Suppl. 3), 463–473. doi: 10.1289/ehp.00108s3463

**Conflict of Interest Statement:** The author VM was employed by company LifeGlimmer GmbH. This author contributed to the data analysis and the writing of the manuscript.

All other 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 © 2017 Benis, Kar, Martins dos Santos, Smits, Schokker and Suarez-Diez. 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.

# Formal Modeling of mTOR Associated Biological Regulatory Network Reveals Novel Therapeutic Strategy for the Treatment of Cancer

Zurah Bibi 1 , Jamil Ahmad1 \*, Amnah Siddiqa1 , Rehan Z. Paracha1 , Tariq Saeed1 , Amjad Ali 2 , Hussnain Ahmed Janjua2 , Shakir Ullah3 , Emna Ben Abdallah4 and Olivier Roux 4

*1 Research Centre for Modeling and Simulation, National University of Sciences and Technology, Islamabad, Pakistan, 2 Atta-Ur-Rahman School of Applied Biosciences, National University of Sciences and Technology, Islamabad, Pakistan, 3 School of Business, Stratford University, Falls Church, VA, United States, 4 IRCCyN UMR Centre National de la Recherche Scientifique 6597, BP 92101, Nantes, France*

### Edited by:

*Kai Breuhahn, Heidelberg University, Germany*

### Reviewed by:

*Katharina Beuke, Lille University of Science and Technology, France Melanie Boerries, Deutsches Krebsforschungszentrum (DKFZ), Germany*

\*Correspondence:

*Jamil Ahmad<EMAIL_ADDRESS>
### Specialty section:

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

Received: *24 April 2017* Accepted: *30 May 2017* Published: *13 June 2017*

### Citation:

*Bibi Z, Ahmad J, Siddiqa A, Paracha RZ, Saeed T, Ali A, Janjua HA, Ullah S, Ben Abdallah E and Roux O (2017) Formal Modeling of mTOR Associated Biological Regulatory Network Reveals Novel Therapeutic Strategy for the Treatment of Cancer. Front. Physiol. 8:416. doi: 10.3389/fphys.2017.00416* Cellular homeostasis is a continuous phenomenon that if compromised can lead to several disorders including cancer. There is a need to understand the dynamics of cellular proliferation to get deeper insights into the prevalence of cancer. Mechanistic Target of Rapamycin (mTOR) is implicated as the central regulator of the metabolic pathway involved in growth whereas its two distinct complexes mTORC1 and mTORC2 perform particular functions in cellular propagation. To date, mTORC1 is a well defined therapeutic target to inhibit uncontrolled cell division, while the role of mTORC2 is not well characterized. Therefore, the current study is designed to understand the signaling dynamics of mTOR and its partner proteins such as PI3K, PTEN, mTORC2, PKB (Akt), mTORC1, and FOXO. For this purpose, a qualitative model of mTOR-associated Biological Regulatory Network (BRN) is constructed to predict its regulatory behaviors which may not be predictable otherwise. The depleted expression of PTEN and FOXO along with the overexpression of PI3K, mTORC2, mTORC1 and Akt is predicted as a stable steady state which is in accordance with their observed expression levels in the progression of various cancers. The qualitative model also predicts the homeostasis of all the entities in the form of qualitative cycles. The significant qualitative (discrete) cycle is identified by analyzing betweenness centralities of the qualitative (discrete) states. This cycle is further refined as a linear hybrid automaton model with the production (activation) and degradation (inhibition) time delays in order to analyze the real-time constraints for its existence. The analysis of the hybrid model provides a formal proof that during homeostasis the inhibition time delay of Akt is less than the inhibition time delay of mTORC2. In conclusion, our observations characterize that in homeostasis Akt is degraded with a faster rate than mTORC2 which suggests that the inhibition of Akt along with the activation of mTORC2 may be a better therapeutic strategy for the treatment of cancer.

Keywords: mTOR signaling pathway, SMBioNet, Biological regulatory networks (BRNs), René Thomas, Qualitative modeling, Model checking, Cancer

# INTRODUCTION

Cells need continuous supply of resources that maintain intracellular energy and require nutrient levels contributing to macromolecular biosynthesis and serving as an upstream regulator of cell size and growth rate (Schmelzle and Hall, 2000; Wullschleger et al., 2006; Avruch et al., 2006; Sonntag et al., 2012). Compromised growth homeostasis can lead to several diseases including metabolic disorders, aging and cancer (Zoncu et al., 2011). A serine/threonine protein kinase mTOR, a member of the phosphatidylinositol kinase related kinases (PIKKs) family (Schmelzle and Hall, 2000), acts as a central regulator of homeostasis during growth and starvation (Zoncu et al., 2011). Recent studies have shown that dysregulation in mTOR signaling could lead to cancer and other pathologies (Menon and Manning, 2008). In these studies, abnormally elevated levels of mTOR have been linked with several human cancers including prostate, pancreas, liver, breast, colorectal, urinary tract, and female reproductive organs. On the other hand, due to excess nutrients supply, hyperactivation of mTOR has also been implicated to cause diabetes (Zoncu et al., 2011). Moreover, being the central regulator of growth, mTOR also monitors the process of aging (Zoncu et al., 2011). mTOR functions in the form of two distinct complexes namely mTOR Complex 1 (mTORC1) and 2 (mTORC2) (Wullschleger et al., 2006; Guertin and Sabatini, 2007). These complexes are distinguished by their unique accessory proteins, i.e., raptor and PRAS40 in case of mTORC1 and rictor, Protor and mSin1 in case of mTORC2 (Hara et al., 2002; Kim et al., 2002; Sarbassov et al., 2004). The function of these accessory proteins is to specify their binding with different substrates and regulators (Hara et al., 2002; Kim et al., 2002; Nojima et al., 2003; Schalm et al., 2003; Wullschleger et al., 2005; Sancak et al., 2007; Pearce et al., 2007). Both mTORC1 and mTORC2 share also some components including mLST8 and Deptor that act as positive and negative regulators, respectively (Loewith et al., 2002).

### Signaling of mTOR

The mTOR signaling pathway (**Figure 1**) is initiated by insulin, insulin-like growth factor 1 (IGF1) and Ras (Laplante and Sabatini, 2012) along with others. The binding of insulin with insulin/IGF1 signaling (IIS) receptors causes its autophosphorylation with subsequent recruitment of insulin receptor substrate (IRS) with its cytosolic domain. Activated IRS then activates several downstream effector proteins including PI3K.

The role of PI3K pathway in metabolism and growth is well-established and its dysregulation could result in certain metabolic disorders and cancers (Carracedo and Pandolfi, 2008). In addition to the activation of JNK pathway (Vivanco et al., 2007) that down-regulates PTEN transcription and promotes cellular proliferation by hindering apoptosis, it is also involved in the activation of Akt (Carracedo and Pandolfi, 2008). PI3K phosphorylates and converts phosphatidylinositol (4,5)-bisphosphate (PIP2) to phosphatidylinositol (3,4,5) trisphosphate to (PIP3) (Engelman et al., 2006; Manning and

Cantley, 2007). This event is important for the activation of phosphatidylinositol dependent protein kinase 1 (PDPK1) that eventually stimulates Akt. The activation of Akt is achieved through phosphorylation at two sites i.e., Thr308 and Ser473. Activated PDPK1 phosphorylates Akt at position Thr308 whereas another protein mTORC2 phosphorylates it at Ser473 (Alessi et al., 1997; Sarbassov et al., 2005). Both these phosphorylation events are essential for the complete activation of Akt. Thobe et. al. examined the influence of PI3K pathway kinases on mTORC2 and found PI3K mediated regulation essential for mTORC2 recruitment and further activation Thobe et al. (2017).