In silico models of cancer

Software and Database Resources for Cancer Systems Biology

With the rapid advancement of technologies for high-throughput biological measurement, the construction of tools for the distribution, integration and assessment of this raw information has become increasingly important [49]. Multiple online databases (Table 1) have been created for the storage and distribution of genome-scale data, including transcriptomics [50, 51], regulatory sequences [52, 53], and proteomics [54]. Importantly, these resources provide easily accessible high-throughput data by which cancer models can be constructed and validated. This access enables the global study of disease in a manner not previously possible and will continue to be transformative for systems biology and systems medicine. Models of cancer at different scales of size and complexity often employ a priori biological information as constraints to impose parsimony on model architecture. Other web-based services provide annotation of both individual genetic elements [55–59] and their collective associations within distinct biochemical pathways or functional modules [60–65]. Additionally, published mathematical models and required parameters (e.g., reaction kinetics) can be obtained for diverse biological systems [66–68]. Among the major challenges for both statistical and biochemical cancer modeling is knowledge of protein interactions [69]. Multiple public repositories of human and cancer protein interaction networks are currently available based on yeast-two-hybrid data and other experimental methods [70–73]. Additionally, standardized markup languages for biochemical and signaling networks [74–77], as well as software packages [78, 79] have been developed to enable rapid porting of diverse datasets, annotations, and models between different programs and formats. The Systems Biology Markup Language (SBML) in particular is a machine-readable format for representing models of biological processes that has been evolving over the past several years [76]. SBML is becoming more widely adopted, both in the reporting and distribution of models, and in the import and export capabilities of software toolboxes.

The application of advanced analytical techniques to the examination of high-throughput experimental data will continue to facilitate assembly of cancer models with increasing biological precision. However, several practical challenges remain to be addressed which impede the successful integration of data acquired from different experimental sources, or analyzed through different analytical methods. For example, a recent report found that curated protein interaction datasets integrating multiple literature sources can be highly error-prone [80]. The fusion of data across multiple studies, performed at different times or by different laboratory groups, poses a substantial challenge for systems-level analyses that rely so strongly upon large and information-rich datasets. The Minimum Information for Biological and Biomedical Investigations (MIBBI) project is an example of ongoing efforts focused on developing and sharing ‘minimum information checklists’ for experimental data. Such checklists place more stringent guidelines on the quality and content of published data [81]. Many journals and funding agencies already require authors reporting microarray data to comply with the Minimum information About a Microarray Experiment (MIAME) checklist [82]. Other minimal information guidelines which are gaining support in their respective communities include Minimum Information About a Proteomics Experiment (MIAPE) [83], Core Information for Metabolomics Reporting (CIMR), and Minimal Information About a RNAi Experiment (MIARE). The ability to estimate and compare the quality of different data integration schemes in the context of protein interactomics has also been the subject of recent investigation [84, 85]. As these challenges continue to be addressed going forward, they promise the development of ever more powerful and biologically accurate in silico representations of human cancer.

Gene Expression Models for Cancer Diagnosis, Stratification, and Prognosis

The discrete mutational events found in the cancer genome and epigenome substantially modulate transcriptional profiles within cancer cells. Models based on these perturbed gene expression signatures can be applied for the diagnosis of disease subtypes, as well as for stratification of tumor grade or prediction of therapeutic response (Figure 3). These models apply deterministic or probabilistic rules to model the molecular context through which a cell or tissue progresses to malignancy. Different transcriptional classifiers have been developed for the discrimination of cancer type, subtype, and grade [89–91], including hierarchical clustering [92, 93], k-means clustering [94], support vector machines [95], artificial neural networks [96, 97], and classifiers based on the relative expression of gene pairs [98–101]. Transcriptomic signatures, which represent oncogenic processes as a function of gene expression levels, have also been applied to model relapse and overall survival in diverse cancers, including non-small-cell lung cancer and pediatric leukemia [102–105]. Additionally, transcriptional classifiers have been used to predict tumor response to chemotherapeutic agents [106–108], which represents one of the most important clinical determinations in oncology practice.

