The liver operates in a gluconeogenic mode during the shortterm fasting trajectory in the present study. In this state, the liver takes up amino acids, lactate, and glycerol to produce glucose and urea. The liver also takes up non-esterified fatty acids to produce ketone bodies. We constrained the uptake rates of amino acids, fatty acids, lactate, and glycerol, using values reported in the literature from in vivo measurements in rats undergoing short-term fasting (**Supplementary Table S2**).

### Application of Transcriptionally Inferred Metabolic Biomarker Response (TIMBR) Algorithm

Transcriptionally inferred metabolic biomarker response (TIMBR) is a recently published method developed for predicting changes in extracellular metabolites due to gene expression changes under defined physiological operating conditions by integrating those changes into genome-scale network reconstructions (see Blais et al., 2017 for details). In the present study, we applied TIMBR to predict metabolite changes during a 5–6-h window of short-term fasting, where gene expression changes have little influence on metabolic state (Ikeda et al., 2014), in contrast to the changes in the central carbon metabolism fluxes. TIMBR calculates the global network demand required for producing a metabolite (Xmet) by minimizing the weighted sum of fluxes across all reactions for each condition and metabolite, while satisfying the steady-state mass balance and a defined optimal fraction of maximum network production flux capability (νopt) to produce a metabolite as shown below:

$$\begin{aligned} X\_{\text{met}} &= \min \sum |\nu| \\ \text{s.t.} &\text{ : } \nu\_{\text{X}} \ge \nu\_{\text{opt}}; \nu\_{\text{lb}} < \nu < \nu\_{\text{ub}} \text{ ; } \text{S} \cdot \nu = \text{0} \end{aligned} \tag{1}$$

where ν is a vector of reaction fluxes and S is the stoichiometric matrix. We included boundary conditions for uptake and secretion rates into the algorithm by fixing the respective lower (νlb) and upper bounds (νub) of the metabolite exchange reactions (νex), as shown in Eq. (2). Similarly, we integrated measurements from 13C-labeled tracer studies for some of the central carbon metabolism fluxes into the TIMBR algorithm by constraining the lower and upper bounds of the respective reactions in the model (νmfa) (Eq. 3).

$$\nu\_{\rm lb} \prec \nu\_{\rm ex} \prec \nu\_{\rm ub} \tag{2}$$

$$\nu\_{\rm lb} \prec \nu\_{\rm rnfa} \prec \nu\_{\rm ub} \tag{3}$$

Using this method, we determined the relative production scores for all metabolites (Xraw) from 5 to 7 h (X5−7) and 10 to 13 h (X10−13) time points (Eq. 4), and then calculated the TIMBR production scores (Xs) as the z-transformed scores across all exchangeable metabolites (Eq. 5).

$$X\_{\text{raw}} = \frac{X\_{5-7} - X\_{10-13}}{X\_{5-7} + X\_{10-13}} \tag{4}$$

$$X\_{\rm s} = \frac{X\_{\rm raw} - \mu}{\sigma} \tag{5}$$

**Figure 1** shows the workflow for the application of the TIMBR algorithm (adapted from Pannala et al., 2018). We performed the model computations in MATLAB R2017b using the linear programming solver provided in the GNU Linear Programming Kit. We refer the reader to the original publication for detailed descriptions of the TIMBR algorithm and the corresponding computer codes (Blais et al., 2017).

### RESULTS AND DISCUSSION

### Liver Glucose Production and Glycogenolysis Fluxes Decrease With Fasting Duration

During fasting, the liver produces glucose by synthesizing it from glycerol, lactate, and amino acids, as well as by breaking down glycogen. **Figure 2** shows a schematic of the

metabolism to those changes.

liver glucose production pathways, which include reactions of glycogenolysis, gluconeogenesis, and the tricarboxylic acid cycle. The aforementioned fluxes are collectively termed central carbon fluxes. The flux values through individual reactions at 10 and 13 h of fasting (**Figure 3**, Studies 1–3) were measured by stable isotope tracer studies, and those at 5–7 h of fasting (**Figure 3**, Est. 5–7 h) were compiled from the literature under conditions similar to our studies. In all studies considered for flux values at 5–7 h, food was withdrawn at the beginning of the light cycle. To reduce the influence of potential confounding factors, we first obtained absolute flux of liver glucose production from Rossetti et al.'s (1993) study conducted in 322 g male Sprague-Dawley rats [standard error (SE) = 7 g, n = 35] fed standard chow under conscious unrestrained conditions. The fractional contribution of glycogen to liver glucose production (48%) at 5–7 h was reported to be invariant to rat strain, body weight, state of anesthesia, and measurement technique (Rossetti et al., 1993; Neese et al., 1995; Peroni et al., 1997; Sena et al., 2007; Jin et al., 2013). The remaining 52% of glucose output came from glycerol, and lactate and amino acids (Rossetti et al., 1993; Neese et al., 1995; Peroni et al., 1997; Sena et al., 2007; Jin et al., 2013). The reported range of glycerol contribution was 15–19% and that of lactate and amino acids was 37–41% (Peroni et al., 1997; Sena et al., 2007; Jin et al., 2013) in various rat strains and a wide range of body weights. We selected fractional contributions of glycerol and lactate from the study of Jin et al. (2013) where they used 324 g male Sprague-Dawley rats (SE = 4 g, n = 9). **Table 2** shows the fractional contributions of various precursors to liver glucose output at 5–7 h and **Figure 3** shows the absolute flux values.

Overall glucose output progressively decreased by 30% from 5 to 7 h until 13 h of fasting. Much of this reduction was due to a decrease in the flux of glycogenolysis, whose fractional contribution to glucose output decreased from 48% at 5 h to 2.3% at 13 h of fasting (**Table 2**). Thus, the contributions of the remaining precursors—glycerol, lactate, and amino acids—to glucose output remained nearly constant as absolute values but increased as fractions of glucose output. As a result, the absolute fluxes through the reactions downstream of glycogen breakdown (PYGL in **Figure 2**), beginning with glucose-6-phosphate isomerase (GPI in **Figure 2**) and ending in the tricarboxylic acid cycle at succinate dehydrogenase (SDH in **Figure 2**), were nearly

equal in magnitude at 10 and 13 h of fasting but higher than the values at 5–7 h of fasting (**Figure 3**).

The major conclusions from the central carbon flux data (**Figure 3**) were that glycogenolysis and overall glucose output decline with fasting duration. A key observation was that the glycogenolysis flux was almost completely depleted after 13 h of fasting. The flux analysis assumption that liver metabolism operated in a pseudo-steady state at 5–7 h and 10–13 h is consistent with numerous observations reported in the literature (McGarry et al., 1973; Rossetti et al., 1993). The 5–7-h time interval represented the end of an early post-absorptive period where glycogen breakdown contributed to half of the liver glucose output—which was followed by a steep decline in glycogenolysis and a steep increase in ketogenesis plateauing at the 10–13 h time interval. Although the absolute flux of gluconeogenesis from glycerol was nearly equal at all time points, the flux of gluconeogenesis from lactate and amino acids was higher at the 10–13-h time interval, which indicated the coupling of liver metabolism to extra-hepatic sources of precursors for gluconeogenesis after longer fasting durations. Finally, a key approximation in the central carbon flux analysis was that the

liver provided all of the glucose output. Although the kidney is also known to contribute to overall gluconeogenesis, its contribution is important only at fasting durations beyond 24 h (Mithieux et al., 2006). Together with previous evidence, our data suggest the presence of distinct metabolic states after 5–7 h and 10–13 h of fasting.

Abbreviated reaction names on the y-axis follow their definitions in the legend

# Metabolite Changes Observed During Short-Term Fasting

Plasma metabolites changed after short-term fasting (**Table 3**). Given the similarity in liver central carbon fluxes, we treated the 5-h (Studies 1 and 2) and 7-h (Study 3) fasting durations as early time points, and the 10-h (Studies 1 and 2) and 13-h (Study 3)

for Figure 2.

TABLE 2 | Fractional contributions of metabolic precursors glycogen, glycerol, and lactate and amino acids to liver glucose production at varying durations of fasting.

The fractional contribution of lactate alone to liver glucose production is shown in parentheses adjacent to that of lactate and amino acids in total. a (Rossetti et al., 1993; Jin et al., 2013). bThis study. c (Lopez et al., 1998).

TABLE 3 | Observed changes in metabolites between early (5–7 h) and late (10–13 h) time intervals, experimentally measured in the plasma, and in the subset that is represented in the rat metabolic network model as exchangeable between the hepatocyte and plasma.

durations as later time points for determining metabolite foldchange values and their statistical significance. Of the 884 metabolites observed across the three studies, 198 changed significantly (p < 0.01). Of these, 39 metabolites were represented in the rat metabolic network model (iRno) as exchangeable between liver cells and the extracellular space or plasma. We compared our model predictions for the direction of change with fasting to those for the 39 metabolites, 33 of which showed an increase and 6 of which showed a decrease.

We also compared the significant changes in plasma metabolites observed in the present study to those reported in the literature on short-term fasting in the rat (McGarry et al., 1973; Ho, 1976; Brass and Hoppel, 1978; Palou et al., 1981; Kotal et al., 1996; Ikeda et al., 2014). In terms of major metabolite pathways, most of the changes reported in the literature were in agreement with those found in our study (**Table 4**). Important changes indicative of fasting were a reduction in glucose and phospholipids, and an elevation of ketone bodies, fatty acyl carnitines, corticosterone, and choline. Furthermore, key liver-specific metabolite changes observed here and in the literature were the elevation of primary and secondary bile acids, and the elevation of bile pigments bilirubin and biliverdin. **Supplementary Table S3** provides detailed lists of those metabolites and the entire summary of statistical analysis of all metabolites.

Reports on large-scale data on plasma metabolite changes during a short-term fast, the number of biological replicates required to resolve them, and their sensitivity to the type of vehicle administered, do not exist in the literature. The number of metabolites measured in Studies 1, 2, and 3 were 569, 645, and 633, respectively, where the vehicle administered to the rats was different for each study. The metabolite fold-change values needed to be pooled across the three studies to resolve metabolite changes above experimental and biological noise during shortterm fasting. The sum total of unique metabolites measured in the plasma in all three studies was 824 (**Table 3**), of which 420 were common to all three studies, 183 were common to exactly any two studies, and 221 were observed in exactly any one study. We calculated bootstrapped 99% confidence intervals of the foldchange values of the 420 common metabolites and confirmed that the vehicle was not a significant factor influencing metabolite changes (see **Supplementary Table S3**).

Among the 193 significantly changed metabolites (**Table 3**), 104 (54%) were measured in all three studies, 44 (23%) in exactly any two studies, and 45 (23%) in exactly any one study. Similarly, among the 631 unchanged metabolites, 316 (50%) were measured in all three studies, 139 (22%) in exactly any two

TABLE 4 | Concordance of observed changes in plasma metabolite data with reported changes in the literature due to short-term fasting.

studies, and 176 (27%) in exactly any one study. Taken together, there was no study-wise representation bias in the proportion of metabolites among the changed and unchanged groups, nor was there any differential effect of the vehicle on metabolite changes between studies, ensuring that pooling of metabolite fold-change data across studies was not confounded by known experimental differences between studies.

Of the 216 metabolites represented in iRno as exchangeable metabolites, 163 (76%) were measured in all three studies, which indicated the overall reliability of the data on exchangeable metabolites. Similarly, among the 39 significantly changed metabolites, 35 (90%) were measured in all three studies, which indicated the reliability of the metabolite data against which our model predictions were compared. Of the remaining four, N-carbamoylaspartate was measured in Study 3, acetylcarnitine in Study 1, inosine in Studies 1 and 2, and isocitrate in Studies 1 and 3.

Metabolite pathway annotations showed that lipids, amino acids, and cofactors and vitamins account for 49%, 23%, and 4% of the 824 metabolites, respectively, which indicated that lipid metabolites constituted the single largest category. Among the 193 metabolites that changed significantly, lipid metabolites again constituted the single largest group at 58%. The fraction of significantly changed lipid metabolites among all lipid metabolites was also highest at 28%, when compared to changes in other major pathways (19% or less). These results underscore the significance of lipids during short-term fasting.

