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

Characterization of an in vitro angiogenesis co-culture model by multi parameter high content image analysis

HUVECs form a capillary structure when co-cultured with NHDF in vitro [16]. Unlike the HUVEC Matrigel™ assay, the HUVEC/NHDF cord formation co-culture allows us to study cell-cell interactions during long term (up to 14 days) angiogenesis. The formation of the cord structure is a highly organized morphogenic process. We developed a high content imaging based assay to document and quantitate cord formation using the Cellomics ArrayScan platform. We were able to simultaneously assess dozens of morphological measurements, such as tube length, area, connectivity of the cord, branching, and number of nuclei within the cord. Using this assay system, we tested a clinically approved antiangiogenic compound Sunitinib (SUTENT, SU11248). While 20 ng/ml VEGF significantly increased angiogenesis, 30 nM Sunitinib significantly inhibited angiogenesis (Fig. 1A).

To explore the phenotypic signatures induced by VEGF or Sunitinib, we performed principal component analysis (PCA) and hierarchical clustering (Fig. 1B &1C). Multivariate analyses using 47 morphological measurements clearly distinguished the three treatment groups, i.e. control, VEGF, and Sunitinib. In the PCA analysis, the first three principal components captured more than 80% of the total variance. Examining the PC loadings, we found that among the most discriminative features are the cord length (either connected or unconnected or in general), cord area, the branching segment count. The cord CD31 immunofluorescence intensity is also one of the features that significantly contributed to PC1, likely due to the overall correlation of fluorescence staining to cord mass. Interestingly, the cord width is not a discriminative feature for either VEGF or Sunitinib action. The current multi-parameter, high content analysis allowed for a detailed dissection of the phenotypic characteristics associated with a given factor or compound. Shown in Fig. 1D is the mean value plot of a few representative parameters. Sunitinib significantly inhibited the mass of the cord in 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 HUVECs in the co-culture. Sunitinib treated cultures are less morphologically complex, as the number of connected cords, the branching nodes, and segment measurements are all significantly reduced.

HUVEC cord formation competence is cell passage dependent

It has been well known that HUVECs change their biochemical and cellular behaviors when passaged in vitro [17–20]. We cultured HUVEC cells from a single batch to various passages (P6, P14, P22) and tested them for cord formation competence. We found that they gradually lost cord formation capabilities upon passage (Fig. 2A). While early passage (P6) HUVECs formed very nice capillary structures and responded robustly to VEGF when co-cultured with NHDF, late passage (P22) HUVECs lost cord formation capability and had diminished response to VEGF (Fig. 2A &2B).

Profiling passage dependent angiogenesis support by dermal fibroblast cells

Long term sustained angiogenesis in co-culture is supported by dermal fibroblast feeder layer cells, which are also affected by cell passage. While late passage (P12) adult NHDFs demonstrated poor support for HUVECs to form capillary structures, early passage (P5) adult cells did well (Fig. 3A), suggesting that supportiveness of NHDFs in the angiogenesis assay is indeed cell-passage dependent. An early passage NHDF from a neonatal donor also supported angiogenesis in this assay (Fig. 3A).

Establishment of a gene functional association networkThe fact that cell passage dependent gene expression in both HUVEC and NHDF cells influences angiogenesis supports modeling the angiogenesis network using the DEGs identified by these DNA microarray experiments. To understand the functional association of these genes, we constructed a molecular interactome among the total 2012 differentially expressed genes identified from either NHDFs or HUVECs. The interactome consisted of over 23,000 gene-to-gene connections retrieved from PubMed abstracts using the natural language processing algorithm in PathwayAssist. We then refined the interactome by removing genes whose expression was not detectable by microarray analysis in the co-culture system. The resulted molecular interactome consists of 17226 molecular relationships described by the following terms: binding, expression, molecular transport, protein modification and regulation (see Additional file 3). There were 1201 unique genes in the interactome, among which 887 genes were differentially expressed in the co-culture system. Thus, the interactome covers roughly 74% of the identified DEGs. The interactome displays the typical characteristics of a scale-free network [22] (Fig. 4). Thirty genes with the highest connectivity are listed in Table 3. We further identified 823 angiogenesis-related genes by conducting automated querying of the PubMed abstract database using an in-house developed text mining tool, TargetMiner. The thirty most connected genes in the interactome were all angiogenesis related, accounting for over 33% of the total 17226 connections in the interactome. This data suggests a high degree of relevancy of the interactome to angiogenesis.

