Systems Analysis of Insulin and IGF1 Receptors Networks in Breast Cancer Cells Identifies Commonalities and Divergences in Expression Patterns

Insulin and insulin-like growth factor-1 (IGF1), acting respectively via the insulin (INSR) and IGF1 (IGF1R) receptors, play key developmental and metabolic roles throughout life. In addition, both signaling pathways fulfill important roles in cancer initiation and progression. The present study was aimed at identifying mechanistic differences between INSR and IGF1R using a recently developed bioinformatics tool, the Biological Network Simulator (BioNSi). This application allows to import and merge multiple pathways and interaction information from the KEGG database into a single network representation. The BioNsi network simulation tool allowed us to exploit the availability of gene expression data derived from breast cancer cell lines with specific disruptions of the INSR or IGF1R genes in order to investigate potential differences in protein expression that might be linked to biological attributes of the specific receptor networks. Modeling-generated information was corroborated by experimental and biological assays. BioNSi analyses revealed that the expression of 75 and 71 genes changed during simulation of IGF1R-KD and INSR-KD, compared to control cells, respectively. Out of 16 proteins that BioNSi analysis was based on, validated by Western blotting, nine were shown to be involved in DNA repair, eight in cell cycle checkpoints, six in proliferation, eight in apoptosis, seven in oxidative stress, six in cell migration, two in energy homeostasis, and three in senescence. Taken together, analyses identified a number of commonalities and, most importantly, dissimilarities between the IGF1R and INSR pathways that might help explain the basis for the biological differences between these networks.