### Metabolic Gene Expression Did Not Change Significantly During Short-Term Fasting

Gene expression changes in the liver during short-term fasting in all three studies (**Table 5**) revealed that the transcripts from each study mapped to a similar total number of genes (about 14,000), of which 2,258 were mapped to 2,240 in iRno. Out of the 2,325 genes in iRno, which were annotated with NCBI gene identifiers, 2,240 had 2,258 ENSEMBL gene identifiers that were used to annotate our transcriptomic data, with several genes mapping to than one ENSEMBL identifier. Based on the criteria of a false discovery rate of less than 0.1 and a biological effect size cutoff of 0.6 (corresponding to the 90th percentile), we found no statistically and biologically significant change in the expression of metabolic genes mapping to iRno in Studies 1 and 3 except for 100 genes in Study 2. Therefore, we did not use any differential gene expression-based weights in our implementation

TABLE 5 | Summary of gene expression changes with fasting.

of the TIMBR algorithm to predict plasma metabolite changes. **Supplementary Table S4** shows the results of the gene expression analysis.

### Liver Metabolism Accounts for 64% of Plasma Metabolite Changes

We integrated liver central carbon flux data, as well as known physiological flux bounds for metabolite exchange fluxes at early (after 5–7 h of fasting) and late (after 10–13 h of fasting) time points, with iRno using the TIMBR algorithm. We then used the TIMBR algorithm to compute a TIMBR score, whose positive or negative sign indicated the tendency of a metabolite to increase or decrease in the plasma, respectively, owing to changes in the liver metabolic network demand induced by fasting. The TIMBR predictions agreed overall with the metabolite changes observed here; TIMBR scores accurately predicted five out of six depressed, and 20 out of 33 elevated metabolites (**Table 6**). A summary of the 39 metabolites, their observed log2(foldchange) values, and corresponding TIMBR scores (**Figure 4**) revealed an overall accuracy of 64% for predicting any changes, and accuracies of 61% and 83% for predicting elevated and depressed metabolites, respectively. The probability that 64% or higher prediction accuracy could be achieved by chance was calculated to be 0.054, using the exact binomial test. Therefore, our network model of liver metabolism could account for 64% of plasma metabolite changes (increase or decrease) that were represented in the model, during short-term fasting.

The results in **Figure 4**, organized by metabolite pathways, revealed three major pathways represented in our data set: amino acids (8 metabolites), cofactors and vitamins (7 metabolites), and lipids (18 metabolites). The model accuracy in predicting metabolite changes for these three major pathways was 75% for amino acids, 42% for cofactors and vitamins, and 78% for lipids, providing estimates of both the reliability of the network model and the hepatic origin of metabolite changes in the pathways. In particular, the model achieved 100% accuracy in predicting the elevation of five primary and secondary bile acids (under lipids in **Figure 4**), and two bile pigments (under cofactors and vitamins), which are specific to the liver.

### Computational Model Assumptions, Limitations, and Interpretation of Predictions

The rat metabolic network model, iRno, currently the most comprehensive genome-scale model of rat metabolism, instantiated with physiological flux bounds pertinent to the

TABLE 6 | Concordance of TIMBR predictions with observed directions of change in metabolite data.

liver, was tested for satisfying defined liver-specific metabolic functionalities (Blais et al., 2017). The implicit assumption in our model was that overall liver metabolism could be represented by a single network with a representative set of physiological boundary conditions. This assumption seemed to contradict the known metabolic differences in hepatocytes between perivenous and periportal regions in the liver (Thurman et al., 1986). Despite not representing those different kinds of hepatocytes in our model, the overall satisfaction of liver metabolic tasks attested to a sufficient representation of liver metabolic functions originating in both regions. Additionally, the physiological flux bounds and central carbon fluxes employed to constrain the model did not include any metabolic heterogeneity. Finally, a key assumption in analyzing the model was that the network maintained a steady state, which was reasonable given the known metabolic flux conditions at 5–7 h and 10–13 h.

A limitation of our modeling analysis was the restricted coverage of metabolites exchanged between the plasma and liver cells. Additional curation of iRno, which included addition of exchange fluxes to improve network coverage of plasma metabolites, was limited by the paucity of literature evidence on the exchangeability of those metabolites. Consequently, the fraction of lipid metabolites among the 216 exchangeable metabolites (37%) was lower than that of the overall data set (49%). However, the fraction of lipid metabolites among the 39 significantly changing metabolites was higher at 46%, which is consistent with the trend in lipid metabolite fractions observed in the overall data set. Therefore, metabolite changes mapped to the network model are not biased by their limited coverage.

The measured changes in the circulating metabolites in plasma reflected the fasting response of the whole body. Our modeling effort sought to investigate plasma metabolite changes that can be associated with changes in liver metabolism under shortterm fasting conditions where the primary observation was a decrease in the hormonally regulated flux of liver glycogenolysis and no significant transcriptomic changes of liver enzymes (Lin and Accili, 2011). Our metabolic network analysis was made liver specific and relevant to liver metabolism by the flux constraints.

We used the in vivo central carbon fluxes derived from our tracerinfusion studies under short-term fasting conditions coupled with literature data from several studies during short-term fasting that sets the overall metabolite uptake and secretion fluxes of the liver (Lopez et al., 1998; Jin et al., 2013). This analysis assumed that the bulk of the glucose production flux captured by the in vivo metabolic flux analysis was of hepatic origin under these conditions (Hasenour et al., 2015). Thus, even though the measured metabolite changes were reflective of the overall systemic response, our computational analysis estimated those changes that were in concordance with a hepatic origin. To assess the impact of liver transcriptomic changes, we repeated our implementation of the TIMBR method using all of the transcriptomic changes regardless of their statistical significance and found that the predicted directions of metabolite changes were unaltered from those shown in **Figure 4** (see **Supplementary Figures S1–S3**).

Finally, the estimated model accuracy in predicting bile acids and bile pigments (100%, p = 0.004, subset of lipids), lipids (78%, p = 0.03), and amino acids (75%, p = 0.29) demonstrated the capability of the model to describe liver metabolic functions, and provided estimates of contributions of liver metabolism that agreed with metabolite changes observed in those pathways. In particular, lipid metabolite changes emerged as indicators of changes in liver metabolism, which were characterized both experimentally and computationally with sufficient statistical significance.

### CONCLUSION

Liver glycogenolysis became vanishingly small over the course of a short-term fast of 13 h, which resulted in a decline in the overall liver glucose output from 5 h until 13 h. Metabolites in plasma during this period showed changes known to be associated with short-term fasting, whereas liver gene expression did not change significantly. Finally, our computational analysis showed that two-thirds of the metabolite changes in plasma between 5–7 h and 10–13 h of fasting could be explained by central carbon flux changes in the liver without significant changes in gene expression.

### DATA AVAILABILITY

RNA-seq datasets generated for this study can be found in the Gene Expression Omnibus repository (accession numbers GSE123935, GSE124004, and GSE123987). Metabolomics data sets generated for this study are provided in **Supplementary Table S3**.

### AUTHOR CONTRIBUTIONS

KV conceived the short-term fasting analysis, analyzed RNAseq and metabolomics data, computed model predictions, and wrote the manuscript. VP curated the rat metabolic network model. MW and MR performed metabolic flux calculations based on isotope labeling measurements. SE performed all of the animal studies, including catheterization surgeries. IT processed plasma samples for metabolic flux analysis. TO collected and analyzed all of the blood samples. RP contributed to RNA extraction from tissue and purifications. JR helped to conceive and supervise the study, and helped to edit the manuscript. MS conceived the study, supervised and carried out the experiments on rats to generate the raw data, and helped to write the manuscript. JY conceived the study, supervised the metabolic flux analysis, and helped to write the manuscript. AW conceived and supervised the study, and helped to edit and write the final manuscript.

## FUNDING

This study was supported by the United States Army Medical Research and Materiel Command, Fort Detrick, MD, United States, as part of the United States Army's Network Science Initiative and under Contract No. W81XWH-14-C-0058 (to JY). Vanderbilt University Medical Center's Vanderbilt Technologies for Advanced Genomics (VANTAGE Core) was supported in part by CTSA grant (5UL1 RR024975-03), the Vanderbilt Ingram Cancer Center (P30 CA68485), the Vanderbilt Vision Center (P30 EY08126), and NIH/NCRR (G20 RR030956).

### ACKNOWLEDGMENTS

The authors gratefully acknowledge the assistance of Drs. Shivendra Tewari, Francisco Vital-Lopez, and Srinivas Laxminarayan for helpful discussions at various stages of this study. The authors also gratefully acknowledge Dr. Tatsuya Oyama contributions in editing the manuscript. The Vanderbilt University Medical Center's VANTAGE Core provided the genome-wide RNA sequencing data; Metabolon Inc. provided the global metabolic profiling data and some technical assistance for this work. The opinions and assertions contained herein are the private views of the authors and are not to be construed as official or as reflecting the views of the United States Army, the United States Department of Defense, or the Henry M. Jackson Foundation for the Advancement of Military Medicine, Inc. This paper has been approved for public release with unlimited distribution.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2019.00161/full#supplementary-material

FIGURE S1 | Binary heat map of TIMBR scores of significantly changed exchangeable metabolites in plasma represented in iRno compared with measured fold-change values, computed using both liver physiological flux changes and gene expression changes in Study 1 and grouped by major biochemical pathways: amino acid, carbohydrate, cofactors, and vitamins, TCA cycle, lipid, nucleotide, and peptide. The values in the left-hand side column (data) are measured log2 (fold change) values of metabolites, grouped as depressed

(black background), or elevated (white background) metabolites. The values in the right-hand side column are the computed TIMBR scores whose negative (black background) or positive (white background) sign indicates a predicted tendency of the metabolite to be depressed or elevated in plasma.

FIGURE S2 | Binary heat map of TIMBR scores of significantly changed exchangeable metabolites in plasma represented in iRno compared with measured fold-change values, computed using both liver physiological flux changes and gene expression changes in Study 2 and grouped by major biochemical pathways: amino acid, carbohydrate, cofactors, and vitamins, TCA cycle, lipid, nucleotide, and peptide. The values in the left-hand side column (data) are measured log2 (fold change) values of metabolites, grouped as depressed (black background), or elevated (white background) metabolites. The values in the right-hand side column are the computed TIMBR scores whose negative (black background) or positive (white background) sign indicates a predicted tendency of the metabolite to be depressed or elevated in plasma.

FIGURE S3 | Binary heat map of TIMBR scores of significantly changed exchangeable metabolites in plasma represented in iRno compared with measured fold-change values, computed using both liver physiological flux changes and gene expression changes in Study 3 and grouped by major biochemical pathways: amino acid, carbohydrate, cofactors, and vitamins, TCA cycle, lipid, nucleotide, and peptide. The values in the left-hand side column (data) are measured log2 (fold change) values of metabolites, grouped as depressed

### REFERENCES

(black background), or elevated (white background) metabolites. The values in the right-hand side column are the computed TIMBR scores whose negative (black background) or positive (white background) sign indicates a predicted tendency of the metabolite to be depressed or elevated in plasma.

TABLE S1 | Rat genome-scale metabolic network reconstruction iRno, with modifications.

TABLE S2 | Flux constraints for iRno during short-term fasting obtained from the literature.

TABLE S3 | Metabolomic imputed data from Studies 1–3, and the summary of statistical analysis of metabolite changes during short-term fasting in all three studies.

TABLE S4 | Results of RNA-seq data analysis, using the Wald and the likelihood ratio tests, to assess gene expression changes during short-term fasting in Studies 1–3.

DATA SHEET S1 | Rat genome-scale metabolic network reconstruction iRno in SBML format.

DATA SHEET S2 | Computer code to analyze metabolomics data.

DATA SHEET S3 | Computer code to reproduce TIMBR analysis results in Figure 1 and Supplementary Figures S1–S3.

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

Copyright © 2019 Vinnakota, Pannala, Wall, Rahim, Estes, Trenary, O'Brien, Printz, Reifman, Shiota, Young and Wallqvist. 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.