Rank	Gene	Connectivity	Percentage
1	TGFB1	397	2
2	TP53	301	4
3	IGF1	290	6
4	MAPK1	270	7
5	FGF2	261	9
6	PDGF	251	10
7	AKT1	248	12
8	SRC	243	13
9	JUN	242	15
10	AGT	231	16
11	SP1	230	17
12	IL6	220	18
13	IL1B	216	20
14	NGFB	207	21
15	VEGF	183	22
16	MYC	173	23
17	MAPK8	169	24
18	HGF	168	25
19	RAF1	162	26
20	MAPK3	156	27
21	CREB1	146	28
22	EP300	141	28
23	EGFR	138	29
24	CDKN1A	122	30
25	RAC1	117	31
26	STAT3	115	31
27	BCL2	109	32
28	PRL	108	33
29	PPARG	103	33
30	NR3C1	101	34

Genes within the regulation network or interactome were ranked by 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 highest connectivity. Gene symbols are those assigned by Human Gene Nomenclature Committee (HGNC).

Subnetwork analysis of oncogene involvement in angiogenesis

The information on which the network was based comes from many cell types. Although we filtered and refined the network by only including genes that are present in the model, validating 17226 individual network edges is challenging. We explored the feasibility of using the present network model to address known biological phenomena as a pathway toward functional validation of the network. Examining the connectivity hierarchy, we found that oncogenes and tumor suppressor genes are among the highest ranked molecules (hub nodes) based on the number of connections in which they are directly engaged. For example, the number of connections that P53, SRC, and MYC directly engaged is 301, 243, and 173 respectively, ranking at the 2^nd, 8^th and 16^th places in terms of connectivity in the interactome. This observation recapitulates the emerging role of oncogenes and tumor suppressor genes in angiogenesis regulation [1, 5, 23, 24]. We focused our analysis on MYC by first dissecting its interactions with other genes in the network. A subnetwork centered on MYC was constructed from the initial interactome (Fig. 5A). It was found that MYC is linked to a variety of angiogenesis regulators with high connectivity within the network, suggesting a central role of MYC in angiogenesis regulation. We hypothesized that inactivation of MYC function would have profound angiogenic defects in this co-culture assay. We tested this hypothesis by attenuating MYC expression in HUVECs using c-MYC siRNA. Cord formation was significantly suppressed in siMYC treated cultures (Fig. 5B &5C). Cord mass, branching pattern, and responsiveness to exogenous VEGF were all inhibited, which is consistent with the network prediction. Knockdown of MYC expression by the siRNA was confirmed in parallel by TaqMan RT-PCR analysis (Fig. 5D). The involvement of TGF-β, the top node in the present angiogenesis network was also validated (data not shown).

Assessment of compound antiangiogenic activity by network perturbation analysis