These models often begin with genome-scale microarray data, and through either computational or combined computational and experimental analyses, derive prognostic classifiers comprised of a smaller number of highly relevant transcripts. The statistical analysis of molecular signatures has potential utility for informing small molecule, radiological, and surgical cancer treatment choices that cannot be filled by standard histopathological or clinical analyses. While these transcriptomic models have not entered widespread clinical use, they embody a promising direct application of cancer models to medicine. Furthermore, such models may be more widely adopted as nucleotide measurement technology itself becomes more common in non-research contexts. Transcriptional models of cancer also provide an accessible characterization of the molecular processes that mediate tumor initiation, progression, and sensitivity to medical interventions.

Pathway Enrichment Models and Biological Discovery

The genomics revolution has prompted substantial efforts towards comprehensively annotating biological functions of and associations among identified genes [58, 86]. Several curated databases have been established for grouping genes into modules of defined activity according to physical or functional association, such as Gene Ontology [61] and the KEGG database [88]. Multiple computational methods have been developed for the purpose of mapping genetic and transcriptional activity into these groups [109, 110]. These tools use diverse numerical methods to determine gene groupings whose activity, as represented by gene expression, is substantially modulated through malignant transformation. For example, Gene Set Enrichment Analysis [111] and related tools have been applied to identify pathway perturbations in human cancers based on transcriptomic data, including tumors of the breast and prostate [112–114]. An alternative approach combines human protein-protein interaction network data with gene expression profiles for the classification of breast and lung cancer samples [7, 115]. These related studies have considered the cumulative activity of discrete pathways, as defined by protein interaction data, to discriminate between phenotypes. This approach yielded results of similar predictive accuracy to the more common individual gene-based classifiers discussed previously, and demonstrates the importance of integrating various data sources for biological modeling. Similar pathway-based analysis methods also have been applied as prognostic indicators of clinical tumor progression – and the most important pathways identified through these computational methods often overlap significantly with those known experimentally to be oncogenic [116–118].

These frameworks have been applied to diverse cancer systems and serve as a robust source of biological discovery at the genomic scale. That said, mapping genetic activity into discrete, human-defined gene sets clearly does not accurately recapitulate the full biological complexity encoded by historical evolutionary constraints. Rather, it is a useful tool that leads to further, more detailed studies. Indeed, a large portion of fundamental eukaryotic cell biology remains poorly characterized. This includes complex events such as non-coding transcription and epigenetic phenomena, but also more basic information such as the specific components and mechanisms of various signaling cascades. As such, pathway annotations themselves are subject to incomplete information and biases in curation. Nonetheless, these techniques, while coarse-grained, enable the simultaneous consideration of many cellular pathways and processes using high throughput data.

Network Inference at the Pathway Level

This scale of challenge to construct quantitatively predictive network models at the genome-scale in humans is not likely to be met in the near future [127], but more modest inference problems focusing on pathway-based subnetworks continue to prove useful. These studies examine different forms of biological data to infer causal relationships between genes at the transcriptional or signaling levels. For example, a multiparametric flow cytometry approach recapitulated a cell surface signaling pathway in the activation of human CD4+ leukocytes using a Bayesian network inference procedure [124]. Other studies have examined apoptotic signaling networks following experimental perturbation with pro-apoptotic TNF and other agents, using an innovative method in which downstream molecular perturbations are mapped into concordant ‘principal signaling components’ [128, 129]. These studies have shown that a small number of canonical signaling patterns can accurately describe even a complex cellular response such as apoptosis. Furthermore, by comparing signaling in different cell types with differing phenotypic responses, it was found that differential upstream signaling – at the level of kinase activation, for example – most strongly mediates cell-specific stimulus response, but that different signals are propagated through common intracellular networks.

