As a proof of concept analysis using our CMN model, we chose to use the data from a recent study utilizing mouse xenograft models to study acquired resistance of human lung adenocarcinomas to angiogenesis inhibitors (Cascone et al., 2011
). The authors in the aforementioned study injected H1975 cells (human non-small cell lung cancer cell line that is sensitive to bevacizumab) into mice and after a while started treating them with a vehicle control or bevacizumab. Despite initial response, majority of tumors developed resistance to the drug. Gene expression microarray analyses were then performed on control and bevacizumab-resistant tumors (n
=3 for each) using Illumina mouse-specific (WG-6 v2) and human-specific (WG-6 v3) arrays. Mouse-specific arrays were used to measure stromal gene expression and the human array was used to measure the tumor gene expression.
In order to identify community-wide molecular networks involved in the acquired resistance to bevacizumab, we first performed t-tests on each gene in the stromal and tumor datasets to score the extent of their change in resistant versus control cases. The resultant t-values (the value of t-statistic from t-tests) for tumor and stromal gene expression data had different distributions; stromal cells had a wider t-value distribution (not shown), consistent with the original study's report that stromal gene expression changes were more pronounced than that in tumor data. Therefore, we transformed the two t-value distributions using quantile normalization to obtain identical distributions for tumor and stromal t-values. These t-values (two sets, one for stromal and one for tumor cells) were used as node weights for our NetWalk-based network analyses ().
Fig. 2. Workflow of the method in this study. First, t-tests were performed for each gene to score the difference of their expression between control and bevacizumab-resistant samples in tumor and stroma datasets, respectively. CMN was created for two cell types, (more ...)
3.1 NetWalk analysis of CMN
First, we constructed a CMN model of two cell types: tumor (T
) and stroma (S
). This CMN is defined as
is the adjacency matrix of this CMN, T
is the adjacency matrix for tumor and S
is the adjacency matrix for stromal intracellular networks (note that T
are identical, we denoted them separately for clarity), rt:s
is the adjacency matrix for intercellular interactions between tumor and stromal cells defined by rules in ‘Methods’ section. This model is composed of 26 512 genes (13 256 genes/cell type) connected by 375 570 interactions (all intra- and intercellular interactions). By overlaying both t
-value distributions onto the respective nodes on this CMN, we performed NetWalk analysis to obtain numeric scores of relevance for each interaction in this network (i.e. EF values). First, in order to understand the performance of our approach with respect to identifying intra- and intercellular networks, we tested the composition of increasing numbers of EF values for intra- and intercellular interactions. A shows the compositions of draws of 50–5000 highest EF values (A). Relative compositions of the individual intracellular interactions are similar, with intercellular interactions comprising <10% of all interactions, as expected. Therefore, NetWalk does not bias the scoring of interactions to one or the other cell type.
Fig. 3. NetWalk analysis of CMN. (A) Fraction of intracellular interactions in cells 1 and 2 (i.e. tumor and stroma) and intercellular interactions within the indicated number of highest scoring EF values from the CMN of two cell types. (B) Boxplot of t-value (more ...)
We have shown that networks extracted from EF values are more faithful to the underlying data than other network extraction algorithms (i.e. high scoring networks contain genes with high data values and vice versa) (Komurov et al., 2010
). We wanted to test if the EF values from NetWalk analysis of the CMN are also faithful to the underlying data distribution from the two cell types. B shows the t
-value distributions of genes from each cell type within the networks of 50–5000 highest EF values. The boxplots show that even in the network of 5000 highest scoring interactions from this NetWalk analysis, the distribution of t
-values of genes in both cell types is strictly at the higher end of the individual t
-value distributions. We have also conducted additional randomization simulations to verify the strict coherence of NetWalk output within the CMN context to the input data (see Supplementary Fig. S2
and legend). This indicates that the CMN subnetworks extracted from NetWalk analysis are highly related to the input data, so that high scoring networks are composed of genes with high t
-values, rather than unrelated genes. This feature of NetWalk analysis increases confidence in the relevance of the extracted networks to the underlying data.
3.2 Analysis of individual cell networks
First, we analyzed the individual intracellular networks by visualizing networks corresponding to highest EF values from each cell type. An individual analysis of the network from highest EF values in tumor cells shows an upregulation of networks involved in survival signaling centered around phosphoinositide 3 kinase (PI3K), chemokine signaling, TGF-β and apoptosis signaling, among others (C). PI3K in this network seems to be activated by multiple mechanisms, including Toll-like receptor (TLR2), insulin receptor substrate (IRS2) and G-protein signaling. Stromal cells, however, are characterized by upregulation of different networks, primarily in the extracellular matrix (ECM) remodeling, epidermal growth factor receptor (EGFR) signaling, oxidative energy metabolism and others (D). Upregulation of EGFR in stromal cells was also reported in the original study, suggesting that our approach, in addition to revealing novel associations, can reproduce previous findings. Although revealing potentially important information on the respective intracellular network interactions involved in acquired drug resistance, individual network analyses do not allow for analyses of interactions between these two networks.
3.3 Analysis of the community networks of tumor and stromal cells
In order to identify community-wide molecular networks involved in drug resistance in this experimental model, we visualized the network corresponding to highest scoring 400 EF values (A). As would be expected from the total percentage of intercellular interactions in the CMN (~6%), most of the network is composed of intracellular interactions of the respective cell types. There are however a few potentially important interactions occurring between stromal and tumor cell gene products (highlighted in A). This subnetwork contains both physical and metabolic intercellular interactions.
Fig. 4. CMN analysis of tumor and stromal t-values. (A) Highest scoring interactions from the CMN analysis of tumor and stromal t-values. Nodes are colored by the respective t-values of tumor and stromal cells as shown in the color key. Interactions are colored (more ...) 3.3.1 Physical interactions
Some notable examples for physical interactions include HMGB1 secreted by stromal cells acting on the TLR2 and PLG (plasminogen) in tumor cells to activate them (B), NID1–LAMC2 (nidogen and laminin C2) involved in the ECM remodeling, fibroblast growth factor signaling (tumor FGF8–stromal FGFR2), tumor cell bone morphogenetic protein (BMP) proteins acting on the BMP receptor on stromal cells, chemokine signaling between stromal and tumor cells, and neuropeptide signaling (neuropeptide (NPY) from stromal cells activating NPY5R, the receptor, on tumor cells). Upregulation of FGFR2 in stromal cells was reported in the original study of Cascone et al., again showing that our analysis can capture findings from traditional approaches. However in addition, this analysis shows multiple potentially important intercellular physical interactions associated with drug resistance. Importantly, most of these intercellular physical interactions have been described in the literature to play roles in different aspects of cancer progression.
3.3.2 Metabolic interactions
In addition to physical interactions, the subnetwork in A contains several metabolic intercellular interactions. A particularly intriguing set of interactions is shown separately in C, which contains several enzymes in tumor and stromal cells involved in energy metabolism. A sketch of the reactions performed by these enzymes in the stromal and tumor cells is shown in D. Stromal cells have a significant upregulation (based on t-value) of the enzymes lactate dehydrogenase B (LDHB), dihydrolipoamide dehydrogenase (DLD), glutamic pyruvate transaminase (GPT2), pyruvate carboxylase (PC) and asparagine synthetase (ASNS), while tumor cells have upregulated lactate dehydrogenase C (LDHC) and asparaginase like 1 (ASRGL1). While other lactate dehydrogenase enzymes primarily catalyze conversion of pyruvate to lactate as the terminal reaction of glycolysis, LDHB mainly converts lactate back to pyruvate. DLD is a component of several enzyme complexes, including the pyruvate dehydrogenase complex, which catalyzes conversion of pyruvate to acetyl CoA, which enters the tricarboxylic acid cycle (TCA). Similarly, PC catalyzes conversion of pyruvate to oxaloacetate, which also enters the TCA cycle. Enzyme GPT2 in stromal cells catalyzes the reversible reaction pyruvate + glutamate ↔ alanine + 2-oxoglutarate, which is one of the important reactions in gluconeogenesis. ASNS and ASRGL1 catalyze, respectively, asparagine synthesis producing glutamate from glutamine, and breakdown of asparagine to aspartate producing NH3.
It is obvious from this picture, that stromal cells have a significant upregulation of the enzymes involved in oxidative energy generation (LDHB, DLD, PC) and gluconeogenesis (GPT2), which is also evident from the significant upregulation of the mitochondrial electron transport chain components (A). Tumor cells, however, seem to mainly perform glycolysis, as evidenced from preferential expression of LDHC and glycolytic enzymes hexokinase I (HK1), GPI and the glucose transporter SLC2A10 (see the highlighted subnetwork 9 in A). Three metabolites, lactate, asparagine and aspartate mediate the metabolic interactions between the tumor and stromal cells in this scenario. While the functional significance of the interaction of ASNS in stromal cells and ASRGL1 in tumor cells potentially involving an aspartate/asparagine shuttle may not be immediately clear, the interaction of LDHC with LDHB and other stromal enzymes involved in oxidative energy metabolism is particularly intriguing. Based on the above observations, a hypothesis emerges where tumor cells primarily consume glucose by glycolysis and secreting lactate, which is then uptaken by stromal cells and oxidized through mitochondrial TCA cycle and/or used for gluconeogenesis to produce glucose or other building blocks (e.g. by GPT2) that may be used again by tumor cells. This hypothesis is not far-fetched, as a similar metabolic symbiosis involving lactate shuttling between stromal and tumor cells has been reported and has been implicated with increased tumorigenicity in several cancer models (Koukourakis et al., 2006
; Pavlides et al., 2009
; Sonveaux et al., 2008
). Therefore, the metabolic interactions identified in this study may underlie bona fide
cross-talk mechanisms between tumor and stromal cells that play roles in meeting the metabolic demand of the tumorigenic and drug resistance phenotypes.