We further explored the application of network analysis to understand compound polypharmacology. Since kinases have increasingly become effective targets for anti-angiogenic drug discovery [25–27], we assessed kinase contributions to angiogenesis regulation in the co-culture model by network analysis. First, 39 kinases were mapped to the interactome with a total of 525 directly-connected genes, accounting for 44% of the total number of genes in the refined interactome. Any mapped kinase is engaged in the rest of the network through multi-step connections. Accordingly, we developed a "network walking" algorithm to model the connection cascade of kinases within the network. The methods are described in the method section and shown in a schematic diagram in Fig. 6. Using this algorithm, the characteristic connection curves of the 39 mapped kinases were determined and are shown in Fig. 7A. All kinases were able to connect to the rest of the interactome (except for singletons), and the connectivity plateaus at the connection of 1040 within the interactome. This property was used to assess the potential maximal perturbation of the interactome if a given kinase is inhibited. While every connected molecule plateaus at the same number of connections, each arrives at a different rate, which can be assessed by the slope of the connectivity curve. Since the spread of activity along the network can not reach 100 percent penetration at each step, and the penetration is expected to decay dramatically at later steps due to multiple inputs and outputs at each node, we thus defined a perturbation index calculated from the data collected by the network walking algorithm using a weighted sum method as described in the method section.

Using network perturbation analysis to assess drug phenotypic activity, we tested three multi-targeted anti-angiogenic kinase (MAK) inhibitors: Sunitinib, PTK787, and LY2401401. The analysis was performed based on the number and the identities of the kinases that are inhibited by a given compound at its respective potency as described in the method section. Sunitinib and PTK787 have been previously tested in clinical trials [25], while LY2401401, is a novel Lilly MAK inhibitor with a distinct kinase activity profile. LY2401401 (Fig. 7B), an imidazopyridine, is an orally bioavailable, broad-spectrum inhibitor of the VEGF and PDGF receptors. Additionally, Flt3, Kit, Tie-2 and the Eph family of receptors are also sensitive to this small molecule inhibitor (SMI). To differentiate these MAK compounds, it is important to determine if the unique kinase profile of LY2401401 is advantageous to achieve anti-angiogenic efficacy. We simulated the potential network perturbation by the three molecules based on the integrated perturbation index of every kinase that a given SMI is active against. The resulting theoretical dose response curves were used to estimate the simulated EC50 of network perturbation (Fig. 7C). The compounds were ranked by their simulated activities as LY2401401>Sunitinib>PTK787, with the EC50 values estimated to be 7.8 nM, 13.6 nM, and 133.7 nM, respectively. The network analysis clearly differentiates these three MAK molecules. Interestingly, the network analysis predicted that Sunitinib would be far more potent than PTK-787, although their activities against KDR are rather close (49 nM, and 58 nM, respectively for PTK-787 and Sunitinib).

We then determined the actual angiogenesis activity of each compound by testing its dose dependent response directly in the cord formation assay. They were ranked in the same order as inferred by the network simulation (Fig. 7D). The observed EC50 values from multiple independent experiments were determined to be 9.2 nM (± 9.1 nM s.d.; n = 3), 39.0 nM (± 9.2 nM s.d.; n = 4), and 140.4 nM (± 53.7 nM s.d.; n = 2) for LY2401401, Sunitinib, and PTK787, respectively (Fig. 7D). The actual compound activity in the cord formation bioassay was consistent with the predicted values based on the network perturbation analyses.

Reagents and Cells

Recombinant VEGF was purchased from R&D Systems (Minneapolis, Minnesota). Human umbilical vein endothelial cells (HUVECs) and Normal Human Dermal Fibroblast (NHDF) cells and culture medium were purchased from Cambrex (Walkersville, Maryland). Cells were maintained in EGM with 10% FBS and FGM-2, respectively for HUVECs and NHDFs, according to manufacture's instructions. In vitro angiogenesis co-culture AngioKit and Optimized medium were purchased from TCScellworks (Buckingham, UK).

In vitro angiogenesis assay