# Analysis of Genes Involved in Ulcerative Colitis Activity and Tumorigenesis Through Systematic Mining of Gene Co-expression Networks

Wanting Shi1,2† , Rongjun Zou3† , Minglei Yang4† , Lei Mai1 , Jiangnan Ren2 , Jialing Wen5,6 , Zhaoshi Liu1,2 and Renxu Lai1,2 \*

1 Department of Gastroenterology, Fifth Affiliated Hospital, Sun Yat-sen University, Zhuhai, China, 2 Digestive Endoscopy Center, Fifth Affiliated Hospital, Sun Yat-sen University, Zhuhai, China, 3 Department of Cardiovascular Surgery, Sun Yat-sen Memorial Hospital, Sun Yat-sen University, Guangzhou, China, 4 Department of Genetics, Zhongshan School of Medicine, Sun Yat-sen University, Guangzhou, China, 5 Guangdong Institute of Gastroenterology, Guangdong, China, 6 Guangdong Provincial Key Laboratory of Colorectal and Pelvic Floor Diseases, The Sixth Affiliated Hospital, Sun Yat-sen University, Guangzhou, China

### Edited by:

Steven Dooley, Universität Heidelberg, Germany

### Reviewed by:

Shikha Prasad, Northwestern University, United States Ingrid Arijs, KU Leuven, Belgium

### \*Correspondence:

Renxu Lai<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: 26 June 2018 Accepted: 09 May 2019 Published: 31 May 2019

### Citation:

Shi W, Zou R, Yang M, Mai L, Ren J, Wen J, Liu Z and Lai R (2019) Analysis of Genes Involved in Ulcerative Colitis Activity and Tumorigenesis Through Systematic Mining of Gene Co-expression Networks. Front. Physiol. 10:662. doi: 10.3389/fphys.2019.00662 Ulcerative colitis (UC) is an idiopathic, chronic inflammatory disorder of the colon, characterized by continuous mucosal inflammation. Recently, some studies have considered it as part of an inflammatory bowel disease-based global network. Herein, with the aim of identifying the underlying potential genetic mechanisms involved in the development of UC, multiple algorithms for weighted correlation network analysis (WGCNA), principal component analysis (PCA), and linear models for microarray data algorithm (LIMMA) were used to identify the hub genes. The map of platelet activation, ligand-receptor interaction, calcium signaling pathway, and cAMP signaling pathway showed significant links with UC development, and the hub genes CCR7, CXCL10, CXCL9, IDO1, MMP9, and VCAM1, which are associated with immune dysregulation and tumorigenesis in biological function, were found by multiple powerful bioinformatics methods. Analysis of The Cancer Genome Atlas (TCGA) also showed that the low expression of CCR7, CXCL10, CXCL9, and MMP9 may be correlated with a poor prognosis of overall survival (OS) in colorectal cancer (CRC) patients (all p < 0.05), while no significance detected in both of IDO1 and VCAM1. In addition, low expression of CCR7, CXCL10, CXCL9, MMP9, and IDO1 may be associated with a poor prognosis in recurrence free survival (RFS) time (all p < 0.05), but no significant difference was identified in VCAM1. Moreover, the NFKB1, FLI1, and STAT1 with the highest enrichment score were detected as the master regulators of hub genes. In summary, these results indicated the central role of the hub genes of CCR7, CXCL10, CXCL9, IDO1, VCAM1, and MMP9, in response to UC progression, as well as the development of UC to CRC, thus shedding light on the molecular mechanisms involved and assisting with drug target validation.

Keywords: ulcerative colitis, colorectal cancer, pathway enrichment, molecular mechanism, bioinformatics analysis

# INTRODUCTION

fphys-10-00662 May 29, 2019 Time: 19:42 # 2

Ulcerative colitis (UC) is a global, progressive and complex disease, the incidence of which is still growing, according to large-scale epidemiological statistics studies (Ng et al., 2018). Epidemiological reports shows that the highest annual incidence of UC was 24.3 per 100,000 person-years in Europe, 6.3 per 100,000 person-years in Asia and the Middle East, and 19.2 per 100,000 person-years in North America; and adjusted prevalences have exceeded 0.3% in many countries, especially in Europe and North America (Molodecky et al., 2012; Ng et al., 2018). A further concern is that the incidence of UC and the widespread use of therapeutic agents are associated with an increased risk of cancer (Biancone et al., 2015). Regarding the pathogenesis of UC, most of the emerging evidence supports the concept of an "inflammatory bowel disease (IBD) interactome," that is, UC is considered as part of a global disease network, with a complex interplay between host genetics, immunity, and environmental factors (Dabritz and Menheniott, 2014). According to this model, gene–environment interactions have pivotal roles in UC progression and mediate UC-related comorbidities and complications, including colitis-associated cancer (Dabritz and Menheniott, 2014). In the past two decades, novel genotyping and sequencing technologies, including RNA expression profiling, DNA methylation profiling, single-cell DNA analysis, chromatin immunoprecipitation sequencing, and RNA sequencing, have launched the era of genetic diseases; so far, 242 susceptibility loci and over 50 hub genes have been discovered in relation to IBD and various phenotypes (Mirkov et al., 2017). Moon et al. performed deep resequencing of UC-associated genes, showing that genetic variants of rs10035653 in C5orf55, rs41417449 in BTNL2, rs3117099 in HCG23, rs7192 in HLA-DRA, rs3744246 in ORMDL3, and rs713669 in IL17REL were significant (Moon et al., 2018). Hong et al. (2018) performed a trans-ethnic meta-analysis based on Asian IBD patients and subsequently identified three novel susceptibility loci at MYO10- BASP1, PPP2R3C/KIAA0391/PSMA6/NFKB1A, and LRRK1; as well as four previously known loci at NCF4, TSPAN32, CIITA, and VANGL2. Similarly, Peters et al. (2017) presented a predictive model of immune-related genes and further analyzed the functional and regulatory annotations based on genomewide association studies. Consequently, a driver set including DOCK2, GPSM3, AIF1, NCKAP1L, and DOK3 was selected, representing a high predictive efficiency in the integrated circuits of genetics, molecular, and clinical traits of IBD (Peters et al., 2017). The candidate biomarkers identification of UC activity and tumorigenesis in prior studies were presented in **Table 1**.

In summary, these novels genotyping and sequencing technologies and validated hub gene or susceptibility loci not only confer new regulators of pathophysiology, but also open a new horizon to find drug targets and redefine the disease's regulatory framework. However, these genetic variants combined only explain one in four cases of UC (Uhlig and Muise, 2017). The results also suggest that: (1) The genetic variants considered as personal pathogenic components cannot be isolated in the gene-environment network; (2) These hub genes show better statistical significance while loss of the functional and regulatory TABLE 1 | The candidate biomarkers identification of ulcerative colitis activity and tumorigenesis in prior studies.

annotations or may play an important part in protein-protein interaction (PPI) networks without statistical power; (3) some of the hub gene information may have been missed, owing to low abundance or small fold change (FC); and (4) diseasebased co-expression network analysis may further improve the mining efficiency beyond classical methods (Uhlig and Muise, 2017). Based on the above notes, we may apply linear models for microarray data power differential expression analyses (LIMMA), weighted correlation network analysis (WGCNA), and principal component analysis (PCA) to explore the hub gene regulatory network, using high-throughput gene expression arrays in UC, to further elucidate the molecular mechanisms of gene–environment interactions.

### MATERIALS AND METHODS

### Materials

Raw expression microarray array (CEL data) from the GSE13367, GSE38713, GSE16879, GSE48958, GSE75214, GSE4183, GSE37283, and GSE31106 datasets were downloaded from Gene Expression Omnibus1 (Barrett et al., 2013). Probe annotations and platform information were generated by matching with the GPL6244 (HuGene-1\_0-st) Affymetrix Human Gene 1.0 ST Array (Affymetrix, Santa Clara, CA, United States).

In this study, we analyzed the patients with colitis exclusively, no other IBD cases included. Here, GSE48958 and GSE75214 are matching with the GPL6244 (HuGene-1\_0-st) Affymetrix Human Gene 1.0 ST Array, while GSE13367, GSE38713, GSE16879, GSE37283, and GSE4183 are pairing with the GPL570 (HG-U133\_Plus\_2) Affymetrix Human Genome U133 Plus 2.0 Array and the GSE31106 is in line with GPL1261 (Mouse430\_2) Affymetrix Mouse Genome 430 2.0 Array. Total RNA extracted from mucosal biopsies was used to analyze mRNA expression via Affymetrix arrays, and corresponding grouping

1http://ncbi.nlm.nih.gov/geo/

information from each sample was subsequently pooled for further correlation analysis. Statistical analysis was performed with the R (version 3.3.2).

### Data Processing

fphys-10-00662 May 29, 2019 Time: 19:42 # 3

To remove bias and variability (resulting from heterogeneity and latent variables) from the high-throughput data for the different microarrays, the "ComBat" function in the SVA package was used to directly adjust the batch effects and latent variables (Leek et al., 2012). Subsequently, all of the microarray raw data analyzed using bioinformatics methods, including background correction, quantile normalization, and probe summarization of the expression values (Irizarry et al., 2003; Ritchie et al., 2015). Some advanced algorithms were used, including: (1) robust multi-array average for background-adjusted, normalized, and log-transformed probe expression values; (2) k-nearest-neighbor for displacing missing values of probes; (3) the t-test in the "LIMMA" package to identify differentially expressed genes (DEGs) in mucosal biopsy specimens from the comparative analysis among normal, UC, adenoma, and colorectal cancer (CRC) for GSE4183; and (4) the Benjamini–Hochberg method to adjust p-values and thus calculate the false discovery rate and FC (Irizarry et al., 2003; Ritchie et al., 2015). Gene expression values with | log2FC| > 1.5 and adjusted p-value < 0.05 were used to define DEGs. The co-annotated genes (a total of 16,653 genes) between GPL570 and GPL6244 platform were selected for further co-expression network analysis. The analysis strategy is presented in **Figure 1**.

### Weighted Co-expression Network Construction and Module Detection

The advantages of co-expression network analysis include the ability to integrate external information and avoid information loss in the case of low-abundance or small-FC genes. Systemslevel insight gives WGCNA an edge over other approaches (Langfelder and Horvath, 2008). Therefore, we carried out a systems-level analysis based on WGCNA. The analysis involved the following processes: (1) identifying the appropriate sample basing on the flash-Clust method; (2) selecting a "soft" threshold using the scale-free topology criterion; (3) identifying coexpression modules by employed the dynamic hybrid cut method; (4) relating the co-expression modules to sample traits based on the gene significance (GS) measures, which are defined as the statistical significance of the difference between the gene profile and the sample trait; and (5) accessing the interactions and connectivity of eigengenes among different co-expression modules by the topological overlap matrix method (Langfelder and Horvath, 2008).

# PPI Networks and Functional Enrichment Analysis

We accessed gene biological knowledge, protein functional associations, and PPIs with respect to genetic function, using a web-based analytic tool. The analysis flowchart as flowing that: (1) the gene ontology (GO) functions enrichment was extracted from the DAVID database2 (Huang et al., 2007) for annotation, visualization, and integrated discovery bioinformatics resources; GO terms for which p < 0.05 were considered to be significantly enriched in the gene modules of interest; and (2) the network of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was identified form Metascape database3 (Zhou et al., 2019), and p < 0.05 and enrichment score >1.0 was set as the cut-off criteria; (3) after carried out for the genes enriched in KEGG pathway of immunity, inflammation and tumorigenesis for the interesting gene modules, we subsequently constructed the PPI biological networks based on the STRING online database (V10.54 ) with the nodes association confidence score >0.4 (Szklarczyk et al., 2017). In addition, the Cytoscape software (V3.5.15 ) was used to visualize and evaluate interactions and identifying the hub gene in functional networks (Szklarczyk et al., 2017). The top 10 highest-degree nodes were defined as functional hub genes in the PPI network.

## Identification of Candidate Biomarkers Involving in UC-Associated Carcinogenesis

