Grantee Research Project Results

2012 Progress Report: Model of Toxicant Response in Engineered Liver

EPA Grant Number: R834998
Title: Model of Toxicant Response in Engineered Liver
Investigators: Rajagopalan, Padmavathy , Murali, T. M. , Ehrich, Marion
Institution: Virginia Polytechnic Institute and State University
Current Institution: Virginia Tech
EPA Project Officer:
Project Period: June 1, 2011 through May 30, 2015
Project Period Covered by this Report: June 1, 2012 through May 30,2013
Project Amount: $750,000
RFA: Computational Toxicology: Biologically-Based Multi-Scale Modeling (2010) RFA Text |  Recipients Lists
Research Category: Chemical Safety for Sustainability

Objective:

Objective 1: Use 3D liver mimics as a platform for testing the hepatotoxicity of prototypic hepatotoxicants. 

Objective 2: Compute biological process linkage networks that anchor phenotypes triggered by hepatotoxicants to gene expression profiles. 

Objective 3: Discover the effects on hepatic cells of combinations of the toxicants examined in Objectives 1 and 2. Prioritize toxicity pathways for further proteomic and metabolomic analyses.

Our overall approach has the following components. We are assembling 3D liver mimics composed of hepatocytes and liver sinusoidal endothelial cells, separated by a nanoscale polyelectrolyte layer that mimics the Space of Disse. We are measuring the phenotypes triggered by two prototypic hepatotoxicants and their mixtures in the liver mimics. By analyzing the phenotypes, we plan to select specific concentrations of individual toxicants and mixtures for further analysis using DNA microarrays. We are integrating the transcriptional data with comprehensive molecular interaction networks and functional annotation databases to compute biological process linkage networks. Based on a careful comparative analysis of the linkage networks, we plan to prioritize proteins and metabolites for detailed experimental study.

We expect that this approach will establish a combined experimental and computational pipeline for toxicity testing and risk assessment centered on 3D liver mimics and biological process linkage networks. Our approach allows perturbations in toxicity pathways to be traced across multiple scales: from regulatory proteins and transcriptional perturbations to communication between different hepatic cell types. Unique aspects of our approach include modeling the cell’s response to combinations of toxic chemicals and in more than one hepatic cell type.

Progress Summary:

Objective 1: The design of organotypic liver models and investigating the hepatotoxicity of acetaminophen (APAP).

 1A. 3D Organotypic Liver Models:

We have been refining the liver models developed by our group.  Our current model incorporates primary rat hepatocytes, liver sinusoidal endothelial cells (LSECs), Kupffer cells and a detachable polymeric Space of Disse (pSD).  In this liver model the nanoscale polymeric membrane mimics the Space of Disse, and we have incorporated more than two hepatic cell types (Figure 1).  It is of significance to note that this liver model simultaneously maintains their phenotypes and cellular ratios.   In this study, detachable, nanoscale PEMs whose properties were tuned to those of the Space of Disse played a critical role in maintaining a physical barrier between hepatic parenchymal and non-parenchymal cells. Remarkably only the 3D organotypic hepatic model simultaneously exhibited proliferation of all cell types while maintaining cell ratios observed in vivo. We hypothesize that this polymeric interface promoted heterotypic cellular interactions via soluble molecules.  The key characteristics of this liver model are the following

1. Maintenance of phenotype of all three cell types up to 16 days.
2. Enhanced hepatic function (urea and albumin production).
3. In vivo like trends upon addition of KCs (manifested by the reduction of CYP1A1 activity upon addition of KCs).
4. Proliferation of all cell types while maintaining cellular ratios that are almost identical to those found in vivo. This is a very encouraging result since it demonstrates that these liver models can serve as good physiological models.
5. Gene expression data that validate the experimental trends of proliferation and CYP activity.

We highlight our results from the gene expression analysis in Section 1B of this report since they are relevant to this EPA-funded project. 

Figure 1. A schematic of the assembly of 3D organotypic liver models. 