In vitro HUVEC cord formation co-culture was performed using the AngioKit commercially available from TCScellworks (Buckingham, UK) adapted to a 24-well format. To test HUVEC and NHDF passage dependence effects, the co-culture was done in house adapted to a 96-well format. Briefly, HUVEC cells were co-cultured onto human dermal fibroblast feeder layers in Optimised medium (TCScellworks, Buckingham, UK) for a total of about 12 days. Treatment started the day when the plates were received; usually three to four days after initial plating. Medium with or without growth factors or compounds was replenished every 2 to 3 days. Cells were fixed with 70% cold ethanol at the end of the culture and processed for anti-CD31/Alexa 488 immunofluorescence. To visualize feeder layer fibroblasts, some cultures were stained using a polyclonal antibody against Vimentin (Abcam, Cambridge, Massachusetts) along with an Alexa 568 conjugated secondary antibody to examine the feeder cell integrity. Nuclear staining was performed using Hoechst staining (Invitrogen Molecular Probe, Carlsbad, California).

Immunofluorescence stained co-culture plates were scanned using ArrayScan VTI (Cellomics, Pittsburgh, Pennsylvania) and the multi-parameter tube formation Bioapplication was used for quantification. Morphological measurements included cord length, width, area, percent connected cord, branching node and segment count, nuclei per cord, angiogenic index (expressed as the percentage of area covered by connected tube), etc. General cytotoxicity was monitored by examination of the nuclear staining of the feed layers.

Attenuation of gene expression by siRNA

HUVECs were trypsinized by 0.05% Trypsin/EDTA (Invitrogen, Carlsbad, CA), and 1 × 10⁶ cells were nucleofected with siMYC (200 pmol each, Ambion, Austin, TX) using a nucleofector device (Amaxa, Gaithersburg, MD) adopting optimized A-034 protocols for HUVECs according to the 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 formation assay. Cells were also plated in parallel into 6-well plates and collected in Trizol (Invitrogen, Carlsbad, CA) 48 h or 72 h later for total RNA extraction.

Real Time RT-PCR

Total RNAs were extracted from Trizol and cleaned by RNeasy mini kit (Qiagen, Chatsworth, CA). Genomic DNA contamination was removed by DNA-free kit (Ambion, Austin, TX), and cDNA synthesis was done by using a high capacity cDNA archive kit (Applied Biosystems, Forster City, CA). The cDNAs were used as templates for real-time PCR with a universal PCR master mix (Applied Biosystems, Foster City, CA). Real time PCR was performed on a ABI 7900 HT. The average values were normalized to 18s rRNA. Delta CT method was applied to calculate the relative gene expression level [40].

RNA isolation and microarray experiments

RNA isolation and microarray studies were performed as described previously [41]. Briefly, HUVECs and NHDFs with different passages were cultured to about 80% confluency. RNA was isolated using Trizol™ (Life Technologies, Inc. Invitrogen, Carlsbad, California) according to manufacture's instructions, and followed by cleanup using RNeasy spin columns (Qiagen, Valencia, California). RNA labelling was performed on the Onyx robotic system from MWG-Biotech and Aviso using the standard labelling protocol. Biotin labelled cRNA was fragmented and hybridized to human whole genome U133 plus 2 GeneChip (Affymetrix, Santa Clara, California). Chip processing, imagine capturing, and raw data analysis were performed using the Affymetrix Microarray Suite MAS with default parameters.

Statistical analysis to identify differentially expressed genes

Hybridized microarrays were analyzed by the Affymetrix microarray analysis software microarray suite 5 (MAS 5). The chip signals were normalized using the default method in Affymetrix MAS5, but setting the target at 2%-trimmed mean to 1500 instead of 500. HUVEC were analyzed with four passages. We designated the cell passages 5, 8, 12 and 22 as 1, 2, 3 and 4 respectively. For each probeset, the signal data were fitted to a linear regression model of time to identify consistent gene expression changes over cell passage. Three batches of the NHDF cells were analyzed, including early passage neonatal NHDF, early passage adult NHDF, and late passage adult NHDF cells. To reduce a possible batch to batch effect, the MAS 5 signals for each chip were re-normalized by Local Polynomial Regression Fitting (loess) approach using one neonatal sample as the baseline on a log scale. The re-normalized signal data were fitted to an ANOVA model to identify differential gene expression changes between the different fibroblast cells.