Additionally, raw data of GSE4183 was used to analysis the co-DEGs involved in UC, adenoma, and CRC. Here, we overlapped the co-DEGs and the PPI network's functional hub genes, which constructed by WGCNA key modules and identified as important parts in response to the pathway of UC immunity-inflammation and tumorigenesis, to detect the UC-associated carcinogenesis in hub genes. Additionally, the suitable dataset of GSE37283, including the expression profiling of UC with neoplasia, UC and normal mucosa samples, was used to validate the UC-associated carcinogenesis biomarkers; as well as the mouse dataset of GSE31106 involved in the multistep process of "inflammationdysplasia-cancer." The human and mouse genes were matched by Gene database6 (Brown et al., 2015).

The Cancer Genome Atlas (TCGA) colon adenocarcinoma normalized gene expression value (fragments per kilobase of exon model per million reads mapped, FPKM) were downloaded from the "TCGA biolinks" package (Colaprico et al., 2016). Subsequently, the FPKM data transformed into transcripts per kilobase million (TPM; Li et al., 2010), a comparable data type, which used to apply the survival analysis. The sample and corresponding clinical features were included in further survival and DEGs analysis.

# Investigating the Functional Role and Transcription Factor of Hub Genes

Importantly, the DAVID7 and Metascape database were used to explore the GO terms and KEGG pathway enrichment analysis of candidate targets, respectively. The enrichment cut-off criteria

2http://david.abcc.ncifcrf.gov/ 3http://metascape.org/

6https://www.ncbi.nlm.nih.gov/gene/

4http://string-db.org/

5http://cytoscape.org/

7http://david.abcc.ncifcrf.gov/

keep the same with the chapter and Section "PPI Networks and Functional Enrichment Analysis." Subsequently, to identify the transcription factor (TF) of the hub genes, the plug in iRegulon for Cytoscape software was applied, with the parameters were set to: (1) minimum identity between orthologous genes = 0.05; (2) maximum FDR for motif similarity = 0.001; and (3) normalized enrichment score (NES) ≥ 3.0 (Janky et al., 2014). Here, the top three regulators with the highest NES value were detected to construct the regulatory network involved in UC-associated carcinogenesis process.

### RESULTS

### Data Processing

The normal and UC mucosa without additional treatment from GPL6244 (GSE48958 and GSE75214) and GPL570 (GSE13367, GSE38713, and GSE16879) platform were selected in their entirety for further analysis, including the 58 normal, 55 UC inactive and 170 UC active samples (**Supplementary Table S1**). After merging the co-annotated genes, 16,653 genes were retained in further analysis (**Figure 2A** and **Supplementary Table S1**). The PCA of co-annotated genes in response to pre- and postcorrect the batch effects were showing in **Figure 2B**, which presenting a significant distinction between control, UC inactive and UC active samples.

# Construction of Co-expression Network and Gene Modules

After sample cluster analysis, the 283 samples with 16,653 gene variables were divided into 13 clusters (MEblack, MEblue, MEbrown, MEgreen, MEgreenyellow, MEgrey, MEmagenta, MEpink, MEpurple, MEred, MEtan, MEturquoise, and MEyellow; **Figure 2C**), and no samples removed in this process (**Supplementary Figure S1A**).

Following the WGCNA: (1) when the critical parameter of the power value was 12, the scale independence was up to 0.8 and had a higher mean connectivity (**Supplementary Figure S1B**); (2) two key modules and the relationship with the clinical traits were detected (Colitis: MEbrown Pearson coefficient = 0.94, p = 7E-134, MEgreen Pearson coefficient = 0.75, p = 2E-51; UC

FIGURE 2 | Co-expression modules construction and identify interesting modules of ulcerative colitis and activity. (A) Venn diagram showing the overlap of the co-annotated genes in GPL6244 and GPL570 platform. (B) The principal component analysis (PCA) for co-annotated genes among the various microarrays in response to pre- and post-adjusting of the batch effects of the status of ulcerative colitis activity. (C) Construction of co-expression modules based on a dynamic branch-cutting method. (D) The relationship between the co-expression modules and clinical traits. Red represents a positive correlation, and the green represents a negative correlation. (E) The connectivity of eigengenes. Red represents a positive correlation, and blue represents a negative correlation. (F) The PCA for interesting module genes in response to the status of ulcerative colitis activity. (G) The gene significance (GS) and module membership (MM) analysis of interesting modules in response to ulcerative colitis and activity.

active: MEbrown Pearson coefficient = 0.69, p = 2E-41; MEgreen Pearson coefficient = 0.86, p = 1E-85; **Figure 2D**); (3) the result of interaction analysis among co-expression modules suggested a high degree of independence among different module genes, such that the heatmap showed no significant interaction among module genes (**Supplementary Figure S1C**); and (4) The connectivity of eigengenes in different modules allowed us to identify three clusters identified, and the eigengenes of different modules within the same cluster showed significant connectivity, whereas there was no difference among different clusters' modules (**Figure 2E**). The interesting module gene list was presented in **Supplementary Table S2**.

Additionally, in related to UC activity, two-dimensional PCA results also showing satisfactory connectivity and distinguish ability of MEgreen and MEbrown module genes in response to UC and activity (**Figure 2F**; MEgreen: first principal component: 64.9%, second principal component: 7.2%; MEbrown: first principal component: 69.8%, second principal component: 23.9%). And, as shown in **Figure 2G**, the GS analysis results showing a tight correlation between the gene and the trait of colitis (MEbrown: Pearson coefficient = −0.41, p = 1.1E-23; MEgreen: Pearson coefficient = 0.99, p < 1E-200), as well as the trait of UC active (MEbrown: Pearson coefficient = −0.58, p = 7.2E-91; MEgreen: Pearson coefficient = 0.93, p = 2.6E-133). Similarly, the MEbrown and MEgreen module genes also presented a significant contribution to the module membership (MM; **Figure 2G**).

# Functional Enrichment Analysis and PPI Networks Construction

The MEgreen and MEbrown modules were assessed for further functional enrichment, consisting of GO term enrichment analysis of module genes of interest. Regarding GO terms enrichment, the MEgreen module was mainly enriched in GO: 0007041∼lysosomal transport (5 genes enriched; p = 1.11E-05), GO: 0098609∼cell-cell adhesion (13 genes enriched; p = 1.68E-04), and GO: 0043254∼regulation of protein complex assembly (4 genes enriched; p = 6.43E-04). Genes in the MEbrown module were predominantly enriched in GO: 0071805∼ion transmembrane transport (21 genes enriched; p = 3.98E-11), GO: 0042391∼regulation of membrane potential (13 genes enriched; p = 4.80E-07), and GO: 0034765∼regulation of ion transmembrane transport (18 genes enriched; p = 3.95E-09). These results are illustrated in **Figure 3A**. Regarding the KEGG pathway enrichment, the MEgreen module genes were significantly enriched in viral carcinogenesis (15 genes enriched; enrichment score = 5.01; p = 3.51E-05), proteoglycans in cancer (10 genes enriched; enrichment score = 3.97; p = 1.01E-03), and platelet activation (8 genes enriched; enrichment score = 4.91; p = 1.46E-03). However, genes in the MEbrown module were significantly enriched in ligand-receptor interaction (22 genes enriched; enrichment score = 3.55; p = 3.41E-07), calcium signaling pathway (13 genes enriched; enrichment score = 3.21; p = 1.72E-03), and adenosine 30 , 50 -cyclic monophosphate (cAMP) signaling pathway (12 genes enriched; enrichment score = 2.72; p = 2.36E-03).

Additionally, the **Figure 3B** illustrating a part of the visible pathway that tightly correlated with cancer from the KEGG pathway network (**Table 2**). These cancer-correlated pathways were related to immunity-inflammation response and tumorigenesis. These results also illustrated in **Table 2**. After submitting the genes enriched in cancer-correlated pathways to the STRING database, 50 and 48 PPI nodes were obtained for the MEbrown and MEgreen modules, respectively, with a confidence threshold greater than 0.4. After analyzed by Cytoscape software as an undirected method, the top 10 highest connectivity nodes of each PPI network were considered to be central agents. The PPI network of interesting modules was presented in **Figure 3B**.

# Identification of Candidate Biomarkers Involving in UC-Associated Carcinogenesis

What's more, 184 DEGs were obtained in the comparison of UC and CRC (160 DEGs down-regulated and 24 DEGs upregulated in CRC samples), and 344 DEGs were identified in the comparison of UC and adenoma in GSE4183 (332 DEGs down-regulated and 12 DEGs up-regulated in adenoma samples) (**Supplementary Table S3**). After overlapped, 106 co-DEGs were selected. Subsequently, we've further identified the same genes between co-DEGs and central agents of each PPI network, which have been selected as UC-related tumorigenesis genes. And 6 (CXCL10, VCAM1, CXCL9, MMP9, IDO1, and CCR7) out of the 106 co-DEGs remained after selection (**Figure 2B** and **Supplementary Table S5**). We also found a statistical difference in the gene expression levels of these genes between healthy individuals and UC patients (**Figure 4A** and **Supplementary Table S4**), as well as the expression levels between healthy individuals, UC, adenoma, and CRC patients of GSE4183 dataset (**Figure 4B** and **Supplementary Table S5**).

After validated by GSE37283 and GSE31106 datasets, the six hub gene expression levels in phases of UC were significantly increased in compared with the normal sample in both of the human and mouse's colonic mucosa (all the p < 0.05; **Supplementary Table S6**). Additionally, in comparison with UC mucosa, the expression level of six hub genes was decreased in phases of adenocarcinoma in human's colonic mucosa (all the p < 0.05; **Figures 4C,D**). Additionally, in comparison with normal samples, the expression level of CCR7, CXCL10, IDO1, and MMP9 were increased in phases of adenocarcinoma's tissue, while the expression level of CXCL9 and VCAM1 were decreased. These results are shown in **Figure 4** and **Supplementary Table S6**.

To extend our findings, the gene expression levels in CRC and para-cancerous tissues were compared based on TCGA database. Consequently, the gene differential expression level with regrading to hub genes was constructed; and all of the hub genes shown a significant difference in expression level between cancer and para-cancer tissues (all the p < 0.05; **Figure 5A**). Importantly, the Kaplan–Meier survival curves indicated that a higher expression level of CXCL10 (hazard ratio = 0.63; p = 0.035), CXCL9 (hazard ratio = 0.63; p = 0.037), MMP9 (Hazard Ratio = 0.61; p = 0.023), and CCR7 (Hazard Ratio = 0.59;

p = 0.013) were significantly associated with the poor prognosis for CRC patients; although the difference in overall survival (OS) between high and low expression of IDO1 (Hazard Ratio = 0.78; p = 0.27) and VCAM1 (Hazard Ratio = 0.74; p = 0.215) were not significant. Additionally, there was a clear tendency for lower expression of CXCL10 (Hazard Ratio = 0.38; p < 0.001), CXCL9

(Hazard Ratio = 0.11; p = 0.004), MMP9 (Hazard Ratio = 0.59; p = 0.046), IDO1 (Hazard Ratio = 0.53; p = 0.012) and CCR7 (Hazard Ratio = 0.61; p = 0.051) to be associated with a better prognosis in the recurrence free survival (RFS) time. This suggests that, to some extent, the effect of CXCL10, CXCL9, MMP9, IDO1, and CCR7 overexpression on early survival time resulted in a decrease in the survival rate. And these results are shown in **Figure 6**.

# Investigating the Functional Role and TF of Hub Genes

To further understand how the hub genes were correlated with UC-associated carcinogenesis, we applied DAVID and Metascape online database to explore the biological function. The results of GO term enrichment indicated that the GO:0032496∼response to lipopolysaccharide (Enriched genes: CCR7, CXCL9, IDO1, CXCL10; p = 1.87E-06), GO:0030816∼regulation of cAMP metabolic process (Enriched genes: CXCL9, CXCL10; p = 1.25E-03), and GO:0006954∼inflammatory/immune response (Enriched genes: CCR7, CXCL9, CXCL10; p = 2.11E-03) were mainly enriched, while pathway of Ecb04668: TNF signaling pathway (Enriched genes: CXCL10, VCAM1, MMP9; p = 1.36E-03), Ecb04062: Chemokine signaling pathway (enriched genes: CCR7, CXCL9, CXCL10; p = 3.76E-03), and Ecb04060: Cytokinecytokine receptor interaction (enriched genes: CCR7, CXCL9, CXCL10; p = 4.64E-03) (**Figure 5B**).