A recent examination of receptor tyrosine kinases found that the complex and overlapping crosstalk involved in such signal transduction can be explained substantially by linear combinations for docking affinities for downstream proteins [130]. This study represents a systems-level characterization of an important cellular system that is too complex to be fully explicable using traditional ‘univariate’ biological modeling. Another recent study considered the cross-talk between phosphorylation pathways to predict the effects of inhibitor treatment on cancer cell migration phenotype [131]. This study characterized off-target and cross-talk effects of kinase inhibitor treatment to predict the effects of combinatoric treatments to forestall cell migration – an important malignant phenotype in cancer. Using different experimental and computational tools, these pathway-centric inference procedures have revealed important insights into diverse systems-level events in cancer that can have significant applications to pharmaceutical development and target selection.

These network inference studies incorporate a priori information in the design of the experiment itself – for example, in guiding the selection of what to measure (genes, proteins, and/or phosphorylation events) and the planning of perturbation experiments to informatively test the system. This model-guided approach constitutes a reduction and focusing of network search space that increases the fidelity of network inference. These methods are not easily extendable to the inference of larger and more complex networks, as the larger number of variables introduces a much larger amount of possible interactions, with associated experimental noise and statistical uncertainty. As larger networks are examined, the task of ‘model identifiability’ becomes an increasingly refractory but critically important challenge.

Network Inference at the Genome Scale

The inference of biological networks at the genome scale involves several challenges for eukaryotic systems, and thus most network inference methods have been developed initially using unicellular organisms. One example is the ‘Inferelator’ computational suite that has been used to infer a genome-scale transcription regulatory network in the unicellular archaeon Halobacterium [132, 133]. This study constructed a network model of transcription able to determine the level of expression of each gene given the expression state of transcriptional regulators. This model was able to predict quantitatively the response of the organism to novel stimuli that were not used to train the model itself. This type of quantitative inferred genome-scale transcriptional regulatory network in humans remains a difficult and ongoing challenge.

Nonetheless, efforts have been made to craft network-based statistical models of human cells, including in breast cancer and lymphoma, in which the architecture of regulatory networks for a portion of the human genome is characterized [134–137]. These models have applied Bayesian network inference methods [138] to discriminate physical or functional interactions between several thousand genes. The related probabilistic Boolean network formalism [139] has been used to construct regulatory networks in the context of several cancer systems, including glioma [140]. One study of melanoma adapted Boolean network analysis to include the use of ‘seed’ genes, in which groups of genes believed to interact were used as a subnetwork starting point which was then iteratively grown to incorporate other related genes [141]. Many network inference methods appear robust when validated on synthetic data, but exhibit substantially reduced performance when applied to experimental gene expression datasets [142–144]. These results speak to the substantial complexity of gene regulation networks – even in the context of unicellular organisms – and to the requirement for large datasets in network inference tasks [145].

Statistically inferred network models can be used to study the topology of complex cellular systems, and elucidate important genetic interactions and control points. Such statistical networks map the status of regulatory agents into qualitative interactions with cognate targets and discrete network states. Going forward, the acquisition of larger high-throughput datasets encoding different facets of regulatory interactions, combined with innovative methods for their integration and analysis, will enable the construction of increasingly precise numerical network models. The important ability to reliably model the behavior of cancer cells, for example upon response to treatment or genetic evolution, would hold significant implications for the development and administration of cancer treatments.

Stoichiometric Models of Biochemical Reaction Networks

