|Home | About | Journals | Submit | Contact Us | Français|
Simulating cancer behavior across multiple biological scales in space and time, i.e., multiscale cancer modeling, is increasingly being recognized as a powerful tool to refine hypotheses, focus experiments, and enable more accurate predictions. A growing number of examples illustrate the value of this approach in providing quantitative insight on the initiation, progression, and treatment of cancer. In this review, we introduce the most recent and important multiscale cancer modeling works that have successfully established a mechanistic link between different biological scales. Biophysical, biochemical, and biomechanical factors are considered in these models. We also discuss innovative, cutting-edge modeling methods that are moving predictive multiscale cancer modeling toward clinical application. Furthermore, because the development of multiscale cancer models requires a new level of collaboration among scientists from a variety of fields such as biology, medicine, physics, mathematics, engineering, and computer science, an innovative Web-based infrastructure is needed to support this growing community.
Cancer is a class of diseases characterized by out-of-control cell growth and tissue invasion. Because cancer initiation seems to depend on a series of genetic mutations affecting intrinsic cellular programs, to date, the vast majority of cancer research has focused on the identification and characterization of these genetic and molecular properties of cancer cells themselves (1). However, tumors are also heterogeneous cellular entities whose growth is dependent upon dynamical interactions among the cancer cells themselves, and between cells and the constantly changing microenvironment (2). For example, such interactions include signaling through cell adhesion molecules such as cadherins and integrins (3) and differential cell responses to growth factors and other external signals (4). All of these interactive processes act together to control cell phenotypic behaviors such as proliferation, apoptosis, and migration. There is increasing consensus that these dynamic interactions cannot be investigated purely by using biological experiments, since experimental complexity usually restricts the accessible spatial and temporal scales of observations. Consequently, as of late, there is a transition within the cancer research community to treat and study cancer as a systems disease.
Cancer systems biology is a newly emerging field that uses an interdisciplinary approach to provide the systemic understanding of cancer initiation and progression by investigating how individual components interact to give rise to the function and behavior of the cancerous system as a whole (5). It seeks to decipher emergent behavior rather than to focus only on the system’s constituents’ properties. In addition to conventional biological and medical experiments, a systems approach often involves mathematical and computational modeling to interpret and integrate the massive amount of data that experimentalists are currently uncovering, especially in molecular and cell biology (6). In silico cancer models are necessarily simplified, yet we argue can provide adequate representations of a particular cancer phenomenon to be investigated. In fact, such a theoretical approach has been increasingly recognized as having the capability 1) to simulate experimental procedures and to optimize and predict clinical therapies and outcome, and 2) to test and refine medical hypotheses (7). The development of a successful in silico cancer model is a long-term process, expected to be iteratively conducted, with available experimental data used to guide the model design, and to verify and validate model results.
Most current computational cancer models focus on cancer behaviors that occur at a single biological scale (7), and the decades of dedicated efforts by cancer modelers have made scale-specific models not only possible, but advanced enough to be practical for applications in oncology (8). However, because cancer growth and invasion indeed span multiple scales, using a scale-specific model is insufficient to uncover cross-scale mechanisms let alone to render predictions about the clinical outcome of the disease system as a whole (9). Modeling cancer across different biological scales has recently begun to play a more important role in moving the field of cancer systems biology towards clinical implementation. In the context of biology and physiology at large, a model is considered to be “multiscale” if it spans two or more different spatial scales and/or includes processes that occur at two or more temporal scales. Because a multiscale cancer model has to quantify parameters on, and relationships between biological processes that occur at different scales, the complexity of model development is significantly increased.
In this review, we present representative works that have investigated key questions on cancer progression, invasion, angiogenesis and metastasis using a multiscale modeling approach. Systems approaches emphasize the integration and coordination between computational modeling and experimental efforts. Hence, attention will be given to the evaluation of the models’ capacity to account for in vitro, in vivo, or clinical data, and more importantly, an analysis of how the computational findings drive new hypotheses for further experimental investigation.
Multiscale cancer modelers already have a wealth of useful, mostly scale-specific resources to refer to or base their innovative work on, but they also face the enormous challenge of developing more realistic and more accurate predictive models. The fundamental reason is that, when considering the increasing number of components at multiple scales (whether in time or space), more model parameters and the relationships between them will have to be defined, quantified, and frequently adjusted according to data from the literature, from experiments or clinics. Additionally, defining the linkage between different scales poses a significant barrier to model development; in many cases, the link is bidirectional, meaning that higher-and lower-level variables, parameters, and functions characterizing the model are influenced by each other. Furthermore, it has proven challenging to adaptively interface discrete (individual-based) and continuum (population-based) models while 1) ensuring mathematical consistency between the models, and 2) adequately conserving mass and momentum when switching between the models; indeed, these are continuing challenges in cutting-edge scientific computing (10-11). Generally, when developing a multiscale cancer model, particularly if with a clinical application in mind, it is imperative to maximize the accuracy and predictive power of the model while at the same time constraining the size of the model as much as possible in order to produce tractable results within a reasonable time frame.
Various definitions of biological scales have been created in different life sciences fields (see (12) for an excellent review). We focus our discussion of multiscale on four main biological spatial scales: atomic, molecular, microscopic, and macroscopic (Figure 1). While each of these spatial scales may have multiple temporal scales, biological processes that take place at a lower-scale generally happen much faster than those at a higher-scale, i.e., spatial and temporal scales tend to vary together, with slower to faster at the temporal scale corresponding to smaller to larger at the spatial scale. In the following, we briefly explain what processes occur at each of these four spatial scales and the main methods employed to simulate these processes, all from the cancer modeling perspective and practice. We note here that, as also commented in (12), model design and development should take into account the nature of the system being modeled rather than simply follow an invariable hierarchical structure. That is, depending on modeling purpose, it is no violation to merge two of the levels discussed here or to split one or more of them into distinct sub-levels.
This scale is used to study the structure and dynamic properties of proteins, peptides, and lipids, as well as their dependence on the features of the environment or on ligand binding (13). The most common modeling method used at this scale is molecular dynamics (MD) simulation, where atoms and molecules are allowed to interact for a period of time. Atomic-scale models deal with length scales in the order of nm and time scales of ns.
Models at this scale do not represent the molecular dynamics of individual proteins, but represent an average of the properties of a population of proteins. Cell signaling mechanisms, the natural regulators of biological systems, are usually investigated at this scale (14). Analysis of this scale constitutes an intensely active field of biomedical research and has the potential to provide new therapeutic targets to combat disease. Signal transduction starts with the binding of extracellular molecules (ligands) to cell-surface receptors, and ends with a change in cell function. Most of the current modeling efforts focus on this scale, adding insights into quantifying signal-response relationships and signaling events that control cellular responses (15). Ordinary differential equations (ODEs) are often used to represent biochemical reactions contained in a signaling pathway. Molecular-scale models deal with length scales in the order of nm~μm and time scales of μs~s.
This scale is also referred to as the tissue or multicellular scale, and in our definition, it also includes the cellular scale (i.e., single cell behaviors and properties). Individual cells are contained in a selectively permeable cell membrane (16). Models at this scale must suitably describe the malignant transformation of normal cells, associated alterations of cell–cell and cell-matrix interactions, the heterogeneous tumor environment and the element of tumor heterogeneity. These models usually use partial differential equations (PDEs) or agent-based modeling (ABM) rather than ODEs to simulate these factors and processes. The simulation run time can increase substantially if individual cell behaviors are investigated in fine detail. Tissue-scale models deal with length scales in the order of μm~mm and time scales of min~hour.
Models at this scale focus on the dynamics of the gross tumor behavior including morphology, shape, extent of vascularization, and invasion, under different environmental conditions (17). Microscopic details of tissue structure are averaged over short spatial scales to produce a description of the macroscopic-level tissue properties. At this scale, because the number of cells in the model is sufficiently large, it becomes possible and sometimes necessary to treat some or all of the cells as a single continuum. This in turn allows for cell and substrate transport to be modeled with conservation laws for spatiotemporally-varying densities (i.e., PDEs), rather than keeping track of individual cell activities. Models generally consider cell responses to gradient fields of diverse origins, such as concentration gradients of diffusible or non-diffusible molecules as well as strain and stress gradients generated by the growing tumor mass. Macroscopic-scale models deal with length scales on the order of mm~cm and time scales of day~ear.
Since lower-level processes are much faster, it is reasonable to assume that they are in quasi-equilibrium with the slower, higher-level processes, i.e., lower-level processes can be included at the higher level via e.g. constitutive equations or force fields (18). For example, if a reaction occurs on a very fast time scale, then we can assume that the chemicals in that reaction are in equilibrium. This assumption eliminates one of the differential equations from the integrated system, making its solution more straightforward and less computationally intensive, while still maintaining the accuracy of the model. This type of multiscale modeling, where lower-level processes (small spatial scales, fast dynamics) are coupled to higher-level processes (large spatial scales, slow dynamics), has received the most attention in the current quantitative cancer research field, and is particularly useful when developing multiscale models because the system of governing equations is generally large.
Modeling cancer behaviors can involve techniques that are discrete, continuum, or hybrid, i.e. the integration of both [see Summary box]. Continuum models are capable of capturing larger-scale volumetric tumor growth dynamics (which are also accessible to conventional clinical imaging modalities) at a comparatively lesser computational cost (19). Continuum descriptions of tumor growth are also to some degree supported by fundamental physical principles, and thus benefit from the knowledge gained in this field (20). Despite these advantages, the averaging over space realized in continuum formalisms often cannot fully account for the diversity of cellular and sub-cellular dynamical features as well as genetic or epigenetic regulatory mechanisms exhibited at the individual cell level. In other words, a continuum technique is often a lesser choice when exploring tumor heterogeneity when the cell properties vary over small spatiotemporal scales, which is an inherent feature of cancer cells (21). Moreover, continuum models cannot easily describe processes where individual cell effects are important or dominate, such as the epithelial–mesenchymal transition (EMT) process (22). Alternatively, discrete models are suitable to address these shortcomings, since they operate at the scale of individual cells or cell clusters. Discrete models can easily incorporate biological rules (based on biomedical data or data-driven assumptions), such as those defining cell-cell and cell-matrix interactions involved in both chemotaxis and haptotaxis. However, discrete techniques also have drawbacks. The most serious one is the large computational demand when modeling each cell in fine detail, which limits the model to a relatively small number of cells. As a result, a typical discrete model is usually designed with a sub-millimeter or lower domain size (23). Furthermore, some effects (e.g., tissue biomechanics) are best described at the macroscopic scale with continuum techniques (11). For these reasons, cancer modelers are increasingly turning to hybrid techniques that combine the benefits of both continuum and discrete descriptions. In fact, in current discrete-based cancer models, extracellular factors are often modeled as continuous quantities, thereby rendering the models hybrid in nature (9, 23).
Regardless of which techniques are used to develop a cancer model, all of them are abstractions of reality. All hypotheses on which the models are based will eventually prove to be incomplete in one way or another (27). However, models have the advantage of being quantitative and interactive rather than solely descriptive. Thus, modeling and experimentation in cancer research can and should work cooperatively, supplement, and promote each other. In the following, we highlight some of the most recent and representative multiscale modeling works that demonstrate the importance and necessity of this approach in current cancer research.
Signaling pathways induced by growth factors and mediated by receptors are responsible for transmitting information into the cell to regulate a diverse array of cellular processes including proliferation, migration, differentiation, and apoptosis. Molecular alterations (mutation, overexpression, or underexpression of genes or proteins) of these pathways contribute to cancer initiation and progression. It has been found that differential phosphorylation rates are associated with different phosphorylation docking sites in receptors, resulting in differing downstream signaling response. Therefore, understanding the effects of a receptor’s different docking sites on preferential signaling will help differentiate signaling characteristics of mutant cell lines derived from cancer patients. In this section, we introduce such a modeling approach linking atomic and molecular scales that can be used to investigate the effects of point mutations on preferential downstream signaling response.
Epidermal growth factor receptor (EGFR) is frequently overexpressed and mutated in a variety of cancers, and as a result, small molecule tyrosine kinase inhibitors for EGFR are of significant interest as anti-cancer drugs (28). It has been reported that cytosolic signaling proteins bind to different phospho-tyrosine docking sites of EGFR, which may cause differential patterns of EGFR downstream signaling (29). Hence, quantifying the difference in wild-type and mutant EGFR signaling and accounting for differential signaling through different phospho-tyrosines docking sites of EGFR is useful in identifying the role and significance of drug sensitizing mutations of EGFR in cancers.
Radhakrishnan and co-workers have been developing multiscale models across atomic and molecular scales to understand the mechanisms of how altered signaling induced by point-mutations in EGFR leads to the onset of oncogenic transformations. In models presented in (30-31), they mathematically implemented a wild-type and a mutant EGFR activation mechanism, which differ in that the wild-type receptor tyrosine kinase (RTK) initiates phosphorylation of C-terminal tail substrate tyrosines only as a dimer, whereas the mutants can initiate phosphorylation as a monomer and as a dimer. In order to translate these context-specific phospho-tyrosine mechanisms into differences in the downstream response, they introduced a branched signaling approach through EGFR. As shown in Figure 2, two parallel phosphorylation pathways are implemented, corresponding to tyrosine 1068 (Y1068) and tyrosine 1173 (Y1173), where phosphorylated Y1068 (pY1068) binds only to Gab-1 and Grb2, while phosphorylated Y1173 (pY1173) binds only to Shc. Accordingly, phosphorylation events at the Y1068 and Y1173 sites transduce signals through different signaling routes. This approach allows for transcribing the effects of somatic mutations (L834R and Del) in the EGFR to different downstream responses in phosphorylated ERK (ERK-p) and Akt (Akt-p) in the EGF-induced signaling pathway (31-32). Using these models, the authors predicted that the ratio of Akt-p/ERK-p increases in mutants compared to wild-type (33), and found that the perturbation of the phosphotyrosine kinetics of Y1068 and Y1173 through mutations is directly responsible for the differential signaling leading to preferential Akt activation. They also investigated the inhibitory effects of small molecule tyrosine kinase inhibitors on EGFR phosphorylation and downstream ERK and Akt activation.
Another major finding of their models is that the clinically identified mutations of EGFR kinase induce network fragility in the stabilizing interactions of the inactive kinase conformation. This in turn provides a persistent stimulus for kinase activation even in the absence of any growth factor. Moreover, parameters that drive network hypersensitivity through the enhancement of phosphorylated ERK and Akt levels show a striking correlation with observed mutations of specific proteins in oncogenic cell lines as well as the observed mechanisms of drug resistance to EGFR inhibition. Therefore, the authors suggested that cascading mechanisms of network hypersensitivity and fragility may enable molecular-level perturbations (clinical mutations) to induce oncogenic transformations and mechanisms of drug resistance (31-32). As a specific example, they described a possible mechanism for preferential Akt activation in non-small cell lung cancer (NSCLC) harboring EGFR activating mutations. This preferential Akt activation makes these mutant cell lines conducive to oncogenic addiction, and this pathway addiction mechanism also results in a remarkable sensitivity to EGFR kinase inhibition. These structural studies on kinase activation can be used to forecast the mutation landscape associated with other ErbB family members that may have not yet been reported (34).
Cancer models across molecular and microscopic scales are important and necessary due to (at least) two characteristics of cancer. One is that cancer is widely viewed as a disease involving irreversible genomic changes affecting intrinsic cellular programs (1), and thus it is difficult for wet-lab cancer researchers to trust a cancer model which misses correlations of molecular-level alterations with cancer cell properties. The other is that cancer is a context-dependent disease (35): depending on different microenvironments, cancer cells and the emergent tumor exhibit different phenotypes and behaviors. Cancer models introduced in this section are discrete-based models because they need to function on a single cell level, and most of them use ABM technique. ABM is of particular interest to cancer modelers because it helps to address the role of diversity in cell populations and also within each individual cell, which in turn enables us to explore how cancer growth and invasion patterns (due to cell proliferation and migration) emerge as a result of individual dynamics, including cell-cell and cell-environment interactions and intracellular signaling of individual cells. In the following, we introduce the development of some of the most recent molecular-microscopic cancer models.
Deisboeck and coworkers have been working extensively on the development of molecular-microscopic ABMs to simulate tumor properties within both brain tumors and NSCLC. Their models aim to quantify the relationship between extracellular stimuli, intracellular signaling dynamics, and multicellular tumor growth and expansion. At the molecular level are signaling pathways induced by growth factors and mediated by growth factor receptors. Similar to general pathway analysis studies, the pathway models are implemented using ODEs. At the microscopic level is a biochemical microenvironment constructed to represent a virtual tissue in a two dimensional (2D) or three dimensional (3D) domain. Heterogeneous environments are attained by distributing external diffusive chemical cues (such as growth factors, glucose, and oxygen tension) throughout the microenvironment. Throughout the simulation, the concentrations of the chemical cues are continuously diffused and updated at a fixed rate based on PDEs. Each cell carries a self-maintained signaling pathway or networked pathway combination, meaning that cells will differ in signaling profiles and thus will exhibit different phenotypic behaviors as the simulation progresses.
A molecular-level module, i.e., an EGFR signaling pathway, was first incorporated into a 2D brain tumor model to study and describe how context-dependent single-cell activities potentially affect the dynamics of the entire tumor system (36). In this model, an experimentally-supported molecularly-driven cellular phenotypic decision algorithm was established to link molecular changes to cell phenotype determinations. It is noteworthy that, for processing phenotype transitions, a series of subsequent modeling works all adopt this algorithm directly, or a variety of it. We briefly demonstrate how this algorithm determines the cell’s migration fate. Experimental studies have shown that the transient acceleration of accumulating phospholipase Cγ (PLCγ) levels leads to cell migration (37), and thus the rate of change of PLCγ can be used to determine the cellular migration decision. In this algorithm, the potential for any individual tumor cell to migrate was assessed by evaluating Mn[(PLCγ)]=[d(PLCγ)/dt]n, where d(PLCγ)/dt is the change in concentration of PLCγ over time, t, for a cell with an ID, n; if Mn exceeds a pre-specified threshold, then the cell becomes eligible to migrate; otherwise it can proliferate or become quiescent (i.e. it neither proliferates nor migrates yet remains viable). Note that a cell additionally has to meet other microenvironmental requirements, such as sufficient local nutrient conditions and available adjacent space, in order to proliferate or migrate successfully. In a follow-up study (38), increasing the EGFR density per cell was found to lead to an acceleration of the entire tumor system’s spatiotemporal expansion dynamics in accordance with experimental data (39). To simulate brain tumor growth in a more realistic microenvironment, a 3D model was developed (40), where a simplified cell-cycle description based on (41) and a more complicated extracellular matrix representation were implemented. Simulation results indicated that over time, proliferative and migratory cell populations oscillate and have a direct effect on the entire spatiotemporal tumor expansion pattern. Using the 3D model (40), an element of genetic instability was added, to investigate how heterogeneity impacts brain tumor progression patterns (42). The extended model found that cell clones with higher EGFR density are comprised of a larger migratory fraction and smaller proliferative and quiescent fractions, a result that agrees well with reported experimental data (43).
EGFR signaling also plays a vital role in the pathogenesis and progression of NSCLC (44). Deisboeck, Wang and colleagues attempted to look into how the methods applied to modeling brain tumors could be applied to the case of NSCLC. They also explored the possibility of using a cross-scale approach to determine critical molecular parameters that have a significant impact on tumor outcome at the microscopic level (45). A multiscale NSCLC model was first developed, with a revised EGF-induced EGFR-mediated signaling pathway implemented at the molecular level and a 2D heterogeneous biochemical tumor growth environment implemented at the microscopic level (46). This model revised the previous cellular phenotypic decision algorithm used in the brain tumor models by adding another decision molecule, ERK, in determining the cell’s proliferation fate. This was also based on an experimental study which reported that the transient acceleration of accumulating ERK concentration levels triggers cell replication (47). Figure 3a demonstrates this essential algorithm. With this 2D NSCLC model, they found that a minimal increase in EGF concentration can temporarily abolish the proliferative phenotype in the cancer cell closest to the nutrient source. More recently, Wang et al. (48) presented a 3D model in which both EGF and transforming growth factor beta (TGFβ) as well as their interplay were taken into account. This model was used to investigate how the effects of individual and combined changes in EGF and TGFβ concentrations at the molecular level alter tumor growth dynamics on the multicellular level, i.e., tumor volume and expansion rate. As shown in Figure 3b, when EGF and TGFβ concentrations were jointly varied asynchronously, a particular region of tumor system stability, generated by unique pairs of EGF and TGFβ concentration variations, was discovered. The simulated tumor system became sensitive to external variations in EGF and/or TGFβ when they occurred outside this robust region. The expansion rate for the standard simulation (with all kinetic parameters set to their reference values) was 2.07 μm/hr, which is in very good agreement with both computational modeling (49) and experimental studies (50).
The importance of the tumor microenvironment is currently of great interest to both the biological and the modeling communities (51). To investigate the impact of the microenvironment on the tumor growth dynamics, a series of hybrid ABMs were developed by Gerlee and Anderson. Each cell in their models was equipped with a microenvironment response network modeled using a feed-forward artificial neural network. In the neural network, microenvironmental variables such as local oxygen concentration, glucose concentration and extracellular matrix (ECM) gradient were represented as the input layer, regulatory genes as the hidden layer, and the response for cell phenotypes as the output nodes. This way, a genotype-to-phenotype mapping with the use of a neural network approach was established, in order to determine the actions of the cell based on its genotype, the microenvironment in which it resides, and their interactions. Furthermore, the neural network is subject to mutations when the cells divide. This means that the behavior of cancer cells can change from one generation to the next, implying that the models have the capability to capture the evolutionary dynamics of tumor growth. In (52), it was revealed that tumors grown in low oxygen concentrations exhibited branched morphologies, and the oxygen concentration influenced the evolutionary dynamics. A subsequent extension of the model involving the effect of the ECM and anaerobic metabolism was used to examine the emergence of a glycolytic phenotype and the influence of the tissue oxygen concentration and ECM density on the dynamics of the model (53). Simulation results showed that this glycolytic phenotype was most likely to occur in low oxygen concentrations and within a dense ECM. Moreover, it was observed that, while a low oxygen concentration results in branched tumor morphology, increased ECM density gave rise to more compact tumors with less fingering morphology. Figure 4 shows a series of simulation results with this model. More recently, these authors investigated the impact of the microenvironment on the emergence of a motile invasive phenotype in an evolving tumor population using the same approach (54). The model focused on haptotaxis, i.e., directed cell movement along ECM gradients in the tissue, which is known to be the dominant mechanisms in tumor invasion (55). Results showed that in identical simulation conditions, the tumor’s evolutionary dynamics converge to a proliferating or migratory phenotype, which suggests that the introduction of cell motility into the model changes the shape of the fitness landscape on which the cancer cell population evolves.
E-cadherin mediates cell-cell adhesion and plays a critical role in the formation and maintenance of cell contact, and loss of E-cadherin-mediated adhesion is a key feature of the EMT. Only recently, cancer modelers began to explore linking the intracellular signaling dynamics-related cell-cell adhesion and the extracellular consequences in invasive tumors. A molecular-microscopic multiscale agent-based lattice-free model was first developed for this purpose by taking into account the intracellular dynamics of the E-cadherin and β-catenin interactions and the physical forces on the cells (56). The model focused on a simplified β-catenin pathway which captures the key features of the cell adhesion process. Simulation results showed that down-regulation of β-catenin can be mainly driven by cell-cell contacts, and EMT can be achieved and reversed depending on the regulation of soluble β-catenin by local contacts. The intra- and intercellular protein interactions that govern cell–cell adhesion combined with cellular physical properties are also the driving forces of an essential mechanism that a cancer cell uses to attach to the endothelial wall, i.e., transendothelial migration (TEM) (57). Thus, with necessary extensions and modifications to the previous model (e.g., the addition of the Src pathway, which plays a principal role in TEM), the same modeling group examined the influence of different protein pathways in the achievement of TEM (58). Four cancer cell genotypes that differ in the adhesion protein pathways were considered. The genotypes were characterized by their capacity of creating N-cadherin-mediated bonds with the tunica intima and by their capacity of inducing a detachment of the endothelial–endothelial bonds by Src activity. Their results indicate that the slowest migration was found in the case when both N-cadherin and Scr were knocked out, while the fastest case occurred when both N-cadherin and Scr remained active.
Molecular-microscopic modeling has enjoyed good success in predicting the local behavior of cancer on scales of hundreds of microns to millimeters, and has proved a valuable tool for investigating the links between our current wealth of experimental molecular and cellular biology data and complex, emergent behavior of large multicellular systems. However, this approach is too computationally intense to simulate full tumors and their surrounding tissues. Furthermore, some important aspects of tumor biology, such as mechanical stresses, are best characterized by macroscopic, continuum models (10-11). Hence, there is need to dynamically combine microscopic models to describe important cell-scale phenomena (e.g., EMT) with the efficiencies of continuum models.
Current microscopic-macroscopic models can be roughly categorized (with increasing levels of multiscalarity) as 1) composite hybrid models, where a microscopic model of one phenomenon is coupled with a macroscopic model of another, 2) continuum models with functional parameters, where continuum models express microscopic effects through heterogeneous, time-dependent parameters, and 3) adaptive hybrid models, where one or more phenomenon is represented using both microscopic and macroscopic descriptions, with a physically-motivated means of adaptively choosing the appropriate characterization. We focus here on representative examples of these three approaches; interested readers may also consider recent reviews (10) and books (11, 59).
To date, most hybrid microscopic-macroscopic models have been composite, where a discrete cell-scale representation of one phenomenon is combined with a continuum tissue-scale model of another. Some of the most important examples have been coupled models of tumor growth (using a continuum model) with angiogenesis (using a discrete model) (60-65). Zheng et al. (60) used a sharp interface continuum model of tumor growth, where the moving tumor boundary is modeled as the boundary between two incompressible fluids, with a discrete cellular automaton model of angiogenic sprout tip motion first introduced by McDougall, Chaplain, and Anderson (66-70). In the viable rim of the tumor region (where substrate—oxygen, glucose, and growth factors—levels were higher than a threshold value), the tumor volume increased with a proliferation rate proportional to the substrate level. In the necrotic core (where the substrate level fell below the threshold), volume was lost due to cell lysis and subsequent fluid flux. Cell mechanics were modeled by introducing a mechanical pressure P that was related to the local tissue velocity v by Darcy’s law (v = −μP), and cell-cell adhesion was modeled as a surface tension along the tumor boundary. The boundary of the necrotic region released angiogenesis-promoting factors that diffused out of the tumor, reached the sprout tips, guided their chemotaxis and branching behavior, and ultimately determined the vascular topology. Simultaneously, the authors solved for the local substrate concentration throughout the tumor and surrounding tissue using a quasi-static reaction-diffusion equation. The models were further coupled by releasing substrate from the neovasculature and by restricting this nutrient source whenever the tumor pressure exceeded a threshold value. This model was able to reproduce tumor morphological instability caused by substrate heterogeneity, as well as biased tumor growth along the neovasculature (Figure 5a).
More recently (64), Macklin et al. coupled a refined sharp interface model by Macklin and Lowengrub (71-74) with a more advanced discrete model of angiogenesis, i.e., dynamic adaptive tumour-induced angiogenesis (DATIA), by McDougall, Chaplain, and Anderson (70). In the DATIA model, blood flow (including the effect of hematocrit) is solved throughout the vasculature as a coupled system of Poiseuille flow equations in individual segments of the neovasculature. The vessel segment radii vary according to the balance of wall shear stresses, vessel fluid pressure, and vessel elasticity—this, in turn, affects transport throughout the vasculature. The refined tumor model introduced a hypoxic region, where substrate levels were too low to promote cell proliferation but not so low as to induce necrosis. This region, rather than the necrotic core, released angiogenesis-promoting factors that drove chemotaxis in the angiogenesis model. The proliferating rim released matrix metalloproteinases (MMPs) that degraded the ECM in and near the tumor; this, in turn, affected the development of the vasculature due to haptotaxis of the sprout tips, as well as the evolution of the tumor by altering the distribution of the mechanical pressure. As an advance over the work by Zheng et al. (60), oxygen was released by the neovasculature at a rate proportional to the hematocrit level (to model oxygen transport by red blood cells). Furthermore, the vessel radius was determined by the balance of vessel pressure, wall shear stress, vessel elasticity, and the mechanical tumor pressure, leading to tumor vessel collapse as an emergent phenomenon (Figure 5b). This work was able to capture the complex dynamics of 1) hypoxia-induced angiogenesis, 2) increased oxygen supply by the vasculature following anastomosis of the vessels, 3) rapid subsequent tumor growth leading to increased proliferation-induced mechanical pressure, followed by 4) new regions of hypoxia due to vessel collapse (64). The results led the authors to hypothesize repeated hypoxia-angiogenesis-growth cycles due to the combined effects of tumor shape instabilities, mechanical collapse of the neovasculature, and heterogeneous oxygen distributions. Similar recent work by Cristini, Frieboes, Lowengrub, Wise and co-workers combined a more general, multiphase model of tumor growth with a “free-swimming”, off-lattice angiogenesis model, and observed continued hypoxia and angiogenesis (61-63, 65).
In this approach, a continuum model is endowed with spatiotemporally-varying parameters that encapsulate smaller-scale biophysics, effectively including these as constitutive relations. Early notable examples include nonlinear oxygen consumption and proliferation terms in various tumor growth models. Ward and King used nonlinear, Michaelis-Menten-type uptake and proliferation terms to model the relationship between oxygen availability, proliferation, and observations that cell growth tends to saturate due to independent limiting factors even when supplied unlimited growth substrates (75-76). Gatenby, Smallbone, and other colleagues more recently included similar such relationships between proliferation and substrate (oxygen and glucose) uptake, based upon the analysis of molecular-scale metabolism modeling (77-79). Macklin et al. recently observed Michaelis-Menten population dynamics as an emergent phenomenon from an ABM that varied the quiescent-to-proliferative phenotypic transition probability with oxygen, and validated the result by comparing against breast cancer patient data (Figure 6) (80-81).
Based in part upon earlier parameter studies (60, 73, 82-83), Macklin and co-workers formulated a phenomenological relationship between the local ECM density E and a tumor’s mechanical response to pressure gradients by varying the mobility coefficient in Darcy’s law as μ = μ0/(1 + bE), where b is constant (10, 61, 64). Thus, the tumor tissue’s mechanical response (μ) decreases as the ECM density increases, modeling 1) the increased mechanical resistance of a denser ECM, and 2) the greater number of cell-ECM integrin bonds that must be broken for cell motion (73). This functional relationship can be derived through a careful upscaling analysis (e.g., as in (10) and later expanded in (84)) of the adhesive, repulsive, and frictional forces in mechanistic, force-based agent models (e.g., (49, 56, 58, 80, 85)). Similar models have been used by Macklin et al. in nonlinear chemotaxis and haptotaxis coefficient functions (10, 64), and other groups have used a nonlocal, integral term to model cell-scale adhesion effects in continuum tumor invasion models (86-88). Indeed, such functional relationships between microenvironmental quantities and averaged cell phenotypic behavior could be incorporated in functional, spatiotemporally varying model parameters by analyzing or numerically sampling appropriate molecular-microscopic models.
This approach has been applied to multiscale chemotherapy modeling. Sinek, Frieboes, Cristini and coworkers developed a cell-scale compartmental model of chemotherapeutic drug (doxorubicin) transport within a cell’s cytoplasm and nucleus, which they combined with a mechanistic model of DNA-drug adduct formation, repair, and apoptosis (89-91). Their model accounted for cell cycling effects by varying the probability of apoptosis with the substrate level via a Hill-type function. This model, in turn, informed a spatiotemporally-varying apoptosis parameter in a continuum model of tumor growth, which, when combined with reaction-diffusion equations for substrate and drug transport, was very successful in predicting 3D, whole-tumor response to doxorubicin administered via the (neo)vasculature. As shown in Figure 7, cell monolayer drug response data was used to calibrate the model, and subsequently used to predict 3D tumor response. As in our previous examples, the functional parameters are informed by upscaling mechanistic microscale models. Indeed, the Hill-type substrate dependence (which was imposed a priori in (89-91)) can be predicted as the emergent behavior of an ABM of Macklin et al (see Figure 6) (80-81, 85).
In the fullest realization of the “multiple hierarchies” paradigm discussed earlier, discrete and continuum representations of the same phenomenon are applied dynamically and simultaneously, along with a biophysical rule to appropriately select the correct representation throughout the spatiotemporal modeling domain. This fully hybrid approach is a ongoing, cutting-edge topic in the scientific computing community; until recently, this approach has largely been computationally infeasible, and only recently have tools begun to emerge that can enforce mathematical and biophysical consistency between such simultaneous representations (10).
In recent work, Kim, Stolarska, and Othmer coupled a viscoelastic model of a tumor’s quiescent and necrotic regions with an ABM of tumor cells in the outer proliferative rim (of thickness ~100-200 μm), allowing for more detailed treatment of molecular and cellular biology within the viable rim while making use of the computational efficiency of a continuum model for the majority of the millimeter-sized tumor (Figure 8a) (92-93). The discrete cellular model, which updated earlier work by Dallon and Othmer to include apoptosis and proliferation (94), described the cells as deformable, viscoelastic ellipsoids subject to adhesive, repulsive, and drag forces, as well as growth-induced stresses. The model included a careful transfer of forces from the discrete cells onto the continuum model in the quiescent region by interpolation of the discrete cells’ forces onto the nearest nodes of a triangular mesh; mass could be transferred between the discrete and continuum models as cells changed between the quiescent and proliferative states using a least-squares projection algorithm. This work is characteristic of one main approach to hybrid modeling, where the discrete and continuum representations are applied in pre-specified regions.
More recently, Wise, Lowengrub, Cristini and colleagues combined a 3-D multiphase continuum model (where each “voxel” of tissue is modeled as a mixture of fluid, ECM, and various cell populations, and a Cahn-Hilliard type equation to regulates the mechanics of the mixture) with a lattice-free, discrete model of migrating tumor cells (10, 62-63, 65). As in previous hybrid models, the discrete and continuum representations exchange forces with one another. As an advance over earlier hybrid modeling, their approach could dynamically convert between the discrete and continuum representations while conserving mass and momentum, using criteria such as cell density and hypoxia to trigger the conversion. Mutation dynamics were modeled as transfers between the multiphase mixture cell components. The group investigated the epithelial-mesenchymal transition in hypoxic glioma (brain tumor) cells, and observed such behavior as single-file like strands of motile, palisading cells moving from the hypoxic regions towards high-oxygen environments—an effect observed in clinical histopathology (Figure 8b-c). The model was also able to treat the formation of satellite tumors in regions where discrete cells aggregated sufficiently to satisfy the continuum hypothesis (e.g., in normoxic regions) (Figure 8d).
We outlined the different ways that recent models link lower-scale modules with those at higher-scale, in the interest of developing more realistic and predictive models of cancer. We emphasize the viewpoint that mathematical and computational cancer models, whether scale-specific or multiscale, are most efficient when they are as complex as necessary, yet as simple as possible (95). The integration of an ever growing body of experimental and clinical data will demand even broader ranges of expertise and knowledge thus will require the formation of cross-disciplinary and multi-institutional groups focusing on particular fine-grained building blocks or functional modules (96). In the following, we discuss several of these daunting challenges and conclude with future directions of multiscale cancer modeling.
It is widely accepted that the more quantitative experimental data are available to build and constrain the parameter values, the more likely models are to accurately describe observed behaviors (97). Parameters of scale-specific cancer models are already difficult to quantify, not to mention those for multiscale cancer models for which parameters have to be defined at different scales. Not all of these parameters are experimentally available or even measurable with current techniques hence some have to be estimated by comparing model results to experiments. There are two classes of parameter estimation techniques that are frequently used in systems biology: local optimization (e.g., direct-search, and gradient-based methods) and global optimization (e.g., simulated annealing, branch and bound, and evolutionary algorithm). Each method has its advantages and disadvantages, and interested readers should refer to (98-100) for a comprehensive survey of parameter estimation methods.
Due to the complexity of cancer models and the quality of the data, over-fitting is a common problem in parameter estimation (101). When a model has many estimated parameters and the corresponding data are insufficient, over-fitting can lead to some unwarranted conclusions. To narrow the range of the parameter space, constraint-based modeling, which involves a multi-step procedure and accounts for both quantitative and qualitative data at the same time, may be an alternative choice (102). However, cancer modelers should be aware that, given the current underlying difficulty in estimating parameter values, many models may remain limited in the mechanistic insight they provide and in their capability to predict system dynamics in unforeseen conditions. To constrain these models and thereby improve predictive power, modelers must integrate insights from numerical parameter studies, direct experimental data from collaborators, and any information that can be gleaned from the experimental and theoretical biology literature (which itself may require reinterpretation with a mathematical point of view) (81). The hybrid modeling approach may also be useful for integrating patient data across a variety of spatial scales, where calibration to each datum occurs at the appropriate scale, and the information is subsequently allowed to propagate throughout the multiscale framework (11, 81).
Significant efforts have been made to seek mathematical approaches and computational resources that will enable cancer models to be run in a reasonable period of time (103). It is generally accepted that if a model works on higher spatial and temporal resolution, it will need higher compute power and thus longer run time to accomplish computing tasks. Discrete or discrete-based hybrid models are more seriously affected by the compute intensity constraints because they are generally too detailed to simulate over a long period of time, particularly in large, 3D domains. One solution is to take advantage of very powerful, massively parallel supercomputers and develop numerical algorithms that can run on these machines. However, doing so may not resolve all the difficulties in handling the enormous amount of experimental and clinical data; it is also not practical in many current clinical settings thus would have limited usage.
As discussed, adaptive hybrid modeling has the potential to save computational time while maintaining the predictive power of the model (see Section 2.2). Other methods for solving the compute intensity issue are also available, but are all at an experimental stage. For instance, the equation-free approach developed by Kevrekidis and co-workers (104-106) leverages the spatiotemporal scale separation to allow for significant gains in computational efficiency by alternating short bursts of appropriately initialized microscopic simulations with accelerated result processing at the macroscopic, continuum scale. Furthermore, methods from other modeling communities, such as the Heterogeneous Multiscale Method (HMM, e.g., (107-108)), can provide useful insight into efficient numerical methods that may be incorporated into the development of multiscale cancer models as well. By drawing on the strengths of these potentially very useful methods and integrating them into the next generation of multiscale cancer modeling, we may be able to produce more comprehensive and computationally efficient models to simulate tumor progression and predict treatment impact.
Data sharing is an active topic in the field of systems biology (95). We also advocate that, in order to avoid duplication of effort and thus waste of resources, achievements of publicly funded work including data and models should be made accessible to the community in a form where others can review and reproduce the modeling results, and then reuse and revise these resources in future works. In the following, we focus on two aspects of data sharing and model reuse.
First, the data and models should be presented in standardized formats with clearly stated dependencies and problems. This will improve the reusability of models and clarify their limitations. In systems biology at large, there have been a few standards established, such as SBML for biochemical pathways, CellML for biophysical models, and FieldML for spatial fields (109). These standards define how models are structured, how the mathematical equations are encoded, and how units are defined. Establishing these standards has been a significant step towards achieving a robust foundation for the modeling community; such standardization efforts, including amending existing or creating new ontologies, can facilitate translation of multiscale models into and eventual adoption of these in silico tools in clinics.
Second, to store and exchange modeling tools, the cancer modeling community needs digital repositories of models provided by researchers from academic, governmental and corporate laboratories. Web services are needed to connect and integrate these currently often incompatible models and tools for data acquisition, storage, analysis and presentation in order to bridge the integration gap (18). In brief, cancer modelers need to efficiently share data, information and knowledge across geographic and organizational boundaries within the context of distributed, multi-disciplinary and multi-organizational collaborative teams, while at the same time effectively protecting data ownership and data security.
Many web-based systems biology communities or databases have been created over the past years, and some of them are developed based on semantic technology (see (110) for a review). Among them, the Center for the Development of a Virtual Tumor (CViT) (111) is building an ever-growing community of researchers around the world dedicated to cancer modeling. It currently provides cancer researchers a community-driven web platform with functions including wikis, blogs, forums, member profiles, and RSS-based news updates (112). It also develops NIH/NCI-caBIG® compliant infrastructure tools to facilitate interaction among its contributing scientists. Members of CViT come together online on a regular basis to discuss cutting-edge literature on CViT’s online forums which separates CViT from other model repositories such as BioModels (113). Additionally, the social networking aspect of the site allows research teams to collaborate from anywhere around the world in a workflow designed specifically for the cancer modeling community. Recently, CViT released the Digital Model Repository (DMR), an innovative web platform for the exchange of cancer models (114). Cancer investigators can upload their models, experimental data, and simulation results and publish them to other DMR users of their choice who will then be able to access those files. The DMR also implements an innovative elicensing workflow and is built using semantic web technologies through the Resource Description Framework (RDF; a standard metadata model for describing the relationship between two objects (115)). The newest feature of the DMR is the Computational Model Execution Framework (CMEF). Its aim is to allow executable files to be uploaded as part of a model. The owner of the model can indicate variable parameters which can be defined for each run, and specify details of the model, including programming language and runtime computing environment. Others can either download the model and run it on their own computers, or simply execute the model using computing resources provided by the DMR. Altogether, CViT and its semantic services help to advance the cancer modeling field in facilitating web-based multidisciplinary cancer research.
Approaches provided by systems biology are beginning to impact medicine (116). In treating cancers, the integration of computational and experimental techniques may improve our abilities to design more efficient cancer therapies. Corresponding applications include biomarker validation, development of more accurate diagnostic tests, and targeted drug discovery - all geared towards the optimization of individualized cancer therapy.
The conventional population-based approach for therapeutic developments in clinics still relies on a series of randomized clinical trials aimed at searching for favorable yet averaged treatment outcome (117). However, patient responses to a particular drug or therapy are known to fall into a more or less wide range that deviates from this averaged behavior. Multiscale cancer modeling may eventually help to explain not only why some therapies fail while others prove to be effective in controlling tumor progression, but also why a particular treatment works only in a fraction of patients. Models will have to incorporate patient-specific experimental and clinical data at all biological levels, including, e.g., genomic, proteomic, anatomical, physiological, and pathological data. Training the model on a patient’s data will yield a more accurate description of the specific kinetics of disease progression. Hence, this approach should provide a higher predictive power than that achieved with pooled data only and thus should be able to guide personalized treatment strategies more accurately to improve outcome (117). Moreover, multiscale cancer models can be used to help characterize a patient’s specific biomarkers (ranging from critical pathway proteins to phenotypic profiles and imaging patterns) and then compare these ‘cross-scale signatures’ with pooled data from conventional clinical practice or from a trial with the candidate drug being applied, in an effort to simulate the specific response of the patient. Introducing multiscale cancer modeling to medicine has the potential to facilitate the breakthrough of personalized medicine, and eventually maximize advances in science and technology for the benefit of cancer patients by helping select or optimize preventative and therapeutic patient care.
In summary, multiscale cancer modeling is a most promising, innovative research area that constitutes a critical driver for the field of integrative cancer systems biology. Challenges to the success of this approach arise as a result of our still limited understanding of the complex, dynamic nature of cancers, the often constrained access to appropriate experimental and clinical data, the difficulties in validating models against these data, and the challenges involved in communicating and sharing modeling methods amongst the field’s multiple stakeholders. However, by drawing on the collaborative effort and expertise of scientists from different disciplines and the continuing development of advanced, innovative computational and mathematical methods, we believe that multiscale cancer modeling will reach its full potential in guiding targeted experimental research, in enabling patient-specific predictions and thus in accelerating personalized medicine.
Discrete modeling can explicitly represent individual cells in space and time, and track and update their internal states according to a pre-defined set of biological and biophysical rules. This approach is particularly useful for studying carcinogenesis, genetic instability, natural selection and cell-cell and cell-matrix interaction mechanisms. The dynamics of discrete cancer cells can be investigated with lattice-based or lattice-free methods, where the former builds a grid with a finite number of dimensions in which cells live, while the latter describes cell actions in arbitrary locations and their interactions in arbitrary directions. Because discrete modeling is based on a series of rules applied to each cell, it is possible to translate detailed biological findings into rules for the model (24). However, the computational demand increases rapidly with the number of cells modeled, limiting these models in the spatial and temporal scales they can represent.
Continuum modeling describes tumor tissue as a continuum medium rather than working at the resolution of individual cells, and thus is a good choice for modeling larger-scale systems. This approach draws upon principles from continuum mechanics to describe model variables as continuous fields mostly by means of partial differential or integro-differential equations. Common continuum model variables (e.g., cell volume fractions, density, and cell substrate concentrations, e.g. nutrient, oxygen, and growth factors) are somewhat easier to obtain, analyze, and control, compared with those in the discrete case (25). Although these models can characterize global properties of gross tumor growth and invasion at the tissue or higher-scales, they cannot be used to examine individual cell dynamics and discrete events, such as EMT, since in a process that small, changes to a cell or a set of cells can move a nonlinear cancer system to a different state (26). This may be important when studying the effects of genetic, cellular, and microenvironment characteristics on overall tumor behavior.
Hybrid modeling attempts to integrate and draw on the strengths of both continuum and discrete descriptions. Different definitions on hybrid modeling exist in the filed (9-10, 23), but in this article, hybrid models are roughly divided into two categories: composite and adaptive hybrid modeling. In composite hybrid models, individual cells are treated discretely but interact with other chemical and mechanical continuum fields. Under appropriate formalisms, such models are able to couple different scales impacted by the growth process with biophysical, biochemical, and biomechanical information passed between scales. In adaptive hybrid models, both discrete and continuum representations of cells are chosen dynamically and adaptively where appropriate, e.g., discrete modeling for EMT and continuum modeling for the tumor bulk. Thus, the adaptive hybrid modeling can achieve discretely high resolution wherever and whenever necessary, while at the same time reducing the compute intensity as much as possible to support scalability of the approach to clinically relevant levels. Together, although continuum and discrete approaches have each provided important insights into cancer-related processes occurring at particular spatial and temporal scales, the complexity of cancer and the interactions among the cells and their complex tumor environment call for a multiscale continuum-discrete (hybrid) approach (10). Hybrid models have the potential to couple biological phenomena from the molecular and cellular scales to the tumor scale.
This work has been supported in part by the National Institutes of Health (NIH) grant CA 113004 and by the Harvard-MIT (HST) Athinoula A. Martinos Center for Biomedical Imaging and the Department of Radiology at Massachusetts General Hospital. Vittorio Cristini acknowledges the NIH for support through the Physical Sciences in Oncology Center grants 1U54CA143837 and 1U54CA143907. He also acknowledges the NIH Integrative Cancer Biology Program for support under grant 1U54CA149196, and the National Science Foundation (NSF) for support under grant DMS-0818104.
DISCLOSURE STATEMENT The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.