1B. Gene Expression Analysis:

Gene expression data were obtained on day 12 and analyzed for hepatocyte monolayers (HM), collagen sandwich (CS), co-cultures of hepatocytes and LSECS (2DHL), co-cultures of hepatocytes, LSECs and Kupffer cells (2DHLK), a 3D culture of hepatocytes+pSD+LSECs (3DHL) and a 3D culture of hepatocytes +pSD+LSECs+KCs (3DHLK) . All samples passed the quality controls imposed by the Simpleaffy package (1). The samples were normalized using the GeneChip RMA (GCRMA) method (2). For each of the contrasts functional enrichment was performed on the normalized data using the Gene Set Enrichment Analysis (GSEA) package (Figure 2).

To obtain information on hepatocyte proliferation and enhanced function in 3DHL and 3DHLK, genome-wide gene expression data were obtained from hepatocytes in 3DHL and 3DHLK cultures and analyzed using GSEA.  The gene set “Proteinaceous extracellular matrix” (GO: Gene Ontology 0005578) was highly up-regulated (p-value 0, rank 1). Leading edge genes annotated to this GO term include over seven types of collagen, nidogen 2, and laminins.  These proteins are typically found in basement membranes (3).  For example, nidogen 2 is known to link collagen and laminin molecules.  Although basement membranes are not usually associated with adult liver, they are expressed during hepatic regeneration and liver development (3, 4).  The Reactome pathway “NCAM1 interactions” (p-value 1.7 x 10^-3, rank 27) was also up-regulated. The neural cell adhesion molecule, NCAM1 is a surface glycoprotein belonging to the immunoglobulin family that plays a role in liver development (5). While NCAM1 is itself not highly up-regulated in the comparison between 3DHL and 3DHLK, the leading edge of this pathway includes several types of collagens, many of which are also members of the leading edge of “Proteinaceous extracellular matrix” gene set.

The gene set corresponding to the GO term “Cell migration” (GO:0016477, p-value 6.8 x 10^-4, rank 22) is highly up-regulated in hepatocytes in 3DHLK. Of particular interest in this leading edge was the gene cadherin 13 (CDH13 or T-cadherin). Although this protein is shown to be very weakly expressed in normal hepatocytes (6), a recent study has shown that hepatocellular functions are enhanced in the presence of CDH13 (7).  Western blots for CDH13 demonstrated that its precursor protein (MW 130kDa) was observed in 3DHLK but only weakly expressed in HM, CS and 2DHL cultures (Figure 3B). 

The Reactome pathway for Phase 1 functionalization of compounds is up-regulated with p-value 2.7 x 10^-4 and rank 15.  Phase 1 CYP enzymes are known to catalyze the detoxification of a wide range of compounds and they exhibit broad substrate specificity.   Although the leading edge of this gene set contains several members of the cytochrome P450 family, CYP1A1 was not a member of this leading edge. In fact, it was among the highly down-regulated genes in the comparison of 3DHLK to 3DHL cultures, corroborating our earlier observation that the presence of KCs decreased CYP1A1 activity.

Figure 2:  Gene expression analysis of hepatocytes. The 30 most differentially expressed gene sets identified by GSEA. The table shows the results for all pairs of culture conditions that were compared. Yellow highlights indicate gene sets discussed in the text. 

Figure 3. (A) Statistics and leading edges of gene sets discussed in the text. (B) Western immunoblotting for CDH13.

1C. Investigating APAP hepatotoxicity:

Over this past year, we have conducted multiple experiments on APAP exposure to in vitro cultures of hepatocytes.  In Figure 4, we outline the pathway by which APAP confers hepatotoxicity.   We have conducted measurements on hepatocyte monolayers (HMs), collagen sandwich (CS) and a 3D liver model comprised of primary rat hepatocytes, a polymeric Space of Disse, and primary liver sinusoidal liver endothelial cells (LSECs).  

Figure 4.  Pathway for APAP-related hepatotoxicity. 

1D. Expression of CYP2E1:  

