Jennifer S Park3, Walter K Schlage1, Brian P Frushour3, Marja Talikka2, Gary Toedter3, Stephan Gebel1, Renee Deehan3, Emilija Veljkovic2, Jurjen W Westra3, Michael J Peck2, Stephanie Boue2, Ulrike Kogel2, Ignacio Gonzalez-Suarez2, Arnd Hengstermann1, Manuel C Peitsch2, and Julia Hoeng2*
Received date: December 19, 2012; Accepted date: January 31, 2013; Published date: February 04, 2013
Citation: Park JS, Schlage WK, Frushour BP, Talikka M, et al. (2013) Construction of a Computable Network Model of Tissue Repair and Angiogenesis in the Lung. J Clinic Toxicol S12:002. doi: 10.4172/2161-0495.S12-002
Copyright: © 2013 Park JS, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Clinical Toxicology
We have recently developed methodologies to quantify the biological impact of exposure to environmental toxicants via the use of computable network models and network scoring methods. In this study, we extend our collection of lung pathophysiology specific network models to tissue repair and angiogenesis. It is important to understand the molecular mechanisms of wound healing which, if unresolved, could eventually progress to irreversible disease. The Tissue Repair and Angiogenesis (TRAG) Network consists of nine modular subnetworks that describe the following processes: hypoxia-inducible factor 1 alpha (HIF1A) signaling, sprouting and tubulogenesis, Vascular Endothelial Growth Factor (VEGF)-mediated angiogenesis, growth factor-mediated angiogenesis, immune regulation of angiogenesis, immune regulation of tissue repair, cell migration, differentiation of progenitor cells and fibrosis. We used a data-driven approach to augment the initial literature-based network, and to evaluate a portion of the network using two independent gene expression data sets. This approach increases the confidence in the network’s ability to accurately describe tissue repair processes. The TRAG Network serves as a valuable research tool for assessing the biological impact of exposure to environmental insults and in understanding the initial molecular events that may lead to disease. The TRAG network, which consists of 666 nodes linked by 1215 relationships or edges covering 1371 PubMed IDs, is expressed in the Biological Expression Language and made available in computer readable formats including XGMML and .xls.
Computable; Network model; Tissue repair; Angiogenesis
BEL: Biological Expression Language; RCR: Reverse Causal Reasoning; ADAM: A Disintegrin and Metalloproteinase Domain; AKT1: Protein kinase B; AP1: Activator protein 1; COPD: Chronic Obstructive Pulmonary Disease; CTNNB1: Beta-catenin; DKK1: Dickkopf1; EGR1: Early Growth Response 1; EMT: Epithelial– Mesenchymal Transition; ETS2: v-etserythroblastosis virus E26 oncogene homolog 2; FRAP1: FK506 binding protein 12-rapamycin associated protein 1; HIF1A: Hypoxia-Inducible Factor 1; HYPs: Hypothesis; IFNG: Interferon gamma; IL: Interleukin; ITGB6: Integrin Beta 6; MAPK: Mitogen-Activated Protein Kinase; MMP: Matrix Metalloproteinase; mTOR: Mammalian Target of Rapamycin; PI3K: Phosphoinositide 3-kinase; PPARA: Peroxisome Proliferator-Activated Receptor Alpha; PPARG: Peroxisome Proliferator-Activated Receptor Gamma; RCR: Reverse Causal Reasoning; SHH: Sonic Hedgehog Homolog; SMAD: Mothers Against Decapentaplegic Homolog; SP1: Specificity Protein 1; SRC: Proto-oncogene tyrosine-protein kinase Src; STAT: Signal Transducers and Activators of Transcription; TGFB1: Transforming Growth Factor beta 1; TGFBR1: Transforming Growth Factor beta receptor 1; TNF: Tumor Necrosis Factor; TRAG: Tissue Repair and Angiogenesis; VEGF: Vascular Endothelial Growth Factor; WNT3A: Wingless-type MMTV integration site family, member 3A; YY1: Transcriptional repressor protein yin yang 1
The pulmonary system serves as an interface between an organism and the gaseous environment. Consequently, the pulmonary epithelium is exposed to a variety of environmental insults, including complex aerosols resulting from combustion such as vehicle exhaust and cigarette smoke. These aerosols contain a number of toxicants, including oxidants, in both gaseous and particulate phases. It is well known that exposure to these substances will increase the risk of developing lung disease, including Chronic Obstructive Pulmonary Disease (COPD) and lung cancer. With the ultimate aim to quantify the dose-dependent biological impact of these types of exposures we have recently developed methods to quantify the perturbation of biological networks [1,2]. These methods require two types of input: firstly, genome-scale experimental data such as transcriptomics and secondly, the networks models for which a Network Perturbation Amplitude score is to be computed [3,4]. As these methods are geared towards quantifying the biological impact of inhaled substances, the input networks should adequately describe the mechanisms that are perturbed by these types of exposures. Our collection of published computable networks currently covers cell proliferation , and cellular stress . We are currently working on additional network models on DNA damage and repair, autophagy, cell death and senescence, as well as on inflammatory processes in healthy lung tissue. The aim of this study is to expand our collection of lung pathophysiology specific computable network models and thereby cover a broader range of disease-relevant mechanisms. This will further improve our ability to quantify the biological impact of environmental insults. We thus set out to construct and evaluate a network model describing processes involved in tissue repair and angiogenesis in the context of non-diseased pulmonary cells and tissues.
The pulmonary system will generally respond to environmental insults by mounting an inflammatory response which will result in resolution, chronic inflammation or fibrosis . The ultimate outcome is influenced by the nature of the insult and the lung tissue in which the inflammation has occurred . Angiogenesis is an important component of the repair mechanisms that occurs in following tissue damage . It is a complex process comprised of a number of mechanisms, including sprouting and tubulogenesis mediated by growth factors (including vascular endothelial growth factor (VEGF), a major driver of angiogenesis), Hypoxia-Inducible Factor 1 (HIF1A) regulation, and immune processes that regulate angiogenesis [8-10]. Angiogenesis has a complex interplay with inflammation, and hypoxic inflammatory tissue promotes angiogenesis through activation of hypoxic mechanisms as well as cytokines that have both inflammatory and angiogenic effects. VEGF plays a role not only in angiogenesis, but cell proliferation and apoptosis as well, therefore angiogenesis has far-reaching effects when activated. Tissue repair is also a multifaceted process influenced by immune processes including cytokines and chemokines that induce cell migration and adhesion as well as differentiation of progenitor cells, contributing to wound healing. As a final step of wound healing, fibrosis and Epithelial–Mesenchymal Transition (EMT) mechanisms lead to closure and repair of the wound [11-14]. Interactions of many different cell types including fibroblasts, epithelial cells and immune cells orchestrate this elaborate interplay of processes to result in tissue repair.
In this study, our objective was to capture these tissue repair and angiogenesis mechanisms in modular subnetworks comprising the TRAG Network. Tissue repair subnetworks were evaluated by applying two independent molecular profiling data sets with measured tissue repair endpoints available in the literature (i.e., bleomycin-induced inflammation and mechanical injury), providing confidence in the ability of the tissue repair subnetworks to accurately describe specific tissue repair mechanisms induced by the different perturbations. Angiogenesis subnetworks were not evaluated due to the lack of an appropriate data set with measured angiogenic endpoints that could be used to confirm inclusion of relevant angiogenesis mechanisms.
The TRAG Network was constructed using the same iterative process used to create previously published network models [3,4]. The network was populated with nodes and edges from two main sources: firstly, prior knowledge described in the scientific literature and which is now referenced in the network models and secondly, results obtained from the computational analysis of two relevant transcriptomic profiling data sets via Reverse Causal Reasoning (Figure 1).
Figure 1: Workflow used to construct and evaluate the TRAG Network.
The literature-based network was initially constructed from causal relationships extracted from relevant scientific literature following the definition of network boundaries. The network was then augmented with additional nodes derived from RCR analysis of two transcriptomic data sets with perturbations known to stimulate tissue repair and/or angiogenesis. Manual review and refinement resulted in the final TRAG Network.
Boundary criteria that define the network content were identified and applied to the network construction. Starting with a list of nodes identified by surveying published literature in the areas of pulmonary tissue repair and angiogenesis, causal relationships were identified describing the mechanistic links between these nodes with literature support for both normal lung and vascular cell types. For angiogenesis, the search was restricted to the context of non-diseased and particularly non-tumor lung tissue. In the case of tissue repair, the focus was on wound healing and acute injury (e.g., bleomycin, mechanical stress, etc.) that lead to tissue repair. Disease-based mechanisms such as emphysema, asthma, COPD and bacterial/viral infection were excluded in order to represent the “normal” healing processes in otherwise healthy tissue. Likewise, embryonic contexts were also excluded as they describe organogenesis and development, which involves processes that may fundamentally differ from tissue repair.
In cases where the relevant experiments have not been published in these contexts, relationships derived from non-lung contexts using cell types similar to those from normal lung (e.g., fibroblasts, epithelial cells, endothelial cells, etc.) were used. “Canonical” pathways that are well-known in the literature were also included in the TRAG Network even if literature support explicitly demonstrating the presence of the pathway in normal lung or cardiovascular tissues was not found. For direct and proximal connections such as a kinase phosphorylating a residue on a target or protein-protein interactions, evidence from cellfree in vitro systems were also used when normal lung or cardiovascular tissues were not available. Lastly, human contexts were favored, but relationships derived from rodent (specifically mouse and rat) contexts were also included and homologized.
Using the established boundary criteria, a literature network model was created by compiling causal relationships extracted from the Selventa knowledgebase.
The Selventa knowledgebase, previously described in [3,4] is a unified collection of over 1.5 million elements/nodes of biological knowledge and 7 million connecting edges. Edges are relationships between nodes, and may be either causal or non-causal. The knowledgebase is derived from peer-reviewed scientific literature describing findings in human, mouse and rat species. The edges in the knowledgebase are coded using Biological Expression Language (BEL), a language that allows for the semantic representation of life science relationships in a computable format . BEL enables the development of computable pathway models comprised of cause and effect relationships. A context-specific collection of causal edges is referred to as either the human, mouse or rat Knowledge Assembly Model (KAM) . Using these KAMs, the final network models were homologized to human, mouse and rat.
When causal edges did not exist in the knowledgebase to support network building, they were identified and manually curated from peer-reviewed literature into the knowledgebase. The literature network model was then augmented with additional nodes derived from the computational analysis of molecular profiling data sets relevant to tissue repair or angiogenesis using RCR (Figure 1).
Analysis of transcriptomic data sets
Data sets were prioritized according to three main criteria: 1) whether the biological process relevant to the TRAG Network was induced in non-diseased cell types found in normal lung, 2) whether phenotypic endpoint data were available to provide experimental evidence for the biological relevance of the transcriptomic data set, and 3) the statistical rigor of the design of transcriptomic profiling experiments. Four publicly available gene expression data sets were used in the construction or evaluation of the TRAG Network (GSE11341 Hypoxia ; GSE25640 Bleomycin ; GSE18800 Bleomycin ; GSE5372 Mechanical Injury ) (Table 1). Data sets were downloaded from Gene Expression Omnibus (GEO; http://www. ncbi.nlm.nih.gov/ geo). The analysis of the data sets was performed as previously described [3,4]. Probe sets were considered to have statistically-significant changed expression levels in a specific comparison if they had an FDR adjusted p-value of lower than 0.05 and an absolute fold change greater than 1.3. Probe set differences were considered significant only if the average expression intensity was above 150. Gene expression changes that met these criteria were called “State Changes” and have the directional qualities of “increased” or “decreased”, (upregulated or downregulated, respectively) in response to the experimental condition.
|Network Building Data Sets||Network Evaluation Data Sets|
|Data Set ID||GSE25640||GSE11341||GSE18800||GSE5372|
|Pub Med ID||21602491||18469115||19966781||17164391|
|Context||In vivo||In vitro||In vivo||In vivo|
|Cell Type||Whole Lung||Lung Endothelial cells||Whole lung||Airway epithelium from non-smokers|
|Insult||Bleomycin 7.5 U/kg||Hypoxia 1% O2||Bleomycin 1 mg/kg||Brushing to denude epithelium|
|Timepoint(s)||21 days post treatment||48 hours||14 days post treatment||7 days post treatment|
|Control||Saline Day 21||Normoxia 21% O2||Saline Day 0||Day 0|
|Measured Endpoint||Hydroxyproline, BAL inflammatory cells, Lung T-cells||GREM1 protein, CXCR7 gene expression||Hydroxyproline, at Day 21; TGFβ1 and BAL inflammatory cells at Day 7 and 21||Partially differentiated epithelial layer, <1% of sample is inflammatory cells|
|# State Changes||2672||639||3425||1101|
|# of HYPs||501||292||521||265|
Table 1: Data sets analyzed by RCR for network augmentation and evaluation.
Reverse causal reasoning
Reverse Causal Reasoning (RCR) analysis of the transcriptomic data sets was used to aid in the selection of nodes for the TRAG Network. RCR interrogates a literature-based KAM to identify upstream controllers of the mRNA State Changes (gene expression changes) observed in the data set. For the hypoxia and mechanical injury data sets, a human KAM was used, while a mouse KAM was used for the bleomycin data sets. The potential upstream controllers identified by RCR are called “hypotheses” (HYPs), as they are statistically significant potential explanations for the observed RNA State Changes. A node can be predicted as a HYP via RCR (RCR-capable) when there are at least four gene expression changes causally downstream of the node in the Selventa knowledgebase.
Each HYP is scored according to two probabilistic scoring metrics, richness and concordance. Richness is the probability that the number of observed RNA State Changes connected to a given HYP could have occurred by chance alone, calculated using the hypergeometric distribution. Concordance is the probability that the number of observed RNA State Changes that match the direction of the HYP (e.g., increased or decreased activity or abundance of a node) could have occurred by chance alone, calculated using a binomial distribution. HYPs meeting both richness and concordance p-value cutoffs of 0.1 were considered to be statistically (although not necessarily biologically) significant. For the purposes of network model construction, only HYPs with p<0.05 were used for enhancement of the network.
RCR-derived HYPs were included as new nodes in the TRAG Network if they had literature support (determined by manual curation, “HYP vetting”) for a mechanistic role in the process and context of interest, thereby uncovering relevant nodes and edges that were not identified during the construction of the literature network. RCRbased augmentation of the TRAG Network was performed using two transcriptomic data sets (one for tissue repair and one for angiogenesis), referred to as “building” data sets (Table 1). The tissue repair set selected was the GSE25640 Bleomycin data set , which examined the effect of bleomycin on mouse lung after 21 days of exposure. The angiogenesis data set selected was the GSE11341 Hypoxia data set , which examined the effect of hypoxia on human lung microvascular endothelial cells at 48 hours. In the “HYP vetting” step, a total of 414 HYPs derived from RCR were evaluated from these two data sets, 49 of which were considered biologically plausible in the context of previous literature reports.
As a final step in the construction of the TRAG Network, the nodes and edges were manually reviewed and refined to include relevant context and pathways (e.g., by additional specific literature curation), producing the final TRAG Network model.
Two transcriptomic data sets were selected to evaluate the final TRAG Network tissue repair subnetwork models: GSE18800 Bleomycin data set , which examined the effect of bleomycin treatment on mouse lung at day 14, and GSE5372 Mechanical Injury data set , which examined the effect of mechanical injury due to brushing on human airway epithelial cells in non-smokers after 7 days. Both data sets were relevant for tissue repair and independent from those used for network construction, and were selected using the same criteria described for the building data sets.
Network structure and content
Here we report the construction of a network model focused on tissue repair and angiogenesis processes. The TRAG Network is provided in .XGMML format for human (Supplementary file 1), mouse (Supplementary file 2) and rat (Supplementary file 3), and .xls format (Supplementary file 4), and can be viewed using freely available visualization software such as Cytoscape . The human TRAG Network contains 666 unique nodes and 1215 unique edges, supported by 1371 PubMed-referenced literature citations (Table 2). 306 of the 666 nodes in the human TRAG Network are RCR-capable (see Materials and Methods) .
|Nodes (RCR-capable)||Edges (Causal)|
|TRAG Network||666 (306)||1215 (828)|
|Tissue Repair subnetworks|
|Cell Migration and Adhesion in Wound Healing||201 (70)||317 (191)|
|Differentiation of Progenitor Cell in Wound Healing||38 (13)||44 (13)|
|Fibrosis and EMT||157 (97)||293 (203)|
|Immune Regulation of Tissue Repair||127 (83)||254 (201)|
|Growth Factor-Mediated Angiogenesis||99 (63)||135 (78)|
|HIF1A Regulation and Signaling||45 (18)||48 (32)|
|Immune Regulation of Angiogenesis||50 (31)||54 (34)|
|Sprouting and Tubulogenesis||82 (36)||93 (60)|
|VEGF-Mediated Angiogenesis||45 (21)||61 (43)|
Table 2: TRAG Network statistics. For the full TRAG Network and for each TRAG subnetwork, the total number of nodes is given for the human network versions, along with the number of RCRcapable nodes in parentheses (a node known to causally regulate at least four downstream gene expressions in the knowledgebase, making it capable of being predicted as a HYP by RCR). Additionally, the total number of edges is shown, with the number of causal edges in parentheses. Note that the sum of the totals for the individual subnetworks is greater than the total number of unique nodes/edges in the full TRAG Network due to overlap of nodes/edges amongst the subnetworks.
The network was constructed using a modular design of nine discrete subnetworks that, when merged, form the integrated TRAG Network. Tissue repair processes were represented by four subnetworks and verified using evaluation data sets, and angiogenic processes were represented by five subnetworks but remain unverified due to the unavailability of appropriate evaluation data sets with measured angiogenic endpoints (Table 2). These subnetworks were organized to contain pathways that describe particular processes within angiogenesis and tissue repair. These pathways are made of nodes which represent biological entities such as protein abundances, mRNA expression levels and protein activities, or biological processes (e.g., EMT). Some nodes can be predicted as increased or decreased HYPs by RCR based on the gene expression changes in a data set (see Materials and Methods). A node in one subnetwork can also belong to another subnetwork and regulate nodes in other subnetworks, resulting in a highly interconnected model when the subnetworks are combined to yield the final TRAG Network.
Data sets from rodent bleomycin instillation and mechanical injury of the human airway, with measured tissue repair endpoints were used to evaluate the corresponding subnetworks. Bleomycin is a chemotherapeutic whose major side effect is pulmonary fibrosis and impaired lung function. Because bleomycin is used widely in animal models of fibrosis, including investigation of the molecular pathways leading to its inflammatory and fibrotic effects, we chose to use this perturbation to construct and evaluate the tissue repair subnetworks. RCR analysis of the two evaluation data sets resulted in the prediction of 521 HYPs for the GSE18800 Bleomycin data set and 265 HYPs for the GSE5372 Mechanical Injury data set. The results of this evaluation are shown in Table 3.
|Tissue Repair Subnetworks||# RCR-capable nodes||GSE18800 Bleomycin||GSE5372 Mechanical Injury|
|#consistent HYPs (%)||#inconsistent HYPs (%)||#consistent HYPs (%)||#inconsistent HYPs (%)|
|Cell Migration and Adhesion in Wound Healing||70||24 (34%)||5 (7%)||15 (21%)||2 (3%)|
|Differentiation of Progenitor Cell in Wound Healing||12||1 (8%)||1 (8%)||2 (17%)||0 (0%)|
|Fibrosis and EMT||97||43 (44%)||5 (5%)||24 (25%)||6 (6%)|
|Immune Regulation of Tissue Repair||83||40 (48%)||5 (6%)||15 (18%)||4 (5%)|
Table 3: Summary of results obtained with evaluation data sets. Shown are the numbers of RCR-capable nodes in each of the three subnetworks used for evaluation, as well as the number of HYPs matching these nodes that are predicted by RCR in the GSE18800 Bleomycin and GSE5372 Mechanical Injury data sets. HYPs may be predicted in directions either consistent or inconsistent with their regulation of the process in the subnetworks, relative to an increase in the process described by the subnetwork. Percentages are calculated based on the total number of RCR-capable nodes in each subnetwork.
Up to 48% of the RCR-capable nodes in a subnetwork were predicted by RCR in consistent directions (relative to an increase in the process described by the subnetwork) based on the GSE18800 Bleomycin data set, indicating that the networks contain tissue repair mechanisms relevant to bleomycin injury. Up to 25% of the RCR-capable nodes in a subnetwork were predicted by RCR in consistent directions based on the GSE5372 Mechanical Injury data set, indicating that although the tissue repair subnetworks contain mechanisms relevant to mechanical injury, there were less mechanisms from the mechanical injury data set present in the subnetworks compared to the bleomycin data set. This could be explained by the fact that there were less total mechanisms predicted in the mechanical injury data set compared to the bleomycin data set. In all cases, no more than 8% of the RCR-capable nodes were predicted in inconsistent directions. Because the Differentiation of Progenitor Cells in Wound Healing subnetwork largely consisted of gene expression nodes that could not be predicted as HYPs, this resulted in the low percentage of predicted nodes. The specific HYPs representing pathways in the Cell Migration and Spreading in Wound Healing, Fibrosis and EMT, and Immune Regulation of Tissue Repair subnetworks are elaborated in greater detail in the subsequent Results sections. Overall, the HYPs predicted by RCR in these data sets indicate activation of pathways that are in line with previous literature findings as well as the endpoints measured in the data set.
Fibrosis and EMT pathways predicted by RCR
GSE18800 Bleomycin Data Set: RCR predicted activation of many of the pathways described in the Fibrosis and EMT subnetwork and known in bleomycin literature. In addition, prediction of fibrotic pathways was in line with the fibrotic endpoints measured in the GSE18800 mouse study. 44% of the RCR-capable nodes in the Fibrosis and EMT subnetwork were predicted in a direction consistent with increased fibrosis (Table 3). Evidence for major fibrotic pathways including Tgf-beta, angiotensin II, β-catenin, Pi3k and Hedgehog was supported in the data set (Figures 2 and 3). Increased Tgfb1, its receptor Tgfbr1 and downstream effectors Smad1, Smad3 and Smad4 transcription factors were predicted to be increased in the data set. Tgfb1 and Smad family members have been widely studied in bleomycin-induced injury, and prediction of these HYPs is not only consistent with the known role of Tgfb1 in driving bleomycin-induced fibrosis, but also with increased Tgfb1 protein measured in the data set [18,21-28]. PPAR agonists have been shown to decrease fibrosis by inhibition of TGFB1 [29,30], and decreased Pparg activity along with increased Tgfb1 predicted after bleomycin treatment is consistent with this finding. Additionally, it has been shown that Tgfb1 inhibits Pparg, captured in the Fibrosis and EMT subnetwork. Pparg-deficient mice showed enhanced Bleomycin-induced skin fibrosis after 2-6 weeks, and enhanced sensitivity to Tgfb1 . Angiotensin II has been shown to promote bleomycin-induced lung fibrosis via Agtr1a, Tgfb1 and p38 Mapk family members [32-34], all of which are pathways predicted in the bleomycin data set. In addition to angiotensin II, increased activity of angiotensin-converting enzyme and angiotensin II receptor, Agtr1a was predicted in the bleomycin data set. Angiotensin II has been shown to promote fibrosis through Tgfb1 along with downstream signaling mediators Src and p38 Mapk family, both of which were predicted to have increased kinase activity in this data set [32-34]. In addition to regulating p38 Mapk family and Tgfb1, angiotensin can also regulate Pi3k/Akt, predicted to be increased in the data set and shown to regulate bleomycin-induced fibroblast proliferation and collagen production . The activation of mTOR/was indicated by a decreased prediction of the Sirolimus HYP, anmTOR/inhibitor. Because many experiments studying the mTOR pathway use Sirolimus (rapamycin), these experiments have been curated in the knowledgebase to the Sirolimus HYP. In addition to mTOR, increased beta-catenin (Ctnnb1) and beta-catenin activator Wnt3a were predicted to be increased, and WNT inhibitor Dkk1 was predicted to be decreased. This is in line with published results, which have shown that downstream of Akt-signaling, the bleomycin-induced extracellular matrix synthesis in fibrosis is mediated by mTOR , and beta-catenin .
Figure 2: Fibrosis and EMT HYPs predicted in the GSE18800 Bleomycin data set.
HYPs predicted by RCR in the GSE18800 Bleomycin data set that were present in the Fibrosis and EMT subnetwork are listed in the right column with the predicted direction. Yellow indicates increased HYPs, blue indicates decreased HYPs. Direction of HYP indicating increased fibrosis is indicated in the left column.
Figure 3: A subset of fibrosis pathways predicted in bleomycin and mechanical injury data sets.
HYPs predicted in the bleomycin data set (left) and mechanical injury data set (right) painted on a subset of the Fibrosis and EMT subnetwork. Yellow indicates increased HYP, blue indicates decreased HYP, and red indicates increased gene expression State Change. See Supplementary Files 6 and 7 for complete subnetwork painted with HYPs predicted in GSE18800 and GSE5372, respectively.
In addition to the major fibrotic pathways described above, increased Hedgehog signaling, which regulates EMT, was supported by HYPs for increased Shh and downstream effector Gli1 transcriptional activity, as well as a HYP for decreased Ptch1, a negative regulator of the Hedgehog pathway. Finally, the activity of multiple transcription factors known to regulate extracellular matrix gene expression were predicted to be increased, including Egr1, Sp1, Ets2, Yy1 and Ap-1 Complex. These transcription factors represent final downstream elements of multiple signaling pathways resulting in the regulation of Extracellular Matrix (ECM) gene expression. The prediction of increased activity and abundances of these transcription factors along with the measurement of increased ECM gene expression indicate that both upstream and downstream signaling components that regulate fibrosis are activated with bleomycin (Figure 3).
GSE5372 Mechanical Injury data set: 23% of the RCR-capable nodes in the Fibrosis and EMT subnetwork were predicted in the mechanical injury human data set. Angiotensin II, PI3K/AKT and beta-catenin pathways were predicted by RCR and have been shown to regulate wound healing processes [38-40] (Figures 3 and 4). Angiotensin signaling has been shown to play a role in ventilator-induced lung injury, with increased angiotensin II levels. Additionally, decreased lung injury upon ACE inhibition or angiotensin receptor antagonist treatment has been shown [41-43]. AKT has been shown to contribute toepithelial cell wound-healing response and ventilator-induced lung injury [44,45]. Mechanical ventilation increased β-catenin protein expression in rat lung . In addition, PPARG was predicted to be increased in the mechanical injury data set, which has also been shown to be activated during wound healing . TGFB1 was predicted, but there is no evidence of TGFB receptor or SMAD activity. Finally, many transcription factors that regulate extracellular matrix gene expression were predicted to be increased such as the transcriptional activities of SP1, ETS2, YY1 and AP-1 and the protein abundance of EGR1, as well as increased ECM gene expression, indicating that not only upstream but also downstream signaling components that contribute to fibrosis are activated with mechanical injury.
Figure 4: Fibrosis and EMT HYPs predicted in the GSE5372 Mechanical injury data set.
HYPs predicted by RCR in the GSE5372 Mechanical injury data set that were present in the Fibrosis and EMT subnetwork are listed in the right column with the predicted direction. Yellow indicates increased HYPs, blue indicates decreased HYPs. Direction of HYP indicating increased fibrosis is indicated in the left column.
Cell Migration and Adhesion Pathways Predicted By RCR
GSE18800 Bleomycin data set: In the Bleomycin data set, RCR predicted 34% of the RCR-capable nodes in the Cell Migration and Adhesion subnetwork, indicating evidence for cytoskeletal, integrin and MMP pathways. Pi3k activates Rhoa, Rac1 and Ilk [48-50], key regulators of cytoskeletal rearrangement and cell adhesion important for migration, which were predicted as HYPs (Figures 5 and 6). Ilk increases actin cytoskeletal organization and can promote focal adhesion formation to promote cell attachment [51,52]. Mmp activation was also predicted in response to bleomycin with increased predictions in Mmp3 and Adam17. Mmp and Adam proteins promote migration by matrix degradation and shedding to allow cell movement. Mmp3 is a mediator of pulmonary fibrosis, and Mmp3-null mice were protected against bleomycin-induced pulmonary fibrosis . Bleomycin is also known to induce Adam17 activity and protein expression in mouse lung , and ADAM17 promotes FGF7/FGFR2b-dependent cell migration through the EGFR/ERK1/2 pathway in human epithelial and endothelial cells . Fgf7, Fgf2 and Erk Family activity were also predicted to be increased in the bleomycin data set. In addition to Fgf2 and Fgf7 growth factors, Pdgf was predicted to be increased. These growth factors can stimulate PI3K and MAPK to activate cell migration and adhesion, which were also predicted in the data set [56,57].
Figure 5: Cell migration and adhesion HYPs predicted in the GSE18800 Bleomycin data set.
HYPs predicted by RCR in the GSE18800 Bleomycin data set that were present in the Cell Migration and Adhesion subnetwork are listed in the right column with the predicted direction. Yellow indicates increased HYPs, blue indicates decreased HYPs. Direction of HYP indicating increased migration is indicated in the left column.
Figure 6: A subset of cell migration and adhesion pathways predicted in bleomycin and mechanical injury data sets.
HYPs predicted in the bleomycin data set (left) and mechanical injury data set (right) painted on a subset of the Cell Migration and Adhesion subnetwork. Yellow indicates increased HYP, blue indicates decreased HYP. See Supplementary Files 8 and 9 for complete subnetwork painted with HYPs predicted in GSE18800 and GSE5372, respectively.
GSE5372 Mechanical injury data set: In the mechanical injury data set, RCR predicted 21% of the RCR-capable nodes in the Cell migration and adhesion subnetwork (Figure 7). Many of the same cell migration pathways predicted in the bleomycin data set were predicted in the mechanical injury data set, including increased FGF2, PI3K, AKT and RHOA activities as well as ECM proteins collagen I and FN (Figure 6). In addition to RHOA, downstream effectors ROCK1 was predicted to be increased, which regulates actin polymerization essential for adhesion and migration. Growth factor activin A was also predicted to be increased, which can stimulate migration through phosphorylation of focal adhesion proteins .
Figure 7: Cell Migration and Adhesion HYPs predicted in the GSE5372 Mechanical injury data set.
HYPs predicted by RCR in the GSE5372 Mechanical injury data set that were present in the Cell Migration and Adhesion subnetwork are listed in the right column with the predicted direction. Yellow indicates increased HYPs, blue indicates decreased HYPs. Direction of HYP indicating increased fibrosis is indicated in the left column.
Immune regulation of tissue repair pathways predicted by RCR
GSE18800 Bleomycin data set: In the bleomycin evaluation data set, 48% of the RCR-capable nodes of the Immune Regulation of Tissue Repair subnetwork were predicted as HYPs in a direction consistent with an increased immune response (Table 3). The major immune pathways identified included IL6/STAT3, TNF/NFKB, TH2 cytokines, PPARα, and macrophage activation (Figures 8 and 9).
Figure 8: Immune HYPs predicted in the GSE18800 Bleomycin data set.
HYPs predicted by RCR in the GSE18800 Bleomycin data set that were present in the Immune Regulation of Tissue Repair subnetwork are listed in the right column with the predicted direction. Yellow indicates increased HYPs, blue indicates decreased HYPs. Direction of HYP indicating increased migration is indicated in the left column.
Figure 9: A subset of fibrosis pathways predicted in bleomycin and mechanical injury data sets.
HYPs predicted in the bleomycin data set (left) and mechanical injury data set (right) painted on a subset of the Immune Regulation of Tissue Repair subnetwork. Yellow indicates increased HYP, blue indicates decreased HYP, and red indicates increased gene expression State Change. See Supplementary Files 10 and 11 for complete subnetwork painted with HYPs predicted in GSE18800 and GSE5372, respectively.
The IL6/STAT3 pathway includes IL6 which activates the transcriptional activity of STAT3, which is inhibited by SOCS3. The IL6/ STAT3 pathway has been shown to stimulate collagen transcription and accumulation in bleomycin-treated mice. An increase in Stat3 activity mediated by Il6 was associated with the presence of inflammatory infiltrates in the lung, and a more severe bleomycin-induced fibrotic response .
TNF/NFkB is a key regulator of the immune response, including inflammation . NFkB activity was predicted in the bleomycin data set, which is activated by Tlr4 and Tnf through the activity of its receptor, Tnfrsf1a, also predicted HYPs. Bleomycin is known to induce NFkB activity shortly after treatment of endothelial cell cultures, and this activity continues for up to 21 days following in vivo treatment of mice [61-63]. TNFα is known to have a central role in the regulation of lung inflammation and fibrosis .
Previous work has shown that bleomycin administration increases the mRNA expression of TH2-associated cytokines (Il4, Il10, and Il13) [65,66]. Increase of Il5 and Il13 was predicted in the bleomycin treatment data set are consistent with the literature reported results.
The PPARα pathwaywas predicted to be decreased in the bleomycin data set, consisting of the following predictions: the transcriptional activity of Ppara, Ppara protein abundance and the transcriptional activity of Ncor1. In a study comparing bleomycin-induced fibrosis between WT and PPARα KO mice, the KO mice showed enhanced fibrosis and Tnf and Il1b activity, indicating that PPARα can regulate these inflammatory cytokines in an anti-inflammatory context , consistent with the increased predictions in Tnf and Il1b in this data set.
Macrophage activation, and the catalytic activities of thrombin (F2), Csf1, Csf2 and Csf3 and transcription factor PU.1 (Spi1) were predicted as HYPs in this data set, reflecting the increased inflammatory activity of macrophages and macrophage-associated growth factors. The predicted HYPs play key roles in macrophage activity and differentiation. The continuous expression of Spi1 on myeloid precursor cells commits the cell to macrophage differentiation . Thrombin induces the chemotaxis and the proliferation of macrophages , and the Csf proteins are involved in macrophage or neutrophil proliferation and chemotaxis . These results are in line with increased bronchoalveolar lavage macrophage numbers measured in the GSE18800 study .
GSE5372 Mechanical injury data set: In the mechanical injury evaluation data set, 15% of the RCR-capable nodes were predicted as HYPs. As in the bleomycin data set, the pathways in the Immune Regulation of Tissue Repair subnetwork were largely predicted as increased but with many fewer HYPs. TNF, IL6, STAT3 and TH2 cytokine IL4 were predicted in the mechanical injury data set, along with cytokines that point to macrophage activation (Figure 10). The differences in the predictions between the two data sets may be due to biological differences between bleomycin and mechanical injury, as explained in the Discussion.
Figure 10: Immune HYPs predicted in the GSE5372 Mechanical injury data set.
HYPs predicted by RCR in the GSE5372 Mechanical injury data set that were present in the Immune Regulation of Tissue Repair subnetwork are listed in the right column with the predicted direction. Yellow indicates increased HYPs, blue indicates decreased HYPs. Direction of HYP indicating increased fibrosis is indicated in the left column.
The TRAG Network builds upon our previous work aimed at constructing network models for lung-relevant processes through the development of an additional network describing the molecular pathways of tissue repair and angiogenesis. These networks contain causal relationships displayed in BEL, a biological expression language that represents life science relationships in a computable format. This enables quantitative network perturbation amplitude scoring and thereby to compare different toxicants for their biological impact [1,2]. Additionally, the TRAG Network, which covers non-cancerous angiogenesis and reversible tissue repair processes, addresses the initial steps leading to disease onset. Thereby the network lays the foundation to construct networks that describe the onset of lung diseases which involve irregular wound healing and irreversible tissue remodeling.
Tissue repair subnetworks were evaluated for their ability to reflect known biological signaling by using two transcriptomic data sets from experiments where wound healing leading to tissue repair was induced following experimental insult. Evaluation of the angiogenesis subnetworks within the TRAG Network was not performed due to the lack of data sets with measured angiogenic phenotypes. Evaluation of angiogenesis subnetworks can be undertaken in the future when transcriptomic profiling data sets fulfilling our acceptance criteria are available. The evaluation of the tissue repair subnetworks confirmed that the subnetworks contain pathways predicted by RCR that are involved in fibrosis and EMT, immune regulation of tissue repair, and migration and cell adhesion in wound healing. The subnetworks reflect a priori expectations for tissue repair based on previous literature reports as well as measured phenotypic endpoints in the data sets.
GSE18800 Bleomycin data set validates relevance of tissue repair subnetworks
Bleomycin is a chemotherapeutic agent, the use of which can induce pulmonary fibrosis in patients. Due to the ability of bleomycin to induce fibrosis, administration of bleomycin to rodents has become an extensively studied model of lung fibrosis, fibroblast activation and EMT transition [71,72]. The GSE18800 bleomycinstudy used for network evaluation measured increased Tgf-betapathway activation and fibrotic endpoints . RCR-predicted HYPs present in the Fibrosis and EMT subnetwork includeTgf-beta, angiotensin II Pi3k/Akt, betacatenin and Hedgehog pathways, many of which are known to be induced by bleomycin. Increased ECM production was also predicted including fibronectin and Ctgf as well as multiple transcription factors that regulate ECM genes. In addition to predicting bleomycin- induced pathways described in the literature, RCR predicted pathways not previously studied in the context of bleomycin-induced injury. Wnt3a, an activator of beta-catenin, and Dkk1, inhibitor of Wnt3a, were predicted to be increased and decreased, respectively. Although Wnt3a is known to stimulate beta-catenin in mouse fibroblasts [73,74] and Dkk1 overexpression ameliorates skin fibrosis , Wnt3a and Dkk1 have not been shown to be involved specifically in bleomycininduced lung fibrosis. Therefore, RCR predictions can potentially offer more mechanistic detail to explain bleomycin-induced beta-catenin signaling. Additional experimental investigation is a future next step to confirm this prediction.
In addition to fibrosis, bleomycin is known to induce the migration of fibroblasts and inflammatory cells [76,77]. Within the Cell Migration and Adhesion in Wound Healing subnetwork, Mapk, Pi3k, and Rho kinase family member activities were predicted to be increased in the bleomycin data set. Although these pathways are known to regulate cytoskeletal rearrangements important in cell migration and Mmp signaling [78-80], they are also involved in many other biological processes. Mmp, Adam and integrin pathways give additional support for the activation of cell migration pathways [55,81], and were also predicted in the bleomycin data set. Mmp3 and Adam17 in particular are known from the literature to be increased with bleomycin [53,54]. Some HYPs were predicted in directions inconsistent with increased migration (decreased Itgb6, Adam10 and Mmp12). This discrepancy is likely due to a lack of suitable existing context-specific knowledge of the genes downstream of these HYPs in the public domain. MMPs and integrins have been largely studied in the context of cancer, which can reflect different biology than non-cancerous tissue, and could be one reason for the decreased predictions.
Bleomycin has been shown to induce an initial inflammatory response that gradually decreases with time, during which the initiation of gene expression associated with fibrosis between days 9 and 14 in rodents occurs . The bleomycin gene expression data sets evaluated in this study were obtained at day 14 post-exposure, and inflammatory endpoints measured in the study related to the data set to monitor the effects of bleomycin on the mouse model support the prediction of many inflammatory HYPs identified in the present work . These HYPs include Il6, Stat3, Ifng, Tnf, NFkB and many inflammatory cytokines. An increase in the total number of cells obtained by bronchoalveolar lavage in this study was seen at days 7 and 21 postexposure. Macrophage numbers in the lavage fluid were also increased at these time points. Alveolar macrophages have a key role in lung inflammation and resolution , and the activation of macrophages, indicated by associated growth factors Csf1 and Csf2 predicted in the bleomycin data set, supports this role. RCR successfully predicts mechanisms in that reflect measured phenotypes, however a caveat is that bleomycin phenotypes were measured on day 7 and day 21 while the transcriptomic data is from day 14. Therefore, we are making the assumption that the signaling was the same for day 14 as it was on day 7 and day 21.
GSE5372 mechanical Injury data set validates relevance of tissue repair subnetworks
In the GSE5372 data set, mechanical injury was performed by bronchial brushing; the gene expression data used for analysis in our study was obtained at day 7 following brushing. The GSE5372 study found that the injured area was completely covered by a partially redifferentiated epithelial layer after day 7, the expression data was dominated by cell cycle genes, and less than 1% of the sample consisted of inflammatory cells . This suggests that overt inflammation had largely subsided by the day 7 time point , and that the epithelial cells had migrated to the wound site and were actively proliferating to cover the wound. These measured endpoints indicate that at day 7, the mechanical injury model includes a wound healing component that has passed the initial inflammatory stages. Support for this view comes from the predicted increase in IL4, decrease in IFNG and increase in anti-inflammatory PPARs, PPARA and PPARG. IL4 has been identified as being important in wound closure  and suppressing an inflammatory response , and mechanical injury has been shown to induce increased IL4 and less IFNG production . In addition, PPARA and PPARG have a role in wound healing, contributing to the resolution of inflammation after injury [87,88]. PI3K and RHOA were predicted as HYPs within the Cell Migration and Adhesion subnetworkfor the mechanical injury data set; however, these molecules are also involved in fibrosis in addition to other biological processes, and additional factors related to cell migration (ADAM, integrin and MMP signaling) were not predicted as HYPS in the mechanical injury data set, in contrast to the bleomycin dataset. This indicates that the cell migration mechanisms at day 7 have subsided, and wound healing processes to complete wound closure such as cell proliferation and secretion of ECM may be taking place. Indeed, an increase in Collagen I is predicted along with wound healing mechanisms in the Fibrosis and EMT subnetwork leading to fibroblast activation and ECM secretion, indicated by predicted increase of angiotensin II, PI3K/AKT and betacatenin. RCR predicts wound healing mechanisms 7 days after wound injury, which is in line with a partially differentiated epithelial layer measured on day 7 indicating wound healing. In addition, inflammatory mechanisms are not strongly predicted by RCR which is in line with the small amount of inflammatory cells measured in the study on day 7. The lack of signaling however is not as convincing as having a positive measurement of inflammation perhaps 1 day after wounding and gene expression collected to confirm that RCR could successfully predict the presence of inflammation shortly after wounding.
Differences between evaluation data set predictions validates the applicability of the TRAG Network to multiple tissue injury data sets
RCR was able to predict common elements in tissue repair pathways in bleomycin and mechanical injury. RCR could also distinguish signaling that was unique to each perturbation, highlighting the ability of the subnetworks to apply to different types of tissue injury contexts. Bleomycininjury displayed a strong fibrotic and inflammatory signal with some evidence for cell migration, whereas mechanical injury had induced wound healing pathways, indicating that inflammation and cell migration had subsided by day 7. These predictions matched the measured endpoints in both studies.
RCR predicted shared fibrosis pathways including angiotensin, beta-catenin and PI3K and various transcription factors that regulate ECM genes. Increased collagen and fibronectin gene expression also indicated that ECM secretion was a component of the wound healing processes in both data sets. Tgfb, a major regulator of ECM secretion and fibrosis, was predicted in the bleomycin data set, matching literature expectations and measured endpoints. However, the mechanical injury data set had little support for TGFB/SMAD pathway activation, and predicted to be increased rather than decreased PPARG. Mechanical injury via bronchial brushing would not be expected to induce fibrosis, but an activation of pathways leading to extracellular matrix production for wound healing may occur. Increased PPARG is predicted in response to mechanical ventilation, in agreement with the absence of a strong prediction in TGFB/SMAD signaling which would inhibit PPARG . It has been shown that PPARG inhibitor (GW9662)- treated animals have significant degree of lung injury following 8 hours of mechanical ventilation, indicating a role for PPARG in normal wound healing . PPARG’s role in normal wound healing processes is also illustrated by the observation of increased PPARG expression in neointimal hyperplasia after balloon injury and in cultured smooth muscle cells after scraping injury . Therefore, the difference in PPARG and TGFB pathways between bleomycin and mechanical injury may be explained by the differences in the bleomycin-induced TGFBMechanism driven fibrotic response versus a less fibrotic wound healing response resulting in ECM secretion to repair the damaged tissue.
HYPs that are involved in cell migration such as ADAMs and MMPs were predicted in the bleomycin data set but not the mechanical injury data set. This finding is in line with measured endpoints in the mechanical injury data set indicating that by day 7, the epithelium was fully covered and cells had completed migration into the wound site.
Finally, the bleomycin data set had an increased immune signal compared to the mechanical injury data set, in line with measured endpoints in each of the data sets. NFkB, STAT1 and STAT6 activities along with multiple cytokines were predicted in the bleomycin data set, but not predicted in the mechanical injury data set. NFkB activity, a critical driver of inflammation induced by bleomycin, has been shown to be sustained for 14 days following treatment [91,92]. In contrast, mechanical injury has been shown to induce NFkB activity within minutes of injury [93,94], and in wound closure (as seen in the mechanical injury data set), NFkB expression occurs at early stages along the margins of the wound . By seven days, the mechanical injury study indicated that wound closure was completed and as a result, NFkB expression would no longer be increased. This supports the lack of a strong immune response prediction in the mechanical injury data set and the view that this inflammatory stage had passed, unlike in the bleomycin study.
These two evaluation data sets illustrate the application of the tissue repair subnetworks to multiple models of injury and highlight the ability of RCR to predict HYPs unique to these perturbations. Additionally, the TRAG Network describes non-disease, reversible tissue repair pathways that, if disturbed, may lead to disease initiation, and could be used as a substrate to bridge toward construction of networks describing lungrelevant diseases, providing an additional tool to help in the elucidation of disease-driving mechanisms.
The objective of this study was to build and evaluate a network model describing molecular mechanisms of pulmonary tissue repair and angiogenesis in order to understand how these molecular events can lead to disease initiation, and to use this model with our computational methodologies to quantify the biological impact of different types of exposures. We built the TRAG network, composed of nine subnetworks, and evaluated the tissue repair subnetworks with two independent data sets using RCR, a computational approach that predicts upstream controllers of gene expression changes observed in a data set. Results from RCR analysis of these tissue repair data sets validated the content of the tissue repair subnetworks, confirming that they contained tissue repair pathways relevant to two types of wounding perturbations. These results indicate that the tissue repair subnetworks can be used to assess the impact of various types of tissue damage, and the signaling perturbations during the repair processes. Moreover, unique signaling events can be predicted that reflect phenotypic endpoints measured in the study.
We would like to acknowledge Michael J. Maria for project management and support with preparation of this manuscript.
Selventa and PMI authors performed this work under a joint research collaboration funded by PMI.