|Home | About | Journals | Submit | Contact Us | Français|
Multiscale modeling is increasingly being recognized as a promising research area in computational cancer systems biology. Here, exemplified by two pioneering studies, we attempt to explain why and how such a multiscale approach paired with an innovative cross-scale analytical technique can be useful in identifying high-value molecular therapeutic targets. This novel, integrated approach has the potential to offer a more effective in silico framework for target discovery and represents an important technical step towards systems medicine.
During the past decade, technical developments in cell and tissue analysis and in molecular diagnostics have been significant, but progress in drug development and clinical implementation has not been able to satisfy market needs [Turk 2006]. Less than 10% of new compounds that enter clinical trials ultimately reach the market with many more fail in the preclinical stages of development [Lalonde et al. 2007]. Given the low productivity and escalating costs of drug development, it is critical for the pharmaceutical industry to explore more creative, efficient, and cost-effective methods to develop new therapeutic compounds [Couzin 2002]. In cancer drug development, the discovery of small molecules targeting specific oncogenic pathways has revolutionized anti-cancer therapy. The term ‘targeted therapy’ refers to drugs that specifically act on well-defined pathway components or proteins to impair abnormal cell proliferation, inhibit tumor blood vessel development, or promote the specific death of cancer cells through a process known as apoptosis. However, targeted therapy often fails to demonstrate survival benefit in cancer patient populations due to the acquisition of drug resistance [Petricoin et al. 2005] which, to some extent, reflects the fact that cancer is such a complex disease that only focusing on particular genetic and molecular properties of the cancer cells may not be sufficient.
Cancer systems biology is a newly emerging, multi-disciplinary field in current oncology research that focuses on the systemic understanding of cancer initiation and progression, by investigating how individual subcellular components interact to generate the whole tumor system [Deisboeck et al. 2007]. Other than being based only on conventional wet lab experiments, a systems approach often also involves computational modeling to simulate experiments and to predict and optimize therapies. In drug development, an important goal is the selection of a dose and a dosing regimen that achieves the targeted clinical effect while minimizing undesirable adverse effects [Nelander et al. 2008]. Computational modeling can be very helpful in this area, providing a means to study the dose- and time-response behaviors of safety and efficacy conditions, in a controlled setting and at a fraction of the cost of conventional experimental studies. Thus far, computational modeling has been primarily used to assess the effects of drugs on the activities of large-scale signaling networks [Araujo et al. 2007]. Such a focus on the molecular level is useful for identifying which parts of a signaling network give rise to robustness and which parts are more fragile, and this analysis may lead to new ways to rationally design multi-target or multi-component therapies. However, cancers are heterogeneous cellular entities whose growth is dependent upon dynamical interactions between genetically altered cells, and their surrounding microenvironment changes dynamically [Bissell and Radisky 2001]. Even extrinsic environmental conditions alone, independent of genetic mutations, can induce the carcinogenic transformation of cells [Postovit et al. 2006]. Hence, computational cancer models should take into account both signaling proteins or networks and different biochemical and biophysical contexts, i.e., environmental conditions, when identifying potential therapeutic drug targets. The technique used to develop this type of multiscale cancer models integrates tumor properties across multiple spatial and temporal scales and then investigates the emergence of larger scale phenomena from the interaction of units at smaller scales [Tracqui 2009]. Further discussion of multiscale cancer models developed in other cancer research related fields is beyond the scope of this review and interested readers should refer to Sanga et al. ; Wang and Deisboeck  and Zhang et al. .
The identification and validation of for the overall outcome critical sites of the cancer system is an essential first step in therapeutic target discovery and thus drug development [Sioud 2007]. To move a step further, i.e., to identify cross-scale context-dependent therapeutic targets, it is necessary to assess the influence of molecular-level parameters on tumor growth behaviors at a higher level. There are two preconditions to meet in order to perform such a task: (1) developing a reliable multiscale cancer model and (2) finding a suitable cross-scale target evaluation method. In this review, by taking the most recent developments contributed by our laboratory as examples, we introduce how a multiscale cancer model can be applied to the discovery of therapeutic targets. We also discuss some of the current challenges of this approach and future directions for continuing research.
Because its aim is to identify potential molecular targets that once hit lead to control if not inhibition of cancer growth, the in silico model put forward should involve molecular events that give rise to higher-level behaviors. Only a small number of multiscale models incorporate an underlying molecular module to investigate molecular dynamics and their effects on cancer growth and invasion (see [Alarcon et al. 2006; Anderson et al. 2006; Gerlee and Anderson 2009; Ramis-Conde et al. 2009; Ramis-Conde et al. 2008; Smallbone et al. 2008] for representative examples). None of these models however attempts to provide insights into drug target discovery. We have been exploring the possibility of using a cross-scale approach to determine critical molecular parameters that have a significant impact on tumor phenotype or ‘outcome’ at the microscopic level [Wang et al. 2009; Wang et al. Submitted; Wang et al. 2008b]. These critical parameters have the potential to being developed into therapeutic targets, and hence warrant further experimental examination. In the following, we will first introduce the development of the two essential parts, i.e. multiscale cancer models and related cross-scale analysis methods, respectively. We then present some preliminary results that we have gained in identifying cross-scale molecular targets.
We use agent-based modeling (ABM) to simulate molecular signaling dynamics within a cell and the interactions of cancer cells with each other and with their local biochemical environment in order to predict higher-level emergent tumor phenomena [Wang et al. 2009; Wang et al. 2007]. In our models, each cell is represented by a software agent, and this direct one-to-one mapping between real and virtual cells in terms of parameter acquisition from experiments and in silico model validation is an advantage of the model design.
At the molecular level, we introduce signaling pathways into the multiscale framework. The pathways are induced by growth factors, mediated by growth factor receptors, and terminated with extracellular signal-regulated kinase (ERK) signaling. In Wang et al.  we implemented an epidermal growth factor receptor (EGFR) pathway, and in Wang et al.  we extended it to a combined pathway with transforming growth factor β (TGFβ) signaling. Both pathways are designed specific to Non-Small Cell Lung Cancer (NSCLC). Pathway dynamics are regulated by material balance and kinetic equations, as well as by reaction rates that are dependent on the changes of species concentrations over time, and the pathway models are implemented using ordinary differential equations (ODEs). At the microscopic level, a biochemical microenvironment is constructed to represent a virtual lung tissue in a two dimensional (2D) [Wang et al. 2007] or three dimensional (3D) [Wang et al. 2009] space. Heterogeneous environments are attained by distributing external diffusive chemical cues (such as growth factors, glucose, and oxygen tension) throughout the microenvironmental lattice. In our specific model setup, a nutrient source maintains the highest weight for each of the aforementioned cues, and hence it is the most attractive location for the chemotactically acting tumor cells. Throughout the simulation, the concentrations of the chemical cues are continuously updated at a fixed rate. These diffusion processes are implemented based on partial differential equations (PDEs). Each cell carries (a) self-maintained signaling pathway(s), meaning that cells will differ in signaling profiles and thus will exhibit different phenotypic behaviors as the simulation progresses.
Here, our approach to establish the link between molecular and microscopic scales uses an experimentally supported molecularly-driven cellular phenotypic ‘decision’ algorithm (Figure 1). In particular, PLCγ and ERK (two downstream signaling molecules of EGFR) are used to determine the emergence of two important phenotypic traits: migration and proliferation. Experimental studies have shown that the transient acceleration of accumulating PLCγ levels leads to cell migration [Dittmar et al. 2002], while that of ERK leads to cell replication [Santos et al. 2007]. Therefore, in our model, the rate of change of PLCγ determines the cellular migration decision, and the rate of change of ERK dictates the cellular proliferation fate. If a cell decides to migrate or proliferate, it will search for an appropriate neighborhood to move to or for its offspring to occupy. It should be emphasized that tumor growth and invasion patterns due to cell proliferation and migration are neither pre-defined nor intuitive: rather, they emerge as a result of the intracellular signaling of individual cells and the dynamic cellular interactions within the framework of the biochemical microenvironment.
Sensitivity analysis has been adopted as a useful tool by the pharmaceutical industry to complement traditional hypothesis-based approaches for identifying effective drug targets [McDermott and Settleman 2009]. It is a method that systematically identifies specific perturbations in system inputs that have significant effects on system responses, and this approach is especially helpful when it is not possible or practical to conduct numerous experiments on the living system itself [van Riel 2006]. When identifying drug targets using cell signaling models, typical biological measurements of system output include signaling duration, amplification, and activation patterns [Materi and Wishart 2007], i.e., all at the molecular level; however, since the system investigated, i.e. cancer, is multiscale in nature, system output assessment in silico has to capture tumor growth behavior at a higher level. For example, in our approach, the cross-scale analysis evaluates how changes occurring at the molecular level percolate across the scales of the model and affect tumor growth behaviors at the microscopic level, which can be any microscopic tumor property, e.g., cell density, tumor volume, morphology, or motility speed.
Sensitivity analysis methods fall into two categories: local sensitivity analysis (LSA) and global sensitivity analysis (GSA). LSA (one-at-a-time parameter variation) evaluates changes in system outputs with respect to each parameter at a particular point, resulting in a very narrow coverage of parameter space. GSA, on the other hand, addresses model behavior over a wide range of parameter operating conditions [Frey and Patil 2002; Helton et al. 2006]. Local and global sensitivity analysis methods may be separately appropriate for different systems depending on the purpose of their implementation. While we have recently started to apply GSA, local sensitivity analysis is a reasonable method to begin with because its implementation is relatively simpler than that of the global. In a typical local sensitivity analysis, a sensitivity coefficient is used as an index to evaluate the effects of perturbations of individual parameter values on the overall system outcome. The coefficient is defined as follows [Rabitz et al. 1983]:
where p represents the parameter that is varied in a simulation and M, the response of the system; M0 is obtained by setting all parameters to their reference values, and thus (Mi−M0) is the change in M due to the change in p, i.e., (pi−p0). In our case, M now corresponds to a microscopic tumor response, while p corresponds to any of the individual, molecular input parameters (e.g., pathway component concentrations, kinetic reaction rates). The bigger the absolute value of the sensitivity coefficient, the more sensitive the given parameter.
By utilizing the multiscale cancer modeling platforms and the cross-scale analysis technique introduced above, we have been able to contribute to the discovery of molecular drug targets with (in silico) proven cross-scale impact. We present here two pioneering studies, which should give the reader a clearer understanding of how to integrate the cross-scale discovery process into his/her experimental and/or theoretical works.
The first example [Wang et al. 2008b] applies this analysis technique to a 2D NSCLC model [Wang et al. 2007]. Figure 2 schematically illustrates how the evaluation process identifies a cross-scale therapeutic target. The signaling pathway consists of 20 components downstream of EGF stimulation (see Wang et al. 2007 for detail). Input parameters are the initial concentrations of the seven implemented pathway components. Only one parameter is varied at a time when running a simulation, and all other parameters are held fixed at their reference values (i.e., LSA). The model output is the spatiotemporal expansion rate of the tumor, a microscopic level behavior. In the 2D model, a simulation run is terminated when the first cell reaches the nutrient source, and thus cells travel the same distance to the source under all simulations (each corresponding to an individual parameter variation). Hence, we use the number of simulation steps as a measure for the tumor expansion rate, with a greater number of time steps until termination indicating a slower tumor expansion rate in a given simulation. The analysis finds that three signaling molecules, PKC, MEK, and ERK, are critical to the stability of the model in influencing the tumor expansion rate, the result of which is in good agreement with experimental and pharmaceutical studies [Allen et al. 2003; Green et al. 2006; Han et al. 2005]. Moreover, each of these three critical parameters has a critical interval, and if a variation of concentration occurred in that interval, then varying these parameter values had a strong influence on the tumor expansion rate. Another interesting finding is that the sensitive parameters exhibit a similar changing pattern in that a relatively larger increase or decrease in their reference values (obtained from the literature) result in a lesser influence on the multicellular performance of the system.
It is a common phenomenon that a parameter that is significant in determining one specific model output may not be significant for others. For example, in a mitogen activated protein kinase (MAPK) pathway study, MEK dephosphorylation had a significant impact on the duration and integrated output, but not the amplitude, of ERK activation [Hornberg et al. 2005]. Hence, an evaluation function that concurrently accounts for two or more model outputs can be highly useful. In the second example [Wang et al. Submitted], we therefore propose a therapeutic index (TI) function and apply it to a 3D NSCLC model [Wang et al. 2009] where, at the molecular level, both EGF and TGFβ can trigger downstream ERK signaling through different routes. The purpose of this study is to identify the molecular parameters that are most effective in suppressing overall tumor growth, by evaluating two tumor growth indices at the microscopic scale: tumor volume and expansion rate. Indeed, the TI function attempts to find a high-value target that reduces the ability of cancer cells to grow (and/or causes them to die) and at the same time, diminishes cancer cell motility (i.e., reduces invasion and contains metastasis) as much as possible, in order to minimize overall tumor expansion. The TI function is defined as follows [Wang et al. Submitted]:
where x and y represent the two tumor outcomes: simulation steps (as a measure for tumor expansion rate) and the final number of alive cells (as a measure for tumor volume), respectively; xi and yi indicate their corresponding values for the ith simulation run; x0 and y0 are their corresponding values for the simulation when all model parameters are set to their reference values. Parameters that produce high TI values merit particular attention since, generally, they will produce favorable anti-tumor outcomes and hence will be valuable therapeutic targets. There are exceptions to this case (e.g., (xi>x0 and yi>y0) or (xi<x0 and yi<y0) may lead to a large TI value, but do not produce a desirable therapeutic outcome) and interested readers should refer to the original publication for a more detailed description of the model and results. Because targeted combinatory therapies may well play an increasingly vital role for cancer treatment in the future, we seek to gain insights into the influence of both individual and combined parameter perturbations on the two tumor outcomes. We simulated two types of treatment strategies with relevance for the pharmaceutical and biotechnology industries: parameter inhibition and amplification. Analysis finds that MEK is the top therapeutic target under the parameter inhibition strategy, while ERK emereges as the most valuable target in an amplification regimen. However, most combined parameter perturbations were not an improvement over individual parameter inhibition or amplification treatments, and this discouraging finding is in agreement with the sobering outcomes of current clinical trials of multi-targeted therapy in NSCLC treatment [Felip et al. 2007].
The emergence of systems biology is now beginning to influence health care, through its early applications to drug discovery and its vision for personalized medicine [Popel and Hunter 2009]. In discovering therapeutic targets for cancer treatment, it becomes imperative that more efforts be directed towards a quantitative approach to the analysis of therapeutic success and failure. It is our understanding that our cross-scale analysis technique can help to identify true, effective therapeutic targets, because it accounts for both molecular and microenvironment factors of cancer, as well as for interactions among these factors. A target identified using a signaling pathway in silico model alone may have little relevance for what happens in vivo and later in clinics across the scales of interest. Such a potentially ‘false’-positive target may prove ineffective in an integrated, multiscale model early on. Indeed, it is increasingly agreed that considering the full biological context of therapeutic targets and moving beyond individual genes and proteins can increase the productivity of drug discovery [Cho et al. 2006]. Multiscale cancer modeling paired with a cross-scale target evaluation technique provides such a new, effective framework for therapeutic target discovery. Moreover, our multiscale models may also facilitate the tracking of molecular components as the entire tumor system evolves (see [Wang et al. 2009] for an example). Hence, this integrated approach, once the input data are personalized, can also be used to identify potential molecular biomarkers, which may help to (1) determine which patients are likely to respond to a particular agent or regimen, and (2) indicate disease status and treatment progress in a particular patient. The combination of targeted drugs is considered to be a promising strategy to improve treatment efficacy, reduce off-target effects and/or prevent evolution of drug resistance [Borisy et al. 2003; Chou 2006; Keith et al. 2005]. It is intuitive to assume that the simultaneous inhibition of a cascade of some key regulating proteins could enhance the efficacy of the therapy, compared to inhibiting a single protein in isolation. In fact, clinical trials are currently in progress, to evaluate the efficacy of various drug combinations, including monoclonal antibodies and tyrosine kinase inhibitors [Guarino et al. 2009], anti-receptor therapy with EGFR downstream signaling inhibitors [Milton et al. 2007], and angiogenesis inhibitors [Morabito et al. 2009]. However, these clinical trails indicate that there is little or only modest benefit from using combination therapy. From our systemic cross-scale analysis (the second example shown above), most of the combined parameter perturbations (inhibitions, amplifications, or a combination of both) indeed did not show improvement when compared to individual parameter inhibition or amplification treatments. More work will have to be done on this strategy in the future and, as demonstrated here, multiscale models can add value by generating reasonable, testable hypotheses as to which targets to hit when and in what combination.
While most promising, the field however faces a number of formidable technical challenges. These include, for the multiscale cancer modeling part, dealing with data heterogeneity (as, for the foreseeable future, model parameters are likely obtained from different sources, disparate experimental settings, and a variety of cell lines), performing model parameter verification and model validation, and addressing the computational intensity issue when a model becomes more complicated (as more pathways and a richer microenvironmental milieu will need to be incorporated) [Hunter and Borg 2003; Walker and Southgate 2009]. As for the latter, we further note here that current models entirely or partly (i.e., hybrid models) developed with discrete modeling techniques can only handle a relatively small amount of cell populations, compared to pure continuum models, because they are too detailed to simulate over a long period of time, particularly in large, 3D domains. Hence, a more innovative way of thinking about modeling is required to solve this problem. New modeling methods, such as multiscale, multi-resolution modeling [Wang and Deisboeck 2008], iterative down-scaling and up-scaling processes [Kevrekidis and Samaey 2009], and heterogeneous multiscaling [Ren and E 2005] are exciting technical frontiers that require further examination. For the cross-scale analysis part, challenges include dealing with model parameter uncertainty, handling multiple outputs at the same time (we only have provided a pilot study in this review [Wang et al. Submitted]), and taking into account multiple sensitivity analysis methods. The last point is crucial because each individual sensitivity analysis method has its advantages and disadvantages. In other words, there exists no single ultimate solution that best fits all types of applications [Frey and Patil 2002; Helton et al. 2006].
Despite these challenges, further advances in the technique of multiscale cancer modeling and related cross-scale analytic methods have the potential to add much needed insights as to why some drug treatments fail while others prove to be effective in controlling tumor expansion. In the future, one can imagine developing functional modules at different scales, specified to a particular cancer type and tailored to individual patients, thus using multiscale in silico modeling to advance the field of personalized predictive medicine.
This work has been supported in part by 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.