The most basic mathematical representation of a biochemical reaction network is a stoichiometric model. Stoichiometric models describe the interconversion of biomolecules purely in terms of the number of reactants and products participating in each reaction. The generation of stoichiometric and analysis of their properties is a well-established process [146–148], and genome-scale models of metabolism have been completed for a diverse range of organisms, including the prokaryote E. coli [149], archaea [150], and eukaryotes such as S. cerevisiae [151] and the protozoan L. major [152]. Most important for cancer research is the recent reconstruction of known human metabolism at the genome-scale [153, 154]. Methods have also been developed for reconstructing signaling networks [155–157], transcriptional and translational networks [146], and regulatory networks [158, 159]; these models are fundamentally analogous to reconstructed metabolic networks (Figure 5). The reconstruction of a biochemical reaction network results in a database of stoichiometric equations that can be represented mathematically to form the foundation of a genome-scale, computable model. Computational tools for constraint-based analysis are then used to interrogate the properties of the reconstructed network in silico, and to facilitate model-driven validation and refinement [160]. Physico-chemical and environmental constraints under which the network operates are applied in the form of balances, including mass, energy, and charge, and bounds, such as flux capacities and thermodynamic constraints. The statement of constraints defines a solution space comprising all non-excluded network states, thereby describing possible functions or allowable phenotypes. These methods are now being adapted for modeling human systems in greater detail.

Constraint-based analysis of biochemical reaction networks has been applied to a number of human systems. Using the reconstruction of the human mitochondrial metabolic network [161], linear programming and random sampling were applied to identify candidate steady states of the network under normal, diabetic, ischemic, and dietetic conditions [162]. In a related study, Monte Carlo sampling of flux spaces was used to study the effects of enzymopathies on the human erythrocyte metabolic network [163]. A stoichiometric formalism was also used to reconstruct reactions comprising the JAK-STAT signaling network in the human B-cell [164]. In this study, a topological constraint-based approach was employed to evaluate network crosstalk, redundancy in signaling inputs and outputs, and to delineate correlated reaction sets or systemic modules. The recent completion of a global reconstruction of the human metabolic network [153] represents a significant milestone in human systems biology – and clears a path forward for reconstructing metabolic models for all of the 200 cell types in the human body and their modified forms in various cancers. This reconstruction comprises 1496 genes and 3748 reactions divided into 88 metabolic pathways. In addition to the typical network capabilities determined by constraint-based modeling, the initial genome-scale reconstruction has enabled analysis of relationships between network topology and human metabolic diseases [165]. An independently reconstructed human metabolic network [154] – comprising 2322 and 2823 reactions among 70 different pathways – was also used to demonstrate the potential of systems modeling in human metabolism to aid in drug discovery [166]. An important effort going forward will be to resolve the discrepancies between different versions of the human reconstructions, establishing a consensus global reconstruction of human metabolism.

Modeling Metabolism in Human Cancer

The global human metabolic reconstruction provides a basis for the known set of metabolic reactions catalyzed by human proteins. However, the utility of these models for cancer research going forward depends upon overcoming several challenges. First, further refinement of the global human metabolic map is essential to increase its accuracy. Second, each of the approximately 200 cell types in the human body exhibits only a portion of the full metabolic capability contained in the genome [167] (Table 2). The high percentage of undetermined activities for metabolic enzymes in human tissues clearly shows how much more we have to learn about even this very well studied cellular process. Effectively representing which portions of the global human metabolic network are active in any given cell type, and at what level, is thus of critical importance. Cancers in particular are known to exhibit diverse metabolic phenotypes compared to their progenitor cells, typically with an increased rate of overall metabolic activity to support their increased growth and the highest metabolic activity observed in the most aggressive malignancies [168, 169]. Multiple other hallmarks of cancer including angiogenesis, metastasis, evasion of apoptosis, and avoidance of immune detection have been previously linked to human tumor metabolism [170, 171]. Metabolic targets have also been used in cancer chemotherapy [172–174]. For these reasons, metabolic networks in human cancer have the potential to be a rich focus area for systems modeling going forward.

Kinetic Models of Biochemical Reaction Networks

