Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Developing and applying a gene functional association network for anti-angiogenic kinase inhibitor activity assessment in an angiogenesis co-culture model

Developing and applying a gene functional association network for anti-angiogenic kinase... Background: Tumor angiogenesis is a highly regulated process involving intercellular communication as well as the interactions of multiple downstream signal transduction pathways. Disrupting one or even a few angiogenesis pathways is often insufficient to achieve sustained therapeutic benefits due to the complexity of angiogenesis. Targeting multiple angiogenic pathways has been increasingly recognized as a viable strategy. However, translation of the polypharmacology of a given compound to its antiangiogenic efficacy remains a major technical challenge. Developing a global functional association network among angiogenesis-related genes is much needed to facilitate holistic understanding of angiogenesis and to aid the development of more effective anti-angiogenesis therapeutics. Results: We constructed a comprehensive gene functional association network or interactome by transcript profiling an in vitro angiogenesis model, in which human umbilical vein endothelial cells (HUVECs) formed capillary structures when co-cultured with normal human dermal fibroblasts (NHDFs). HUVEC competence and NHDF supportiveness of cord formation were found to be highly cell-passage dependent. An enrichment test of Biological Processes (BP) of differentially expressed genes (DEG) revealed that angiogenesis related BP categories significantly changed with cell passages. Built upon 2012 DEGs identified from two microarray studies, the resulting interactome captured 17226 functional gene associations and displayed characteristics of a scale-free network. The interactome includes the involvement of oncogenes and tumor suppressor genes in angiogenesis. We developed a network walking algorithm to extract connectivity information from the interactome and applied it to simulate the level of network perturbation by three multi-targeted anti-angiogenic kinase inhibitors. Simulated network perturbation correlated with observed anti- angiogenesis activity in a cord formation bioassay. Conclusion: We established a comprehensive gene functional association network to model in vitro angiogenesis regulation. The present study provided a proof-of-concept pilot of applying network perturbation analysis to drug phenotypic activity assessment. Page 1 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 has also been noticed that HUVECs change their cellular Background Angiogenesis, the generation of new blood vessels, plays and molecular properties upon passage in vitro, a phe- an essential role under normal physiological conditions nomenon thought to be due to in vitro cellular senescence as well as during the pathogenesis of many diseases [17-20]. In the present study, we determined that HUVEC including atherosclerosis, macular degeneration, wound competence and NHDF supportiveness for angiogenesis healing, diabetic retinopathy, and human malignancy [1- in this co-culture system are both cell passage dependent. 3]. Remarkably, tumor dormancy is believed to be attrib- Gene Ontology (GO) analysis of differentially expressed uted, at least in part, to the lack of angiogenesis support. genes showed that the cell passage dependent global tran- The transition from an avascular, dormant tumor to an scriptional changes are highly related to angiogenesis. aggressively growing angiogenic cancer is referred to as the This enabled us to construct a comprehensive functional "angiogenic switch" [4,5]. More than 30 years ago, it was association network of differentially expressed genes hypothesized that inhibition of tumor angiogenesis using a natural language processing algorithm. We further would inhibit solid tumor growth [6]. Since then, the developed a "network walking" algorithm to estimate net- expansion of angiogenesis research has resulted in the work perturbation by small molecule kinase inhibitors. identification of various pro- and anti-angiogenic factors, The simulated compound activity via network perturba- such as FGF, VEGF, angiopoietin, endostatin, vasostatin, tion analysis was in good agreement with actual pheno- and neuronal cell axon guidance molecules [2,7-9] and typic activity in the cord formation bioassay. the development of several anti-tumor angiogenesis med- icines that have recently proven efficacious in the clinic Results [1,2,7,10,11]. Characterization of an in vitro angiogenesis co-culture model by multi parameter high content image analysis Malignant tumors often express an array of angiogenic HUVECs form a capillary structure when co-cultured with factors to potentiate angiogenesis and tumor growth [8]. NHDF in vitro [16]. Unlike the HUVEC Matrigel™ assay, Ultimate angiogenic outcomes depend on the dynamic the HUVEC/NHDF cord formation co-culture allows us to equilibrium between positive and negative regulators and study cell-cell interactions during long term (up to 14 the interplay among their signal transduction pathways, days) angiogenesis. The formation of the cord structure is not on a single discrete pathway. Avastin, a monoclonal a highly organized morphogenic process. We developed a antibody that specifically blocks VEGF, is the first high content imaging based assay to document and quan- approved anti-angiogenic therapy. Notably, combination titate cord formation using the Cellomics ArrayScan plat- with chemotherapies has been necessary for its clinical form. We were able to simultaneously assess dozens of efficacy. Moreover, tumors often develop drug resistance morphological measurements, such as tube length, area, in response to anti-VEGF therapy [12]. Hanahan et al connectivity of the cord, branching, and number of nuclei demonstrated one possible mechanism by which such within the cord. Using this assay system, we tested a clin- drug resistance might develop[13]. Inhibition of VEGF ically approved antiangiogenic compound Sunitinib signalling by neutralizing antibody to VEGFR2 (KDR) (SUTENT, SU11248). While 20 ng/ml VEGF significantly induces elevated expression of hypoxia associated proan- increased angiogenesis, 30 nM Sunitinib significantly giogenic factors such as FGF and EphrinA1, which subse- inhibited angiogenesis (Fig. 1A). quently reactivates VEGF independent angiogenesis and tumor growth. Thus, disrupting a single proangiogenic To explore the phenotypic signatures induced by VEGF or pathway by itself is often insufficient to achieve sustained Sunitinib, we performed principal component analysis therapeutic benefits. In this light, it is necessary to explore (PCA) and hierarchical clustering (Fig. 1B &1C). Multivar- global functional association among angiogenesis-related iate analyses using 47 morphological measurements genes rather than focusing on an individual or a few ang- clearly distinguished the three treatment groups, i.e. con- iogenesis factors discretely. trol, VEGF, and Sunitinib. In the PCA analysis, the first three principal components captured more than 80% of Importantly, angiogenesis involves intercellular interac- the total variance. Examining the PC loadings, we found tion among vessel-forming endothelial cells, nonmalig- that among the most discriminative features are the cord nant cells such as the supporting pericytes, immune and length (either connected or unconnected or in general), stromal cells, as well as the malignant tumor cells [1,14]. cord area, the branching segment count. The cord CD31 Therefore, global gene-gene interactions during angiogen- immunofluorescence intensity is also one of the features esis need to be explored in a multi-cell type context. that significantly contributed to PC1, likely due to the Human umbilical vein endothelial cells (HUVECs) are overall correlation of fluorescence staining to cord mass. widely used to study vascular biology [15]. They form a Interestingly, the cord width is not a discriminative fea- lumen bearing capillary structure when co-cultured with ture for either VEGF or Sunitinib action. The current normal human dermal fibroblast (NHDF) cells [16]. It multi-parameter, high content analysis allowed for a Page 2 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 Multi p Figure 1 arameter characterization of in vitro angiogenesis in a co-culture model Multi parameter characterization of in vitro angiogenesis in a co-culture model. (A) HUVECs formed capillary structure (cord) when plated on dermal fibroblast cells over the period of 12–14 days. The angiogenic cord structure was vis- ualized by human CD31 staining (green), fibroblast cells were stained with an intermediate filament marker vimentin (red), while the nucleus of both cell types were visualized with Hoechst staining (blue). Cord formation was greatly enhanced by addi- tion of 20 ng/ml VEGF (VEGF) to the culture, while 30 nM Sunitinib (CMPD) drastically inhibited cord formation as compared to the medium control (CTL). (B). Morphological features of the cord formed in the co-culture were captured and quantitated by high content image analysis using the ArrayScan platform adopting the tube formation bioapplication. 47 measurements were assessed, normalized using the Z-score across different treatment groups and analyzed by Principal Component Analysis (PCA). The first three principal components accounting for more than 80% of the variance were selected for sample scatter plotting. Blue are the control cultures (CTL); yellow, cultures treated with 20 ng/ml VEGF; and red, cultures treated with 30 nM Sunitinib (CMPD) treated samples. (C). The Z-scores of 47 morphological measurements were clustered in Spotfire Deci- sion Site with the correlation coefficient being used as the distance metric. Complete linkage was used as the clustering algo- rithm. Red stands for higher and green for lower than mean (black) expression. (D). Six representative parameters were selected and plotted. Except for the cord width, the cord length, area, branching point, the number of nucleus within the cord, and a composite angiogenesis index (Angio Index) all manifested the anti and pro angiogenic effects from Sunitinib and VEGF, respectively. Data are expressed as Mean ± S.D. detailed dissection of the phenotypic characteristics asso- HUVECs in the co-culture. Sunitinib treated cultures are ciated with a given factor or compound. Shown in Fig. 1D less morphologically complex, as the number of con- is the mean value plot of a few representative parameters. nected cords, the branching nodes, and segment measure- Sunitinib significantly inhibited the mass of the cord in ments are all significantly reduced. length and area, but not width. Moreover, the number of nuclei per tube is significantly reduced, indicating that Sunitinib could have anti-proliferation activity on Page 3 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 HUVEC cord formation competence is cell passage various passages (P6, P14, P22) and tested them for cord dependent formation competence. We found that they gradually lost It has been well known that HUVECs change their bio- cord formation capabilities upon passage (Fig. 2A). While chemical and cellular behaviors when passaged in vitro early passage (P6) HUVECs formed very nice capillary [17-20]. We cultured HUVEC cells from a single batch to structures and responded robustly to VEGF when co-cul- H Figure 2 UVEC passage dependent angiogenic competence in cord formation assay HUVEC passage dependent angiogenic competence in cord formation assay. HUVECs cultured to various passage numbers (P6, early; P14, intermediate; and P22, late) were co-cultured with NHDFs in the cord formation assay. The cord was visualized by human CD31 staining (green) and quantitated. The cultures were treated with or without 20 ng/ml exogenously added VEGF. Representative images are shown in (A). 6 well-level parameters were quantitated and plotted in (B). Please note that HUVECs gradually lost their cord formation competence when passaged in vitro. Data are expressed as Mean ± S.D. (C- D). HUVECs with various passage numbers (P5, P8, P12, and P22) were assayed for basal gene expression using the Affymetrix GeneChip platform. A linear regression method was used to identify differentially expressed genes over cell passage. An enrichment testing was used to identify significantly changed Biological Process (BP) categories associated with HUVEC pas- sage-dependent gene expression. Representative BP categories were shown using heatmap visualization, for genes associated with mitotic cell cycle (C), and cell death (D). Page 4 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 tured with NHDF, late passage (P22) HUVECs lost cord Profiling passage dependent angiogenesis support by formation capability and had diminished response to dermal fibroblast cells Long term sustained angiogenesis in co-culture is sup- VEGF (Fig. 2A &2B). ported by dermal fibroblast feeder layer cells, which are We next profiled gene expression changes associated with also affected by cell passage. While late passage (P12) HUVEC passage. Since the loss of angiogenesis capability adult NHDFs demonstrated poor support for HUVECs to was proportional to the round of cell passage, we applied form capillary structures, early passage (P5) adult cells did a statistical approach to identify genes with regressive well (Fig. 3A), suggesting that supportiveness of NHDFs expression patterns, i.e. their expression either increased in the angiogenesis assay is indeed cell-passage depend- or decreased over passage numbers. We identified 1103 ent. An early passage NHDF from a neonatal donor also differentially expressed (DE) genes which are listed in supported angiogenesis in this assay (Fig. 3A). Supplemental Table 1 (see Additional file 1). In order to understand the biological relevance of the DE genes, we Differential angiogenesis supportiveness by NHDFs in the analyzed the DE genes by using the DAVID functional co-culture system represents a unique opportunity for us annotation system [21]. Table 1 lists Gene Ontology (GO) to dissect factors that may contribute to the complex inter- terms significantly enriched with the DE genes. The most cellular communications during angiogenesis. To assess significant Biological Processes (BP) categories are cell this, we profiled gene expression in the three listed prepa- cycle, proliferation and programmed cell death. Expres- rations of NHDF cells by DNA microarray analysis. sion profiles of the DE genes involved in representative Between the supportive NHDF cells versus the non-sup- biological processes were analyzed by hierarchical cluster- portive ones, we identified 787 non-redundant common ing analysis. The majority of the DE genes involved in DE genes which are listed in Supplemental Table 2 (see mitotic cell cycle and proliferation were down-regulated, Additional file 2). Similarly, we used the DAVID func- whereas the majority of genes related to cell death were tional annotation system to examine the functional rele- up-regulated (Fig. 2C–D), which supports the replicative vance to angiogenesis of these DEGs. Table 2 lists the cell senescence hypothesis. significantly enriched GO terms. Interestingly, blood ves- Table 1: Gene Ontology (GO) terms enriched with differentially expressed genes in HUVEC of different passages 1 2 3 GO GO Term DEG FDR Biological Process cell cycle 93 0.00075 mitotic cell cycle 44 0.00042 M phase of mitotic cell cycle 36 0.00124 programmed cell death 44 0.00717 mitosis 35 0.00511 cell division 37 0.00332 regulation of cell cycle 60 0.01370 cell cycle checkpoint 13 0.01512 cell organization and biogenesis 113 0.05000 cell proliferation 51 0.01164 Molecular Function protein binding 270 0.00051 ATP binding 106 0.00048 nucleotide binding 137 0.00473 catalytic activity 309 0.00144 purine nucleotide binding 119 0.00210 cytoskeletal protein binding 33 0.00545 kinase activity 70 0.03651 binding 545 0.03367 transferase activity 116 0.01551 Cellular Component spindle 14 0.00058 microtubule cytoskeleton 36 0.00182 intracellular 439 0.00138 cytoskeleton 73 0.00558 nucleus 237 0.00615 intracellular membrane-bound organelle 317 0.04183 intracellular organelle 364 0.01793 1. Gene Ontology; 2. Differentially expressed genes; 3. False discovery rate Listed are GO categories identified through enrichment testing of the differentially expressed genes from HUVEC using the DAVID functional annotation system as described in the methods section. Page 5 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 NHDF pa Figure 3 ssage dependent angiogenic supportiveness in cord formation assay NHDF passage dependent angiogenic supportiveness in cord formation assay. P6 NHDFs from a neonatal donor, P5 and P12 NHDFs from an adult donor were tested for their angiogenesis supportiveness in the cord formation assay (A). Cultures were performed either with or without 20 ng/ml exogenously added VEGF. Cords were visualized by human CD31 staining (green). Both early passage neonatal and adult NHDF cells are more supportive for cord formation than the late pas- sage adult NHDFs, with the late passage NHDF cells being non-permissive. (B). Clustering analysis of genes involved in angio- genesis. P6 NHDFs from a neonatal donor (Neo P6), P4 NHDFs from an adult donor (Ad P4), and P13 NHDFs (Ad P13) from the same adult donor were assayed for basal gene expression using the Affymetrix GeneChip platform. 787 nonredundant genes were identified as common differentially expressed genes between the two early passage NHDF and the late passage NHDF cells. Enrichment testing was used to identify significant Biological Process (BP) categories. Representative angiogenesis BP categories are shown in a heatmap visualization. sel development and morphogenesis, angiogenesis and both HUVEC and NHDF cells influences angiogenesis vasculature development were among the most signifi- supports modeling the angiogenesis network using the cantly enriched biological processes, suggesting that these DEGs identified by these DNA microarray experiments. molecular changes were related to cellular potential to To understand the functional association of these genes, support angiogenesis. Representative expression profiles we constructed a molecular interactome among the total of genes involved in blood vessel development and angio- 2012 differentially expressed genes identified from either genesis are shown in Fig. 3B. This result is consistent with NHDFs or HUVECs. The interactome consisted of over the supporting roles of NHDFs in this co-culture system. 23,000 gene-to-gene connections retrieved from PubMed Taken together, passaging NHDF cells in vitro adversely abstracts using the natural language processing algorithm affects cellular capacity to support angiogenesis due to sig- in PathwayAssist. We then refined the interactome by nificantly orchestrated changes of many genes involved in removing genes whose expression was not detectable by a few related biological processes. microarray analysis in the co-culture system. The resulted molecular interactome consists of 17226 molecular rela- Establishment of a gene functional association network- tionships described by the following terms: binding, The fact that cell passage dependent gene expression in expression, molecular transport, protein modification Page 6 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 Table 2: Gene Ontology terms enriched with differentially expressed genes in NHDF cells of different passages 1 2 3 GO GO Term DEG FDR Biological Process blood vessel development 64 0.00700 vasculature development 64 0.00700 actin filament organization 40 0.00850 actin cytoskeleton organization and biogenesis 81 0.00850 actin filament-based process 84 0.00850 regulation of signal transduction 98 0.01020 angiogenesis 56 0.01140 blood vessel morphogenesis 60 0.01140 cell morphogenesis 161 0.02020 Cellular Component extracellular matrix 75 0.00200 proteinaceous extracellular matrix 74 0.00200 actin cytoskeleton 60 0.04800 adherens junction 30 0.06140 cell-matrix junction 25 0.06140 cytoskeleton 132 0.06140 Molecular Function protein binding 864 0.00820 binding 1195 0.01740 cytoskeletal protein binding 95 0.02380 1. Gene Ontology; 2. Differentially expressed genes; 3. False discovery rate. Listed are GO categories identified through enrichment testing of the differentially expressed genes from NHDF cells using the DAVID functional annotation system as described in the methods section. and regulation (see Additional file 3). There were 1201 of connectivity in the interactome. This observation reca- unique genes in the interactome, among which 887 genes pitulates the emerging role of oncogenes and tumor sup- were differentially expressed in the co-culture system. pressor genes in angiogenesis regulation [1,5,23,24]. We Thus, the interactome covers roughly 74% of the identi- focused our analysis on MYC by first dissecting its interac- fied DEGs. The interactome displays the typical character- tions with other genes in the network. A subnetwork cen- istics of a scale-free network [22] (Fig. 4). Thirty genes tered on MYC was constructed from the initial with the highest connectivity are listed in Table 3. We fur- interactome (Fig. 5A). It was found that MYC is linked to ther identified 823 angiogenesis-related genes by conduct- a variety of angiogenesis regulators with high connectivity ing automated querying of the PubMed abstract database within the network, suggesting a central role of MYC in using an in-house developed text mining tool, Target- angiogenesis regulation. We hypothesized that inactiva- Miner. The thirty most connected genes in the interac- tion of MYC function would have profound angiogenic tome were all angiogenesis related, accounting for over defects in this co-culture assay. We tested this hypothesis 33% of the total 17226 connections in the interactome. by attenuating MYC expression in HUVECs using c-MYC This data suggests a high degree of relevancy of the inter- siRNA. Cord formation was significantly suppressed in actome to angiogenesis. siMYC treated cultures (Fig. 5B &5C). Cord mass, branch- ing pattern, and responsiveness to exogenous VEGF were Subnetwork analysis of oncogene involvement in all inhibited, which is consistent with the network predic- angiogenesis tion. Knockdown of MYC expression by the siRNA was The information on which the network was based comes confirmed in parallel by TaqMan RT-PCR analysis (Fig. from many cell types. Although we filtered and refined the 5D). The involvement of TGF-β, the top node in the network by only including genes that are present in the present angiogenesis network was also validated (data not model, validating 17226 individual network edges is chal- shown). lenging. We explored the feasibility of using the present Assessment of compound antiangiogenic activity by network model to address known biological phenomena as a pathway toward functional validation of the network. network perturbation analysis Examining the connectivity hierarchy, we found that We further explored the application of network analysis to oncogenes and tumor suppressor genes are among the understand compound polypharmacology. Since kinases highest ranked molecules (hub nodes) based on the have increasingly become effective targets for anti-ang- number of connections in which they are directly iogenic drug discovery [25-27], we assessed kinase contri- engaged. For example, the number of connections that butions to angiogenesis regulation in the co-culture P53, SRC, and MYC directly engaged is 301, 243, and 173 model by network analysis. First, 39 kinases were mapped nd th th respectively, ranking at the 2 , 8 and 16 places in terms to the interactome with a total of 525 directly-connected Page 7 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 Table 3: The most connected thirty molecules in the interactome Rank Gene Connectivity Percentage 1TGFB1 397 2 2TP53 301 4 3IGF1 290 6 4 MAPK1 270 7 5FGF2 261 9 6PDGF 251 10 7AKT1 248 12 8SRC 243 13 9JUN 242 15 10 AGT 231 16 11 SP1 230 17 12 IL6 220 18 13 IL1B 216 20 14 NGFB 207 21 In HUVEC/NH Figure 4 teractome analysis of DF co-cultur the gene functional association in the e 15 VEGF 183 22 Interactome analysis of the gene functional associa- 16 MYC 173 23 tion in the HUVEC/NHDF co-culture. The interactome 17 MAPK8 169 24 that was built with PathwayAssist from a union of differen- 18 HGF 168 25 tially expressed genes identified in the microarray studies 19 RAF1 162 26 consists of 17226 connections among 1201 genes. The distri- 20 MAPK3 156 27 bution of the number of connections for 1201 genes is 21 CREB1 146 28 sorted by its connectivity, which acts as a typical scale-free 22 EP300 141 28 network. 23 EGFR 138 29 24 CDKN1A 122 30 25 RAC1 117 31 26 STAT3 115 31 the penetration is expected to decay dramatically at later 27 BCL2 109 32 steps due to multiple inputs and outputs at each node, we 28 PRL 108 33 thus defined a perturbation index calculated from the data 29 PPARG 103 33 collected by the network walking algorithm using a 30 NR3C1 101 34 weighted sum method as described in the method section. Genes within the regulation network or interactome were ranked by Using network perturbation analysis to assess drug phe- their connectivity in the network, i.e., the number of genes that they directly engaged with. Listed above are the top 30 genes with the notypic activity, we tested three multi-targeted anti-ang- highest connectivity. Gene symbols are those assigned by Human iogenic kinase (MAK) inhibitors: Sunitinib, PTK787, and Gene Nomenclature Committee (HGNC). LY2401401. The analysis was performed based on the genes, accounting for 44% of the total number of genes in number and the identities of the kinases that are inhibited the refined interactome. Any mapped kinase is engaged in by a given compound at its respective potency as the rest of the network through multi-step connections. described in the method section. Sunitinib and PTK787 Accordingly, we developed a "network walking" algo- have been previously tested in clinical trials [25], while rithm to model the connection cascade of kinases within LY2401401, is a novel Lilly MAK inhibitor with a distinct the network. The methods are described in the method kinase activity profile. LY2401401 (Fig. 7B), an imidazo- section and shown in a schematic diagram in Fig. 6. Using pyridine, is an orally bioavailable, broad-spectrum inhib- this algorithm, the characteristic connection curves of the itor of the VEGF and PDGF receptors. Additionally, Flt3, 39 mapped kinases were determined and are shown in Kit, Tie-2 and the Eph family of receptors are also sensitive Fig. 7A. All kinases were able to connect to the rest of the to this small molecule inhibitor (SMI). To differentiate interactome (except for singletons), and the connectivity these MAK compounds, it is important to determine if the plateaus at the connection of 1040 within the interac- unique kinase profile of LY2401401 is advantageous to tome. This property was used to assess the potential max- achieve anti-angiogenic efficacy. We simulated the poten- imal perturbation of the interactome if a given kinase is tial network perturbation by the three molecules based on inhibited. While every connected molecule plateaus at the the integrated perturbation index of every kinase that a same number of connections, each arrives at a different given SMI is active against. The resulting theoretical dose rate, which can be assessed by the slope of the connectiv- response curves were used to estimate the simulated EC50 ity curve. Since the spread of activity along the network of network perturbation (Fig. 7C). The compounds were can not reach 100 percent penetration at each step, and ranked by their simulated activities as Page 8 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 A Figure 5 nalysis of the involvement of MYC oncogene in angiogenesis Analysis of the involvement of MYC oncogene in angiogenesis. (A). Subnetwork analysis of MYC in the interactome. Differentially expressed genes identified from the HUVEC and NHDF microarray studies were used to build a gene regulation interactome as described in the method section. Genes connected to MYC were extracted from the interactome and visual- ized by building a MYC centric subnetwork. (B). Validation of the involvement of MYC in cord formation by siRNA mediated gene silencing. Expression of MYC in the HUVECs was attenuated by transfecting the cells with siRNA against MYC (siMYC). The siRNA treated HUVECs were co-cultured with regular NHDF cells. siRNA against luciferase (siLUC) was used as the neg- ative control. Knocking down the MYC expression in HUVECs impaired cord formation either in the presence or absence of 20 ng/ml VEGF. (C). Quantitation of cord area and angiogenesis index from the siRNA experiment. (D). Validation of reduced MYC expression in the siMYC treated HUVECs by TaqMan RT-PCR analysis. Data are expressed as Mean ± S.D. LY2401401>Sunitinib>PTK787, with the EC50 values directly in the cord formation assay. They were ranked in estimated to be 7.8 nM, 13.6 nM, and 133.7 nM, respec- the same order as inferred by the network simulation (Fig. tively. The network analysis clearly differentiates these 7D). The observed EC50 values from multiple independ- three MAK molecules. Interestingly, the network analysis ent experiments were determined to be 9.2 nM (± 9.1 nM predicted that Sunitinib would be far more potent than s.d.; n = 3), 39.0 nM (± 9.2 nM s.d.; n = 4), and 140.4 nM PTK-787, although their activities against KDR are rather (± 53.7 nM s.d.; n = 2) for LY2401401, Sunitinib, and close (49 nM, and 58 nM, respectively for PTK-787 and PTK787, respectively (Fig. 7D). The actual compound Sunitinib). activity in the cord formation bioassay was consistent with the predicted values based on the network perturba- We then determined the actual angiogenesis activity of tion analyses. each compound by testing its dose dependent response Page 9 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 Schematic diagram of the net Figure 6 work algorithm to estimate connectivity index (C ) from the interactome Schematic diagram of the network algorithm to estimate connectivity index (C ) from the interactome. Filled th circles represent the i gene (i = 1,2, ..., K), and filled diamonds the other genes reported in the literature to have one or more th functional associations with the i gene. Broken arrows are associations that were not taken into the calculation due to the th th reverse direction along the connection path. L is the number of unique genes connecting to the i gene at the j step, C is the i,j i,j th th th connectivity of the i gene at the j step, and C is the connectivity of the i gene, which is the limit of C . i i,j Discussion which are used by various malignant tumors preferen- Angiogenesis involves various cellular and multi-cellular tially [1,8]. However, the proangiogenic signals used by events such as cell proliferation, survival, differentiation, tumors are not static. It is thus unlikely to find a pleio- migration, branching and sprouting. Such complicated tropic therapy by targeting a single proangiogenic factor. processes are integrated via intrinsic and extrinsic regula- Moreover, tumors develop escape strategies by expressing tory mechanisms. Studying angiogenesis in the context of alternative proangiogenic factors [13]. Unsurprisingly, a global gene association network is critical to better mixed results have been observed in clinical trials with understand angiogenesis and to improve cancer drug tar- anti angiogenesis agents being used as a mono therapy or geting strategies. For instance, through network analysis, in combination therapies [27,29,30]. It has been well peroxisome proliferative-activated receptor δ was recently accepted in the field that targeting a single gene is not suf- identified as a "hub" node in a network that plays critical ficient [31]. However, how to target multiple genes? What roles in mediating a tumor angiogenesis switch in human are the desirable combinations of targets? And how to pancreatic cancers [28]. Over several decades, many achieve efficacy with less rebound possibility? These are proangiogenic factors have been uncovered, subsets of some key issues to be addressed. In this light, a compre- Page 10 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 A Figure 7 ssessing compound phenotypic activity through network perturbation analysis Assessing compound phenotypic activity through network perturbation analysis. (A). Characteristic connection curves of 39 kinases mapped to the interactome. 39 kinases were mapped to the interactome; collectively, they directly con- nect to 525 genes in the network. A network walking algorithm was developed to estimate the spread of the connectivity of each kinase by calculating the total number of nonredundant genes it connected at each step along the connection pathway. The smoothed connectivity curve of each of the 39 kinases was shown, with the maximum reaching 1040 for each kinase. Sum- mation of connectivity for 39 kinases is represented by the broken line. (B). The chemical structure of LY2401401, a novel small molecule and a multi angiogenesis kinase inhibitor tested in the network perturbation experiment. (C). The connectivity index was calculated by an exponentially weighted sum of the connectivity at each stage of each kinase. The dose dependent network perturbation of each small molecule inhibitor was simulated as described in the methods section. The data were imported into Graphpad Prism software to estimate the EC50 values of network perturbation using Sigmoidal dose-response curve fitting. (D). Three small molecule kinase inhibitors were tested for cord formation. They all showed dose dependent inhibition of angiogenesis. Representative image summaries of the dose dependent responses are shown. The cord area data were imported to Graphpad Prism software, and normalized to the control group for EC50 calculation, using the same method as the one that was used to calculate the simulated network perturbation values. Shown in the figure are mean EC50s for each compound, averaged from two to four independent experiments. hensive gene association network governing angiogenesis cells at different passages, from which we identified more is needed to aid rational anti-angiogenesis drug design. than two thousand genes whose expression was passage- dependent. Using a natural language processing algo- In the present study, we first characterized an in vitro ang- rithm, we retrieved a large collection of reported func- iogenesis co-culture model using high content imaging tional associations (binding, regulation, transport and analysis and found that cord formation is highly cell-pas- modification) among identified genes. Analysis of the sage dependent. We then profiled HUVEC and fibroblast interactome indicated that it is functionally related to ang- Page 11 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 iogenesis. Topologically, it displayed a typical characteris- prehensive microarray studies, or directly generated pro- tic of a scale-free network system, that is, 20% of genes tein-protein interaction data. In addition, the complexity account for 80% of functional connections. The interac- of the current interactome, which contains data of differ- tome recapitulated the involvement of oncogenes in ang- ent types, makes systematic filtering in an unbiased way iogenesis regulation, which was subsequently validated difficult and should be considered for future studies. by siRNA mediated gene silencing. We further hypothe- sized that the interactome we constructed is essential and Inferring a regulatory network from DNA microarray data sufficient to support angiogenesis in the co-culture model remains a great challenge despite the fact that numerous system. computational methods have been published [33-39]. Biological processes are complex and involve hundreds of Anti-angiogenic efficacy from a given drug can be thought genes, thus requiring a large amount of diverse data. In of as its ability to disrupt essential functional associations addition, DNA microarray experiments only interrogate in the interactome. Conceivably, inhibitors capable of transcriptionally-regulated genes, so it is inevitable that a affecting more hub nodes should deliver faster and significant portion of false positive and false negative stronger network perturbation, eliciting enhanced anti- interactions in an inferred network will be generated. In angiogenic efficacy. Therefore, we implemented a network the present study, instead of directly inferring a regulatory walking algorithm to calculate connectivity indices in the network, we constructed an interactome from differen- interactome. The connectivity index for the ith gene is a tially expressed genes by applying a literature text mining quantity extracted from the interactome that measures its technique to extract functional associations that had been effectiveness in perturbing the entire angiogenesis interac- observed from experiments; thus we anticipate fewer false tome on inhibition. This is important as a way to address positives. With the shortest path approach, we could the important polypharmacology of anti-angiogenic com- recover functional associations among genes where the pounds. expression is not transcriptionally regulated, thus we anticipate fewer false negatives. For example, 20% of the The present study takes kinase inhibitors as an example. genes constituting our interactome were not differentially Targeting multiple kinases is generally believed to be a expressed either in the HUVEC or NHDF cells, yet they viable antiangiogenic strategy and has been proven clini- ranked rather highly in the connectivity hierarchy. Since cally [25-27,29]. Sunitinib and Nexavar (sorafenib), two the natural language processing of literature text can be small molecule kinase inhibitors, are currently registered error-prone, the interactome we constructed has its own therapies. However, PTK-787 failed to achieve statistical limits. The current approach is limited to address the com- significance in clinical trials [32], although the exact plexity of in vivo cancer angiogenesis which not only mechanism accounting for the failure is not clear. Target- involves multiple cell type and signalling pathway interac- ing the right spectrum of kinases by a single SMI is critical tions, but the interactions are also geometrically and tem- for the development of successful anti-angiogenic drugs porally restricted. For example, anti-angiogenic drug- [31]. Although many kinases are reportedly angiogenesis induced inhibition of vessel formation subsequently related, they may differ in downstream signalling events, causes changes in blood flow and creates hypoxic regions crossreactivity, kinetics, etc. For a given small molecule in the tumor. Such changes are highly heterogeneous and kinase inhibitor, translation of its kinase profile to pheno- can be further compounded by many other human fac- typic activity is challenging. In the present study, we tors. Although HIF is one of the key nodes identified in explored the feasibility of estimating the anti-angiogenic our interactome, applying network analysis to address in activity of a kinase inhibitor via network perturbation vivo angiogenesis is still challenging. Moreover, the cur- analysis. The current heuristic approach addressed the rent interactome is a collection of functional associations number and identity of the kinases that a given com- among participating genes without the inclusion of feed- pound is active against as well as the respective potency. back regulatory mechanisms; it is not suitable for an The simulation was conducted to estimate network per- advanced dynamic or kinetic analysis. For these reasons, turbation by a novel MAK compound LY2401401 along we took advantage of a general statistical property of a with two clinically tested MAK compounds. The estimated biological network, namely, the characteristics of a scale- anti-angiogenic activities for the three kinase inhibitors free network system [22]. The network walking algorithm from the simulation were found to match well with their captures unevenly distributed connectivity among func- actual activities on separate bioassays (Fig. 7). Because the tionally associated genes in the constructed angiogenesis current network contains more than 17,000 edges, it interactome and uses that information to calculate how would be desirable if the network could be further refined much and how fast an external perturbation of the net- by more stringent filtering to potentially increase its pre- work can be induced by a kinase inhibitor. This was then dictive power. This might be done by refining the network used to estimate potential potency of the inhibitor for its using network connections inferred directly through com- anti-angiogenesis effect. A near perfect match between Page 12 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 estimated and measured potencies proves the statistical with or without growth factors or compounds was replen- properties of the interactome are informative for network ished every 2 to 3 days. Cells were fixed with 70% cold perturbation analysis. It would be also very interesting to ethanol at the end of the culture and processed for anti- extend our current network approach to other angiogen- CD31/Alexa 488 immunofluorescence. To visualize esis related data sets, such as the data generated from a feeder layer fibroblasts, some cultures were stained using time course study, as well as applications beyond angio- a polyclonal antibody against Vimentin (Abcam, Cam- genesis in the future. In summary, our current effort bridge, Massachusetts) along with an Alexa 568 conju- defined a comprehensive gene functional association net- gated secondary antibody to examine the feeder cell work and demonstrated the feasibility of assessing drug integrity. Nuclear staining was performed using Hoechst phenotypic activity through network perturbation analy- staining (Invitrogen Molecular Probe, Carlsbad, Califor- sis. Corroborative studies with a larger pool of com- nia). pounds, refined network parameters and an improved modelling approach would be warranted in further inves- Immunofluorescence stained co-culture plates were tigations. scanned using ArrayScan VTI (Cellomics, Pittsburgh, Pennsylvania) and the multi-parameter tube formation Bioapplication was used for quantification. Morphologi- Conclusion We characterized an in vitro angiogenesis co-culture cal measurements included cord length, width, area, per- model via high content imaging analysis. We determined cent connected cord, branching node and segment count, that HUVEC competence and NHDF supportiveness of nuclei per cord, angiogenic index (expressed as the per- cord formation are highly cell-passage dependent. We centage of area covered by connected tube), etc. General built a comprehensive gene association network to model cytotoxicity was monitored by examination of the nuclear angiogenesis regulation via complementary profiling of staining of the feed layers. passage dependent gene expression in two cell types that constituted the co-culture model. We developed a novel Attenuation of gene expression by siRNA algorithm to extract connectivity information from the HUVECs were trypsinized by 0.05% Trypsin/EDTA (Invit- gene functional association network. Using multi-targeted rogen, Carlsbad, CA), and 1 × 10 cells were nucleofected with siMYC (200 pmol each, Ambion, Austin, TX) using a anti-angiogenic kinase inhibitors as an example, our data indicate the feasibility of applying network perturbation nucleofector device (Amaxa, Gaithersburg, MD) adopting analysis to assess drug polypharmacology and phenotypic optimized A-034 protocols for HUVECs according to the activities. instructions from the manufacturer. Transfected cells were appropriately diluted and plated in 96-well plates onto the NHDF feeder cells for the initiation of the cord forma- Methods Reagents and Cells tion assay. Cells were also plated in parallel into 6-well Recombinant VEGF was purchased from R&D Systems plates and collected in Trizol (Invitrogen, Carlsbad, CA) (Minneapolis, Minnesota). Human umbilical vein 48 h or 72 h later for total RNA extraction. endothelial cells (HUVECs) and Normal Human Dermal Fibroblast (NHDF) cells and culture medium were pur- Real Time RT-PCR chased from Cambrex (Walkersville, Maryland). Cells Total RNAs were extracted from Trizol and cleaned by RNeasy mini kit (Qiagen, Chatsworth, CA). Genomic were maintained in EGM with 10% FBS and FGM-2, respectively for HUVECs and NHDFs, according to manu- DNA contamination was removed by DNA-free kit facture's instructions. In vitro angiogenesis co-culture (Ambion, Austin, TX), and cDNA synthesis was done by AngioKit and Optimized medium were purchased from using a high capacity cDNA archive kit (Applied Biosys- TCScellworks (Buckingham, UK). tems, Forster City, CA). The cDNAs were used as templates for real-time PCR with a universal PCR master mix In vitro angiogenesis assay (Applied Biosystems, Foster City, CA). Real time PCR was In vitro HUVEC cord formation co-culture was performed performed on a ABI 7900 HT. The average values were using the AngioKit commercially available from TCScell- normalized to 18s rRNA. Delta CT method was applied to works (Buckingham, UK) adapted to a 24-well format. To calculate the relative gene expression level [40]. test HUVEC and NHDF passage dependence effects, the co-culture was done in house adapted to a 96-well format. RNA isolation and microarray experiments Briefly, HUVEC cells were co-cultured onto human der- RNA isolation and microarray studies were performed as mal fibroblast feeder layers in Optimised medium (TCS- described previously [41]. Briefly, HUVECs and NHDFs cellworks, Buckingham, UK) for a total of about 12 days. with different passages were cultured to about 80% con- Treatment started the day when the plates were received; fluency. RNA was isolated using Trizol™ (Life Technolo- usually three to four days after initial plating. Medium gies, Inc. Invitrogen, Carlsbad, California) according to Page 13 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 manufacture's instructions, and followed by cleanup To identify biological processes that were significantly using RNeasy spin columns (Qiagen, Valencia, Califor- changed with passages, the identified DE genes were ana- nia). RNA labelling was performed on the Onyx robotic lyzed using the DAVID functional annotation system system from MWG-Biotech and Aviso using the standard (National Institute of Allergy and Infectious Diseases, labelling protocol. Biotin labelled cRNA was fragmented NIH) [21]. Significantly changed biological processes and hybridized to human whole genome U133 plus 2 were determined by the enrichment test provided by the GeneChip (Affymetrix, Santa Clara, California). Chip DAVID system. To control for a false positive rate of test- processing, imagine capturing, and raw data analysis were ing the expression change of thousands of genes simulta- performed using the Affymetrix Microarray Suite MAS neously, the false discovery rate (fdrate, or FDR) was with default parameters. estimated using an algorithm derived from Benjamini and Hochberg [42]. Statistical analysis to identify differentially expressed genes Construct angiogenesis related interactome and calculate Hybridized microarrays were analyzed by the Affymetrix connectivity index A tentative interactome was constructed using PathwayAs- microarray analysis software microarray suite 5 (MAS 5). The chip signals were normalized using the default sist from the union of the differentially expressed genes method in Affymetrix MAS5, but setting the target at 2%- identified either from the HUVEC or the NHDF array trimmed mean to 1500 instead of 500. HUVEC were ana- experiment. PathwayAssist was configured to retrieve five lyzed with four passages. We designated the cell passages different types of relationships defined in the software, 5, 8, 12 and 22 as 1, 2, 3 and 4 respectively. For each namely, binding, expression, molecular transport, protein probeset, the signal data were fitted to a linear regression modification and regulation. The software introduced model of time to identify consistent gene expression non-differentially expressed genes into the interactome in changes over cell passage. Three batches of the NHDF cells order to establish the shortest connection pathway were analyzed, including early passage neonatal NHDF, between any two differentially expressed genes in the ini- early passage adult NHDF, and late passage adult NHDF tial input gene list. This is an informative function since it cells. To reduce a possible batch to batch effect, the MAS 5 could recover potentially non-transcriptional regulatory signals for each chip were re-normalized by Local Polyno- relationships between two or more genes. We cleaned up mial Regression Fitting (loess) approach using one neona- this initial interactome by removing genes that were not tal sample as the baseline on a log scale. The re- detected by either array study. In addition, genes that normalized signal data were fitted to an ANOVA model to could not connect to any other genes were called orphans identify differential gene expression changes between the and removed from the interactome. Thus, the final ver- different fibroblast cells. sion of the interactome we constructed can be described as follows: To control for a false positive rate of testing the expression change of thousands of genes simultaneously, the false Let G be the set of K genes in the interactome I, and R be discovery rate (fdrate, or FDR) was estimated using an one of the five relationships: binding, molecular trans- algorithm derived from Benjamini and Hochberg [42]. port, protein modification, inhibition or up-regulation. Probesets with a false discovery rate of 0.2 or smaller were Thus, the interactome I we built can be formally described considered as significant and followed up for further anal- as ysis. We excluded probesets called "absent" in all chips in a data set. These probesets usually have very low expres- I = {(g , g ) ∈ R | g ∈ G; i, j = 1,2, ... K}(1) i i sion signals. In addition, we used expression fold change between different treatment conditions to filter out genes We assumed R is transitive, that is, if gene A regulates gene to make sure the gene expression changes are biologically B, and gene B regulates gene C, it is assumed that gene A meaningful. regulates gene C. Except for binding among the five types of relationships, molecular transport, protein modifica- PCA analysis was performed in R using the principal com- tion, up-regulation and inhibition are unidirectional. ponent analysis function with standardized data [43]. The first three principal components with > 80% variations We wrote a "network walking" algorithm (Figure 6) in a were exported into Spotfire (Sommerville, MA) to aid in recursive function with which we calculated the number th th creating sample scatter plots. Hierarchical clustering anal- of unique genes, L , connecting to the i gene at the j ij ysis of DE genes was done in Spotfire with Euclidean dis- step along the connection path in the interactome. We th th tance by the complete linkage method. define connectivity of the i gene at the j step along the connection pathway as Page 14 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 j-1 C = C + L /2 (2) Authors' contributions i,j i,j-1 ij JS conceived the project. YC, TW, JMY, and JS designed the th study. YC, FL, and LY performed the experiments. TW con- Thus, for the i gene in the interactome its connectivity at th tributed to bioinformatics analyses and algorithm devel- the j step in its connection pathway is the connectivity at th opment. H–RQ performed the statistical analysis. TW, the (j-1) step plus weighted number of unique genes at th j-1 TPB, JJS, JMY, and JS interpreted the data. TW, JMY, YC, the j step. We assigned 1/2 as an arbitrary weight to the th and JS wrote the manuscript. All authors read and number of unique genes connecting to the i gene at the th approved the final manuscript. j step, that is, the weight exponentially decreases with the th number of steps away from the i gene. Obviously C is a i,j Additional material non-decreasing function. We further define connectivity th of the i gene C as Additional file 1 CC = lim This table lists the 1103 differentially expressed genes identified in i ij , (3) jN → HUVECs. th Click here for file where C is the plateau of connectivity of the i gene and [http://www.biomedcentral.com/content/supplementary/1471- N the number of steps C takes to arrive at its limit or the i i,j 2164-9-264-S1.xls] plateau. Since every gene in the interactome is connected, the total number of genes a gene connects to is the same, Additional file 2 th that is, after N steps the connectivity of the i gene will This table lists 787 non-redundant common differentially expressed genes arrive at its plateau. However, they differ in the number of identified in NHDF cells. steps N or the rate of the arrival at the plateau. Some genes Click here for file may take more steps than the others to arrive at the plat- [http://www.biomedcentral.com/content/supplementary/1471- 2164-9-264-S2.xls] ueas, depending upon what and how many genes it con- nects, especially at early steps. Schematic elucidation of Additional file 3 the algorithm is shown in Fig. 6. This file lists the detailed information about the interactome. Click here for file Simulation of network perturbation [http://www.biomedcentral.com/content/supplementary/1471- Let S be a small molecule kinase inhibitor, and P be the s, D 2164-9-264-S3.xls] level of network perturbation by the small molecule at a hypothetical dose D. P is estimated as s, D Acknowledgements We thank Mark Uhlik, Jeff Hanson, Shaoyu Chu and Jonathan Lee for help (4) PD =×C /25×IC 0 sD,, i i s with high content image analysis, and Xiaoling Xia for the HUVECs. We i =1 would also like to thank the functional genomics group for microarray anal- where K is the number of kinases in the interactome I that ysis, John Calley for in-house gene annotation, Gregory Donoho, John Cal- can be inhibited by the small molecule S, IC50 is the i,s ley, Patrick McGlynn, Bruce Konicek, and Dan Li for a critical reading of the th IC50 value of the small molecule S against the i kinase. manuscript and Michal Vieth and Laura Bloem for kinase data. We would Their relative IC50 values against various kinases were also like to extend our thanks to Karen Morgan for her assistance in pre- measured at Upsate (UBI) or in house. To simulate the paring the manuscript. dose response curve, we set the theoretical dose D from 0.1 nM to 2 μM. If a term in (4) is more than the maxi- References 1. Folkman J: Angiogenesis: an organizing principle for drug dis- mum number of the connected genes, the term is then set covery? Nat Rev Drug Discov 2007, 6(4):273-286. to the plateau. Summed network connectivity was used to 2. Carmeliet P: Angiogenesis in life, disease and medicine. Nature calculate the percent perturbation with the plateau to be 2005, 438(7070):932-936. 3. Hanahan D, Weinberg RA: The hallmarks of cancer. Cell 2000, set at 100% network perturbation. The resulted simulated 100(1):57-70. dose dependent network perturbation was then imported 4. Hanahan D, Folkman J: Patterns and emerging mechanisms of into the program GraphPad Prism to estimate EC50 of the angiogenic switch during tumorigenesis. Cell 1996, 86(3):353-364. network perturbation by the small molecule S by per- 5. Naumov GN, Akslen LA, Folkman J: Role of angiogenesis in forming Sigmoidal dose-response curve fitting. human tumor dormancy: animal models of the angiogenic switch. Cell Cycle 2006, 5(16):1779-1787. 6. Folkman J: Tumor angiogenesis: therapeutic implications. N Abbreviations Engl J Med 1971, 285(21):1182-1186. DE: differential expression; GO: Gene Ontology; HUVEC: 7. Ferrara N, Kerbel RS: Angiogenesis as a therapeutic target. Nature 2005, 438(7070):967-974. human umbilical vein endothelial cell; MAK: multi-tar- 8. Folkman J, Klagsbrun M: Angiogenic factors. Science 1987, geted anti-angiogenic kinases; NHDF: normal human der- 235(4787):442-447. mal fibroblast; PCA: principle component analysis. Page 15 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 9. Carmeliet P, Tessier-Lavigne M: Common mechanisms of nerve 33. D'Haeseleer P, Liang S, Somogyi R: Genetic network inference: and blood vessel wiring. Nature 2005, 436(7048):193-200. from co-expression clustering to reverse engineering. Bioin- 10. Jain RK: Molecular regulation of vessel maturation. Nat Med formatics 2000, 16(8):707-726. 2003, 9(6):685-693. 34. Hartemink AJ, Gifford DK, Jaakkola TS, Young RA: Combining 11. Kerbel RS: Tumor angiogenesis: past, present and the near location and expression data for principled discovery of future. Carcinogenesis 2000, 21(3):505-515. genetic regulatory network models. Pac Symp Biocomput 12. Kerbel RS: Therapeutic implications of intrinsic or induced 2002:437-449. angiogenic growth factor redundancy in tumors revealed. 35. Pe'er D, Regev A, Elidan G, Friedman N: Inferring subnetworks Cancer Cell 2005, 8(4):269-271. from perturbed expression profiles. Bioinformatics 2001, 17 13. Casanovas O, Hicklin DJ, Bergers G, Hanahan D: Drug resistance Suppl 1:S215-24. by evasion of antiangiogenic targeting of VEGF signaling in 36. Segal E, Shapira M, Regev A, Pe'er D, Botstein D, Koller D, Friedman late-stage pancreatic islet tumors. Cancer Cell 2005, N: Module networks: identifying regulatory modules and 8(4):299-309. their condition-specific regulators from gene expression 14. Pietras K, Hanahan D: A multitargeted, metronomic, and max- data. Nat Genet 2003, 34(2):166-176. imum-tolerated dose "chemo-switch" regimen is antiang- 37. Iossifov I, Krauthammer M, Friedman C, Hatzivassiloglou V, Bader JS, iogenic, producing objective responses and survival benefit White KP, Rzhetsky A: Probabilistic inference of molecular net- in a mouse model of cancer. J Clin Oncol 2005, 23(5):939-952. works from noisy data sources. Bioinformatics 2004, 15. Jaffe EA, Nachman RL, Becker CG, Minick CR: Culture of human 20(8):1205-1213. endothelial cells derived from umbilical veins. Identification 38. Benson M, Breitling R: Network theory to understand microar- by morphologic and immunologic criteria. J Clin Invest 1973, ray studies of complex diseases. Curr Mol Med 2006, 52(11):2745-2756. 6(6):695-701. 16. Bishop ET, Bell GT, Bloor S, Broom IJ, Hendry NF, Wheatley DN: An 39. Thomas R, Paredes CJ, Mehrotra S, Hatzimanikatis V, Papoutsakis ET: in vitro model of angiogenesis: basic features. Angiogenesis A model-based optimization framework for the inference of 1999, 3(4):335-344. regulatory interactions using time-course DNA microarray 17. Neubert K, Haberland A, Kruse I, Wirth M, Schimke I: The ratio of expression data. BMC Bioinformatics 2007, 8(1):228. formation of prostacyclin/thromboxane A2 in HUVEC 40. Shou J, Soriano R, Hayward SW, Cunha GR, Williams PM, Gao WQ: decreased in each subsequent passage. Prostaglandins 1997, Expression profiling of a human cell line model of prostatic 54(1):447-462. cancer reveals a direct involvement of interferon signaling in 18. Freedman DA, Folkman J: Maintenance of G1 checkpoint con- prostate tumor progression. Proc Natl Acad Sci U S A 2002, trols in telomerase-immortalized endothelial cells. Cell Cycle 99(5):2830-2835. 2004, 3(6):811-816. 41. Shou J, Bull CM, Li L, Qian HR, Wei T, Luo S, Perkins D, Solenberg PJ, 19. Freedman DA, Folkman J: CDK2 translational down-regulation Tan SL, Chen XY, Roehm NW, Wolos JA, Onyia JE: Identification during endothelial senescence. Exp Cell Res 2005, of blood biomarkers of rheumatoid arthritis by transcript 307(1):118-130. profiling of peripheral blood mononuclear cells from the rat 20. Erusalimsky JD, Kurz DJ: Endothelial cell senescence. Handb Exp collagen-induced arthritis model. Arthritis Res Ther 2006, Pharmacol 2006:213-248. 8(1):R28. 21. Sherman BT, Huang da W, Tan Q, Guo Y, Bour S, Liu D, Stephens R, 42. Benjamini Y, Hochberg: Controlling the false discovery rate:a Baseler MW, Lane HC, Lempicki RA: DAVID Knowledgebase: a practical and powerful approach to multiple testing. The Jour- gene-centered database integrating heterogeneous gene nal of Royal Statistical Society 1995, 57(1):289-300. annotation resources to facilitate high-throughput gene 43. Venables W, Ripley B: Modern Applied Statistics with S-PLUS. functional analysis. BMC Bioinformatics 2007, 8:426. New York, Springer-Verlag; 1997. 22. Jeong H, Tombor B, Albert R, Oltvai ZN, Barabasi AL: The large- scale organization of metabolic networks. Nature 2000, 407(6804):651-654. 23. Folkman J: Tumor suppression by p53 is mediated in part by the antiangiogenic activity of endostatin and tumstatin. Sci STKE 2006, 2006(354):pe35. 24. Giuriato S, Ryeom S, Fan AC, Bachireddy P, Lynch RC, Rioth MJ, van Riggelen J, Kopelman AM, Passegue E, Tang F, Folkman J, Felsher DW: Sustained regression of tumors upon MYC inactivation requires p53 or thrombospondin-1 to reverse the angiogenic switch. Proc Natl Acad Sci U S A 2006, 103(44):16266-16271. 25. Baka S, Clamp AR, Jayson GC: A review of the latest clinical compounds to inhibit VEGF in pathological angiogenesis. Expert Opin Ther Targets 2006, 10(6):867-876. 26. Cabebe E, Wakelee H: Sunitinib: a newly approved small-mole- cule inhibitor of angiogenesis. Drugs Today (Barc) 2006, 42(6):387-398. 27. Zakarija A, Soff G: Update on angiogenesis inhibitors. Curr Opin Oncol 2005, 17(6):578-583. 28. Abdollahi A, Schwager C, Kleeff J, Esposito I, Domhan S, Peschke P, Publish with Bio Med Central and every Hauser K, Hahnfeldt P, Hlatky L, Debus J, Peters JM, Friess H, Folk- scientist can read your work free of charge man J, Huber PE: Transcriptional network governing the ang- iogenic switch in human pancreatic cancer. Proc Natl Acad Sci "BioMed Central will be the most significant development for U S A 2007, 104(31):12890-12895. disseminating the results of biomedical researc h in our lifetime." 29. Sartore-Bianchi A, Ricotta R, Cerea G, Maugeri MR, Siena S: Ration- Sir Paul Nurse, Cancer Research UK ale and clinical results of multi-target treatments in oncol- ogy. Int J Biol Markers 2007, 22(1 Suppl 4):S77-S87. Your research papers will be: 30. Meyerhardt JA, Mayer RJ: Systemic therapy for colorectal can- available free of charge to the entire biomedical community cer. N Engl J Med 2005, 352(5):476-487. 31. Zimmermann GR, Lehar J, Keith CT: Multi-target therapeutics: peer reviewed and published immediately upon acceptance when the whole is greater than the sum of the parts. Drug Dis- cited in PubMed and archived on PubMed Central cov Today 2007, 12(1-2):34-42. 32. Marx J: Cancer. Encouraging results for second-generation yours — you keep the copyright antiangiogenesis drugs. Science 2005, 308(5726):1248-1249. BioMedcentral Submit your manuscript here: http://www.biomedcentral.com/info/publishing_adv.asp Page 16 of 16 (page number not for citation purposes) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png BMC Genomics Springer Journals