APAP is metabolized by CYP2E1 enzymes.  The expression of this enzyme is known to decrease in in vitro rat hepatocyte cultures.  For this reason, we used

Western immunoblotting to determine the levels of this enzyme on days 2 and 4 in culture.  Hepatocytes were lysed in each of these cultures and levels of CYP2E1 (in the absence of an inducer) was determined (Figure 5).  By day 4, the level of CYP2E1 goes down in the HM culture but not in CS or the 3D liver models. 

Figure 5: Western immunoblotting for CYP2E1. 

1D. Dose response for APAP:  

A dose response study was conducted on HM and CS cultures with APAP dose ranging from 10-40 mM.  We had to go to higher doses since rat hepatocytes are not very sensitive to APAP in comparison to mouse hepatocytes.  Based on our results (Figure 6), when APAP was administered on day 4 at a concentration of 40 mM it elicited the best response (greatest cell death). By day 12, the HMs were no longer sensitive to APAP and the CS culture showed a small decrease in viability only at 40mM.  In the case of the 3D liver models, even on day 12, there is a decrease in viability when the APAP dose was either 20 or 40 mM. Cell viability was assessed by a combination of MTT assays and fluorescent staining methods. 

 

Figure 6:  Hepatocyte viability upon APAP exposure.  All data were compared to control (untreated) cultures.   

1E. Changes in glutathione upon APAP exposure:  

When hepatocytes are exposed to APAP, an intermediate compound (NAPQI) is formed that results in toxicity and a significant reduction in glutathione (Figure 4).  We measured the depletion in glutathione upon administering 40mM of APAP to cultures and compared to controls (cultures that were not treated with APAP).  A larger decrease was detected in HM and CS cultures (Figure 7) in comparison to a 3D liver model comprised of hepatocytes, pSD and LSECs.  Based on the levels of CYP2E1 that are very similar in CS and the 3D liver model, we hypothesize that the LSECs may be imparting a cytoprotective effect.  We will conduct additional studies into the effects of APAP on LSECs alone and when they are in close proximity to hepatocytes. 

Figure 7:  Glutathione depletion in hepatocyte cultures upon exposure to 40 mM of APAP.  Data is presented as percent glutathione remaining relative to the control (untreated culture).  

Objective 2: Predicting Signaling Pathways

In our last progress report, we presented a new method for computing Biological Process Networks (BPNs) that seeks to ease the difficulty of interpretation of massive systems biology datasets. In the past year, we have developed a novel approach to predict signaling pathways solely from a list of receptors and transcription factors.

Background. As is well known, cells are constantly bathed in a sea of external signals requiring cellular response. An intricate collection of molecular interactions guides the propagation of external signals through the cell and into the nucleus, thereby eliciting a response in the cell's transcriptional profile. These signaling pathways are essential for the proper functioning of the cell. Identifying and understanding the signaling pathway that responds to an external signal is a major focus of systems biology.  

Signaling pathways are commonly presented in the literature using highly stylized images. While these diagrams are informative, they are primarily visual but cannot be used for computation. A number of repositories including KEGG (8), NetPath (9), and Reactome (10) maintain manuallycurated collections of interactions for human signaling pathways.  These representations provide a more detailed description of the pathway, listing a variety of different types of proteinprotein interactions involved in the pathway (e.g., phosphorylation, physical, translocation, etc.). In our work, we represent a pathway as the collection of proteins through which the signal propagates and the interactions that are activated in response to a bound ligand. We convert each pathway to a directed graph 𝑃= (𝑉_!, 𝐸_!) with a set 𝑉_! indicating pathway proteins and a set of (directed or undirected) interactions 𝐸_!.   