Detailed dynamic models are needed to accurately quantify many important molecular processes, such as feedback and feedforward regulation, competitive inhibition, post-translational modification, and transcriptional regulation. At scales much smaller than the genome, biochemical reaction networks have been used as the basis for forming dynamic differential equation models with the addition of kinetic rate constants [155, 175, 176]. Many of the fundamental considerations and challenges associated with constructing kinetic models of cell signaling pathways have been previously reviewed in detail [177, 178]. Briefly, these models are most commonly formulated as a set of ordinary differential equations (ODEs). In these ODE-based pathway models, the rates of production and consumption of individual biomolecular species is described in terms of mass action kinetics (elementary reactions with forward and reverse rate constants). As the full kinetic detail for reactions in biological pathways is often unknown, simplifying assumptions such as the Michaelis-Menten approximation or time-scale separation can be applied. Kinetic models are advantageous in that they do not employ a steady-state assumption as constraint-based models typically do, and thus can simulate detailed dynamic behavior. Dynamic models are also able to simultaneously account for both the concentrations of compounds and flux through reactions. Capturing dynamics however, requires more parameters, and thus these models are more data intensive to create, as well as more prone to overfitting when parameters must be estimated. Nevertheless, smaller scale dynamic models have been used successfully to study in mechanistic detail key pathways related to human cancer. Most of the pathways and molecular components listed below do not function independently to promote oncogenesis, but rather are often connected through signaling crosstalk, feedback mechanisms, and other forms of up- or downstream regulation. These examples are not an exhaustive list, but serve to demonstrate the benefits of detailed kinetic modeling of biological pathways.

Aberrant activity of the transcription factor NF-κB has been linked to oncogenesis, tumor progression, and resistance to chemotherapy. A computational ODE model was employed to determine the role of inhibitor of NF-κB kinase (IκB) isoforms in the temporal control of NF-κB [155]. Analysis of the model revealed that IκBα is responsible for strong negative feedback and fast turn-off of the NF-κB response to IKK stimulation, while IκBβ and IκBε reduce oscillations in the signaling module and stabilize NF-κB response during longer stimulation. The same model was expanded and used to study dynamics of IκB–NF-κB signaling module when cells were stimulated by lipopolysaccharide (LPS) via Toll-like receptor 4 (TLR4) [175, 179]. Modeling of two TLR4-dependent signaling pathways suggested that one pathway required the expression tumor necrosis factor-α (TNFα) to activate NF-κB.

In another study, the negative feedback loop between the tumor suppressor p53 and the Mdm2 oncogene was examined experimentally and computationally [180]. Different mathematical models of the system were analyzed to draw connections between predicted and experimentally observed behavior. The models indicated that low-frequency noise in protein production rates was an important source of variability in oscillations of the p53-Mdm2 system. Another model was built for the single-cell response of p53 to radiation-induced DNA damage [181]. In agreement with experiments, simulations showed that by assuming stochasticity in the initial number of double-strand break sites and the DNA repair process, p53 and Mdm2 exhibit coordinated oscillatory dynamics with radiation stimulation.

Numerous components of the phosphatidylinositol-3-kinase (PI3K)/AKT pathway are targeted by malignant perturbations such as amplification, mutation, and translocation in cancer patients. Activation of the PI3K pathway results in aberrant cell growth and survival, metastatic competence, and therapy resistance. An excellent review of the challenges, opportunities, and progress in modeling and exploiting the PI3K pathway can be found in Hennessy et al., 2005 [182]. For example, inhibitors of the mTOR pathway have been demonstrated to directly target PI3K signaling. However, due to the intricate role of mTOR in regulation of PI3K/AKT signaling, detailed understanding of the pathway derived from computational modeling is needed to determine the consequences of inhibition of particular components.

A number of models have been used to examine mechanisms that govern mitogen-activated protein kinase (MAPK) pathway dynamics. MAPK pathways are involved in the activation of many cancer-promoting genes and respond to a number of different cancer-associated stimuli, such as epidermal growth factor (EGF). A computational model of general MAPK signaling pathways was used to study the role of negative feedback mechanisms in generating complete signal adaptation (i.e., desensitization to stimulus) [183]. Another model was used to simulate EGF-induced MAPK activation [184]. Analysis of model simulations provided quantitative evidence that feedback inhibition of the MAPK cascade is the most important factor in determining the duration of cascade activation. It was also found that differences in EGF and NGF-stimulated signaling could be accounted for by differential feedback regulation of the MAPK cascade. Similarly, a model was constructed to examine the signal-response relationships between EGF binding and the activation of downstream proteins in the cascade [185]. Model predictions showed that EGF-induced responses are remarkably stable over a 100-fold range of ligand concentration.