Developing and applying a gene functional association network for anti-angiogenic kinase inhibitor activity assessment in an angiogenesis co-culture model

Download PDF
16 pages

Loading next page...
 
/lp/springer-journals/developing-and-applying-a-gene-functional-association-network-for-anti-mtMJPmcTAh

References (46)

Publisher
Springer Journals
Copyright
2008 Chen et al; licensee BioMed Central Ltd.
eISSN
1471-2164
DOI
10.1186/1471-2164-9-264
Publisher site
See Article on Publisher Site

Abstract

Background: Tumor angiogenesis is a highly regulated process involving intercellular communication as well as the interactions of multiple downstream signal transduction pathways. Disrupting one or even a few angiogenesis pathways is often insufficient to achieve sustained therapeutic benefits due to the complexity of angiogenesis. Targeting multiple angiogenic pathways has been increasingly recognized as a viable strategy. However, translation of the polypharmacology of a given compound to its antiangiogenic efficacy remains a major technical challenge. Developing a global functional association network among angiogenesis-related genes is much needed to facilitate holistic understanding of angiogenesis and to aid the development of more effective anti-angiogenesis therapeutics. Results: We constructed a comprehensive gene functional association network or interactome by transcript profiling an in vitro angiogenesis model, in which human umbilical vein endothelial cells (HUVECs) formed capillary structures when co-cultured with normal human dermal fibroblasts (NHDFs). HUVEC competence and NHDF supportiveness of cord formation were found to be highly cell-passage dependent. An enrichment test of Biological Processes (BP) of differentially expressed genes (DEG) revealed that angiogenesis related BP categories significantly changed with cell passages. Built upon 2012 DEGs identified from two microarray studies, the resulting interactome captured 17226 functional gene associations and displayed characteristics of a scale-free network. The interactome includes the involvement of oncogenes and tumor suppressor genes in angiogenesis. We developed a network walking algorithm to extract connectivity information from the interactome and applied it to simulate the level of network perturbation by three multi-targeted anti-angiogenic kinase inhibitors. Simulated network perturbation correlated with observed anti- angiogenesis activity in a cord formation bioassay. Conclusion: We established a comprehensive gene functional association network to model in vitro angiogenesis regulation. The present study provided a proof-of-concept pilot of applying network perturbation analysis to drug phenotypic activity assessment. Page 1 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 has also been noticed that HUVECs change their cellular Background Angiogenesis, the generation of new blood vessels, plays and molecular properties upon passage in vitro, a phe- an essential role under normal physiological conditions nomenon thought to be due to in vitro cellular senescence as well as during the pathogenesis of many diseases [17-20]. In the present study, we determined that HUVEC including atherosclerosis, macular degeneration, wound competence and NHDF supportiveness for angiogenesis healing, diabetic retinopathy, and human malignancy [1- in this co-culture system are both cell passage dependent. 3]. Remarkably, tumor dormancy is believed to be attrib- Gene Ontology (GO) analysis of differentially expressed uted, at least in part, to the lack of angiogenesis support. genes showed that the cell passage dependent global tran- The transition from an avascular, dormant tumor to an scriptional changes are highly related to angiogenesis. aggressively growing angiogenic cancer is referred to as the This enabled us to construct a comprehensive functional "angiogenic switch" [4,5]. More than 30 years ago, it was association network of differentially expressed genes hypothesized that inhibition of tumor angiogenesis using a natural language processing algorithm. We further would inhibit solid tumor growth [6]. Since then, the developed a "network walking" algorithm to estimate net- expansion of angiogenesis research has resulted in the work perturbation by small molecule kinase inhibitors. identification of various pro- and anti-angiogenic factors, The simulated compound activity via network perturba- such as FGF, VEGF, angiopoietin, endostatin, vasostatin, tion analysis was in good agreement with actual pheno- and neuronal cell axon guidance molecules [2,7-9] and typic activity in the cord formation bioassay. the development of several anti-tumor angiogenesis med- icines that have recently proven efficacious in the clinic Results [1,2,7,10,11]. Characterization of an in vitro angiogenesis co-culture model by multi parameter high content image analysis Malignant tumors often express an array of angiogenic HUVECs form a capillary structure when co-cultured with factors to potentiate angiogenesis and tumor growth [8]. NHDF in vitro [16]. Unlike the HUVEC Matrigel™ assay, Ultimate angiogenic outcomes depend on the dynamic the HUVEC/NHDF cord formation co-culture allows us to equilibrium between positive and negative regulators and study cell-cell interactions during long term (up to 14 the interplay among their signal transduction pathways, days) angiogenesis. The formation of the cord structure is not on a single discrete pathway. Avastin, a monoclonal a highly organized morphogenic process. We developed a antibody that specifically blocks VEGF, is the first high content imaging based assay to document and quan- approved anti-angiogenic therapy. Notably, combination titate cord formation using the Cellomics ArrayScan plat- with chemotherapies has been necessary for its clinical form. We were able to simultaneously assess dozens of efficacy. Moreover, tumors often develop drug resistance morphological measurements, such as tube length, area, in response to anti-VEGF therapy [12]. Hanahan et al connectivity of the cord, branching, and number of nuclei demonstrated one possible mechanism by which such within the cord. Using this assay system, we tested a clin- drug resistance might develop[13]. Inhibition of VEGF ically approved antiangiogenic compound Sunitinib signalling by neutralizing antibody to VEGFR2 (KDR) (SUTENT, SU11248). While 20 ng/ml VEGF significantly induces elevated expression of hypoxia associated proan- increased angiogenesis, 30 nM Sunitinib significantly giogenic factors such as FGF and EphrinA1, which subse- inhibited angiogenesis (Fig. 1A). quently reactivates VEGF independent angiogenesis and tumor growth. Thus, disrupting a single proangiogenic To explore the phenotypic signatures induced by VEGF or pathway by itself is often insufficient to achieve sustained Sunitinib, we performed principal component analysis therapeutic benefits. In this light, it is necessary to explore (PCA) and hierarchical clustering (Fig. 1B &1C). Multivar- global functional association among angiogenesis-related iate analyses using 47 morphological measurements genes rather than focusing on an individual or a few ang- clearly distinguished the three treatment groups, i.e. con- iogenesis factors discretely. trol, VEGF, and Sunitinib. In the PCA analysis, the first three principal components captured more than 80% of Importantly, angiogenesis involves intercellular interac- the total variance. Examining the PC loadings, we found tion among vessel-forming endothelial cells, nonmalig- that among the most discriminative features are the cord nant cells such as the supporting pericytes, immune and length (either connected or unconnected or in general), stromal cells, as well as the malignant tumor cells [1,14]. cord area, the branching segment count. The cord CD31 Therefore, global gene-gene interactions during angiogen- immunofluorescence intensity is also one of the features esis need to be explored in a multi-cell type context. that significantly contributed to PC1, likely due to the Human umbilical vein endothelial cells (HUVECs) are overall correlation of fluorescence staining to cord mass. widely used to study vascular biology [15]. They form a Interestingly, the cord width is not a discriminative fea- lumen bearing capillary structure when co-cultured with ture for either VEGF or Sunitinib action. The current normal human dermal fibroblast (NHDF) cells [16]. It multi-parameter, high content analysis allowed for a Page 2 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 Multi p Figure 1 arameter characterization of in vitro angiogenesis in a co-culture model Multi parameter characterization of in vitro angiogenesis in a co-culture model. (A) HUVECs formed capillary structure (cord) when plated on dermal fibroblast cells over the period of 12–14 days. The angiogenic cord structure was vis- ualized by human CD31 staining (green), fibroblast cells were stained with an intermediate filament marker vimentin (red), while the nucleus of both cell types were visualized with Hoechst staining (blue). Cord formation was greatly enhanced by addi- tion of 20 ng/ml VEGF (VEGF) to the culture, while 30 nM Sunitinib (CMPD) drastically inhibited cord formation as compared to the medium control (CTL). (B). Morphological features of the cord formed in the co-culture were captured and quantitated by high content image analysis using the ArrayScan platform adopting the tube formation bioapplication. 47 measurements were assessed, normalized using the Z-score across different treatment groups and analyzed by Principal Component Analysis (PCA). The first three principal components accounting for more than 80% of the variance were selected for sample scatter plotting. Blue are the control cultures (CTL); yellow, cultures treated with 20 ng/ml VEGF; and red, cultures treated with 30 nM Sunitinib (CMPD) treated samples. (C). The Z-scores of 47 morphological measurements were clustered in Spotfire Deci- sion Site with the correlation coefficient being used as the distance metric. Complete linkage was used as the clustering algo- rithm. Red stands for higher and green for lower than mean (black) expression. (D). Six representative parameters were selected and plotted. Except for the cord width, the cord length, area, branching point, the number of nucleus within the cord, and a composite angiogenesis index (Angio Index) all manifested the anti and pro angiogenic effects from Sunitinib and VEGF, respectively. Data are expressed as Mean ± S.D. detailed dissection of the phenotypic characteristics asso- HUVECs in the co-culture. Sunitinib treated cultures are ciated with a given factor or compound. Shown in Fig. 1D less morphologically complex, as the number of con- is the mean value plot of a few representative parameters. nected cords, the branching nodes, and segment measure- Sunitinib significantly inhibited the mass of the cord in ments are all significantly reduced. length and area, but not width. Moreover, the number of nuclei per tube is significantly reduced, indicating that Sunitinib could have anti-proliferation activity on Page 3 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 HUVEC cord formation competence is cell passage various passages (P6, P14, P22) and tested them for cord dependent formation competence. We found that they gradually lost It has been well known that HUVECs change their bio- cord formation capabilities upon passage (Fig. 2A). While chemical and cellular behaviors when passaged in vitro early passage (P6) HUVECs formed very nice capillary [17-20]. We cultured HUVEC cells from a single batch to structures and responded robustly to VEGF when co-cul- H Figure 2 UVEC passage dependent angiogenic competence in cord formation assay HUVEC passage dependent angiogenic competence in cord formation assay. HUVECs cultured to various passage numbers (P6, early; P14, intermediate; and P22, late) were co-cultured with NHDFs in the cord formation assay. The cord was visualized by human CD31 staining (green) and quantitated. The cultures were treated with or without 20 ng/ml exogenously added VEGF. Representative images are shown in (A). 6 well-level parameters were quantitated and plotted in (B). Please note that HUVECs gradually lost their cord formation competence when passaged in vitro. Data are expressed as Mean ± S.D. (C- D). HUVECs with various passage numbers (P5, P8, P12, and P22) were assayed for basal gene expression using the Affymetrix GeneChip platform. A linear regression method was used to identify differentially expressed genes over cell passage. An enrichment testing was used to identify significantly changed Biological Process (BP) categories associated with HUVEC pas- sage-dependent gene expression. Representative BP categories were shown using heatmap visualization, for genes associated with mitotic cell cycle (C), and cell death (D). Page 4 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 tured with NHDF, late passage (P22) HUVECs lost cord Profiling passage dependent angiogenesis support by formation capability and had diminished response to dermal fibroblast cells Long term sustained angiogenesis in co-culture is sup- VEGF (Fig. 2A &2B). ported by dermal fibroblast feeder layer cells, which are We next profiled gene expression changes associated with also affected by cell passage. While late passage (P12) HUVEC passage. Since the loss of angiogenesis capability adult NHDFs demonstrated poor support for HUVECs to was proportional to the round of cell passage, we applied form capillary structures, early passage (P5) adult cells did a statistical approach to identify genes with regressive well (Fig. 3A), suggesting that supportiveness of NHDFs expression patterns, i.e. their expression either increased in the angiogenesis assay is indeed cell-passage depend- or decreased over passage numbers. We identified 1103 ent. An early passage NHDF from a neonatal donor also differentially expressed (DE) genes which are listed in supported angiogenesis in this assay (Fig. 3A). Supplemental Table 1 (see Additional file 1). In order to understand the biological relevance of the DE genes, we Differential angiogenesis supportiveness by NHDFs in the analyzed the DE genes by using the DAVID functional co-culture system represents a unique opportunity for us annotation system [21]. Table 1 lists Gene Ontology (GO) to dissect factors that may contribute to the complex inter- terms significantly enriched with the DE genes. The most cellular communications during angiogenesis. To assess significant Biological Processes (BP) categories are cell this, we profiled gene expression in the three listed prepa- cycle, proliferation and programmed cell death. Expres- rations of NHDF cells by DNA microarray analysis. sion profiles of the DE genes involved in representative Between the supportive NHDF cells versus the non-sup- biological processes were analyzed by hierarchical cluster- portive ones, we identified 787 non-redundant common ing analysis. The majority of the DE genes involved in DE genes which are listed in Supplemental Table 2 (see mitotic cell cycle and proliferation were down-regulated, Additional file 2). Similarly, we used the DAVID func- whereas the majority of genes related to cell death were tional annotation system to examine the functional rele- up-regulated (Fig. 2C–D), which supports the replicative vance to angiogenesis of these DEGs. Table 2 lists the cell senescence hypothesis. significantly enriched GO terms. Interestingly, blood ves- Table 1: Gene Ontology (GO) terms enriched with differentially expressed genes in HUVEC of different passages 1 2 3 GO GO Term DEG FDR Biological Process cell cycle 93 0.00075 mitotic cell cycle 44 0.00042 M phase of mitotic cell cycle 36 0.00124 programmed cell death 44 0.00717 mitosis 35 0.00511 cell division 37 0.00332 regulation of cell cycle 60 0.01370 cell cycle checkpoint 13 0.01512 cell organization and biogenesis 113 0.05000 cell proliferation 51 0.01164 Molecular Function protein binding 270 0.00051 ATP binding 106 0.00048 nucleotide binding 137 0.00473 catalytic activity 309 0.00144 purine nucleotide binding 119 0.00210 cytoskeletal protein binding 33 0.00545 kinase activity 70 0.03651 binding 545 0.03367 transferase activity 116 0.01551 Cellular Component spindle 14 0.00058 microtubule cytoskeleton 36 0.00182 intracellular 439 0.00138 cytoskeleton 73 0.00558 nucleus 237 0.00615 intracellular membrane-bound organelle 317 0.04183 intracellular organelle 364 0.01793 1. Gene Ontology; 2. Differentially expressed genes; 3. False discovery rate Listed are GO categories identified through enrichment testing of the differentially expressed genes from HUVEC using the DAVID functional annotation system as described in the methods section. Page 5 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 NHDF pa Figure 3 ssage dependent angiogenic supportiveness in cord formation assay NHDF passage dependent angiogenic supportiveness in cord formation assay. P6 NHDFs from a neonatal donor, P5 and P12 NHDFs from an adult donor were tested for their angiogenesis supportiveness in the cord formation assay (A). Cultures were performed either with or without 20 ng/ml exogenously added VEGF. Cords were visualized by human CD31 staining (green). Both early passage neonatal and adult NHDF cells are more supportive for cord formation than the late pas- sage adult NHDFs, with the late passage NHDF cells being non-permissive. (B). Clustering analysis of genes involved in angio- genesis. P6 NHDFs from a neonatal donor (Neo P6), P4 NHDFs from an adult donor (Ad P4), and P13 NHDFs (Ad P13) from the same adult donor were assayed for basal gene expression using the Affymetrix GeneChip platform. 787 nonredundant genes were identified as common differentially expressed genes between the two early passage NHDF and the late passage NHDF cells. Enrichment testing was used to identify significant Biological Process (BP) categories. Representative angiogenesis BP categories are shown in a heatmap visualization. sel development and morphogenesis, angiogenesis and both HUVEC and NHDF cells influences angiogenesis vasculature development were among the most signifi- supports modeling the angiogenesis network using the cantly enriched biological processes, suggesting that these DEGs identified by these DNA microarray experiments. molecular changes were related to cellular potential to To understand the functional association of these genes, support angiogenesis. Representative expression profiles we constructed a molecular interactome among the total of genes involved in blood vessel development and angio- 2012 differentially expressed genes identified from either genesis are shown in Fig. 3B. This result is consistent with NHDFs or HUVECs. The interactome consisted of over the supporting roles of NHDFs in this co-culture system. 23,000 gene-to-gene connections retrieved from PubMed Taken together, passaging NHDF cells in vitro adversely abstracts using the natural language processing algorithm affects cellular capacity to support angiogenesis due to sig- in PathwayAssist. We then refined the interactome by nificantly orchestrated changes of many genes involved in removing genes whose expression was not detectable by a few related biological processes. microarray analysis in the co-culture system. The resulted molecular interactome consists of 17226 molecular rela- Establishment of a gene functional association network- tionships described by the following terms: binding, The fact that cell passage dependent gene expression in expression, molecular transport, protein modification Page 6 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 Table 2: Gene Ontology terms enriched with differentially expressed genes in NHDF cells of different passages 1 2 3 GO GO Term DEG FDR Biological Process blood vessel development 64 0.00700 vasculature development 64 0.00700 actin filament organization 40 0.00850 actin cytoskeleton organization and biogenesis 81 0.00850 actin filament-based process 84 0.00850 regulation of signal transduction 98 0.01020 angiogenesis 56 0.01140 blood vessel morphogenesis 60 0.01140 cell morphogenesis 161 0.02020 Cellular Component extracellular matrix 75 0.00200 proteinaceous extracellular matrix 74 0.00200 actin cytoskeleton 60 0.04800 adherens junction 30 0.06140 cell-matrix junction 25 0.06140 cytoskeleton 132 0.06140 Molecular Function protein binding 864 0.00820 binding 1195 0.01740 cytoskeletal protein binding 95 0.02380 1. Gene Ontology; 2. Differentially expressed genes; 3. False discovery rate. Listed are GO categories identified through enrichment testing of the differentially expressed genes from NHDF cells using the DAVID functional annotation system as described in the methods section. and regulation (see Additional file 3). There were 1201 of connectivity in the interactome. This observation reca- unique genes in the interactome, among which 887 genes pitulates the emerging role of oncogenes and tumor sup- were differentially expressed in the co-culture system. pressor genes in angiogenesis regulation [1,5,23,24]. We Thus, the interactome covers roughly 74% of the identi- focused our analysis on MYC by first dissecting its interac- fied DEGs. The interactome displays the typical character- tions with other genes in the network. A subnetwork cen- istics of a scale-free network [22] (Fig. 4). Thirty genes tered on MYC was constructed from the initial with the highest connectivity are listed in Table 3. We fur- interactome (Fig. 5A). It was found that MYC is linked to ther identified 823 angiogenesis-related genes by conduct- a variety of angiogenesis regulators with high connectivity ing automated querying of the PubMed abstract database within the network, suggesting a central role of MYC in using an in-house developed text mining tool, Target- angiogenesis regulation. We hypothesized that inactiva- Miner. The thirty most connected genes in the interac- tion of MYC function would have profound angiogenic tome were all angiogenesis related, accounting for over defects in this co-culture assay. We tested this hypothesis 33% of the total 17226 connections in the interactome. by attenuating MYC expression in HUVECs using c-MYC This data suggests a high degree of relevancy of the inter- siRNA. Cord formation was significantly suppressed in actome to angiogenesis. siMYC treated cultures (Fig. 5B &5C). Cord mass, branch- ing pattern, and responsiveness to exogenous VEGF were Subnetwork analysis of oncogene involvement in all inhibited, which is consistent with the network predic- angiogenesis tion. Knockdown of MYC expression by the siRNA was The information on which the network was based comes confirmed in parallel by TaqMan RT-PCR analysis (Fig. from many cell types. Although we filtered and refined the 5D). The involvement of TGF-β, the top node in the network by only including genes that are present in the present angiogenesis network was also validated (data not model, validating 17226 individual network edges is chal- shown). lenging. We explored the feasibility of using the present Assessment of compound antiangiogenic activity by network model to address known biological phenomena as a pathway toward functional validation of the network. network perturbation analysis Examining the connectivity hierarchy, we found that We further explored the application of network analysis to oncogenes and tumor suppressor genes are among the understand compound polypharmacology. Since kinases highest ranked molecules (hub nodes) based on the have increasingly become effective targets for anti-ang- number of connections in which they are directly iogenic drug discovery [25-27], we assessed kinase contri- engaged. For example, the number of connections that butions to angiogenesis regulation in the co-culture P53, SRC, and MYC directly engaged is 301, 243, and 173 model by network analysis. First, 39 kinases were mapped nd th th respectively, ranking at the 2 , 8 and 16 places in terms to the interactome with a total of 525 directly-connected Page 7 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 Table 3: The most connected thirty molecules in the interactome Rank Gene Connectivity Percentage 1TGFB1 397 2 2TP53 301 4 3IGF1 290 6 4 MAPK1 270 7 5FGF2 261 9 6PDGF 251 10 7AKT1 248 12 8SRC 243 13 9JUN 242 15 10 AGT 231 16 11 SP1 230 17 12 IL6 220 18 13 IL1B 216 20 14 NGFB 207 21 In HUVEC/NH Figure 4 teractome analysis of DF co-cultur the gene functional association in the e 15 VEGF 183 22 Interactome analysis of the gene functional associa- 16 MYC 173 23 tion in the HUVEC/NHDF co-culture. The interactome 17 MAPK8 169 24 that was built with PathwayAssist from a union of differen- 18 HGF 168 25 tially expressed genes identified in the microarray studies 19 RAF1 162 26 consists of 17226 connections among 1201 genes. The distri- 20 MAPK3 156 27 bution of the number of connections for 1201 genes is 21 CREB1 146 28 sorted by its connectivity, which acts as a typical scale-free 22 EP300 141 28 network. 23 EGFR 138 29 24 CDKN1A 122 30 25 RAC1 117 31 26 STAT3 115 31 the penetration is expected to decay dramatically at later 27 BCL2 109 32 steps due to multiple inputs and outputs at each node, we 28 PRL 108 33 thus defined a perturbation index calculated from the data 29 PPARG 103 33 collected by the network walking algorithm using a 30 NR3C1 101 34 weighted sum method as described in the method section. Genes within the regulation network or interactome were ranked by Using network perturbation analysis to assess drug phe- their connectivity in the network, i.e., the number of genes that they directly engaged with. Listed above are the top 30 genes with the notypic activity, we tested three multi-targeted anti-ang- highest connectivity. Gene symbols are those assigned by Human iogenic kinase (MAK) inhibitors: Sunitinib, PTK787, and Gene Nomenclature Committee (HGNC). LY2401401. The analysis was performed based on the genes, accounting for 44% of the total number of genes in number and the identities of the kinases that are inhibited the refined interactome. Any mapped kinase is engaged in by a given compound at its respective potency as the rest of the network through multi-step connections. described in the method section. Sunitinib and PTK787 Accordingly, we developed a "network walking" algo- have been previously tested in clinical trials [25], while rithm to model the connection cascade of kinases within LY2401401, is a novel Lilly MAK inhibitor with a distinct the network. The methods are described in the method kinase activity profile. LY2401401 (Fig. 7B), an imidazo- section and shown in a schematic diagram in Fig. 6. Using pyridine, is an orally bioavailable, broad-spectrum inhib- this algorithm, the characteristic connection curves of the itor of the VEGF and PDGF receptors. Additionally, Flt3, 39 mapped kinases were determined and are shown in Kit, Tie-2 and the Eph family of receptors are also sensitive Fig. 7A. All kinases were able to connect to the rest of the to this small molecule inhibitor (SMI). To differentiate interactome (except for singletons), and the connectivity these MAK compounds, it is important to determine if the plateaus at the connection of 1040 within the interac- unique kinase profile of LY2401401 is advantageous to tome. This property was used to assess the potential max- achieve anti-angiogenic efficacy. We simulated the poten- imal perturbation of the interactome if a given kinase is tial network perturbation by the three molecules based on inhibited. While every connected molecule plateaus at the the integrated perturbation index of every kinase that a same number of connections, each arrives at a different given SMI is active against. The resulting theoretical dose rate, which can be assessed by the slope of the connectiv- response curves were used to estimate the simulated EC50 ity curve. Since the spread of activity along the network of network perturbation (Fig. 7C). The compounds were can not reach 100 percent penetration at each step, and ranked by their simulated activities as Page 8 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 A Figure 5 nalysis of the involvement of MYC oncogene in angiogenesis Analysis of the involvement of MYC oncogene in angiogenesis. (A). Subnetwork analysis of MYC in the interactome. Differentially expressed genes identified from the HUVEC and NHDF microarray studies were used to build a gene regulation interactome as described in the method section. Genes connected to MYC were extracted from the interactome and visual- ized by building a MYC centric subnetwork. (B). Validation of the involvement of MYC in cord formation by siRNA mediated gene silencing. Expression of MYC in the HUVECs was attenuated by transfecting the cells with siRNA against MYC (siMYC). The siRNA treated HUVECs were co-cultured with regular NHDF cells. siRNA against luciferase (siLUC) was used as the neg- ative control. Knocking down the MYC expression in HUVECs impaired cord formation either in the presence or absence of 20 ng/ml VEGF. (C). Quantitation of cord area and angiogenesis index from the siRNA experiment. (D). Validation of reduced MYC expression in the siMYC treated HUVECs by TaqMan RT-PCR analysis. Data are expressed as Mean ± S.D. LY2401401>Sunitinib>PTK787, with the EC50 values directly in the cord formation assay. They were ranked in estimated to be 7.8 nM, 13.6 nM, and 133.7 nM, respec- the same order as inferred by the network simulation (Fig. tively. The network analysis clearly differentiates these 7D). The observed EC50 values from multiple independ- three MAK molecules. Interestingly, the network analysis ent experiments were determined to be 9.2 nM (± 9.1 nM predicted that Sunitinib would be far more potent than s.d.; n = 3), 39.0 nM (± 9.2 nM s.d.; n = 4), and 140.4 nM PTK-787, although their activities against KDR are rather (± 53.7 nM s.d.; n = 2) for LY2401401, Sunitinib, and close (49 nM, and 58 nM, respectively for PTK-787 and PTK787, respectively (Fig. 7D). The actual compound Sunitinib). activity in the cord formation bioassay was consistent with the predicted values based on the network perturba- We then determined the actual angiogenesis activity of tion analyses. each compound by testing its dose dependent response Page 9 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 Schematic diagram of the net Figure 6 work algorithm to estimate connectivity index (C ) from the interactome Schematic diagram of the network algorithm to estimate connectivity index (C ) from the interactome. Filled th circles represent the i gene (i = 1,2, ..., K), and filled diamonds the other genes reported in the literature to have one or more th functional associations with the i gene. Broken arrows are associations that were not taken into the calculation due to the th th reverse direction along the connection path. L is the number of unique genes connecting to the i gene at the j step, C is the i,j i,j th th th connectivity of the i gene at the j step, and C is the connectivity of the i gene, which is the limit of C . i i,j Discussion which are used by various malignant tumors preferen- Angiogenesis involves various cellular and multi-cellular tially [1,8]. However, the proangiogenic signals used by events such as cell proliferation, survival, differentiation, tumors are not static. It is thus unlikely to find a pleio- migration, branching and sprouting. Such complicated tropic therapy by targeting a single proangiogenic factor. processes are integrated via intrinsic and extrinsic regula- Moreover, tumors develop escape strategies by expressing tory mechanisms. Studying angiogenesis in the context of alternative proangiogenic factors [13]. Unsurprisingly, a global gene association network is critical to better mixed results have been observed in clinical trials with understand angiogenesis and to improve cancer drug tar- anti angiogenesis agents being used as a mono therapy or geting strategies. For instance, through network analysis, in combination therapies [27,29,30]. It has been well peroxisome proliferative-activated receptor δ was recently accepted in the field that targeting a single gene is not suf- identified as a "hub" node in a network that plays critical ficient [31]. However, how to target multiple genes? What roles in mediating a tumor angiogenesis switch in human are the desirable combinations of targets? And how to pancreatic cancers [28]. Over several decades, many achieve efficacy with less rebound possibility? These are proangiogenic factors have been uncovered, subsets of some key issues to be addressed. In this light, a compre- Page 10 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 A Figure 7 ssessing compound phenotypic activity through network perturbation analysis Assessing compound phenotypic activity through network perturbation analysis. (A). Characteristic connection curves of 39 kinases mapped to the interactome. 39 kinases were mapped to the interactome; collectively, they directly con- nect to 525 genes in the network. A network walking algorithm was developed to estimate the spread of the connectivity of each kinase by calculating the total number of nonredundant genes it connected at each step along the connection pathway. The smoothed connectivity curve of each of the 39 kinases was shown, with the maximum reaching 1040 for each kinase. Sum- mation of connectivity for 39 kinases is represented by the broken line. (B). The chemical structure of LY2401401, a novel small molecule and a multi angiogenesis kinase inhibitor tested in the network perturbation experiment. (C). The connectivity index was calculated by an exponentially weighted sum of the connectivity at each stage of each kinase. The dose dependent network perturbation of each small molecule inhibitor was simulated as described in the methods section. The data were imported into Graphpad Prism software to estimate the EC50 values of network perturbation using Sigmoidal dose-response curve fitting. (D). Three small molecule kinase inhibitors were tested for cord formation. They all showed dose dependent inhibition of angiogenesis. Representative image summaries of the dose dependent responses are shown. The cord area data were imported to Graphpad Prism software, and normalized to the control group for EC50 calculation, using the same method as the one that was used to calculate the simulated network perturbation values. Shown in the figure are mean EC50s for each compound, averaged from two to four independent experiments. hensive gene association network governing angiogenesis cells at different passages, from which we identified more is needed to aid rational anti-angiogenesis drug design. than two thousand genes whose expression was passage- dependent. Using a natural language processing algo- In the present study, we first characterized an in vitro ang- rithm, we retrieved a large collection of reported func- iogenesis co-culture model using high content imaging tional associations (binding, regulation, transport and analysis and found that cord formation is highly cell-pas- modification) among identified genes. Analysis of the sage dependent. We then profiled HUVEC and fibroblast interactome indicated that it is functionally related to ang- Page 11 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 iogenesis. Topologically, it displayed a typical characteris- prehensive microarray studies, or directly generated pro- tic of a scale-free network system, that is, 20% of genes tein-protein interaction data. In addition, the complexity account for 80% of functional connections. The interac- of the current interactome, which contains data of differ- tome recapitulated the involvement of oncogenes in ang- ent types, makes systematic filtering in an unbiased way iogenesis regulation, which was subsequently validated difficult and should be considered for future studies. by siRNA mediated gene silencing. We further hypothe- sized that the interactome we constructed is essential and Inferring a regulatory network from DNA microarray data sufficient to support angiogenesis in the co-culture model remains a great challenge despite the fact that numerous system. computational methods have been published [33-39]. Biological processes are complex and involve hundreds of Anti-angiogenic efficacy from a given drug can be thought genes, thus requiring a large amount of diverse data. In of as its ability to disrupt essential functional associations addition, DNA microarray experiments only interrogate in the interactome. Conceivably, inhibitors capable of transcriptionally-regulated genes, so it is inevitable that a affecting more hub nodes should deliver faster and significant portion of false positive and false negative stronger network perturbation, eliciting enhanced anti- interactions in an inferred network will be generated. In angiogenic efficacy. Therefore, we implemented a network the present study, instead of directly inferring a regulatory walking algorithm to calculate connectivity indices in the network, we constructed an interactome from differen- interactome. The connectivity index for the ith gene is a tially expressed genes by applying a literature text mining quantity extracted from the interactome that measures its technique to extract functional associations that had been effectiveness in perturbing the entire angiogenesis interac- observed from experiments; thus we anticipate fewer false tome on inhibition. This is important as a way to address positives. With the shortest path approach, we could the important polypharmacology of anti-angiogenic com- recover functional associations among genes where the pounds. expression is not transcriptionally regulated, thus we anticipate fewer false negatives. For example, 20% of the The present study takes kinase inhibitors as an example. genes constituting our interactome were not differentially Targeting multiple kinases is generally believed to be a expressed either in the HUVEC or NHDF cells, yet they viable antiangiogenic strategy and has been proven clini- ranked rather highly in the connectivity hierarchy. Since cally [25-27,29]. Sunitinib and Nexavar (sorafenib), two the natural language processing of literature text can be small molecule kinase inhibitors, are currently registered error-prone, the interactome we constructed has its own therapies. However, PTK-787 failed to achieve statistical limits. The current approach is limited to address the com- significance in clinical trials [32], although the exact plexity of in vivo cancer angiogenesis which not only mechanism accounting for the failure is not clear. Target- involves multiple cell type and signalling pathway interac- ing the right spectrum of kinases by a single SMI is critical tions, but the interactions are also geometrically and tem- for the development of successful anti-angiogenic drugs porally restricted. For example, anti-angiogenic drug- [31]. Although many kinases are reportedly angiogenesis induced inhibition of vessel formation subsequently related, they may differ in downstream signalling events, causes changes in blood flow and creates hypoxic regions crossreactivity, kinetics, etc. For a given small molecule in the tumor. Such changes are highly heterogeneous and kinase inhibitor, translation of its kinase profile to pheno- can be further compounded by many other human fac- typic activity is challenging. In the present study, we tors. Although HIF is one of the key nodes identified in explored the feasibility of estimating the anti-angiogenic our interactome, applying network analysis to address in activity of a kinase inhibitor via network perturbation vivo angiogenesis is still challenging. Moreover, the cur- analysis. The current heuristic approach addressed the rent interactome is a collection of functional associations number and identity of the kinases that a given com- among participating genes without the inclusion of feed- pound is active against as well as the respective potency. back regulatory mechanisms; it is not suitable for an The simulation was conducted to estimate network per- advanced dynamic or kinetic analysis. For these reasons, turbation by a novel MAK compound LY2401401 along we took advantage of a general statistical property of a with two clinically tested MAK compounds. The estimated biological network, namely, the characteristics of a scale- anti-angiogenic activities for the three kinase inhibitors free network system [22]. The network walking algorithm from the simulation were found to match well with their captures unevenly distributed connectivity among func- actual activities on separate bioassays (Fig. 7). Because the tionally associated genes in the constructed angiogenesis current network contains more than 17,000 edges, it interactome and uses that information to calculate how would be desirable if the network could be further refined much and how fast an external perturbation of the net- by more stringent filtering to potentially increase its pre- work can be induced by a kinase inhibitor. This was then dictive power. This might be done by refining the network used to estimate potential potency of the inhibitor for its using network connections inferred directly through com- anti-angiogenesis effect. A near perfect match between Page 12 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 estimated and measured potencies proves the statistical with or without growth factors or compounds was replen- properties of the interactome are informative for network ished every 2 to 3 days. Cells were fixed with 70% cold perturbation analysis. It would be also very interesting to ethanol at the end of the culture and processed for anti- extend our current network approach to other angiogen- CD31/Alexa 488 immunofluorescence. To visualize esis related data sets, such as the data generated from a feeder layer fibroblasts, some cultures were stained using time course study, as well as applications beyond angio- a polyclonal antibody against Vimentin (Abcam, Cam- genesis in the future. In summary, our current effort bridge, Massachusetts) along with an Alexa 568 conju- defined a comprehensive gene functional association net- gated secondary antibody to examine the feeder cell work and demonstrated the feasibility of assessing drug integrity. Nuclear staining was performed using Hoechst phenotypic activity through network perturbation analy- staining (Invitrogen Molecular Probe, Carlsbad, Califor- sis. Corroborative studies with a larger pool of com- nia). pounds, refined network parameters and an improved modelling approach would be warranted in further inves- Immunofluorescence stained co-culture plates were tigations. scanned using ArrayScan VTI (Cellomics, Pittsburgh, Pennsylvania) and the multi-parameter tube formation Bioapplication was used for quantification. Morphologi- Conclusion We characterized an in vitro angiogenesis co-culture cal measurements included cord length, width, area, per- model via high content imaging analysis. We determined cent connected cord, branching node and segment count, that HUVEC competence and NHDF supportiveness of nuclei per cord, angiogenic index (expressed as the per- cord formation are highly cell-passage dependent. We centage of area covered by connected tube), etc. General built a comprehensive gene association network to model cytotoxicity was monitored by examination of the nuclear angiogenesis regulation via complementary profiling of staining of the feed layers. passage dependent gene expression in two cell types that constituted the co-culture model. We developed a novel Attenuation of gene expression by siRNA algorithm to extract connectivity information from the HUVECs were trypsinized by 0.05% Trypsin/EDTA (Invit- gene functional association network. Using multi-targeted rogen, Carlsbad, CA), and 1 × 10 cells were nucleofected with siMYC (200 pmol each, Ambion, Austin, TX) using a anti-angiogenic kinase inhibitors as an example, our data indicate the feasibility of applying network perturbation nucleofector device (Amaxa, Gaithersburg, MD) adopting analysis to assess drug polypharmacology and phenotypic optimized A-034 protocols for HUVECs according to the activities. instructions from the manufacturer. Transfected cells were appropriately diluted and plated in 96-well plates onto the NHDF feeder cells for the initiation of the cord forma- Methods Reagents and Cells tion assay. Cells were also plated in parallel into 6-well Recombinant VEGF was purchased from R&D Systems plates and collected in Trizol (Invitrogen, Carlsbad, CA) (Minneapolis, Minnesota). Human umbilical vein 48 h or 72 h later for total RNA extraction. endothelial cells (HUVECs) and Normal Human Dermal Fibroblast (NHDF) cells and culture medium were pur- Real Time RT-PCR chased from Cambrex (Walkersville, Maryland). Cells Total RNAs were extracted from Trizol and cleaned by RNeasy mini kit (Qiagen, Chatsworth, CA). Genomic were maintained in EGM with 10% FBS and FGM-2, respectively for HUVECs and NHDFs, according to manu- DNA contamination was removed by DNA-free kit facture's instructions. In vitro angiogenesis co-culture (Ambion, Austin, TX), and cDNA synthesis was done by AngioKit and Optimized medium were purchased from using a high capacity cDNA archive kit (Applied Biosys- TCScellworks (Buckingham, UK). tems, Forster City, CA). The cDNAs were used as templates for real-time PCR with a universal PCR master mix In vitro angiogenesis assay (Applied Biosystems, Foster City, CA). Real time PCR was In vitro HUVEC cord formation co-culture was performed performed on a ABI 7900 HT. The average values were using the AngioKit commercially available from TCScell- normalized to 18s rRNA. Delta CT method was applied to works (Buckingham, UK) adapted to a 24-well format. To calculate the relative gene expression level [40]. test HUVEC and NHDF passage dependence effects, the co-culture was done in house adapted to a 96-well format. RNA isolation and microarray experiments Briefly, HUVEC cells were co-cultured onto human der- RNA isolation and microarray studies were performed as mal fibroblast feeder layers in Optimised medium (TCS- described previously [41]. Briefly, HUVECs and NHDFs cellworks, Buckingham, UK) for a total of about 12 days. with different passages were cultured to about 80% con- Treatment started the day when the plates were received; fluency. RNA was isolated using Trizol™ (Life Technolo- usually three to four days after initial plating. Medium gies, Inc. Invitrogen, Carlsbad, California) according to Page 13 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 manufacture's instructions, and followed by cleanup To identify biological processes that were significantly using RNeasy spin columns (Qiagen, Valencia, Califor- changed with passages, the identified DE genes were ana- nia). RNA labelling was performed on the Onyx robotic lyzed using the DAVID functional annotation system system from MWG-Biotech and Aviso using the standard (National Institute of Allergy and Infectious Diseases, labelling protocol. Biotin labelled cRNA was fragmented NIH) [21]. Significantly changed biological processes and hybridized to human whole genome U133 plus 2 were determined by the enrichment test provided by the GeneChip (Affymetrix, Santa Clara, California). Chip DAVID system. To control for a false positive rate of test- processing, imagine capturing, and raw data analysis were ing the expression change of thousands of genes simulta- performed using the Affymetrix Microarray Suite MAS neously, the false discovery rate (fdrate, or FDR) was with default parameters. estimated using an algorithm derived from Benjamini and Hochberg [42]. Statistical analysis to identify differentially expressed genes Construct angiogenesis related interactome and calculate Hybridized microarrays were analyzed by the Affymetrix connectivity index A tentative interactome was constructed using PathwayAs- microarray analysis software microarray suite 5 (MAS 5). The chip signals were normalized using the default sist from the union of the differentially expressed genes method in Affymetrix MAS5, but setting the target at 2%- identified either from the HUVEC or the NHDF array trimmed mean to 1500 instead of 500. HUVEC were ana- experiment. PathwayAssist was configured to retrieve five lyzed with four passages. We designated the cell passages different types of relationships defined in the software, 5, 8, 12 and 22 as 1, 2, 3 and 4 respectively. For each namely, binding, expression, molecular transport, protein probeset, the signal data were fitted to a linear regression modification and regulation. The software introduced model of time to identify consistent gene expression non-differentially expressed genes into the interactome in changes over cell passage. Three batches of the NHDF cells order to establish the shortest connection pathway were analyzed, including early passage neonatal NHDF, between any two differentially expressed genes in the ini- early passage adult NHDF, and late passage adult NHDF tial input gene list. This is an informative function since it cells. To reduce a possible batch to batch effect, the MAS 5 could recover potentially non-transcriptional regulatory signals for each chip were re-normalized by Local Polyno- relationships between two or more genes. We cleaned up mial Regression Fitting (loess) approach using one neona- this initial interactome by removing genes that were not tal sample as the baseline on a log scale. The re- detected by either array study. In addition, genes that normalized signal data were fitted to an ANOVA model to could not connect to any other genes were called orphans identify differential gene expression changes between the and removed from the interactome. Thus, the final ver- different fibroblast cells. sion of the interactome we constructed can be described as follows: To control for a false positive rate of testing the expression change of thousands of genes simultaneously, the false Let G be the set of K genes in the interactome I, and R be discovery rate (fdrate, or FDR) was estimated using an one of the five relationships: binding, molecular trans- algorithm derived from Benjamini and Hochberg [42]. port, protein modification, inhibition or up-regulation. Probesets with a false discovery rate of 0.2 or smaller were Thus, the interactome I we built can be formally described considered as significant and followed up for further anal- as ysis. We excluded probesets called "absent" in all chips in a data set. These probesets usually have very low expres- I = {(g , g ) ∈ R | g ∈ G; i, j = 1,2, ... K}(1) i i sion signals. In addition, we used expression fold change between different treatment conditions to filter out genes We assumed R is transitive, that is, if gene A regulates gene to make sure the gene expression changes are biologically B, and gene B regulates gene C, it is assumed that gene A meaningful. regulates gene C. Except for binding among the five types of relationships, molecular transport, protein modifica- PCA analysis was performed in R using the principal com- tion, up-regulation and inhibition are unidirectional. ponent analysis function with standardized data [43]. The first three principal components with > 80% variations We wrote a "network walking" algorithm (Figure 6) in a were exported into Spotfire (Sommerville, MA) to aid in recursive function with which we calculated the number th th creating sample scatter plots. Hierarchical clustering anal- of unique genes, L , connecting to the i gene at the j ij ysis of DE genes was done in Spotfire with Euclidean dis- step along the connection path in the interactome. We th th tance by the complete linkage method. define connectivity of the i gene at the j step along the connection pathway as Page 14 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 j-1 C = C + L /2 (2) Authors' contributions i,j i,j-1 ij JS conceived the project. YC, TW, JMY, and JS designed the th study. YC, FL, and LY performed the experiments. TW con- Thus, for the i gene in the interactome its connectivity at th tributed to bioinformatics analyses and algorithm devel- the j step in its connection pathway is the connectivity at th opment. H–RQ performed the statistical analysis. TW, the (j-1) step plus weighted number of unique genes at th j-1 TPB, JJS, JMY, and JS interpreted the data. TW, JMY, YC, the j step. We assigned 1/2 as an arbitrary weight to the th and JS wrote the manuscript. All authors read and number of unique genes connecting to the i gene at the th approved the final manuscript. j step, that is, the weight exponentially decreases with the th number of steps away from the i gene. Obviously C is a i,j Additional material non-decreasing function. We further define connectivity th of the i gene C as Additional file 1 CC = lim This table lists the 1103 differentially expressed genes identified in i ij , (3) jN → HUVECs. th Click here for file where C is the plateau of connectivity of the i gene and [http://www.biomedcentral.com/content/supplementary/1471- N the number of steps C takes to arrive at its limit or the i i,j 2164-9-264-S1.xls] plateau. Since every gene in the interactome is connected, the total number of genes a gene connects to is the same, Additional file 2 th that is, after N steps the connectivity of the i gene will This table lists 787 non-redundant common differentially expressed genes arrive at its plateau. However, they differ in the number of identified in NHDF cells. steps N or the rate of the arrival at the plateau. Some genes Click here for file may take more steps than the others to arrive at the plat- [http://www.biomedcentral.com/content/supplementary/1471- 2164-9-264-S2.xls] ueas, depending upon what and how many genes it con- nects, especially at early steps. Schematic elucidation of Additional file 3 the algorithm is shown in Fig. 6. This file lists the detailed information about the interactome. Click here for file Simulation of network perturbation [http://www.biomedcentral.com/content/supplementary/1471- Let S be a small molecule kinase inhibitor, and P be the s, D 2164-9-264-S3.xls] level of network perturbation by the small molecule at a hypothetical dose D. P is estimated as s, D Acknowledgements We thank Mark Uhlik, Jeff Hanson, Shaoyu Chu and Jonathan Lee for help (4) PD =×C /25×IC 0 sD,, i i s with high content image analysis, and Xiaoling Xia for the HUVECs. We i =1 would also like to thank the functional genomics group for microarray anal- where K is the number of kinases in the interactome I that ysis, John Calley for in-house gene annotation, Gregory Donoho, John Cal- can be inhibited by the small molecule S, IC50 is the i,s ley, Patrick McGlynn, Bruce Konicek, and Dan Li for a critical reading of the th IC50 value of the small molecule S against the i kinase. manuscript and Michal Vieth and Laura Bloem for kinase data. We would Their relative IC50 values against various kinases were also like to extend our thanks to Karen Morgan for her assistance in pre- measured at Upsate (UBI) or in house. To simulate the paring the manuscript. dose response curve, we set the theoretical dose D from 0.1 nM to 2 μM. If a term in (4) is more than the maxi- References 1. Folkman J: Angiogenesis: an organizing principle for drug dis- mum number of the connected genes, the term is then set covery? Nat Rev Drug Discov 2007, 6(4):273-286. to the plateau. Summed network connectivity was used to 2. Carmeliet P: Angiogenesis in life, disease and medicine. Nature calculate the percent perturbation with the plateau to be 2005, 438(7070):932-936. 3. Hanahan D, Weinberg RA: The hallmarks of cancer. Cell 2000, set at 100% network perturbation. The resulted simulated 100(1):57-70. dose dependent network perturbation was then imported 4. Hanahan D, Folkman J: Patterns and emerging mechanisms of into the program GraphPad Prism to estimate EC50 of the angiogenic switch during tumorigenesis. Cell 1996, 86(3):353-364. network perturbation by the small molecule S by per- 5. Naumov GN, Akslen LA, Folkman J: Role of angiogenesis in forming Sigmoidal dose-response curve fitting. human tumor dormancy: animal models of the angiogenic switch. Cell Cycle 2006, 5(16):1779-1787. 6. Folkman J: Tumor angiogenesis: therapeutic implications. N Abbreviations Engl J Med 1971, 285(21):1182-1186. DE: differential expression; GO: Gene Ontology; HUVEC: 7. Ferrara N, Kerbel RS: Angiogenesis as a therapeutic target. Nature 2005, 438(7070):967-974. human umbilical vein endothelial cell; MAK: multi-tar- 8. Folkman J, Klagsbrun M: Angiogenic factors. Science 1987, geted anti-angiogenic kinases; NHDF: normal human der- 235(4787):442-447. mal fibroblast; PCA: principle component analysis. Page 15 of 16 (page number not for citation purposes) BMC Genomics 2008, 9:264 http://www.biomedcentral.com/1471-2164/9/264 9. Carmeliet P, Tessier-Lavigne M: Common mechanisms of nerve 33. D'Haeseleer P, Liang S, Somogyi R: Genetic network inference: and blood vessel wiring. Nature 2005, 436(7048):193-200. from co-expression clustering to reverse engineering. Bioin- 10. Jain RK: Molecular regulation of vessel maturation. Nat Med formatics 2000, 16(8):707-726. 2003, 9(6):685-693. 34. Hartemink AJ, Gifford DK, Jaakkola TS, Young RA: Combining 11. Kerbel RS: Tumor angiogenesis: past, present and the near location and expression data for principled discovery of future. Carcinogenesis 2000, 21(3):505-515. genetic regulatory network models. Pac Symp Biocomput 12. Kerbel RS: Therapeutic implications of intrinsic or induced 2002:437-449. angiogenic growth factor redundancy in tumors revealed. 35. Pe'er D, Regev A, Elidan G, Friedman N: Inferring subnetworks Cancer Cell 2005, 8(4):269-271. from perturbed expression profiles. Bioinformatics 2001, 17 13. Casanovas O, Hicklin DJ, Bergers G, Hanahan D: Drug resistance Suppl 1:S215-24. by evasion of antiangiogenic targeting of VEGF signaling in 36. Segal E, Shapira M, Regev A, Pe'er D, Botstein D, Koller D, Friedman late-stage pancreatic islet tumors. Cancer Cell 2005, N: Module networks: identifying regulatory modules and 8(4):299-309. their condition-specific regulators from gene expression 14. Pietras K, Hanahan D: A multitargeted, metronomic, and max- data. Nat Genet 2003, 34(2):166-176. imum-tolerated dose "chemo-switch" regimen is antiang- 37. Iossifov I, Krauthammer M, Friedman C, Hatzivassiloglou V, Bader JS, iogenic, producing objective responses and survival benefit White KP, Rzhetsky A: Probabilistic inference of molecular net- in a mouse model of cancer. J Clin Oncol 2005, 23(5):939-952. works from noisy data sources. Bioinformatics 2004, 15. Jaffe EA, Nachman RL, Becker CG, Minick CR: Culture of human 20(8):1205-1213. endothelial cells derived from umbilical veins. Identification 38. Benson M, Breitling R: Network theory to understand microar- by morphologic and immunologic criteria. J Clin Invest 1973, ray studies of complex diseases. Curr Mol Med 2006, 52(11):2745-2756. 6(6):695-701. 16. Bishop ET, Bell GT, Bloor S, Broom IJ, Hendry NF, Wheatley DN: An 39. Thomas R, Paredes CJ, Mehrotra S, Hatzimanikatis V, Papoutsakis ET: in vitro model of angiogenesis: basic features. Angiogenesis A model-based optimization framework for the inference of 1999, 3(4):335-344. regulatory interactions using time-course DNA microarray 17. Neubert K, Haberland A, Kruse I, Wirth M, Schimke I: The ratio of expression data. BMC Bioinformatics 2007, 8(1):228. formation of prostacyclin/thromboxane A2 in HUVEC 40. Shou J, Soriano R, Hayward SW, Cunha GR, Williams PM, Gao WQ: decreased in each subsequent passage. Prostaglandins 1997, Expression profiling of a human cell line model of prostatic 54(1):447-462. cancer reveals a direct involvement of interferon signaling in 18. Freedman DA, Folkman J: Maintenance of G1 checkpoint con- prostate tumor progression. Proc Natl Acad Sci U S A 2002, trols in telomerase-immortalized endothelial cells. Cell Cycle 99(5):2830-2835. 2004, 3(6):811-816. 41. Shou J, Bull CM, Li L, Qian HR, Wei T, Luo S, Perkins D, Solenberg PJ, 19. Freedman DA, Folkman J: CDK2 translational down-regulation Tan SL, Chen XY, Roehm NW, Wolos JA, Onyia JE: Identification during endothelial senescence. Exp Cell Res 2005, of blood biomarkers of rheumatoid arthritis by transcript 307(1):118-130. profiling of peripheral blood mononuclear cells from the rat 20. Erusalimsky JD, Kurz DJ: Endothelial cell senescence. Handb Exp collagen-induced arthritis model. Arthritis Res Ther 2006, Pharmacol 2006:213-248. 8(1):R28. 21. Sherman BT, Huang da W, Tan Q, Guo Y, Bour S, Liu D, Stephens R, 42. Benjamini Y, Hochberg: Controlling the false discovery rate:a Baseler MW, Lane HC, Lempicki RA: DAVID Knowledgebase: a practical and powerful approach to multiple testing. The Jour- gene-centered database integrating heterogeneous gene nal of Royal Statistical Society 1995, 57(1):289-300. annotation resources to facilitate high-throughput gene 43. Venables W, Ripley B: Modern Applied Statistics with S-PLUS. functional analysis. BMC Bioinformatics 2007, 8:426. New York, Springer-Verlag; 1997. 22. Jeong H, Tombor B, Albert R, Oltvai ZN, Barabasi AL: The large- scale organization of metabolic networks. Nature 2000, 407(6804):651-654. 23. Folkman J: Tumor suppression by p53 is mediated in part by the antiangiogenic activity of endostatin and tumstatin. Sci STKE 2006, 2006(354):pe35. 24. Giuriato S, Ryeom S, Fan AC, Bachireddy P, Lynch RC, Rioth MJ, van Riggelen J, Kopelman AM, Passegue E, Tang F, Folkman J, Felsher DW: Sustained regression of tumors upon MYC inactivation requires p53 or thrombospondin-1 to reverse the angiogenic switch. Proc Natl Acad Sci U S A 2006, 103(44):16266-16271. 25. Baka S, Clamp AR, Jayson GC: A review of the latest clinical compounds to inhibit VEGF in pathological angiogenesis. Expert Opin Ther Targets 2006, 10(6):867-876. 26. Cabebe E, Wakelee H: Sunitinib: a newly approved small-mole- cule inhibitor of angiogenesis. Drugs Today (Barc) 2006, 42(6):387-398. 27. Zakarija A, Soff G: Update on angiogenesis inhibitors. Curr Opin Oncol 2005, 17(6):578-583. 28. Abdollahi A, Schwager C, Kleeff J, Esposito I, Domhan S, Peschke P, Publish with Bio Med Central and every Hauser K, Hahnfeldt P, Hlatky L, Debus J, Peters JM, Friess H, Folk- scientist can read your work free of charge man J, Huber PE: Transcriptional network governing the ang- iogenic switch in human pancreatic cancer. Proc Natl Acad Sci "BioMed Central will be the most significant development for U S A 2007, 104(31):12890-12895. disseminating the results of biomedical researc h in our lifetime." 29. Sartore-Bianchi A, Ricotta R, Cerea G, Maugeri MR, Siena S: Ration- Sir Paul Nurse, Cancer Research UK ale and clinical results of multi-target treatments in oncol- ogy. Int J Biol Markers 2007, 22(1 Suppl 4):S77-S87. Your research papers will be: 30. Meyerhardt JA, Mayer RJ: Systemic therapy for colorectal can- available free of charge to the entire biomedical community cer. N Engl J Med 2005, 352(5):476-487. 31. Zimmermann GR, Lehar J, Keith CT: Multi-target therapeutics: peer reviewed and published immediately upon acceptance when the whole is greater than the sum of the parts. Drug Dis- cited in PubMed and archived on PubMed Central cov Today 2007, 12(1-2):34-42. 32. Marx J: Cancer. Encouraging results for second-generation yours — you keep the copyright antiangiogenesis drugs. Science 2005, 308(5726):1248-1249. BioMedcentral Submit your manuscript here: http://www.biomedcentral.com/info/publishing_adv.asp Page 16 of 16 (page number not for citation purposes)

Journal

BMC Genomics – Springer Journals

Published: Dec 1, 2008

Keywords: life sciences, general; microarrays; proteomics; animal genetics and genomics; microbial genetics and genomics; plant genetics and genomics

There are no references for this article.