Finally, we predicted TFs and found that nuclear factor NFkappa-B1 (NFKB1) (NES = 14.21, target genes = 5, motifs = 30), friend leukemia integration 1 TF (FLI1) (NES = 7.57, target genes = 3, motifs = 5), and signal transducer and activator of transcription 1 (STAT1) (NES = 7.71, target genes = 3, motifs = 2) as the master regulators of the hub genes are involved in UCassociated carcinogenesis (**Figure 5C**).

# DISCUSSION

After adjusting the batch effects, 16,653 co-annotated genes among GPL6244 and GPL570 platform microarray datasets. Subsequently, we included co-annotated genes, some of which were present in low abundance or with small FC, in a further analysis, in which combination with WGCNA could integrate external traits and avoid information loss at a system level. According to the results, both MEbrown and MEgreen appeared to be moderately effective in revealing the UC-based global network. Biologically, following the functional enrichment analysis, the pathways of viral carcinogenesis, proteoglycans in cancer, platelet activation, ligand-receptor interaction, calcium signaling pathway, and cAMP signaling pathway were identified as being significantly associated with UC active. And cancer with highly correlated pathway and enriched genes were selected to construct the PPI network. The most critical genes were CXCL10, VCAM1, CXCL9, MMP9, IDO1, and CCR7, indicating that genetic variability influences susceptibility to the disease global network, and subsequently revealing potential regulatory roles in UC-associated carcinogenesis. Furthermore, these hub genes majorly enriched in tumor necrosis factor (TNF) signaling

pathway, chemokine signaling pathway, and cytokine-cytokine receptor interaction; and potentially regulated by NFKB1, FLI1, and STAT1 in TFs network analysis. Furthermore, a significant association of CCR7, CXCL10, CXCL9, IDO1, and MMP9 with UC-correlated CRC development was identified by integrating gene expression and survival analysis.

Emerging evidence has revealed the central role of geneenvironment interaction in UC-based disease networks. Extrinsic and intrinsic environmental factors may cause chronic or acute inflammation in UC patients. Wang et al. found that calcium signaling pathway contributes to the development of colonic dysmotility in UC and intestinal inflammation, may be related with the calcium-transporting ATPase dysregulation in epithelial cells (Wang et al., 2016). And, additionally, the evidence for the interdependence of platelet abnormalities in UC model and patients, suggesting

that the pathological state of changes in platelet parameters and their activation, may be linked to the inflammatory response and enhanced platelet-leukocyte, and aggregate formation associated with colitis (Senchenkova et al., 2015; Gawronska et al., 2017). Proteoglycans have been found to be critical in the regulation of stem cell through inducing precise and coordinated modulation of key growth factors, resulting in selective mitogen-activated protein kinases (MAPK) and/or another intracellular signaling, demonstrating an aberrant expression of ligand-receptor interaction on immune cells in IBD patients (Elshal et al., 2016; Gawronska et al., 2017). Li et al. demonstrated the multiple proinflammatory signaling pathways and candidate biomarkers, including STAT1, STAT6, and cAMP signaling pathway, in the exacerbation of UC (Li et al., 2012). Boothello et al. (2019) revealed that proteoglycans mediate cancer stem cells induced CRC xenograft's growth in a dose-dependent fashion. Moreover, syndecan-2, a type of proteoglycan, up-regulates MMP-7 expression in colon cancer cells via PKCγ-mediated activation of FAK/ERK signaling (Jang et al., 2017). Therefore, the pathway of the calcium signaling pathway, ligand-receptor interaction, platelet activation, cAMP signaling pathway and the none-cancer pathway involved, may provide insight into the immunological and inflammatory response, and the hypothesis of phospholipid-related barrier defects in the intestinal mucosa offers an opportunity to further understand UC-based pathogenesis.

The association between UC-related chronic inflammation and colon cancer has long been recognized (**Table 1**). According to a systemic review reported by Tatiya-Aphiradee, the pathway immuno-inflammatory response was closely linked to the regulation and maintenance of UC pathogenesis, that directly mediated by dynamic and complex communication between immune cells and cytokines (Tatiya-Aphiradee et al., 2018). Biologically, the dysregulation of antigen recognition, neutrophil chemotaxis, commensal microflora, and epithelial barrier defects may provide insight into the immunological and inflammatory response, it might offer an opportunity to further understand UC-based pathogenesis (Hindryckx et al., 2016). Gene expression profiling by Zhang's group shown the pathways includes PI3K-Akt signaling, cytokine-cytokine receptor interaction and ECM-receptor interaction was significantly associated with the process of colitis-associated carcinogenesis (Zhang et al., 2017). Several potential biomarkers of TNF signaling pathway, including TNF-α, IL-6, IL-1, TGF-β, and IL-10, have been confirmed to be involved in the process of malignant transformation of cells and carcinogenesis (Zhang et al., 2017). The pathway of cytokine-cytokine receptor interaction may also be closely linked to UC-related inflammation and tumorigenesis. Fang et al. (2015) compared IBD microarray datasets and found an important role for cytokine-cytokine receptor pathway dysregulation in both pediatric and mouse models of colitis. In sum, during the procession of intestinal inflammation and carcinogenesis, a variety of immunological and inflammatory signaling events, including the TNF signaling pathway, chemokine signaling pathway, and cytokine-cytokine receptor interaction, are activated and involved in a complex biological process.

Among the candidate biomarkers, the current understanding of the function of CXCL10 and CXCL9 may recruit the leukocytes to inflammation sites. However, a novel report from Shukla's group demonstrated that both CXCL10 and CXCL9 may promote colonic tumorigenesis via promotes the cytokine-mediated mucosal injury and inflammation response (Shukla et al., 2016). Additionally, IDO1 were over-expressed in inflamed and adenoma of the colon, also functioned in promotes colitis-associated tumorigenesis independent of T-cell immune surveillance (Thaker et al., 2013). MMP9 could maintain the microbiota and colonic epithelium mucosal barrier, also correlated with tissue remodeling and carcinogenesis via activates the EGFR signaling pathway (Pujada et al., 2017). The adhesion molecules VCAM-1 and ICAM-1, associated with macrophage infiltration, are directly associated with cell transmigration in inflamed colonic tissue (Wu et al., 2014). In addition, Bernhard et al. revealed that VCAM1 was correlated with different subsets of three immune cells and with high densities of T-cell subpopulations within specific tumor regions in CRC, thus the expression of adhesion molecules also associated with survival prognosis (Mlecnik et al., 2010). What's more, the lymphoid chemokine receptor CCR7 was re-expressed by activated T cells, allowing them to flow from the tissue to the lymph nodes through afferent lymphatics. McNamee's data showed a critical role for CCR7 in orchestrating immune cell traffic (McNamee et al., 2015). The role of chemokines in tumor angiogenesis was achieved in a CCR7-dependent manner through inhibiting Met/ERK/Elk-1/HIF-1α/VEGF-A pathway in CRC (Xu et al., 2018).

Finally, the TFs analysis results shown that NFKB1, FLI1, and STAT1 were significantly predicted in hub gene's regulatory network, correlated with UC-correlated tumorigenesis. Here, STAT1, the first member of signal transducer and activator of transcription (STAT) family, has been involved in cancer suppression, including CRC (Zamanian-Azodi and Rezaei-Tavirani, 2019). Schwiebs et al. (2019) found that STAT1 has involved in the process of tumor immune microenvironment during the crosstalk of "inflammation-to-tumor." NF-kappa-B1 (NF-κB1) signaling is a prominent and widely studied inflammatory signaling cascade in the field of immunology (Eden et al., 2017). Increased transcription of NF-κB is associated with inflammation and angiogenesis. Burkitt proposed that NF-κB1 differentially regulate susceptibility to colitis-associated adenoma development (Burkitt et al., 2015). FLI1, a member of the family of ETS TFs, contains a highly conserved domain that recognizes ETS core consensus sites (GGAA/T; Hollenhorst et al., 2011; Tang et al., 2015). EWS-FLI1 regulates multiple target genes through binding to typical ETS core consensus sites or GGAA microsatellites, then participates in the carcinogenic process (Lessnick and Ladanyi, 2012; Tang et al., 2015). Azuara et al. (2018) have defined FLI1 as a DNA methylation signature that can be distinguished in the early detection of CRC associated with IBD.

In summary, we found that the pathways of platelet activation, ligand-receptor interaction, calcium signaling pathway, and cAMP signaling pathway may play an important role in UC development via multiple physiological and pathophysiological

processes, revealing a potentially attractive therapeutic target for UC-based disease networks. The overlapping results for CXCL10, VCAM1, CXCL9, MMP9, IDO1, and CCR7 were obtained, which are considered to be hub biomarkers involved in UCcorrelated tumorigenesis. Following the expression validation, survival analysis, and functional analysis, our results indicated that the novel biomarkers of CXCL10, VCAM1, CXCL9, MMP9, IDO1, and CCR7 has powerful statistical efficiency and biological function. These genes are also linked to immune dysregulation and inflammation response, and thus provide new insights into the pathogenetic mechanisms of UC development and tumorigenic processes. Finally, our results also subsequently identified that the master regulators of NFKB1, FLI1, and STAT1 have significantly associated with UC activity and carcinogenesis via target the candidate biomarkers.

# AUTHOR CONTRIBUTIONS

WS, RZ, and MY took the responsibility for all aspects of the reliability and freedom from bias of the data presented

### REFERENCES

and their discussed interpretation, and drafting the article. LM, JR, ZL, and JW took the responsibility for statistical analyses and interpretation of data. RL took the responsibility for fulltext evaluation and guidance, and final approval of the version to be submitted.

## FUNDING

This work was supported by National Key Clinical Discipline-National Natural Science Foundation of China (Grant No. 81372142), National Basic Research Program of China (Grant No. 2015CB554001), and Zhuhai science and technology project (Grant No. 20171009E030062).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2019.00662/full#supplementary-material

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

Copyright © 2019 Shi, Zou, Yang, Mai, Ren, Wen, Liu and Lai. 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.

fphys-10-00662 May 29, 2019 Time: 19:42 # 14

# Non-invasive Imaging and Modeling of Liver Regeneration After Partial Hepatectomy

*Sara Zafarnia1† , Anna Mrugalla1† , Anne Rix1 , Dennis Doleschel1 , Felix Gremse1 , Stephanie D. Wolf 2 , Johannes F. Buyel3,4‡ , Ute Albrecht 2 , Johannes G. Bode2 , Fabian Kiessling1 and Wiltrud Lederle1 \**

*1 Institute for Experimental Molecular Imaging, Medical Faculty, RWTH Aachen University, Aachen, Germany, 2 Department of Gastroenterology, Hepatology and Infectious Diseases, Medical Faculty, Heinrich-Heine-University Düsseldorf, Düsseldorf, Germany, 3 Fraunhofer Institute for Molecular Biology and Applied Ecology IME, Aachen, Germany, 4 Institute for Molecular Biotechnology, RWTH Aachen University, Aachen, Germany*

### *Edited by:*

*Andreas Teufel, University of Heidelberg, Germany*

### *Reviewed by:*

*Jun Li, University Medical Center Hamburg-Eppendorf, Germany Martin Meier, Hannover Medical School, Germany*

> *\*Correspondence: Wiltrud Lederle<EMAIL_ADDRESS>
*† These authors have contributed equally to this work*

*‡ orcid.org/0000-0003-2361-143X*

### *Specialty section:*

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

*Received: 30 March 2018 Accepted: 01 July 2019 Published: 17 July 2019*

### *Citation:*

*Zafarnia S, Mrugalla A, Rix A, Doleschel D, Gremse F, Wolf SD, Buyel JF, Albrecht U, Bode JG, Kiessling F and Lederle W (2019) Non-invasive Imaging and Modeling of Liver Regeneration After Partial Hepatectomy. Front. Physiol. 10:904. doi: 10.3389/fphys.2019.00904*