We introduce the following problem of reconstructing signaling pathways from a molecular interaction network. Given a global network of known interactions in an organism (the interactome) and the interactions in a signaling pathway, reconstruct the precise interactions from the signaling pathway when only the receptors and downstream targets of the pathway are known. Defining the problem more formally, let the interactome 𝐺= (𝑉, 𝐸) be a directed network with node set 𝑉 and edge set 𝐸, where each edge (𝑢, 𝑣) is assigned a weight 𝑤_!". The network 𝐺 represents an organism's interactome, which may integrate known interactions from publiclyavailable repositories.  Let 𝑃= (𝑉_!, 𝐸_!) be a subgraph of 𝐺 denoting the collection of interactions involved in a specific signaling pathway, e.g., the Wnt signaling pathway. Let 𝑆   ⊂𝑉_! and 𝑇⊂𝑉_! be disjoint sets of nodes that represent receptors and transcription factors (TFs)  in the pathway, respectively. The set 𝑆 represents the collection of receptor proteins for the pathway that occupy the plasma membrane and initiate a cell's response to an external signal by receiving a signaling ligand (e.g., the Frizzled family of proteins in the Wnt signaling pathway). The set 𝑇 indicates the downstream TFs that are regulated in response to the signal propagation through the pathway. Given 𝐺, 𝑆, and 𝑇, we seek to reconstruct the interactions 𝐸_! in the pathway 𝑃.

Algorithm. Given the graph 𝐺 and receptors source nodes 𝑆, we use the PageRank algorithm to compute per-node visitation probabilities during a biased random walk about the graph. We set the starting probability 𝑠_! for each node 𝑣 as follows: 𝑠_! =1/|𝑆| for  𝑠   ∈𝑆  and 𝑠_! =0 otherwise. Let 𝑤_!" be the weight of the directed edge (𝑢, 𝑣). We normalize the edge weights such that the total weight of the outgoing edges for each node in the graph is 1. We start the process by placing the walker on node 𝑢 with probability 𝑠_! for each node 𝑢∈𝑉.  Now, the walker moves in the network according to the following rules, where 0   ≤𝑞   ≤   1 is a parameter:

• Teleport: With probability 𝑞𝑠_𝑥, she teleports to any node 𝑥∈𝑉, including her current node and its neighbors; the total probability of teleporting from 𝑢 is 𝑞.
• Walk: She can move from her current node 𝑢 to any of 𝑢’s out-neighbors 𝑣∈𝑁_! with probability proportional to (1−𝑞)𝑤_!"; thus, the total probability of walking to some neighbor of 𝑢 is 1−𝑞.

The parameter 𝑞 provides control over how often the walker teleports back to one of the receptors in 𝑆. Under certain conditions that hold for our data, we can prove that this random walk converges to a unique stationary distribution.  In this case, we can compute the stationary visitation probability 𝑝_!for each node 𝑣∈𝑉 by solving the following linear system: 𝑤!"𝑝! 𝑝_!

We can solve the system by inverting an appropriate matrix. In practice, we used the wellknown iterative power method to compute the stationary probabilities. We then computed a flux score 𝑓_!"  for each edge (𝑢, 𝑣) as follows: 𝑤!"𝑝!

where 𝑑_! is the out-degree of node 𝑢. Thus, edges with a high probability of being traversed receive a large edge flux while edges with a low probability of traversal receive a low edge flux. We ranked all the edges in decreasing order of flux values. We call this algorithm PRflux.

Data Sources. We collected direct interactions reported in all NetPath human signaling pathways (9). At the time of writing, NetPath contained 32 signaling pathways. We identified a set of signaling receptors in each pathway to represent the source nodes for that pathway as follows. We collected a list of signal receptors from Almen et al. (11). We computed the intersection of the proteins in each pathway with the receptors from this list, and we manually inspected these lists for any erroneous or missing receptors. Since we sought signaling mechanisms that cascade to downstream transcriptional events, we defined the targets as the set of TFs involved in each pathway. We retrieved the set of human TFs reported in two previously-published studies: i) all TFs reported by Ravasi et al. (12) and ii) high-quality TFs from Vaquerizas et al. (13). We utilized only the largest connected component of each NetPath pathway.  Of the 32 pathways, we retained 18 pathways whose largest connected component contained both receptors and TFs and included at least five but no more than 100 internal nodes (i.e., non-target or non-source); we analyzed these 18 pathways in our results. We constructed a “NetPath union” network that includes the union of edges from all NetPath pathways and contains 1226 nodes and 7171 interactions. 