Another model used dynamic simulations of gene expression and signaling networks to study behavior in human cancer cells [186]. While dynamic interaction networks of cancer systems are challenging to construct and simulate, such models are essential to capture the molecular complexity of oncogenesis. Taken together, these studies provide powerful examples of how building a model that correctly explains data and makes verifiable predictions can lead to deepened insight into biological mechanism and function.

Continuum Models

Extracellular parameters can be represented as continuously distributed variables to mathematically model cell-cell or cell-environment interactions in the context of cancers and the tumor microenvironment. Systems of partial differential equations have been used to simulate the magnitude of interaction between these factors, including the effects of hypoxia on cell cycle progression [192], the impact of mechanical forces on tumor invasiveness [193] and extracellular matrix interactions [194]. Recent studies have examined cell population dynamics within colonic crypts in colorectal cancer [195, 196]. These models consider interactions between stem cells, differentiating cells, and differentiated cell populations to quantitatively predict tissue-level invasion and the growth of tumor mass. Other models have represented solid tumors as a multiphase system of both bound and ‘mobile’ forms [197, 198]. Such ‘mixture’ models consider differential growth and apoptosis rates, as well as mass transfer and regulatory interactions between phases. Alternative models have considered nonlinear and combinatoric effects of multiple factors, including nutrient availability and mechanical parameters [199, 200], and the effects of mutation rate on invasion and metastasis [201]. These numerical systems embody a robust method to incorporate the effects of somatic biological phenomena into computational representations of cancer. Continuum-based models are thus a powerful tool to simulate and characterize interactions between intracellular and extracellular factors in oncogenic processes.

Cellular Automata and Agent-Based Models

Multivariate continuum models are able to represent the effects of several physiological or biochemical events on cancer development. However, in situ, these factors are highly heterogeneous, and interact discontinuously with tumor cells [202, 203]. Cellular automata models represent cancer cells as discrete entities of defined location and scale, interacting with one another and external factors in discrete time intervals according to predefined rules [204]. Agent-based models expand the cellular automata paradigm to include entities of divergent functionalities interacting together in a single spatial representation, including different cell types, genetic elements, and environmental factors [205]. With sensitivity to starting conditions, and the ability to incorporate probabilistic interactions at each time step, these models can exhibit similar stochastic behaviors to those observed in oncogenesis in vivo [206–208]. Phenomena that have been modeled using agent-based models include three-dimensional tumor cell patterning [209], immune system surveillance [210], angiogenesis [211–213], and the kinetics of cell motility [214]. Another recent model integrated diverse parameters such as extracellular signals, blood flow, and tissue degradation to simulate the spatiotemporal formation of tumor vasculature [215].

Increasingly, ‘hybrid’ models have been created which incorporate both continuum and agent-based variables in a modular approach. For example, a recent study considered continuous extracellular biomolecule distributions and discrete cell locations to simulate the interaction between intracellular decision-making processes and malignant growth [216]. Another recent model incorporated a continuous model of a receptor signaling pathway, an intracellular transcriptional regulatory network, cell-cycle kinetics, and three-dimensional cell migration in an integrated, agent-based simulation of solid brain tumor development [217]. The interaction between cellular and microenvironment states have also been considered in a multiscale model that predicts tumor morphology and phenotypic evolution in response to such extracellular pressures [218]. These and other techniques which incorporate multiple, nested scales of interacting biology embody promising paradigms to understand cancer as a cascade of information across levels of size and complexity. This ability to interrogate cancer across multiple biological agents and compartments presents a unique framework to elucidate oncogenic processes, and to evaluate potential therapeutic interventions through digital simulation prior to experimental deployment.