The liver has a unique regenerative capability upon injury or partial resection. The regeneration process comprises a complex interplay between parenchymal and non-parenchymal cells and is tightly regulated at different scales. Thus, we investigated liver regeneration using multi-scale methods by combining non-invasive imaging with immunohistochemical analyses. In this context, non-invasive imaging can provide quantitative data of processes involved in liver regeneration at organ and body scale. We quantitatively measured liver volume recovery after 70% partial hepatectomy (PHx) by micro computed tomography (μCT) and investigated changes in the density of CD68+ macrophages by fluorescence-mediated tomography (FMT) combined with μCT using a newly developed near-infrared fluorescent probe. In addition, angiogenesis and tissueresident macrophages were analyzed by immunohistochemistry. Based on the results, a model describing liver regeneration and the interactions between different cell types was established. *In vivo* analysis of liver volume regeneration over 21 days after PHx by μCT imaging demonstrated that the liver volume rapidly increased after PHx reaching a maximum at day 14 and normalizing until day 21. An increase in CD68+ macrophage density in the liver was detected from day 4 to day 8 by combined FMT-μCT imaging, followed by a decline towards control levels between day 14 and day 21. Immunohistochemistry revealed the highest angiogenic activity at day 4 after PHx that continuously declined thereafter, whereas the density of tissue-resident CD169+ macrophages was not altered. The simulated time courses for volume recovery, angiogenesis and macrophage density reflect the experimental data describing liver regeneration after PHx at organ and tissue scale. In this context, our study highlights the importance of non-invasive imaging for acquiring quantitative organ scale data that enable modeling of liver regeneration.

Keywords: non-invasive imaging, modeling, liver regeneration, partial hepatectomy, macrophages, angiogenesis, FMT-μCT

# INTRODUCTION

The liver is known for its high regenerative potential being able to restore up to 70% of its mass after injury or partial resection (Minuk, 2003). The regeneration process constitutes a complex interplay of various cell types and signaling pathways (Taub, 2004). Depending on the circumstances, two different modes of regeneration are known to be activated. In case of an impaired hepatocyte proliferation, as for instance following severe or chronic liver injury, liver stem cells (also known as oval cells in rodents) become activated as a mechanism of the liver to regenerate and recover its function (Itoh and Miyajima, 2014). In contrast, after partial resection or moderate liver damage, complete liver regeneration is achieved by proliferation of the remaining parenchymal and non-parenchymal cells. Hepatocytes are the first cells to grow and proliferate after partial resection followed by Kupffer cells, biliary epithelial cells, and stellate cells. The process is accompanied by the induction of angiogenesis that is also crucially involved in liver regeneration (Drixler et al., 2002; Uda et al., 2013). Besides other cell types, macrophages have been shown to play a stimulatory role in liver regeneration by producing molecular factors that are pivotal in the regeneration process (Abshagen et al., 2007; Li and Hua, 2017). Depletion of resident macrophages (Kupffer cells) as well as an impaired macrophage recruitment from the periphery and bone marrow results in a delayed regeneration demonstrating the stimulatory function of both resident and infiltrating macrophages (Takeishi et al., 1999; Abshagen et al., 2007; Melgar-Lesmes and Edelman, 2015; Nishiyama et al., 2015). The most frequently used model to study liver regeneration is the model of partial hepatectomy (PHx) described first by Higgins and Anderson (1931) in rats. Mitchell and Willenbring recently developed a modified protocol of a standardized surgical technique for PHx in mice (Mitchell and Willenbring, 2008). In both models, approximately two-thirds of the liver are surgically removed. The regeneration process starts immediately leading to full recovery of liver mass within 7–10 days (Taub, 2004; Nishiyama et al., 2015). The advantage of a surgical model in comparison with toxic injury models is the fact that the regeneration process after PHx is not associated with massive necrosis and necrosis-induced acute inflammation so that all changes observed after PHx can be ascribed to the physiological regeneration process (Michalopoulos, 2010). In a clinical context, this physiological regeneration process becomes important in patients who underwent partial liver resection, in donors and recipients following living-donor liver transplantation and in patients with acute liver failure.

Liver regeneration is a complex process involving the interactions of different cell types on various levels. Thus, we followed an approach using multi-scale methods (Castiglione et al., 2014) including non-invasive imaging and histological analyses in order to investigate liver regeneration after pHx. Non-invasive imaging is a useful tool since it enables a longitudinal and quantitative assessment of morphological, functional, and molecular parameters at the organ and whole body level. Liver volume recovery was measured *via* micro computed tomography (μCT), and the density of CD68+ macrophages was determined by combined fluorescence-mediated tomography and μCT (FMT-μCT) using a newly developed near-infrared fluorescent (NIRF) probe. The *in vivo* results were validated by immunohistochemical analyses of CD68+ and F4/80+ macrophages. At the tissue level, the contribution of tissueresident macrophages and angiogenesis was investigated by additional immunohistochemical analyses. Based on the experimental data, a simple model describing liver regeneration and the interrelation between volume recovery, macrophages, and angiogenesis was generated.

# MATERIALS AND METHODS

### Generation and Purification of the Near-Infrared Fluorescent CD68 Probe

The NIRF probe targeting CD68+ macrophages was generated by coupling an amine-reactive NIR fluorochrome (NHS ester), VivoTag 680 (excitation peak 665 ± 5 nm, emission peak 688 ± 5 nm) (Perkin-Elmer), to a rat anti-mouse CD68 antibody (AbDSerotec) according to manufacturer's instructions. In brief, VivoTag 680 was dissolved in dry dimethyl sulfoxide (Sigma-Aldrich) at a concentration of 10 mg/ml. Prior to the labeling, the buffer of the antibody was exchanged by dialysis into conjugation buffer (50 mM carbonate/bicarbonate buffer, pH 8.5) using Slide-A-Lyzer dialysis cassettes (AbD Serotec) according to the protocol provided by the manufacturer. After buffer exchange, 30 μl of VivoTag 680 was added to the rat antimouse CD68 antibody. Following 1 h of incubation at room temperature protected from the light, the NIRF CD68 probe (approximately 151 kDa) was separated from free fluorescent dye (approximately 1 kDa) and antibody oligomers (larger than 300 kDa) by fast protein liquid chromatography using a Superdex 200 resin in a pre-packed 10/300 GL column (GE Healthcare). Probe concentration was determined using a BCA Protein Assay Kit (Uptima) according to the protocol provided by the manufacturer.

### *In vitro* Binding of the Near-Infrared Fluorescent CD68 Probe

Binding specificity of the NIRF CD68 probe was tested *in vitro* by incubating the macrophage cell line J774A.1 (CLS Cell Lines Service) with the NIRF CD68 probe (10 nM, 2 or 4 h, 37°C). For competitive binding analyses, J774A.1 cells were incubated with 10 nM of the NIRF CD68 probe and a 10-fold molar excess of unlabeled CD68 antibody (100 nM, 2 or 4 h, 37°C). Fluorescent microphotographs were acquired with the Axio Imager M2 (Zeiss) and a high-resolution camera (AxioCamMRm Rev.3; Zeiss) using a Cy5.5 filter and a fixed exposure time. The signal intensities were determined using the software ImageJ 1.47v (W. Rasband, National Institutes of Health). For quantification, the mean fluorescent signal intensity at 695 nm was determined by analyzing five microscopic images per well (*n* = 3 wells per culture condition).

### Animal Studies

All animal experiments were performed according to German legal requirements and animal protection laws and were approved by the Authority for Environment Conservation and Consumer Protection of the State of North Rhine-Westphalia (LANUV).

### Biodistribution and *in vivo* Specificity of the Near-Infrared Fluorescent CD68 Probe

Biodistribution and *in vivo* specificity of the NIRF CD68 probe were analyzed in male C57BL/6 J mice (Charles River) (*n* = 5 per group). The mice were fed a chlorophyll-free diet (ssniff Spezialdiäten GmbH) 7 days before imaging, and the scanning area was depilated prior to the scans. For macrophage imaging, 2.6 μg of the NIRF CD68 probe dissolved in 0.9% w/v NaCl was injected intravenously (i.v.). To analyze the biodistribution, animals were scanned longitudinally immediately before and 1, 3, 6, 12, 24, and 48 h after probe injection by FMT-μCT. *In vivo* specificity of the probe was examined by competitive binding analysis injecting a 5-fold mass excess of unlabeled CD68 antibody (13 μg) intraperitoneally (i.p.) 1 h before i.v. injection of the NIRF CD68 probe (2.6 μg). Directly after the last FMT-μCT measurement, the mice were sacrificed, and the liver was resected and cryoconserved in Tissue-Tek (Sakura) for immunohistochemical analyses.

### *In vivo* Determination of Macrophage Density and Volume Recovery During Liver Regeneration After Partial Hepatectomy

Liver regeneration was analyzed in male C57BL/6 J mice (Charles River) after PHx and sham surgery. One hour before surgery, mice were treated with carprofen [subcutan (s.c.) 5 mg/kg] (Pfizer Animal Health SA). PHx and sham surgeries were performed under sterile conditions as previously described by Mitchell and Willenbring (2008). In brief, mice were anesthetized with isoflurane (2% v/v isoflurane in oxygenenriched air) and positioned on a temperature-controlled pad to regulate body temperature. For PHx, after midline laparotomy, the left lateral lobe was ligated as close to the base of the lobe as possible. The secondary knot was placed above the gall bladder of the median lobe but not closer than 2 mm from the suprahepatic vena cava. The ligated liver lobes were surgically resected. At the end of the surgery, the abdomen was rinsed with saline solution, and the abdominal wall and the skin were sutured separately. For sham surgery, a midline laparotomy was performed with gentle palpation and manipulation of the liver without resection of the liver lobes. Directly after surgery, mice received a s.c. injection of 10 mg/kg enrofloxacin (Bayer). Afterwards, analgesia was continued by s.c. injection of 5 mg/kg carprofen (Pfizer Animal Health SA) once per day for 3 days.

CD68+ macrophages and volume recovery were monitored by FMT-μCT and contrast-enhanced μCT at different time points (4, 8, 14, and 21 days) after PHx and sham surgery. As an additional control, untreated mice were measured. Twelve hours before each FMT-μCT measurement, mice were injected i.v. with 5.7 μg of the NIRF CD68 probe diluted in 0.9% w/v NaCl. In addition, 45 min before the FMT-μCT scans, the mice received an i.v. injection (150 μl) of the contrast agent Imeron 400 (Bracco Imaging) to enable a better segmentation of the liver. Directly after each FMT-μCT measurement, mice were sacrificed, and the liver was resected and cryoconserved in Tissue-Tek (Sakura) for immunohistochemical analyses. The group size was as follows: for assessment of macrophage density: untreated mice: *n* = 10; mice after PHx: *n* = 5 for day 4, *n* = 4 for day 8, *n* = 4 for day 14, *n* = 4 for day 21; shamoperated mice: *n* = 3 for day 4, *n* = 3 for day 8, *n* = 3 for day 14, *n* = 6 for day 21. For volume recovery analysis: *n* = 5 for each time point.

### Imaging Protocols

Three-dimensional (3D) FMT-μCT scans were conducted as described by Kunjachan et al. (2013). For the measurements, mice were anesthetized and positioned in a μCT- and FMT-compatible mouse bed. For anatomical information, mice were scanned in a dual-energy μCT system (TomoScope 30s Duo, CT Imaging GmbH). For biodistribution and *in vivo* competitive binding, analyses scans were performed using the scan protocol SQD-6565-360-29, which acquires 720 projections with 516 × 506 pixels requiring a scanning time of 29 s per subscan. For the assessment of macrophage density and liver volume, the HQD-6565-360-90 protocol was applied, which acquires 720 projections with 1,032 × 1,012 pixels requiring a scanning time of 90 s per subscan. Directly after acquiring the μCT scans, the mouse bed was transferred to the FMT system (FMT 2500LX, PerkinElmer), and FMT scans were performed at 680 nm using 120 excitation positions (3 mm distance). Data fusion and reconstruction of the fluorescence distribution were performed as described (Gremse et al., 2014). Based on the μCT data, organs were manually segmented, and the liver volume and probe concentration in the segmented organs were determined using the Imalytics Preclinical software (ExMI/Gremse-IT) (Gremse et al., 2016). To assess the variability of the organ segmentations, several organs (liver, lung, heart, kidney, and bladder) were segmented in five representative scans by two different persons. The organ volumes of the two analyses correlated strongly (*R*2 = 0.996, *p* < 0.05) between the users. DICE scores, which describe similarity between segmentations, were also high (0.92 ± 0.045), showing good reproducibility of the manual organ segmentation. To further confirm the accuracy of manual organ segmentation, we performed an additional correlation analysis between segmented organ volumes and weights of excised organs. Organ weights (heart, liver, kidneys, spleen, and tumors) and contrast-enhanced μCT scans (*n* = 6) were available from a previous study (Rosenhain et al., 2018). The analysis resulted in a strong correlation between segmented organ volumes from contrast-enhanced μCT scans and the organ weights (*R*2 = 0.976). The slope of the regression line was 0.84, i.e., below 1, which can be explained by the loss of blood during the organ excision and harvesting procedure (thus lower values for the organ weights as compared to the segmented volumes).