To control for a false positive rate of testing the expression change of thousands of genes simultaneously, the false discovery rate (fdrate, or FDR) was estimated using an algorithm derived from Benjamini and Hochberg [42]. Probesets with a false discovery rate of 0.2 or smaller were considered as significant and followed up for further analysis. We excluded probesets called "absent" in all chips in a data set. These probesets usually have very low expression signals. In addition, we used expression fold change between different treatment conditions to filter out genes to make sure the gene expression changes are biologically meaningful.

PCA analysis was performed in R using the principal component analysis function with standardized data [43]. The first three principal components with > 80% variations were exported into Spotfire (Sommerville, MA) to aid in creating sample scatter plots. Hierarchical clustering analysis of DE genes was done in Spotfire with Euclidean distance by the complete linkage method.

To identify biological processes that were significantly changed with passages, the identified DE genes were analyzed using the DAVID functional annotation system (National Institute of Allergy and Infectious Diseases, NIH) [21]. Significantly changed biological processes were determined by the enrichment test provided by the DAVID system. To control for a false positive rate of testing the expression change of thousands of genes simultaneously, the false discovery rate (fdrate, or FDR) was estimated using an algorithm derived from Benjamini and Hochberg [42].

Construct angiogenesis related interactome and calculate connectivity index

A tentative interactome was constructed using PathwayAssist from the union of the differentially expressed genes identified either from the HUVEC or the NHDF array experiment. PathwayAssist was configured to retrieve five different types of relationships defined in the software, namely, binding, expression, molecular transport, protein modification and regulation. The software introduced non-differentially expressed genes into the interactome in order to establish the shortest connection pathway between any two differentially expressed genes in the initial input gene list. This is an informative function since it could recover potentially non-transcriptional regulatory relationships between two or more genes. We cleaned up this initial interactome by removing genes that were not detected by either array study. In addition, genes that could not connect to any other genes were called orphans and removed from the interactome. Thus, the final version of the interactome we constructed can be described as follows:

Let G be the set of K genes in the interactome I, and R be one of the five relationships: binding, molecular transport, protein modification, inhibition or up-regulation. Thus, the interactome I we built can be formally described as

I = {(g _i , g _i ) ∈ R | g ∈ G; i, j = 1,2, ... K} (1)

We assumed R is transitive, that is, if gene A regulates gene B, and gene B regulates gene C, it is assumed that gene A regulates gene C. Except for binding among the five types of relationships, molecular transport, protein modification, up-regulation and inhibition are unidirectional.

We wrote a "network walking" algorithm (Figure 6) in a recursive function with which we calculated the number of unique genes, L _ij , connecting to the i^th gene at the j^th step along the connection path in the interactome. We define connectivity of the i^th gene at the j^th step along the connection pathway as

C_i,j= C_i,j-1+ L _ij /2^j-1 (2)

where C _i is the plateau of connectivity of the i^th gene and N _i the number of steps C_i,jtakes to arrive at its limit or the plateau. Since every gene in the interactome is connected, the total number of genes a gene connects to is the same, that is, after N _i steps the connectivity of the i^th gene will arrive at its plateau. However, they differ in the number of steps N _i or the rate of the arrival at the plateau. Some genes may take more steps than the others to arrive at the platueas, depending upon what and how many genes it connects, especially at early steps. Schematic elucidation of the algorithm is shown in Fig. 6.

Simulation of network perturbation

where K is the number of kinases in the interactome I that can be inhibited by the small molecule S, IC 50_i,sis the IC50 value of the small molecule S against the i^th kinase. Their relative IC50 values against various kinases were measured at Upsate (UBI) or in house. To simulate the 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 maximum number of the connected genes, the term is then set to the plateau. Summed network connectivity was used to calculate the percent perturbation with the plateau to be set at 100% network perturbation. The resulted simulated dose dependent network perturbation was then imported into the program GraphPad Prism to estimate EC50 of network perturbation by the small molecule S by performing Sigmoidal dose-response curve fitting.