Results. We compared PRflux to several algorithms published in the literature: ANAT (14), a technique for connecting sources to targets using compact networks; ResponseNet (15), a flowbased algorithm for the same purpose (abbreviated as RN); and eQED, a linear program-based technique that models networks as electric circuits (16). eQED was originally developed to analyze expression quantitative trait loci (eQTL) data, for which eQED predicted the causal gene involved in a particular change in a gene's expression. eQED computes a score for each edge in the network, and we use these scores to infer candidate edges in a signaling pathway.

We tested the ability of each algorithm to reconstruct each of the 18 NetPath pathways given a background network containing the pathway as a subnetwork and the source and target nodes for the pathway. We tested the algorithms using the union of the NetPath pathways as the background network. For each pathway, we ranked the edges in the network according to each algorithm and computed precision-recall curves based on the ranking. We considered as positive predictions the set of interacting pairs of nodes in that pathway. We defined the negative interactions as all interacting pairs that are not adjacent to any node in the pathway.  Note that the number of negatives was much larger than the number of positives, making this an extremely difficult prediction problem. 

As illustrated in Figure 8, PRflux clearly dominates the other algorithms in the precision recall curves at all values of recall. Remarkably, PRflux attains a precision of approximately 80% even at values of recall as high as 60%. RN computes very small networks; hence, its recall does not reach even 15%. For the same reason, the various versions of ANAT also cannot achieve high values of recall. eQED is competitive with PRflux at values of recall less than 20%. However, its precision drops off dramatically around a recall of 30%. 

Figure 8: Precision-recall plots of four algorithms; we ran ANAT for three values of an internal parameter, 0, 0.25, and 0.5. We aggregated results over all NetPath pathways. The plots for ANAT are single points, since ANAT computes a single network rather than ranking all the edges.

Since signaling pathways are highly incomplete, we expect that many of the predictions made by PRflux may represent valid pathway interactions that simply have not yet been added to the pathway by manual curation. We anticipate that these high-confidence predictions adjacent to the pathway are ideal candidates for further exploration in extending signaling pathways.

Future Activities:

We have conducted experiments on the changes in hepatic function upon addition of APAP.  We have used a prototypic toxicant, APAP, and demonstrated that the liver models assembled in our group are more sensitive to this chemical in comparison to HM and CS cultures.  In year 3 of the project, we plan to investigate toxicity induced by carbon tetrachloride and dichlorotheylene.  We will generate genomewide gene expression data that measure cellular responses.  We seek to apply the MCMC-BPN algorithms (developed in year 1 of the project) and the signaling pathway prediction algorithms to these data in order to summarize the liver mimic's responses to toxicants and to decipher intercellular signaling pathways. We will use these analyses to prioritize experimental validations.

References:

1. Wilson CL, Miller CJ. Simpleaffy: a BioConductor package for Affymetrix Quality Control and data analysis. Bioinformatics 2005;21:3683-3685.
2. Wu Z, Irizarry RA, Gentleman R, Martinez-Murillo F, Spencer F. A Model-Based Background Adjustment for Oligonucleotide Expression Arrays. Journal of the American Statistical Association 2004;99:909-917.
3. Kyriakakis E, Maslova K, Philippova M, Pfaff D, Joshi MB, Buechner SA, Erne P, et al. TCadherin Is an Auxiliary Negative Regulator of EGFR Pathway Activity in Cutaneous Squamous Cell Carcinoma: Impact on Cell Motility. J Invest Dermatol 2012.
4. Arias IM, Boyer JL, Chisari FV, Fausto M, Schachter D, Shafrtx DA. The liver: biology and pathobiology. 4 ed: Lippincott Williams and Wilkins, 2001.
5. Libbrecht L, Cassiman D, Desmet V, Roskams T. Expression of neural cell adhesion molecule in human liver development and in congenital and acquired liver diseases. Histochem Cell Biol 2001;116:233-239.
6. Riou P, Saffroy R, Chenailler C, Franc B, Gentile C, Rubinstein E, Resink T, et al. Expression of T-cadherin in tumor cells influences invasive potential of human hepatocellular carcinoma. FASEB J 2006;20:2291-2301.
7. Khetani SR, Chen AA, Ranscht B, Bhatia SN. T-cadherin modulates hepatocyte functions in vitro. FASEB J 2008;22:3768-3775.
8. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res 2008;36:D480-484.
9. Kandasamy K, Mohan SS, Raju R, Keerthikumar S, Kumar GS, Venugopal AK, Telikicherla D, et al. NetPath: a public resource of curated signal transduction pathways. Genome Biol;11:R3.
10. Joshi-Tope G, Gillespie M, Vastrik I, D'Eustachio P, Schmidt E, de Bono B, Jassal B, et al. Reactome: a knowledgebase of biological pathways. Nucleic Acids Res 2005;33:D428-432.
11. Almen MS, Nordstrom KJ, Fredriksson R, Schioth HB. Mapping the human membrane proteome: a majority of the human membrane proteins can be classified according to function and evolutionary origin. BMC Biol 2009;7:50.
12. Ravasi T, Suzuki H, Cannistraci CV, Katayama S, Bajic VB, Tan K, Akalin A, et al. An atlas of combinatorial transcriptional regulation in mouse and man. Cell 2010;140:744-752.
13. Vaquerizas JM, Kummerfeld SK, Teichmann SA, Luscombe NM. A census of human transcription factors: function, expression and evolution. Nat Rev Genet 2009;10:252-263.
14. Yosef N, Ungar L, Zalckvar E, Kimchi A, Kupiec M, Ruppin E, Sharan R. Toward accurate reconstruction of functional protein networks. Mol Syst Biol 2009;5:248.
15. Yeger-Lotem E, Riva L, Su LJ, Gitler AD, Cashikar AG, King OD, Auluck PK, et al. Bridging high-throughput genetic and transcriptional data reveals cellular responses to alphasynuclein toxicity. Nat Genet 2009;41:316-323.
16. Suthram S, Beyer A, Karp RM, Eldar Y, Ideker T. eQED: an efficient method for interpreting eQTL associations using protein networks. Mol Syst Biol 2008;4:162.

Journal Articles on this Report : 3 Displayed | Download in RIS Format

Publications Views

Other project views:	All 34 publications	14 publications in selected types	All 13 journal articles

Publications

Type	Citation	Project	Document Sources
Journal Article	Larkin AL, Rodrigues RR, Murali TM, Rajagopalan P. Designing a multicellular organotypic 3D liver model with a detachable, nanoscale polymeric Space of Disse. Tissue Engineering Part C: Methods 2013;19(11):875-884.	R834998 (2012) R834998 (2013) R834998 (Final)	Full-text from PubMed Abstract from PubMed Associated PubMed link Abstract: Mary Ann Liebert-Abstract Exit
Journal Article	Lasher CD, Rajagopalan P, Murali TM. Summarizing cellular responses as biological process networks. BMC Systems Biology 2013;7:68.	R834998 (2012) R834998 (2013) R834998 (Final)	Full-text from PubMed Abstract from PubMed Associated PubMed link Full-text: BMC-Full Text PDF Exit Abstract: BMC-Abstract and Full Text HTML Exit
Journal Article	Rajagopalan P, Kasif S, Murali TM. Systems biology characterization of engineered tissues. Annual Review of Biomedical Engineering 2013;15:55-70.	R834998 (2012) R834998 (2013) R834998 (Final)	Abstract from PubMed Full-text: Boston University-Full Text PDF Exit Abstract: Annual Reviews-Abstract Exit

Progress and Final Reports:

Original Abstract

2011 Progress Report

2013 Progress Report

Final Report