### Antibodies

The following primary antibodies were used to stain macrophages/ Kupffer cells: rat anti-mouse CD68 antibody (AbDSerotec), rat anti-mouse F4/80 antibody (AbDSerotec), and rat anti-mouse CD169 antibody (AbDSerotec). To stain endothelial cells, rat anti-mouse CD31 antibody (BD Biosciences) was applied. Goat anti-mouse VEGFR2 antibody (R and D Systems) was used to determine the VEGFR2 density. Secondary IgG antibodies (donkey anti-rat Alexa Fluor 488, donkey anti-rat Cy-3 and donkey anti-goat Cy-3) were obtained from Dianova. Cell nuclei were counterstained with 4′,6-diamidino-2-phenylindole (DAPI; Merck KGaA).

### Indirect Immunohistochemistry

For immunohistochemical analysis, frozen organs were cut into 8-μm slices. Fixation and staining of the cryosections were performed as previously described (Lederle et al., 2010). Per section, five to seven fluorescent microphotographs were acquired with the Axio Imager M2 (Zeiss) and a high-resolution camera (AxioCamMRm Rev.3; Zeiss). The number of CD68+ , F4/80+ , CD169+ and DAPI+ cells per microphotograph was counted manually using the ImageJ 1.47v software (W. Rasband, National Institutes of Health), and the percentage of CD68+ , F4/80+ and CD169+ cells per DAPI+ cells was calculated, respectively. Quantitative analysis of microvessel density and angiogenic activity was done using the AxioVisionRel 4.8 software (Zeiss). The microvessel density was determined by quantifying the CD31+ area fraction, and the angiogenic activity was assessed by determining the proportion of the VEGFR2+ area fraction to the CD31+ area fraction. The group size was as follows: untreated mice: *n* = 10 (CD169: *n* = 5); mice after PHx: *n* = 5 for day 4, *n* = 4 for day 8, *n* = 4 for day 14, *n* = 4 for day 21; sham-operated mice: *n* = 4 for day 4, *n* = 6 (CD169: *n* = 3) for day 8, *n* = 3 for day 14, *n* = 7 for day 21.

# Numeric Modeling

Stimulation and growth of liver cell compartments (hepatocytes, Kupffer cells, macrophages, and endothelial cells) were described by differential equations in a simplified model. Based on data from literature, relative proportions of these cell types in healthy livers were assumed to be 0.8, 0.06, 0.06, and 0.08, respectively (Vekemans and Braet, 2005). To simulate PHx, these values were multiplied by 0.3 at the initial state of simulation to account for the reduced total cell amount. The first derivative of each cell compartment fraction was assumed to depend linearly on the amount, excess, or lack of other cell compartments, with six coefficients describing the strength of the corresponding effect. The parameters *k*HE and *k*ME describe the stimulation and support of vessel growth by a hepatocyte lack and macrophage excess, respectively. *k*HM describes attraction of macrophages by a lack of hepatocytes. *k*HK describes stimulation of Kupffer cell growth by a lack of hepatocytes. *k*EH describes the connection of blood vessel and hepatocyte growth, because the latter is limited by nutrition supply and structural alignment requirements. A parameter *k*homeostasis describes other effects not covered by our simplified model, and it drives the four cell compartments slowly toward the homeostatic situation. Therefore, the derivatives of the four compartments are the weighted sums of the affecting coefficients and expressions for relative or absolute lack as described in the following:

$$dH/dt = k\_{\text{homecostasis}} \left(1 - H/0.8\right) + k\_{\text{HE}} \left(0.8/0.08 - H/E\right)$$

$$dK/dt = k\_{\text{homecostasis}} \left(1 - K/0.06\right) + k\_{\text{HK}} \left(0.8 - H\right)$$

$$dM/dt = k\_{\text{homecostasis}} \left(1 - M/0.06\right) + k\_{\text{HM}} \left(0.8 - H\right)$$

$$dE \wedge dt = k\_{\text{homecostasis}} \left(1 - E/0.08\right) + k\_{\text{HE}} \left(0.8 - H\right)$$

$$+ k\_{\text{ME}} \left(M/\left(H + K + M + E\right) - 0.06\right)$$

The simulation was performed using fourth-order Runge-Kutta integration with the time interval of 60 min over a period of 100 days, resulting in four-cell compartment curves. From these, simulated measurement curves were computed, i.e., volume (sum of all compartments), total macrophages (sum of macrophages and Kupffer cells), and Kupffer cells and angiogenesis (derivative of endothelial growth). The six parameters were iteratively adjusted until simulated and measured curves matched.

# Statistical Analysis

Statistical analysis was performed using Prism 5.0 (GraphPad software). Results are shown as mean ± standard deviation. All statistical analyses were performed using one-way analysis of variance (ANOVA) followed by Bonferroni correction for multiple comparisons (\**p* < 0.05, \*\**p* < 0.01, \*\*\**p* < 0.001).

# RESULTS

# *In vitro* Binding Specificity of the Near-Infrared Fluorescent CD68 Probe

For non-invasive imaging of macrophages by FMT-μCT, a NIRF probe targeting CD68+ macrophages was generated and evaluated *in vitro* and *in vivo*. The binding specificity of the NIRF CD68 probe was evaluated *in vitro* by competitive binding analysis using the macrophage cell line J774A.1. After 2 and 4 h of incubation with the NIRF CD68 probe alone, a strong fluorescent signal in the cells was observed (**Figure 1A**). Incubation of the cells with the NIRF CD68 probe together with a 10-fold excess of unlabeled anti-CD68 antibody resulted in a strongly reduced signal at both time points. Quantification of the fluorescent images revealed an increase in signal intensity from 2 to 4 h and confirmed a significantly lower mean signal intensity in the cells of the competitive binding group as compared to the cells incubated with the NIRF CD68 probe alone (*p* < 0.001) (**Figure 1B**).

### Biodistribution and *in vivo* Specificity of the Near-Infrared Fluorescent CD68 Probe

The biodistribution of the NIRF CD68 probe was analyzed longitudinally in healthy mice by FMT-μCT 1, 3, 6, 12, 24,

and 48 h after injection. Quantitative analysis of the NIRF CD68 probe accumulation in the liver, lung, and kidneys revealed a significantly higher mean concentration in the liver as compared to the kidneys and the lung at all measuring time points (*p* < 0.001 for all time points, respectively) (**Figure 2A**). The concentration in the kidneys and the lung was constantly low without significant changes over time. In the liver, the mean concentration increased after probe injection reaching a maximum concentration at 12 h after injection. Thereafter, the mean probe concentration in the liver declined to a value similar to that observed at 1 h after injection (**Figures 2A,C**; representative 3D rendering of reconstructed FMT-μCT data shown in **Figure 2B**).

The *in vivo* specificity of the NIRF CD68 probe was analyzed in a competitive binding experiment in which a 5-fold excess of unlabeled CD68 antibody was injected 1 h prior to injection of the NIRF CD68 probe (competitive binding group). In the liver, starting from 3 h after probe injection, a lower mean probe concentration was measured in the competitive binding group at each time point compared to the concentration measured after injection of the NIRF CD68 probe alone (control group) (**Figure 3A**; representative transversal FMT-μCT fusion images of competitive binding and control mice are shown in **Figure 3B**). The difference in the mean concentration of NIRF CD68 probe in the liver was significant 12 h after injection.

# Volumetric Recovery of the Liver After Partial Hepatectomy

To investigate liver regrowth after PHx by non-invasive imaging, the volume of the liver was measured *via* contrast-enhanced μCT. Quantitative analysis revealed that >70% of the total liver volume was reached at day 4 and > 80% was regained at day 8 after PHx (**Figure 4A**; representative 3D rendered CT images with segmented organs of an untreated control mouse and mice after PHx are shown in **Figure 4B**). The mean value determined at day 14 was slightly above the mean volume of the liver of untreated control mice. The mean volume at day 21 after PHx was comparable to control values.

### *In vivo* Monitoring of CD68+ Macrophage Density During Liver Regeneration

Macrophages play an important stimulatory role during liver regeneration (Takeishi et al., 1999; Abshagen et al., 2007; Nishiyama et al., 2015). Thus, we investigated the time course of macrophage density after PHx by non-invasive FMT-μCT imaging. CD68+ macrophages were monitored at day 4, 8, 14, and 21 after PHx and sham surgery using the NIRF CD68 probe. As an additional control, untreated mice were measured. Quantitative analysis of probe accumulation revealed a higher mean NIRF CD68 probe concentration in the liver at day 8, 14, and 21 after PHx compared to untreated and sham-operated control mice (**Figure 5A**). The mean concentration was highest at day 8 and 14 after PHx followed by a decline between day 14 and 21 indicating a transient increase in the density of CD68+ macrophages. No major changes in the mean NIRF CD68 probe concentration were measured in the liver of shamoperated mice during the whole observation period. Representative frontal plane FMT-μCT fusion images of an untreated control mouse and mice after PHx are shown in **Figure 5B**.

### Immunohistochemical Analysis of Macrophage Subpopulations During Liver Regeneration

To validate the *in vivo* results, the density of CD68+ macrophages in the liver of untreated and sham-operated control mice and mice after PHx was determined by immunohistochemical analyses. Quantification of the density of CD68+ macrophages confirmed the trend of the *in vivo* findings showing a significant increase in the mean values after PHx until day 8 followed by a decline to levels observed in untreated control animals at day 21 (**Figure 6A**, *p* < 0.001). The mean density in the liver of sham-operated animals remained similar to that of untreated control animals without significant changes over time.

For further validation, we performed an immunohistochemical analysis of cells expressing F4/80, a generic macrophage marker that is independent of CD68 expression. Quantitative analysis of the density of F4/80+ macrophages showed a similar trend over

time as compared to CD68+ macrophages with a significant transient increase and highest mean values at day 8 after PHx (**Figure 6B**, *p* < 0.001). Again, the mean density of F4/80+ in the liver of sham-operated animals was comparable to that of untreated control animals and did not significantly change over time.

An increased density of macrophages during liver regeneration can be a result of tissue-resident Kupffer cell proliferation or caused by an infiltration of macrophages from the blood circulation. Both Kupffer cells and infiltrating macrophages have been shown to play an important role during liver regeneration (Takeishi et al., 1999; Abshagen et al., 2007; Melgar-Lesmes and Edelman, 2015; Nishiyama et al., 2015). To investigate the contribution of tissue-resident macrophages to the increased macrophage density, we analyzed the expression of CD169 by immunohistochemical analysis. In contrast to CD68+ and F4/80+ macrophages, no significant differences in the density of CD169+ macrophages were observed in the liver of untreated and shamoperated control mice and mice after PHx (**Figure 6C**).

# Immunohistochemical Analysis of Angiogenesis During Liver Regeneration

Angiogenesis is a crucial process involved in liver regeneration, and macrophages have been shown to stimulate endothelial cell activation and to regulate vessel growth (Melgar-Lesmes and Edelman, 2015). To investigate the interrelation between macrophage density and angiogenesis, we analyzed the microvessel density and angiogenic activity in the liver of untreated control mice and in mice after PHx and sham surgery using immunohistochemistry.

