|Home | About | Journals | Submit | Contact Us | Français|
Cancer is now appreciated as not only a highly heterogenous pathology with respect to cell type and tissue origin but also as a disease involving dysregulation of multiple pathways governing fundamental cell processes such as death, proliferation, differentiation and migration. Thus, the activities of molecular networks that execute metabolic or cytoskeletal processes, or regulate these by signal transduction, are altered in a complex manner by diverse genetic mutations in concert with the environmental context. A major challenge therefore is how to develop actionable understanding of this multivariate dysregulation, with respect both to how it arises from diverse genetic mutations and to how it may be ameliorated by prospective treatments. While high-throughput experimental platform technologies ranging from genomic sequencing to transcriptomic, proteomic and metabolomic profiling are now commonly used for molecular-level characterization of tumor cells and surrounding tissues, the resulting data sets defy straightforward intuitive interpretation with respect to potential therapeutic targets or the effects of perturbation. In this review article, we will discuss how significant advances can be obtained by applying computational modeling approaches to elucidate the pathways most critically involved in tumor formation and progression, impact of particular mutations on pathway operation, consequences of altered cell behavior in tissue environments and effects of molecular therapeutics.
The advent of genomic sequencing and copy number evaluation has revealed that aspirations for understanding cancer in terms of mutation of some limited number of oncogenes and tumor suppressor genes are unlikely to be fulfilled; instead, even tumors of a particular tissue type bear highly heterogenous sets of defects in dozens of different genes (1–3). Similarly, RNA interference studies show that a large number and wide spectrum of gene products contribute to tumor cell phenotype (4,5). The processes affected closely correspond to defined ‘hallmarks’ of cancer, such as resistance to cell death, extension of replicative potential, enhancement of invasiveness, and escape from immune surveillance, among others (6,7). Thus, the relationship between genomic information per se and malignant disease is more elusive than originally hoped, due to the broadly multivariate nature of the molecular-level changes involved in any given cancer (8). Moreover, even when a specific gene is identified to make a substantial contribution to pathology, that determination does not easily lead to an effective avenue for treatment because of the complex consequences propagated down transcriptional, translational and posttranslational circuits (7,9).
Although the genomic mutations associated with any tumor are many and diverse, a promising organizing principle is emerging: that a recognizable cohort of key pathways governing cell phenotypic behaviors can be identified as pathologically altered by the various underlying genetic defects (10,11). This notion is encouraging for elucidating therapeutic targets because a pathway (or multiple pathways) can be targeted through any of numerous comprised components—whether or not they in fact suffer from genetic defect—rather than being restricted solely to mutated oncogenes or tumor suppressor genes and needing to find a unique approach for each gene product (7). An example of this is the class of mammalian target of rapamycin inhibitors (including temsirolimus and everolimus) that are being found effective for treatment of renal cell carcinoma; the predominant genetic defect in these tumors is loss of Von-Hippel Lindau tumor suppressor (thus reducing degradation of the hypoxia-inducible transcription factors HIF1/2/3) rather than any mutation related to the mammalian target of rapamycin gene itself (12).
A pathway-centric approach remains incomplete, however, because of the intricate cross talk among cell regulatory pathways (13). Indeed, a given molecular component can be identified to be associated with or interact with multiple signaling, transcriptional regulation, metabolic and/or cytoskeletal process pathways (14). Pathways thus cannot properly be considered to operate in isolation of one another, as an alteration of one pathway can lead directly (via protein–protein interactions) or indirectly (via transcriptional/translational influences) to changes in others. Accordingly, cancer—along with other complex diseases such as arthritis and diabetes—is most productively conceived of and strategized for treatment as a dysregulation of a multi-pathway network (15,16). Moreover, these networks connect to components beyond the tumor cells themselves, including other cells in the environment along with the extracellular milieu (17,18). This perspective yields at least two consequences. First, experimental characterization of the pathological dysregulation will need to be multivariate and quantitative. For instance, biomarkers in the form of single molecular components or qualitative component lists will probably be inadequate even for categorization of treatment outcomes. Second, predictive or mechanistic understanding of the pathology will almost certainly elude intuition unaided by computational analysis. Hence, a network perspective on cancer strongly motivates the application of computational modeling approaches (19–21). Methodologies for computational analysis can vary widely depending on the question being posed and the experimental data at hand, ranging from highly abstracted models using correlative regression to highly specified models using differential equations, with network component interaction and logic modeling techniques intermediate to these. A number of reviews have provided discussion of which of these methods are more or less appropriate for employment for various kinds of studies, outlining their respective strengths and weaknesses with respect to different applications (22–25).
Herein, we will present selected examples of recent research contributions that are helping to establish the field of cancer systems biology. These studies demonstrate the unique advances in understanding and prediction that can be gained by integration of computational modeling with quantitative experimental data on molecular and cellular networks. Our examples emphasize systems at the level of dynamic protein operations, such as phosphorylation, because it is this level that most effectively integrates the convolution of genomic information and environmental context. As illustrated in Figure 1, environmental context strongly influences network operations regulating transcription, translation and posttranslational processes, so that restricting experimental measurement information to DNA sequence, mRNA expression or even protein expression will miss important aspects of the molecular components and interactions governing cell and tissue behavior. We organize this presentation into five main categories, denoting the kinds of problems being addressed by the various studies: identifying dysregulated pathways, elucidating consequences of mutations on network activities, integrating network operation into cell behavioral functions, integrating cell behavior into tissue-level processes and predicting effects of molecular interventions.
Historically, dysregulated pathways have been identified in cancers based on reductionist studies surrounding an identified mutation. The shift to examining dysregulated networks requires new techniques to identify pathways within the context of the intact cellular network. An especially compelling avenue for identifying affected pathways is the use of interactomes, which define the molecular interactions in a cell. Although still at an early developmental stage for human cells and tissues, interactomes include protein–protein and protein–DNA associations (26,27) and provide a framework for analysis of empirical data of various types, such as transcriptomic, phosphoproteomic and phenotypic assessments. The biochemical resolution of interactome information has been recently enhanced by directly incorporating phosphoproteomic data to indicate not only physical interactions but also kinase–substrate interactions (28). A further exciting breakthrough toward an increasingly powerful network framework for identifying pathways involved in responses to environmental conditions or in governing phenotypic behaviors is offered by a recent development of a Steiner tree computational algorithm, able to integrate protein–protein and protein–DNA interactomes (29). The initial application of this new methodology has been in yeast, demonstrating success in linking genetic data (e.g. knockouts and knockdowns) with gene expression data and in gleaning new insights concerning network activities characterizing the yeast pheromone response. These techniques promise to improve our ability to discern the complex networks impacted by oncogenic changes in human cells.
Utilizing an internally assembled interactome and a compendium of >200 transcriptional microarray profiles for normal and tumor-related B cells, Mani et al. (30) investigated the pathways dysregulated in human B cells for three kinds of non-Hodgkin's lymphomas (follicular, Burkett's and mantle cell). A mutual information algorithm was applied to identify interactions exhibiting gain-of-correlation or loss-of-correlation with respect to tumor phenotype relative to normal. Dysregulated interactions were defined as those exhibiting mutual information (essentially, correlation) in all samples except for a particular phenotype (loss-of-correlation) or lack of mutual information in all samples except for a particular phenotype (gain-of-correlation). A fundamentally interesting finding was that ~80% of the roughly 65000 network interactions appeared to be unaffected in any of the tumor phenotypes. That is, the interactions demonstrated similar mutual information whether in tumor or normal cellular backgrounds; the authors noted that this implies a network ‘backbone’ that operates consistently across various cellular backgrounds. Nonetheless, hundreds of interactions were discerned to be differentially correlated with each of the particular tumor phenotypes, cutting across multiple diverse pathways. A further important insight obtained was that dysregulated pathway interactions could arise even without mutation of genes explicitly corresponding to the pathway components involved.
Chang et al. (31) pursued an analogous study of the NCI-60 cancer cell line set (32), analyzing a combination of transcriptional profile data with human protein–protein interactome information but using a different computational method called statistical factor analysis. This method undertakes a linear regression of gene expression data with respect to particular gene signature vectors, determining score coefficients relating the strengths of contribution of corresponding vectors to a given set of transcriptional profiles. A major focus of this study was elucidation of how important different effector pathways downstream of RAS are in generating the expression profiles exhibited by certain cancer cell lines. Subsets of cell lines bearing stronger contributions of the extracellular signal-regulated kinase (ERK) effector pathway or the AKT effector pathway were identified and validated, at least in aggregate, by their comparative sensitivities to ERK pathway inhibitors versus AKT pathway inhibitors.
Models that include influence and/or logic aspects of network component interactions can also be constructed, with the advantage of understanding propagation of pathway activities following introduction of stimulatory or inhibitor cues. One example application of this kind of approach was recently contributed by Heiser et al. (33). The authors utilized Pathway Logic formalism to elucidate how the ERK pathway is activated downstream of ErbB family receptors in breast cancer cell lines. A network comprising 286 signaling nodes and 396 interaction rules was constructed via literature curation, and disparities in network topologies across different cell lines were identified from cell line-associated mRNA and/or protein expression data for 191 of these components. The expression levels were discretized into ‘present’ or ‘absent’ categories, thus yielding altered interaction rules among cell lines depending on which components were deemed in one category or the other. Across the 30 breast lines examined, roughly 40 components were found to be disparately expressed (with respect to present versus absent calls), resulting in >50% of the interaction rules varying concomitantly. Clustering the cell lines with respect to network features generated classification associated with pathway interactions, leading to identification of novel points such as the presence of PAK1 in association with strong activation of ERK. Explicit calculations on network information flows can be conducted using dynamic Boolean logic or fuzzy logic algorithms (34–37). To date, these algorithms have been mainly applied to hand-curated networks, but the prospect for being founded on formal interactome databases is clear. Indeed, a new contribution shows how a Boolean logic model appropriate for an actual cell type and particular set of contexts can be derived from combination of empirical data with a prior interactome database (38). Bayesian networks, which describe pathway interactions (whether direct or indirect) in terms of probabilistic influences of certain components upon other components, have been demonstrated useful for understanding operation of cell signaling networks but not yet applied to cancer biology problems (39,40).
We note that one alternative approach to using interactome information is the definition of gene modules based on Gene Ontology categories, which has been successfully applied to discerning types of pharmacological agents effective in killing tumor cells representative of cohorts classified by these modules (9). Another option not reliant on an interactome framework is direct application of mutual information-based algorithms, such as algorithm for the reconstructions of accurate cellular networks (41) to relate transcription factor activity to gene expression profiles. This methodology has been successfully demonstrated in work that ascertained how NOTCH1 and c-MYC work together to regulate growth of T-cell leukemia cells (42).
Oncogenic mutations affect cell behavior by changing the cellular network (13). While these effects include obvious changes to proteins immediately downstream, systems biology studies are beginning to reveal how apparently small changes in the network have broad-reaching effects (43). Here, we will profile recent studies of three of the most common mutated pathways in cancer—p53, the ErbB family of receptors and RAS. While many reports have focused on the identification of biomarkers for tumor detection, our focus is on how the mutation changes the network and how this information could be used to identify new targets or treatment strategies.
In response to DNA damage, a variety of pathways are activated in the cell to arrest cell cycle progress, repair the damage or initiate apoptosis in order to prevent generation of cells with mutations. Accordingly, many tumors have alterations in their DNA damage response pathways. A number of systems biology studies have focused on p53, the tumor suppressor mutated in at least half of all cancers. In single cells, levels of p53 oscillate in response to stress, which has been shown computationally to involve feedback circuits in the DNA damage-signaling network. Batchelor et al. (44) used a combination of dynamical systems modeling and quantitative single-cell experiments to determine that the p53/MDM2-negative feedback loop is critically involved in this phenomenon. Activation of signaling kinases ATM and CHK2 activation downstream of DNA damage drive p53 peaks, with WIP1-mediated feedback maintaining coherence during ensuing cycles. In an ensuing contribution employing stochastic physicochemical modeling, p14ARF was identified as another key participant in the feedback (45). An important further step was provided by Toettcher et al. (46) in connecting these signals to effects on cell cycle control. By integrating the DNA damage-signaling pathways with cell cycle regulatory pathways, the authors were able to investigate the role of p53 in both initial and long-term maintenance of cell cycle arrest. Initial arrest was determined to be p53 independent; however, p53 was necessary to maintain arrest, indicating its key role in preventing accumulation of DNA damage and importance in tumor development. Moreover, a novel prediction of pathological endo-reduplication behavior enabled by defects in the p21 pathways was produced by the model and validated by direct experimental test. Linking model predictions to experimental tests is an important step to further the application of systems biology to the study of oncogenic impact.
The ErbB family of receptor tyrosine kinases and components of the downstream network are frequently mutated in cancer (47). The complexity of the ErbB system, with four receptor isoforms and >12 ligands, makes it an ideal network to analyze by systems biology methods (48). Not surprisingly, this system has been subject of numerous studies applying modeling to quantitative experimental data, with the most recent contributions incorporating multiple members of the receptor family (49,50).
One common perturbation to the ErbB network found in breast, lung and colon cancers is the overexpression of ErbB2 (HER2). A mass-action model accounting for ErbB1–4 dimerization and signaling to ERK and AKT determined that overexpression of ErbB2 shifts the cell to a higher percentage of ErbB1–2 heterodimers in place of ErbB1 homodimers (51). This model suggests that since these heterodimers do not undergo ligand-induced degradation, increase of ErbB2 results in sustained signaling. Utilizing partial least squares regression modeling applied to phosphoproteomic data from human mammary epithelial cells expressing increasing levels of ErbB2, Kumar et al. (52) identified nine phosphorylation events that serve as a ‘network gauge’ to determine the impact of treatment proliferation or migration. The components of this gauge include some of the usual suspects such as PI3K signals and endocytosis proteins, but it is the quantitative integration of these signals, rather than any single event, that enables prediction of cell behavior.
In addition to ErbB2 overexpression, mutations to ErbB1 [epidermal growth factor receptor (EGFR)] are frequently identified in tumors. Phosphoproteomic analysis of cells expressing increasing levels of the constitutively active EGFRvIII mutant demonstrated that the active phosphorylation site of c-MET was highly responsive to EGFRvIII level (53). This suggested that EGFRvIII cross-activated c-MET; indeed, dual inhibition strategies against EGFR and c-MET proved to be more effective in killing cells than single treatments. Identifying cotreatment strategies such as these is one of the promising applications of systems biology. Additionally, modeling can help to explain clinical observations of why certain tumors are more susceptible to treatment by targeted inhibitors. Deletion of exon 19 and the point mutation L858R in EGFR are associated with elevated phosphorylated AKT and sensitivity to the EGFR tyrosine kinase inhibitor gefitinib. Using a mass-action model incorporating receptor internalization and activation of ERK and AKT, slower EGFR internalization was determined to be sufficient to explain the signaling differences observed between wild-type EGFR and these mutants (54). When receptor internalization was quantified, mutant cell lines were indeed found to have slower rates. The model suggests that delayed internalization leads to AKT addiction, which is then inhibited by gefitinib treatment, providing an explanation for the limited benefit for tumors with wild-type EGFR.
Downstream of various receptor tyrosine kinases, members of the RAS family of GTPases are frequently mutated in human cancers. RAS mutations in tumors are primarily point mutations in one isoform and result in insensitivity to GTPase-activating proteins that increase GTP hydrolysis. These mutations result in increased levels of active RAS–GTP in cells, whereas additional non-mutant RAS isoforms are subject to ligand-induced activation (55). The RAS pathway has been modeled downstream of various receptors (50,51,56), and the RAS activation–deactivation cycle has been analyzed quantitatively (55). However, these contributions have not addressed the effects of RAS mutation on downstream effectors, which are not as straightforward as simply enhanced effector activities. Two recent systems biology studies have focused on this issue and have demonstrated that oncogenic RAS impacts multiple targets in complex ways. In the first study, an inducible H-RASG12V construct was activated and rapid induction of two phosphatases, DUSP1 and DUSP6, was observed (57). To examine this behavior, multiple model structures were characterized to see which model structure best explained the experimental data. The best-fit model incorporated ultrasensitive activation of ERK in response to RAS activation, ERK induction of DUSP6 and not DUSP1 and feedback from DUSP6 against phosphorylated ERK. Additionally, we have recently reported the impact of mutations to two different isoforms of RAS, K-RAS and N-RAS, on the response of colon carcinoma cell lines to tumor necrosis factor α (TNFα) (58). Using a quantitative phosphoproteomic data set, we determined that K-RAS-mutant cells have decreased activation of ERK due to a depressed TNFα-induced release of autocrine transforming growth factor-α and that N-RAS-mutant cells do not induce DUSP5 as robustly as either wild-type or K-RAS-mutant cells, resulting in sustained activation of ERK. Studies such as these demonstrate how small perturbations impact the cell network by affecting positive and negative feedback mechanisms.
The process of developing from a single cell with an oncogenic profile to a metastatic cancer is exceedingly complex and multivariate. However, there are several key processes common to most cancers, including excessive proliferation, resistance to apoptosis, angiogenesis and metastasis (6). Systems biology modeling techniques promise to improve our understanding of each of these processes and ultimately, how they interact to drive tumor progression.
Excessive proliferation is perhaps the phenotype most associated with cancer progression. The ErbB-signaling network plays an important role in dysregulated proliferation in many tumors and has been extensively modeled (59). Recent models have highlighted the important multivariate characteristics of the ErbB-signaling network in proliferation as the ErbB signal transduction network interacts with a variety of other signaling pathways. Modeling approaches have proven useful in identifying dominant mechanisms of pathway cross talk. For example, a mass-action kinetic model examining the cross talk between the ErbB network and insulin signaling indicated a role for GAB1 in the synergism between treatment with epidermal growth factor and insulin (60). Additionally, identifying multi-pathway interactions may provide insight into potential new therapeutic targets. Sahin et al. (61) developed a Boolean model from the literature linking ErbB1–3 to phosphorylated retinoblastoma protein (a surrogate for cell cycle progression) in breast cancer resistant to ErbB2-targeted inhibitors. Following refinement, the model was able to predict a variety of new conditions that resulted from knocking down proteins in the network. Importantly, simulations indicated that inhibiting ErbB receptors would be insufficient to halt cell cycle progression in resistant cells and suggested c-MYC as a potential alternative target.
In normal development, apoptosis provides a counterpart to proliferation by removing damaged or unnecessary cells. To evade restriction of tumor growth, cancer cells have devised multiple mechanisms to provide sustained pro-growth stimuli or to counteract apoptosis. Apoptosis can be induced by stress and the mitochondrial pathway (intrinsic) or through activation of death receptors (extrinsic) (62). In the intrinsic pathways, stimuli such as DNA damage activate the apoptotic machinery. To model this process, one approach has been to develop mass-action models that decompose the intrinsic apoptotic pathway into subsections of the molecular network for analysis (63). Using such models to examine the impact of DNA damage on p53 phosphorylation suggested that transient DNA damage leads to a level of phosphorylated p53 that will induce cell cycle arrest, whereas sustained damage will induce apoptosis.
A novel mass-action model of tumor necrosis related apoptosis-inducing ligand-induced (extrinsic) apoptosis and mitochondrial outer membrane permeabilization in single cells has recently been employed to examine several mechanistic questions about the regulation of cell death in the HeLa cancer model (64,65). Model parameters were fit to experimental data from tumor necrosis related apoptosis-inducing ligand dose-response curves and small interfering RNA/overexpression perturbations to the network. A particularly striking result of the model is the impact of compartmentalization—once the steady changes in the cell network lead to mitochondrial compartment opening, there is a drastic release of SMAC down the concentration gradient and an all-or-nothing decision for apoptosis results (64). Additional work demonstrated that XIAP- and proteasome-dependent degradation of effector caspases is required to prevent caspase activity before the cell has committed to apoptosis by mitochondrial outer membrane permeabilization (65). This suggests that if this control is overridden, caspases may be activated and cause damage to target substrates, while the cell is in an ‘undead’ state. The developed model was also used to address the causes of cell-to-cell variation in timing of TRAIL-induced death (66). Using experimental measures of the mean and variation of five key proteins in the model, simulation results closely match the time-to-death variation observed. Combined with other experimental observations, this suggests that cell-to-cell variation is a result of changes in protein concentrations rather than genetic or stochastic mechanisms.
It is important to remember that for a tumor to develop, it is the balance between apoptosis and proliferation that matters. In models of single colonic crypts composed of stem cells, differentiated cells and transit cells, increased proliferation, decreased differentiation or decreased apoptosis all lead to a net increase in cell number (67). Alternatively, it has been argued that resistance to apoptosis may actually decrease the ability for a tumor to expand as it would reduce the number of cell divisions before a tumor reached a critical size and lower the probability of forming the mutant combinations necessary for expansion from that stage (68). At the molecular level, recent experimental studies have demonstrated several mechanisms by which cells respond to both pro-death and pro-proliferative stimuli (69). Using principal components analysis to examine HT-29 apoptosis in response to TNFα, insulin or epidermal growth factor, a series of pro-survival and pro-death autocrine loops were revealed as a result of treatment with TNFα (70). These molecular balancing mechanisms help to explain why cytotoxic ligands are unable to kill all cells, which will be essential to consider when administering targeted chemotherapeutics.
An ultimate challenge is the development of models that incorporate the behavior of the entire tumor. In the early-stages of tumor development, cancer cells are clustered together into an avascular structure, which has been modeled in vitro with multi-cellular spheroids (71) and in silico with mathematical balances of proliferating and dying cells that are coupled to physical constraints such as nutrient diffusion. In general, these models are limited to describing tumor morphology and distribution of necrotic cells (72). More recently, multiscale modeling was used to analyze avascular multi-cellular tumors (73). This model incorporated proliferation, diffusion and consumption of nutrients, diffusion and production of wastes, adhesion and cell–environmental interactions. A notable improvement over previous models that relied on probabilities for cellular decisions was the inclusion of a Boolean logic model for the G1 to S progression, which incorporated molecular components such as transforming growth factor-β, p27, p21 and cyclins. Model simulations closely matched in vitro spheroid morphology, size and cell cycle distribution over a 20 day period (73). Further refinement of detailed models such as this should allow investigators to test the impact of molecular interventions on early-stage tumors.
As a tumor increases in size, the center core becomes necrotic and the tumor needs to develop its own vasculature to continue to grow. Therefore, targeting tumor angiogenesis is an attractive strategy to treat cancer (74). Computational models are beginning to provide tools to address how the tumor microenvironment and growth factor signaling regulate these events. Although numerous models for angiogenesis have been proffered, the complexity of the system presents a challenge regarding how to connect the behavior of the blood vessels to the tumor and how to incorporate molecular control mechanisms. Two recent models link the behavior of the tumor, growing blood vessel and cellular microenvironment. Macklin et al. (75) modeled tumor growth in response to the altered oxygen profile resulting from new vascular growth. In this model, mechanical forces directed blood vessel growth, and simulations demonstrated that when heterogenous oxygen profiles were formed from the new vessels, tumor cell proliferation was heterogenous resulting in invasive tumor morphologies. In the second model, blood vessel growth was modeled as a response to local gradients in angiogenic factors such as VEGF (76). By changing the various rules, these model forms can mimic different environmental or genetic conditions; however, they do not allow for the direct test of molecular interventions on angiogenesis.
As an alternative approach, another recent model linked the cancer cell cycle to sensing of environmental conditions such as oxygen levels, with the assumption that suboptimal oxygen levels result in the production of VEGF (77). In contrast to more phenomenological models, the impact of VEGF on endothelial cell proliferation and migration was determined by the pharmacological Emax model and additional pro- and anti-angiogenic molecules were included. Simulations of the effect of endostatin induced by gene therapy indicated there is a critical rate of production needed; below this rate, longer treatment times were predicted to ‘rebound’ more quickly, whereas above this rate, the time to rebound increased with the duration of treatment (77). Expansion of the molecular components of multiscale models of angiogenesis will be essential to utilize them to determine drug targets and optimal treatment regimens.
An alternative method for a tumor to continue its growth is for cells to leave the primary tumor and implant in other tissues—the metastases that result from this process are responsible for the majority of cancer deaths. It is becoming increasingly apparent that metastasis is more than a random process, with additional mutations required for cancer cells to leave the primary tumor and metastatic cells demonstrating ‘preferences’ for target organs (78). Similar to angiogenesis, models of tumor invasion and metastasis must incorporate multiple scales and are primarily constructed using approximations of cellular signaling. For example, in one virtual tumor model, cells respond to microenvironmental changes (e.g. oxygen level, cell proximity) and decide to proliferate or die (79). At the same time, cells undergo random inheritable mutations that change the ‘phenotype’ of the cell (e.g. the likelihood to proliferate). Changes to the microenvironment to incorporate harsher conditions such as hypoxia or heterogenous matrix led to selection for cells with more aggressive traits and invasion tumor shapes. The impact of the matrix on cancer cell invasion is supported by a recent model of invadopodia, the cellular extensions believed to degrade extracellular matrix as tumor cells invade (80). In this model a single cancer cell sits on top of a matrix, which is modeled from known dimensions and characteristics of extracellular matrix proteins. A series of rules describe how invadopodia invade and interact with the extracellular matrix. The simulations suggested that dense matrix (such as basement membrane) forms an effective barrier to prevent invadopodia penetration and matrix degradation, whereas looser matrices (such as stroma) have gaps sufficient to allow invadopodia penetration and matrix degradation.
In order to metastasize to distant organs, cancer cells must invade into the blood or lymph to be transported to other parts of the body. The first stage in this process is for the cancer cell to invade the endothelial layer. To model this process, a multiscale model of transendothelial migration was developed, incorporating both endothelial cells, cancer cells and the cadherin interactions between the distinct cells (81). In the simulated invasion, cancer cells attach to endothelial cells by N-cadherins, and the endothelial layer is connected by VE-cadherin bonds, which the cancer cell must break in order to squeeze through. The model follows protein concentrations within each cell through a series of ordinary differential equations, cell–cell forces by a modified Hertz model and cell movement according to Langevin equations. As cancer cells contact endothelial cells, a N-cadherin bond forms between the cells, leading to competition for β-catenin and dissolution of the existing VE-cadherin bonds. The observed behaviors are similar to experimental results that noted the loss of N-cadherin slows transendothelial migration (82).
Ultimately, models such as these must be built upon to balance multiple signaling pathways and phenotypic processes in order to be used for identification of drug targets. A novel multiscale model of non-small-cell lung cancer provides an example of how such a process can be implemented. In an early version, ErbB activation of ERK and PLCγ in each cell in a tumor was modeled by a series of ordinary differential equations (83). Phosphorylated ERK and PLCγ are used as readouts for proliferation and migration, respectively, and cells migrated or proliferated in a virtual two-dimensional space. In this early version, overexpression of ErbB1 led to a migration dominant phenotype and accelerated tumor expansion. Building on this initial work, the group expanded the virtual space and tumor to three dimensions and added transforming growth factor-β activation of RAS to each cell's signaling module (84). By comparing the results of simulations with single and cotreatments, the expanded model demonstrated conditions where targeting a single pathway would be ineffective in deterring tumor expansion.
Importantly, computational network models have begun to identify novel targets in regulatory networks with the goal of more effectively treating tumors. Schoeberl et al. developed a mass-action kinetic model of ErbB1–3 incorporating ligand binding (betacellulin and heregulin1-β), receptor dimerization, internalization, degradation, recycling and downstream signaling through PI3K to phosphorylated AKT (85). Sensitivity analysis of the model revealed a dominant role of ErbB3 in the level of phosphorylated AKT, and in silico tests of a monoclonal antibody against ErbB3 demonstrated its effectiveness at decreasing phosphorylated AKT over a range of ErbB1/2 levels. Based on the predictions of their systems biology model, the group developed MM-121, a fully human IgG2 monoclonal antibody against ErbB3, and found that, as in the model, blockade of ErbB3 was effective in blocking phosphorylation of AKT (85). Additionally, when mice implanted with ACHN xenografts were treated with MM-121, tumor growth was slowed even after the end of the treatment period and MM-121 has entered Phase I Trials. This study is particularly provocative as the model suggested targeting the kinase-dead ErbB3 rather than a mutated or overexpressed protein, indicating that excessive growth may be dependent on the multivariate receptor and ligand environment rather than a single oncogenic change.
However, this work—like most systems modeling efforts aimed at predicting effects of therapeutic perturbations of cell regulatory networks—restricts its attention to predicting molecular-level processes (in this particular case, AKT phosphorylation). What is vital, of course, is to predict effects of these perturbations on cell phenotypic functions at the very least; predicting effects of integrated tissue properties would be an even more distant goal. The most difficult problem is connecting molecular-level network activities to cell-level functional behavior, even in absence of therapeutic perturbations. Relational modeling methods, such as partial least squares regression, currently represent one of the most effective approaches for this goal. This approach was tested directly by Kumar et al. (86), who used partial least -squares regression modeling to predict how modulation of key kinase pathways such as ERK, AKT and p38 by pharmacological inhibitors alters migration behavior of ErbB2-overexpressing mammary epithelial cell migration in response to epidermal growth factor and Heregulin. The most important insight gained from this contribution was that while a quantitative combination of five phosphorylation sites on four key proteins could successfully interpret the effects of kinase inhibitors on cell phenotypic behavior across multiple treatment conditions, no individual signal could predict the phenotypic behavior by itself. This finding emphasizes that signal-to-response relationships will in general require multiple signaling pathways to be included in the model and that attempts to predict how cells will behave on the basis of alterations (whether genetic or therapeutic) to a single component or pathway will most probably be in vain.
A major issue to emphasize as we conclude this discussion is the challenge posed by the gap between in vitro studies and in vivo. This, of course, is a crucial problem for the entire cancer biology field, not merely the systems biology approach. We submit that a systems perspective, in which multiple variables are considered integratively in explicit manner, is at least as likely or more so to find some significant success in bridging this gap than a focus restricted to intuitive predictions based on individual component effects. A compelling demonstration of this notion has recently been provided by Jiang et al. (87). These authors showed that assessment of both ATM and p53 together is required in order to move from cell culture experiments to organism studies, with the goal of predicting the efficacy of a DNA-damaging chemotherapeutic agent on destruction of lymphoma tumors in mice and survival of human breast cancer patients. While this particular effort did not incorporate formal mathematical analysis, a schematic logic model (involving DNA-PK and Chk2 along with ATM and p53) was found helpful in explicating the integrative network operation. Extension of this kind of perspective to an increased number of components and wider set of pathways will be facilitated by the kinds of computational approaches we have noted, in order to enable understanding and prediction beyond even skilled intuition. We emphasize the work by Schoeberl et al. (85) again here as a promising example of employing computational modeling for effective prediction of novel and non-intuitive cancer therapeutic drug targets, in this case ErbB3, demonstrating promising translation from in vitro studies to in vivo validation.
National Cancer Institute Integrative Cancer Biology Program (U54-CA112967-03 to D.A.L.); American Cancer Society postdoctoral fellowship (PF-08-026-01-CCG to P.K.K.).
We gratefully acknowledge the helpful comments from J.Pritchard and our anonymous reviewers.
Conflict of Interest Statement: D.A.L. notes his Scientific Advisory Board relationship to Merrimack Pharmaceuticals.