Insulin and IGF1, acting respectively via the INSR and IGF1R, display important metabolic and developmental roles at each period of life (6)(7)(8). While INSR and IGF1R share the majority of their downstream cytoplasmic mediators, clinical, and experimental data indicates that INSR activation leads predominantly to metabolic activities. On the other hand, IGF1R activation precedes growth events (9,10). The IGF2R resembles the mannose 6-phosphate receptor, a membrane protein involved in the recycling of lysosomal enzymes (11). The ability of insulin and IGF1 to bind and activate the opposite receptor has been widely reported in the scientific literature (12). Crossactivation usually takes place at elevated concentrations of the hormone (at least one order of magnitude higher doses than those required by the high affinity-binding ligand to activate its cognate receptor) (13). Thus, high ambient values of insulin can activate the IGF1R, leading to proliferative activities, while excessive IGF1 can stimulate the INSR, leading to metabolic actions (5,9). Pathological dysregulation of the IGF system is linked to a number of conditions, ranging from growth deficits to cancer development (14)(15)(16). The involvement of the INSR and IGF1R in breast cancer has been extensively reported (17)(18)(19)(20). However, the complexity of this hormonal network led to conflicting results regarding the relative impact of the various players (i.e., ligands, receptors, etc.) in malignancy (21)(22)(23). Furthermore, this complexity translated into a disappointingly slow pace in the development of INSR/IGF1R-directed therapies (24)(25)(26). Genome-wide analyses of the mechanisms of action of insulin and IGF1 in breast cancer, as well as identification of the signaling pathways involved, is expected to be of major translational relevance.
Modeling and simulation of regulatory networks became an integral part of biological research (27). Network simulation tools are aimed at elucidating the interactions between genes, proteins and signaling pathways and, furthermore, are designed to provide new insights into complex biological questions. Simulation tools allow the visualization and analysis of the dynamics of local biological networks. The Biological Network Simulator (BioNSi) tool is a computational application that allows to import and merge multiple pathways and interaction information from the KEGG database into a single network representation (http://bionsi.wix.com/bionsi) (28,29). Moreover, BioNSi enables the upload of expression values from highthroughput experiments and simulate the gene or protein interactions acquired from KEGG accordingly. Data generated may facilitate the investigation of the entire dynamic network under different biological conditions (30,31). Given the vital roles of the INSR and IGF1R signaling pathways in breast cancer, and to gain new insight into potential commonalities and divergences in expression patterns between both receptors and their downstream mediators, we employed the BioNSi bioinformatics tool for network simulation. MCF7 breast cancerderived cells with specific disruption of the INSR or IGF1R were used to study gene expression. BioNSi bioinformatics analyses integrated and simulated nine KEGG pathways, including apoptosis, cell cycle regulation, PI3K/Akt signaling, senescence and others, to a single network, based on specific gene expression. Our focus on these biological pathways stems from the fact that an important body of work over the past years revealed critical differences between the INSR and IGF1R pathways in these paths. Modeling-generated information was corroborated by experimental analyses.

Microarray Experiment
MCF7 cells were seeded in 10-cm Petri dishes until reaching 100% confluence. The cells were washed with phosphate-buffered saline (PBS) and total RNA was isolated using the TRIzol reagent (InVitrogen, Waltham, MA, USA), according to the RNA extraction guidelines for Affimetrix GeneChip Assays. Affymetrix R Human Gene 2.1 ST Array Strip (#902114, Santa Clara, CA, USA) was used for gene expression analysis according to manufacturer's instruction manual (www.affymetrix.com).
Partek Genomics Suite was used for analysis (www.partek. com/partek-genomics-suite/). The untreated MCF7 cells gene expression data is available as Supplementary Data Sheet 1.

Cell Cycle Analysis
Cells were seeded in duplicate onto 6-well-plates (10 4 cells/well) for 24 h. Cells were then serum-starved for an additional 24 h and incubated in the presence or absence of IGF1 or insulin (50 ng/ml) for 72 h. After incubation, cells were washed with PBS, trypsinized, and pelleted by centrifugation. The cells were permeabilized with Triton-X100, after which propidium iodide was added. Stained cells were analyzed using a FacsCalibur system (Cytek Development Inc., Fremont, CA, USA).

Cell Senescence Assays
Cells were seeded in 24-well-plates (5 × 10 4 cells/well) for 24 h. Cells were then serum-starved for an additional 24 h and incubated in the presence or absence of etoposide (10 µM) for 48 h. After incubation, cells were washed with PBS, fixed with 2% formaldehyde/0.2% glutaraldehyde for 5 min and incubated overnight at 37 • C in staining solution [40 mM sodium citrate, pH 6.0, 1% 5-bromo-4-chloro-3-indolyl-β-D-galactopyranoside (X-gal), 5 mM potassium ferrocyanide, 5 mM ferricyanide, 150 mM sodium chloride, and 2 mM magnesium chloride] (33). Cultures were examined under phase-contrast microscopy. The development of a blue color in the cytoplasm was detected and photographed using an inverted microscope (Nikon Instruments Inc., Lewisville, TX, USA). Four pictures were taken of each well. β-galactosidase-stained and unstained cells were counted and used to calculate a final average ratio for each well from quadruplicate samples.

Bioinformatics Network Analyses
Nine human KEGG pathways (based on 16 key genes) ( Table 1) were imported to BioNSi tool to form a single network for simulation analysis. The rationale for this choice of genes is the fact that they are regarded as prototypical targets of insulin and/or IGF1 action (5,9,10,16).
The BioNSi global network that includes all of these nine pathways contains 385 nodes (genes and small molecules) that are connected by 631 edges, in addition to gene self-inhibition edges that represent biological degradation (30). KEGG pathway information includes activation and inhibition edges. BioNSi settings were: activation edges were set to +2, inhibition to −1, phosphorylation to +2, dephosphorylation to −1, ubiquitination to −1, binding/association to +2, dissociation to −1, and selfinhibition to −1. In order to simulate the effect of activation and inhibition between nodes over time, untreated MCF7 gene expression microarray results were imported to BioNSi analysis as baseline (GSE145787; described above). The average expression values for replicate control cells were used as initial values for the nodes in the network and were normalized by the BioNSi tool to be between 0 and 9 (grayscale). Three different BioNSi simulation analyses were compared as follows (all simulations were run for 100 steps as shown in Figure 1):

Control simulation, in which IGF1 and insulin initial
expression values were manually set to 20 (very high), instead

Statistical Analyses
The statistical significance of the differences between groups was assessed by Student's t-test (two samples, equal variance). Scanning densitometry analyses were evaluated using the TINA imaging analysis software (http://biochemlabsolutions. com/GelQuantNET.html). Signal intensities of proteins were normalized to the corresponding Hsp70 signals. Data are presented as mean ± SEM of three independent experiments. P < 0.05 or 0.01 were considered statistically significant.

BioNSi Analysis of MCF7-Derived INSR-KD and IGF1R-KD Breast Cancer Cell Lines
BioNSi analyses revealed that the expression of 75 and 71 genes changed during simulation of IGF1R-KD (simulation 2) and INSR-KD (simulation 3), compared to control cells (simulation 1), respectively. As it is confusing to observe the complete network of 385 genes (see cytoscape file in Supplementary Data Sheet 2), Figure 2A shows a reduced network of 68 genes. However, the simulation analyses were performed for the complete network of 385 genes. The reduced network (68 genes) presents 16 genes that were biologically tested

Western Immunoblotting Validation of BioNSi Simulations
In order to validate the BioNSi simulation results, changes in protein expression levels in IGF1R-KD and INSR-KD cells were measured by Western blots and quantified by TINA software from new independent protein samples. IGF1R-KD reduced protein expression of six genes (CCND1, ATM, JAK1, STAT3, SOD2, and CHEK2), compared to INSR-KD and control cells (Figures 5A,D,G,J,N,Q, respectively). However, protein expression of eight genes [HRas, TP53, AKT3, CASP3, MTOR,  (20), and genes are highlighted in turquoise frame. Complexes initial expression is zero by BioNSi default. Activation arrow edges are colored green. Inhibition edges are colored red with flat heads. Simulation analyses were performed on the complete network (385 genes). Node's background was changed after simulations: genes whose expression changed during simulation in both KD cells are colored orange; genes whose expression changed during simulation only in IGF1R-KD cells are colored purple; and genes whose expression changed during simulation only in INSR-KD cells are colored green. (B) Venn diagram of numbers of genes whose expression changed during simulation (specific KD vs. control). Nodes whose expression changed only in IGF1R-KD cells are colored purple (IGF1R, STAT3, GAB1, JAK1, PLCG1, SRC, and diacylglycerol). Genes whose expression changed during simulation only in INSR-KD are colored green (INSR, SH2B2, and CBLC). Overlapping genes whose expression changed during simulation in both KD cells are colored orange. PPKAA1 (AmpK3), p16, and p21] was not markedly changed ( Figures 4A,C Out of 16 proteins that BioNSi analysis was based on, validated by Western blotting, nine were shown to be involved in DNA repair, eight in cell cycle checkpoints, six in proliferation, eight in apoptosis, seven in oxidative stress, six in cell migration, two in energy homeostasis, and three in senescence ( Table 2).

Cell Cycle Analysis
Next, experiments were carried out to characterize IGF1 or insulin effects on cell cycle progression. To this end, IGF1R-KD, INSR-KD and control MCF7 cells were treated with IGF1 or insulin for 72 h, after which flow cytometry was performed on propidium iodide-stained cells. In control cells, both IGF1 and insulin decreased the proportion of cells at the SubG1 phase (from 4.3 to 1.35% or 1.65%, respectively) (Figure 6). In addition, insulin induced a decrease in the portion of cells at S phase (from 14 to 7.95%) and an increase in the portion of cells at the G0/G1 phase (from 49 to 59.55%). In contrast, in IGF1R-KD cells both IGF1 and insulin led to increases in the proportion of cells at the SubG1 phase (from 0.98 to 2.26 or 2.58%, respectively). In addition, insulin elicited a mild increase in the portion of IGF1R-KD cells at G0/G1 phase (from 55.7 to 57.6%) and a decrease in the portion of cells at S phase (from 23.6 to 18.2%). As expected, IGF1 led to a small but significant increase in the portion of cells at G0/G1 phase (from 58.2 to 60.6%) and S phase (from 1.

Senescence Assays
To examine the differential involvement of INSR and IGF1R in acquisition of a senescence phenotype, the ability of etoposide (a DNA-damaging anticancer drug) treatment to activate such a response was examined. Etoposide-treated control and INSR-KD cultures exhibited changes in morphology that are typically associated with senescence. Specifically, cells displayed a flattened aspect, increased size and an altered ratio of nucleus to cytoplasm, with higher granularity compared to untreated cells ( Figure 7A). Interestingly, etoposide treatment of IGF1R-KD cells for 48 h led to a marked decrease in the portion of senescent, compared to control, cells (Figures 7A,B). Thus, etoposide led to a ∼5-fold increase (from 18 to 90%) in β-gal-stained control cells, a ∼1.17fold increase (from 70 to 82%) in β-gal-stained INSR-KD cells, and a 48% decrease (from 68 to 36%) in β-gal-stained IGF1R-KD cells (Figure 7B). Furthermore, levels of the p16 and p21 senescence protein markers were enhanced in etoposide-treated control and INSR-KD cells (Figure 7C). In contrast, the levels of both proteins were diminished in etoposide-treated IGF1R-KD cells. Taken together, data suggest that IGF1R has a critical role in the attainment of senescence. The simulated expression of p16 was minimally changed in both INSR-KD and IGF1R-KD cells with respect to controls (Figure 4N), validated by Western blots (Figure 4M). On the other hand, p21 simulation analysis was not affected by the disruption (Figure 4P), validated by Western blots (Figure 4O).

DISCUSSION
A fundamental question in the area of insulin and IGF1 biology concerns the fact that INSR and IGF1R, even though they share the majority of their downstream cytoplasmic targets and signaling pathways, are yet responsible for mediating distinct physiological roles. Numerous theories have been formulated to explain these discrepancies, including a different tissue distribution of INSR and IGF1R (4), divergent internalization kinetics and subcellular localization of the hormone-receptor complex (34), and dissimilar hormone-receptor affinities (35). Furthermore, a number of substrates that are preferentially activated either by insulin or IGF1 have been identified (36). The present study was aimed at identifying mechanistic differences between INSR and IGF1R using BioNSi, a recently developed bioinformatics tool. This simulation tool allowed us to exploit expression data of breast cancer-derived cell lines with specific disruptions of the INSR or IGF1R genes and their targets in an extended network based on nine KEGG pathways. We thus investigated the differences in protein expression that might be linked to biological attributes of these specific receptor networks.
Using the BioNSi tool, we found 78 genes that were changed between both KDs and control cells. Among them, 76 are downstream genes to INSR and IGF1R. One must note that the BioNSi simulation tool is quite simple in terms of input information available from KEGG. Thus, as the signal is farther down the network, the effect seen may be lost. Therefore, we manually increased the initial expression values of the signal start points (IGF1 and insulin) to 20. At such conditions, the most downstream gene that exhibits a change in expression value during simulations is CHEK2, with 12 edges between it and insulin or IGF1 (the signal molecules) and 11 nodes. During IGF1R-KD simulation, CHEK2 keeps its simulated expression value the same as in the control simulation, up to step 32. Then, its expression value drops to zero and rises up to nine again at step 58. This pattern is a result of the complicated gene network it belongs to. If the initial expression value of IGF1 was lower, such changes at the later steps of the simulation may not have been observed, as the strength of the effect is deteriorating with steps.
In addition to signal expression, one must note that although nine KEGG pathways were selected for this analysis, based on 16 key genes, these represent only a small part of the cellular proteins actually involved in the biological processes investigated. Moreover, BioNSi analysis is based on proteins that are all affecting in a similar manner. As a matter of fact, we know that individual proteins act differently on each other. Hence, it is important to emphasize that the in silico BioNSi simulation is based only on partial available molecular information that might mask some of the in vivo results. Taking into account these limitations, the simulation analysis investigated 16 genes that were validated correctly by Western blotting and biological assays. Not surprising, but a necessary validation, was the observation that KD of both IGF1R and INSR genes resulted in a sharp reduction in the expression of the selected KD gene in both simulation and Western blots (Figure 3).  steps of the simulation, but INSR-KD cells show no difference from control. After 50 steps, differences can be seen between control and both disrupted cell types. This phenomenon can be explained by the differences in time frames in which simulation and biological assays are being measured. Simulation starts from a static point of gene expression measurement and expression values are calculated at each step, until 100 steps. At this time, all expression levels have reached a steady-state [either the maximal (9) or minimal (0) expression level]. At the steady-state level, no additional changes are being made. On the other hand, the biological assays take place at a specific point during the lifetime of the cell and cannot reflect all changes in expression as predicted in simulation. Therefore, the results of simulations and Western blotting must be compared in a careful manner. Thus, Western blot results of a specific protein may reflect simulated expression value seen earlier in simulation for upstream genes and/or later for downstream genes (compare CHEK2 and ATM).
Additional differences between results of simulations as compared to Western blots may be due to the selection of specific pathways. In other words, we did not aim to simulate the entire cell, but rather selected important pathways that are responsible for IGF1R and/or INSR downstream responses. However, the fact that many of the interactors in the network are players in multiple cellular pathways, some of which were not the focus of our simulation, might eventually lead to divergent results between both types of analysis. TP53 is the classical example of this situation. Finally, it is important to emphasize that the present study was conducted on MCF7-derived cells, an estrogen receptor positive/progesterone receptor positive/human epidermal growth factor receptor-2 negative adenocarcinoma line derived from a metastatic site. Caution must be exerted when extrapolating data reported here to other cell types, particularly cells with a negative steroid hormone status. IGF1 has been identified as a key player in cell cycle progression (37)(38)(39)(40). Our analysis revealed a number of differences in the proportion of cells at the different stages between IGF1R-KD and INSR-KD cells. For example, the portion of cells at G2M was 37.2% in untreated INSR-KD cells (expressing predominantly an IGF1R) in comparison to 22.6% in IGF1R-KD cells (expressing mainly an INSR). The fact that IGF1 was able to increase the proportion of cells at SubG1 in IGF1R-KD cells (from 0.98 to 2.26%) may suggest that IGF1 is able to crossactivate the INSR, which is the main receptor in these cells. Alternatively, results may indicate the presence of residual IGF1R in IGF1R-KD cells.
The involvement of the insulin-like axis on senescence (at both cellular and organism levels) has been widely reported in recent years (41)(42)(43)(44)(45). Our results provide evidence that etoposide, a DNA damaging agent, leads to a significant decrease in the proportion of senescent cells in IGF1R-KD, compared to INSR-KD and control, cells. Levels of the p16 and p21 senescence markers were enhanced in etoposide-treated control and INSR-KD cells but markedly decreased in IGF1R-KD cells. These findings suggest that IGF1R expression and action are critical for progression of senescence. Furthermore, results are in agreement with reports showing that IGF1R can induce senescence in fibroblast (46), keratinocytes (47), and endothelial (45) cells. The biological role of the IGF1R in senescence is reflected in the well-documented involvement of the growth hormone-IGF1 axis in longevity. Evidence has accumulated showing that disruption of this hormonal system is correlated with extended lifespan in various animal species, including nematodes, flies, and mouse (48)(49)(50). Female "Laron" mice (GHR/BP −/− ) have a 38% longer lifespan than wild-type animals (51). While the impact of specific mutations on human lifespan are difficult to evaluate, studies have identified functionally significant IGF1R mutations in centenarians (52).
In summary, despite the limitations of the bioinformatics method and its inability to take all known interactions into account, BioNSi simulations constitute an important addition in insulin/IGF1 research. The tool is aimed at predicting protein interaction effects in silico in order to suggest biologically plausible models for experimental evaluation. Our computational analyses, validated by Western blot analyses and biological studies, have identified a number of commonalities and, most importantly, dissimilarities between the IGF1R and INSR pathways that might help explain biological differences between these networks. We believe that comprehensive investigation of the INSR and IGF1R pathways, aided by modern experimental and bioinformatics technologies, might have a translational impact in the area of INSR/IGF1R therapeutic targeting.

DATA AVAILABILITY STATEMENT
The article contains previously unpublished data. Gene expression data has been deposited at the GEO repository. Accession number is GSE145787.

AUTHOR CONTRIBUTIONS
RS, AY, MP-C, and HW conceived of and designed the study. The experimental procedures were performed by RS. The simulation analyses were conducted by AY. The microarray experiments were conducted by TS-L. Additional computational analyses were conducted by AY and MP-C. Statistical analyses were performed by RS. HW and MP-C prepared the manuscript. All authors approved the manuscript.