Microvessel density was not significantly different between mice after PHx, untreated and sham-operated control mice (**Figure 7A**). However, the angiogenic activity, as assessed by the proportion of the VEGFR2+ area fraction to the CD31+ area fraction, was markedly increased on day 4 after PHx followed by a continuous decrease until day 21 (**Figure 7B**). In sham-operated and untreated mice, no significant changes were found over time.

# Modeling of Liver Regeneration After Partial Hepatectomy

The *in vivo* and immunohistochemical results revealed different time courses of volume recovery, macrophage density and endothelial cell activation (angiogenesis). A simplified mathematical model was developed including different cell compartments of the liver (hepatocytes, Kupffer cells, recruited macrophages, and endothelial cells), and the growth and interplay of these compartments after PHx was simulated (**Figures 8A,B**). Based on the resulting cell compartment curves (**Figure 8B**), simulated measurement curves were computed describing liver volume (sum of all compartments), total macrophages (sum of recruited macrophages and Kupffer cells), and Kupffer cells and angiogenesis (derivative of endothelial growth) (**Figure 8C**).

After parameter adjustment, the time courses of the numerical model matched the experimental data of angiogenesis, macrophages, and liver volume obtained by non-invasive imaging and immunohistochemical analyses (**Figures 8C,D**). In detail, there was an early onset of angiogenesis, which was followed by an increase in the overall macrophage density peaking on day 8. The liver volume increased rapidly after PHx reaching levels above healthy liver on day 14 before normalization.

# DISCUSSION

Liver regeneration after injury or partial resection comprises a complex interplay of different cell types and is tightly regulated at various scales (Taub, 2004; Michalopoulos, 2010; Li and Hua, 2017). To address alterations during liver regeneration after PHx at different levels, we used non-invasive imaging in combination with immunohistochemistry and developed a simple mathematical model describing the interrelations between different cell types involved in liver regeneration, angiogenesis and liver volume recovery.

Liver volume recovery was non-invasively monitored by μCT imaging. The measurements revealed that >70% of the total liver volume were regained within 4 days. However, normalization of the liver volume was not reached until day 21 due to an increased volume observed on day 14. The increased volume can be explained by edema formation that sometimes occurs during liver regeneration (Pleskovic et al., 1996).

Since macrophages are known to play an important stimulatory role during liver regeneration, we investigated the density of different macrophage populations after PHx by non-invasive imaging and immunohistochemistry. For non-invasive imaging, we used combined FMT-μCT imaging and generated a NIRF probe targeting CD68. Specific binding of the NIRF probe to CD68+ macrophages was confirmed by competition analyses *in vitro* and *in vivo*. Quantitative *in vivo* FMT-μCT imaging after PHx and sham surgery revealed an increased mean concentration of the CD68 probe in the liver at day 8 and 14 after PHx indicating a transient increase in the density of CD68+ macrophages. Immunohistochemical analyses of CD68+ and F4/80+ macrophages showed the same trend with a significantly higher macrophage density on day 8 after PHx as compared to untreated and sham-operated control mice. However, the immunohistochemical data were more distinct than the results obtained by FMT-μCT imaging. The difference in precision can be explained by the lower spatial resolution of the FMT (about 2 mm). Small changes in the density of macrophages are more difficult to determine by FMT-μCT

than *via* immunofluorescence microscopy. Furthermore, although absorption and scattering of the photons are taken into account in the fluorescence reconstruction algorithms, numerical limitations of the complex diffuse optical behavior can still affect the accuracy of the data, especially in deeper lying organs and in organs with a high blood volume, such as the liver, since blood is the main near-infrared absorber *in vivo* (Gremse et al., 2014). Moreover, the blood pool of the circulating probe and the unspecific hepatic uptake of foreign substances such as probes and contrast agents can lead to an unspecific probe signal. Antibodies are known to have a long blood half-life (Freise and Wu, 2015). Since we used an antibody as targeting molecule, we cannot exclude that the blood pool of the imaging probe has influenced the measurements. To reduce the impact of the blood pool, a different targeting molecule could be chosen, e.g., a nanobody, that has a shorter blood half-life.

To analyze the contribution of tissue-resident Kupffer cells in liver regeneration after PHx, we investigated the expression of CD169 by immunohistochemical analysis. Quantification revealed no significant differences in the density of CD169+ macrophages in the liver of untreated and sham-operated control mice and mice after PHx. Thus, while the overall density of macrophages in the liver increased significantly after PHx, the density of tissue-resident CD169+ Kupffer cells did not. This is in line with previous findings showing that the number of Kupffer cells correlates with liver restoration rate (Melgar-Lesmes and Edelman, 2015). Therefore, the results provide further evidence for the involvement of macrophages recruited from the blood circulation in liver regeneration after PHx (Melgar-Lesmes and Edelman, 2015; Nishiyama et al., 2015; Wen et al., 2015). However, further investigation is needed to unravel the details of resident Kupffer cell and infiltrating macrophage contribution to liver regeneration.

Angiogenesis is an important process involved in liver regeneration after hepatectomy and mutual interactions have been described between hepatocytes, Kupffer cells/macrophages and endothelial cells during liver regeneration (Drixler et al., 2002; Uda et al., 2013; Castiglione et al., 2014). Therefore, we analyzed the angiogenic endothelial cell activity and microvessel density in the liver after PHx by immunohistochemistry. A markedly increased mean angiogenic activity was detected at day 4 after PHx that decreased steadily until day 21 which is in accordance with previously published data (Alizai et al., 2017). Quantitative analysis of the microvessel density showed no significant differences between mice after PHx and control mice. This finding is not in line with results published by other groups that showed an increase in the microvessel density following PHx (Drixler et al., 2002, 2003). The discrepancy can be explained by different quantification methods. While we included larger arterioles and venules in the quantification, they were excluded by the other groups. Larger vessels contribute more to the overall CD31+ area fraction than very small vessels. Thus, changes in the density of these small vessels do not have a major effect on the overall vessel density.

(C) Quantitative analysis of the density of CD169+ mice after PHx. UT, untreated mice.

showed a markedly increased mean activity at day 4 after PHx that decreased steadily until day 21. UT, untreated mice.

The *in vivo* and immunohistochemical results revealed differences in the time courses for volume recovery, macrophage density and angiogenic endothelial cell activity, nevertheless, all three time courses show a progression towards levels of healthy situation over time. To describe the interrelation between volume recovery, macrophage density and angiogenesis occurring at different scales, we developed a numerical model that describes the growth and interplay of the involved liver cell compartments (hepatocytes, Kupffer cells, recruited macrophages, and endothelial cells). A numerical model serves to bridge the gap between hidden parameters (e.g., *k*HE) and observable measurements. The model may contain many direct interactions between cell types which are simple by themselves but result in a complicated situation altogether, which cannot be described by closed formulas but require numerical approaches instead. Our simulated measurements generally reflect the experimental data obtained by non-invasive imaging and immunohistochemical analyses. Differences remain in the earlier increase of the macrophage and Kupffer cell populations in our model as compared to the measurements. The maximum in angiogenic activity precedes the peak of macrophage density and normalization of liver volume. This shows that the model can CD169+

macrophage density,

cells, and angiogenesis showed similar time courses as the experimental data. (D) Diagram showing the time courses for liver volume, CD68+

macrophages, and angiogenic activity of the experimental data acquired in mice after PHx.

describe liver regeneration at organ and tissue scale, and that the model substantially benefits from experimental quantitative non-invasive imaging data. Nevertheless, higher sample numbers would improve the stability and reliability of the model. At tissue scale, different mathematical models with considerable higher complexity than our model have been established that describe and predict specific important processes involved in liver regeneration. Recently, a mathematical model revealed a crucial role of hybrids consisting of hepatocytes and bone marrow cells that trigger proliferation in the regeneration process (Pedone et al., 2017). For liver regeneration after CCl4 damage, Hoehme et al. have established a model that describes structural alignments of hepatocytes to sinusoids as a crucial pre-requisite for regaining the complex microarchitecture (Hoehme et al., 2010). Simple algorithmic models like the one proposed here have the advantage of high robustness and thus the suitability for integrating less quantitative *in vivo* data but also face several limitations. The model does not comprise the full complexity of the interrelations between the hepatocytes, macrophages, and endothelial cells that occur during liver regeneration. In addition, it does not describe causal relationships between the involved cell compartments in a detailed mechanistic way. Furthermore, the model does not take into consideration hepatic stellate cells, resident and monocyte-derived liver macrophages, additional immune cells such as lymphocytes or dendritic cells, the biliary system, different blood vessel compartments, hepatic blood flow or portal vein pressure, or the complex microarchitecture of the liver. Nevertheless, our model links information about liver regeneration and the interaction of different cell compartments (hepatocytes, Kupffer cells/macrophages, endothelial cells) from tissue to organ scale

data. While our model is simplification as explained above, it can be extended to describe further liver cell compartments such as stellate cells, bile duct cells, and additional immune cells beyond macrophages or different macrophage populations. In addition, hepatocyte subpopulations in different activation states such as quiescent, primed, and replicating cells could be included as described by Furchtgott et al. (2009). A numerical model can also be used to extrapolate additional time points. To enable unambiguous parameter estimation, an increase in model complexity should be accompanied by an increase of measurement values. In our study, we used *in vivo* and immunohistochemical analyses for modeling, resulting in mean time curves and therefore one model per group. If longitudinal *in vivo* measurements are used, a model could be applied to analyze individual mice, enabling statistical comparison of kinetic parameters between groups. Either way a numerical model could be used to investigate and explain the effects of genetic modifications, e.g., *csf1-*knock-out resulting in a reduced number of macrophages (Amemiya et al., 2011), macrophage depletion, or extended liver resection (Christ et al., 2017) on different aspects of liver regeneration and such data could be used to refine, extend, or validate our simplified model. In addition, advanced numerical models could be used to get a comprehensive insight into the interrelations between different cells and signaling pathways in chronic liver disease progression or in response to therapeutic interventions (Cook et al., 2015; Schwen et al., 2016).

In summary, based on non-invasive imaging and immunohistochemical analyses, we have established a mathematical model for liver regeneration describing the interrelations between hepatocytes (volume recovery), macrophages, and endothelial cells (angiogenesis) at organ and tissue scale. In this context, non-invasive imaging and suitable probes targeting cell populations such as macrophages are of great value for data acquisition in the course of liver regeneration at organ scale.

### ETHICS STATEMENT

All animal experiments were performed according to German legal requirements and animal protection laws and were approved by the Authority for Environment Conservation and Consumer Protection of the State of North Rhine-Westphalia (LANUV).

### AUTHOR CONTRIBUTIONS

SZ and AM performed the experiments, acquired, and analyzed the data and wrote the manuscript. AR and DD participated

### REFERENCES

in the experiments and data analyses. FG performed the modeling and FMT-μCT reconstruction. SW participated in the *in vivo* experiments. JFB participated in probe generation. UA participated in supervision and study design. JGB participated in study design. FK participated in study design, provided supervisory support, and edited the manuscript. WL designed and coordinated the experiments, provided supervisory support, wrote, and edited the manuscript. All authors read and approved the final version of the paper. SZ and AM contributed equally to this work.

### FUNDING

This study was supported by the BMBF (German Federal Ministry of Education and Research) Projects "The Virtual Liver Network" (Grant No. 0315743 and 0315731) and LiSyM (Grant No. 031L0041), and the German Research Foundation (GR 5027/2-1).

including 3D whole body organ segmentations. *Sci. Data* 5:180294. doi: 10.1038/sdata.2018.294

Wen, Y., Feng, D., Wu, H., Liu, W., Li, H., Wang, F., et al. (2015). Defective initiation of liver regeneration in osteopontin-deficient mice after partial hepatectomy due to insufficient activation of IL-6/Stat3 pathway. *Int. J. Biol. Sci.* 11, 1236–1247. doi: 10.7150/ijbs.12118

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

*Copyright © 2019 Zafarnia, Mrugalla, Rix, Doleschel, Gremse, Wolf, Buyel, Albrecht, Bode, Kiessling and Lederle. 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.*

digital media

of impactful research

article's readership