Models of Cancer Drug Discovery and Response to Therapy

of drug response through the application of statistical modeling and machine learning algorithms (e


Introduction
The discovery of genetic abnormalities together with the fast pace of drug development has led to an expanding number of mutation-tailored treatments for human diseases. This paradigm, in particular, has revolutionized therapeutic strategies for an increasing number of cancers. Many of such molecularly targeted therapies are chosen or designed to block the action of driver oncogenes or some components of their effector pathways. The underlying rationale is that tumor cells sustain a state of pathway dependency, or so-called oncogene addiction, in which cellular growth and survival become highly dependent on the continued activity of the oncogenic pathway (Weinstein and Joe, 2006). This phenomenon has been widely exploited in the treatment of many human cancers. Among the numerous examples are BCR-ABL targeted therapies for BCR-ABL driven leukemia (O'Brien et al., 2003), BRAF and MEK kinase inhibitors for BRAF-mutant melanomas (Chapman et al., 2011), and epidermal growth factor receptor (EGFR) inhibitors for EGFR-mutant lung adenocarcinomas (Mok et al., 2009).
Genomic and molecular profiling of human tumors have enabled the discovery of driver oncogenes and identification of patient populations who would potentially benefit from drugs tailored to interact with oncogenic pathways (Hamburg and Collins, 2010;The Cancer Genome Atlas Research Network et al., 2013). Resistance to therapy due to tumor heterogeneity, however, has been a major challenge that has prevented us from fully realizing the impact of these discoveries (Holohan et al., 2013). Many of the molecularly targeted therapies, even when initially effective in a population of cancer patients, eventually fail to cure disease. A prototypical example of this challenge is observed in the case of BRAF/MEK-targeted therapy in melanoma (Lito et al., 2013;Rambow et al., 2018). Small-molecule inhibitors of BRAF and MEK kinases or their combination benefit a majority of BRAF-mutant melanoma patients. However, they commonly fail to cure disease due to acquired resistance, which evidently arises following an initially incomplete response and adaptation during the early phases of treatment. Therefore, predicting which tumor populations would respond to specific targeted therapies, and preventing or overcoming therapeutic adaptation are likely to be key to durable cancer therapies (Luskin et al., 2018).
Numerous studies in the past decade have been focused on studying the origins of heterogeneity in drug response in tumor cells. Together, they have demonstrated a wide spectrum of phenotypic outcomes following inhibition of oncogenic signaling across tumors of different genetic background, tissue type and differentiation state (Hyman et al., 2015;Sosman et al., 2012;Tsoi et al., 2018). At the cell population level, the spectrum spans a range between "no response" to "complete response," as is evident from short-term drug maximal effect (E max ) measurements in tumor cell lines (Barretina et al., 2012;Basu et al., 2013;Fallahi-Sichani et al., 2013) as well as long-term therapeutic responses in patient-derived xenografts (PDX) and human tumors (Gao et al., 2015;Sosman et al., 2012). At the single-cell level, inhibition of oncogenic signaling leads to a variable mixture of phenotypic outcomes, including (transient or durable) proliferative arrest, apoptosis, changes in differentiation state or cellular senescence (Felsher, 2008). This represents a remarkable diversity in the state of pathway dependency among tumor cells that carry the same mutated oncogene. Such diversity has been a major challenge for the development and use of targeted therapies, arguably the cornerstone of precision cancer medicine. Therefore, among the key approaches to improve the benefits of oncogene-targeted therapies have been efforts to: (i) build a precise understanding of the variable state of oncogene or pathway dependency in tumor cells, (ii) acquire the ability to rationally modulate these states, and (iii) extend molecularly targeted therapies to include actionable nononcogene dependencies (Califano and Alvarez, 2017;Luo et al., 2009). By achieving these goals, we will be able to: (1) accurately predict effective targeted therapies, (2) identify which tumors are likely to respond to such therapies, and (3) predict what combinations of drugs would improve therapeutic outcomes in nonresponsive patients.
Recent advancements in high-throughput molecular profiling combined with systems biology approaches have significantly enhanced our understanding of oncogenesis and the discovery of novel targeted therapies. Furthermore, modeling drug responses in heterogeneous tumor cell populations has developed into an attractive approach toward identifying predictors of sensitivity to treatment, mechanisms of drug resistance, and optimal combinations or dosing schedules to maximize therapeutic response (Michor and Beal, 2015). In this article, we discuss different types of heterogeneity in the state of oncogene and nononcogene dependencies and targeted drug response in tumor cells. We review experimental and modeling approaches used to dissect the dynamic, phenotypic and molecular processes that are associated with such heterogeneities. We will then discuss recent systems biology approaches used to measure, model, predict and connect different types of heterogeneity in tumor cell responses to targeted therapies. A multilevel and integrative understanding of such heterogeneities will be essential for both advancing cancer biology and realizing the promise of precision medicine via promoting the rational development of new drugs or combinations that maximize therapeutic benefit in cancer patients.

Heterogeneity in the State of Pathway Dependency
Numerous lines of evidence now provide support for the contribution of several distinct but not mutually exclusive factors to the observed heterogeneity in the state of pathway dependency in tumor cells. The most conventional are genetic mechanisms. Tumor cells frequently acquire multiple mutations in distinct oncogenes and tumor suppressors. Acquisition of additional mutations may reduce dependency on a single oncogene (Choi et al., 2014). On the other hand, heterogeneity in the state of pathway dependency exists even among genetically homogeneous tumor cells derived from a single clone (Fallahi-Sichani et al., 2017;Shaffer et al., 2017). Recent time-lapse live-cell imaging experiments, followed by single-cell tracking, have shown that genetically homogenous populations of cells respond heterogeneously to inhibitors of oncogenic signaling. While some cells die within the first 1-2 days of drug exposure, some arrest in the G0/G1 phase of the cell cycle, and some progress toward a slowly-dividing phenotype, which may remain stable for weeks (Fallahi-Sichani et al., 2017). Feedback-mediated reactivation of the target pathway or induction of parallel signaling cascades have been recognized as major causes of rapidly-induced adaptive resistance to oncogene-or pathway-targeted therapies (Duncan et al., 2012;Lito et al., 2012). Nevertheless, a number of recent studies have also identified transiently heritable, tissue-specific, epigenetic states linked to transcriptional programs with reduced dependence on the oncogenic activity (Fallahi-Sichani et al., 2017;Roesch et al., 2013;Sharma et al., 2010;Toska et al., 2017). Such epigenetic states dynamically evolve from one state to the next and may ultimately transform to an irreversibly oncogene-independent state (Hata et al., 2016).
Taken together, the observed heterogeneity in the state of pathway dependency may be attributed to an interplay between: (i) nongenetic, feedback-regulated or stochastic changes in the expression and activities of signaling proteins, including the oncogene/oncoprotein itself, or its upstream, downstream and parallel effectors (Cohen et al., 2008;Fallahi-Sichani et al., 2013;Konieczkowski et al., 2018), (ii) epigenetically-regulated (transient or heritable) alterations in gene expression states (Fiziev et al., 2017;Roesch et al., 2010;Sharma et al., 2010), and (iii) the mutational background of the cells (Garnett et al., 2012). All three groups of factors are involved. However, none of them on its own appears to be predictive of how an individual cell will respond to inactivation of an oncogene or suppression of an oncogenic pathway. The reason is that genetic, epigenetic, transcriptional and biochemical information in a cell is processed not in isolation but as part of an intertwined network that determines the functional state of a cell at any time. This network adapts dynamically in response to external (e.g. drug) and internal (e.g. genetic) perturbations. The phenotypic consequence of adaptation is regulated both by deterministic processes (such as genetic, epigenetic and signaling mechanisms) that allow network rewiring, and by stochastic events that increase the variance between otherwise identical cells (Voss and Hager, 2014). A precise understanding of how these processes interact in tumor cells will likely be key to predicting the state of cell dependency on an oncogene or its downstream effectors.

Pharmacogenomic Models of Target and Biomarker Discovery
The extensive amount of heterogeneity among and within tumors has made the prediction of anticancer therapies a challenging task. In the past decade, multiple large-scale studies have taken advantage of a variety of cancer models to establish systematic platforms, identifying predictors of drug sensitivity and resistance in tumor cells. By linking measurements of drug response metrics to the functional complexity of cancer genomes through computational modeling, such pharmacogenomic platforms have provided a powerful guide to the discovery of new and rational cancer therapeutic strategies (Barretina et al., 2012;Basu et al., 2013;Garnett et al., 2012;Haverty et al., 2016).
The most classical and commonly used experimental models in these studies are cancer cell lines, developed through decades of effort by the cancer research community. Despite their historical value, however, conventional cell lines are often considered too simple to reflect the full spectrum of the complexity and heterogeneity observed in human tumors. Therefore, associating measurements of drug sensitivity in cancer cell lines to in vivo responses of human tumors has been a substantial challenge to overcome in pharmacogenomic studies (Grandori and Kemp, 2018). Patient-derived xenograft (PDX) tumors capture both genomic diversity and complexity of the tumor microenvironment. However, practical limitations and the high cost of animal studies, as well as the often long time span required for PDX establishment, have limited the use of these models for highthroughput discoveries to a few studies (Gao et al., 2015). Organoid cultures of patient-derived tumor cells, by contrast, have attracted significant attention in recent pharmacogenomic studies (Grandori and Kemp, 2018;Lee et al., 2018;Pauli et al., 2017). These cancer models could be used in large-scale studies, while preserving the complexity, genotype and phenotype of the tumor from which they have been derived. Furthermore, they could be used in studies with a clinically relevant timeframe, generating results that could be compared with the donor patients' clinical history.
Overall, using tumor models of various type and complexity, high-throughput sequencing and transcriptional profiling together with large-scale genetic and drug screening have uncovered important cancer gene-drug associations. Such pharmacogenomic associations have both led to the discovery of new cancer therapeutics and identified biomarkers that may predict therapeutic response in cancer patients. Recent systems biology efforts have played a key role in developing and improving the predictivity of such pharmacogenomic models. These efforts include: (i) improving quantitative analysis of drug response with dose and time, in order to derive reliable metrics of drug sensitivity, (ii) developing statistical learning and network models that link metrics of drug sensitivity to the activity of individual genes or signaling networks in tumor cells, and (iii) integrating data at different biological levels (e.g. tumor genome, epigenome, transcriptome, proteome and phosphoproteome) and complexity (e.g. tumor microenvironment) in order to improve therapeutic response predictions. We will discuss some of these efforts below.

Quantitative Metrics of Drug Response
Quantitative methods of dose-response analysis have played a central role in the discovery of new cancer drugs and the elucidation of their mechanisms of action. The standard approach to quantify drug response has been relying on the generation of sigmoidal dose-response curves following exposure of cells to a drug over a wide range of concentrations. Conventionally, relative-viability measured at a fixed time-point (typically 72 or 96 h) following drug treatment has been used to extract quantitative metrics such as half-maximum inhibitory concentration (IC 50 ), maximal effect (E max ), Hill slope, and the area under the dose-response curves. These metrics have been routinely used to compare various aspects of drug sensitivity across cancer cell lines and drug types (Barretina et al., 2012;Fallahi-Sichani et al., 2013;Garnett et al., 2012;Haverty et al., 2016).
Despite its wide-spread usage, however, the relative-viability approach in assessing drug response suffers from a fundamental flaw, which is being confounded by variation in cell proliferation rates and assay duration. The reason is that cell count, which is used as a normalization factor in this approach, is nonlinearly time-dependent. As a consequence, the use of relative-viability in pharmacogenomic studies results in false-positive and false-negative associations . Therefore, new generation drug response metrics have recently been introduced to correct for this bias. The nature of these metrics is based on modeling druginduced changes in the net growth rate of the cancer cell population (instead of relative-viability) as a function of drug dose. These metrics include drug-induced proliferation (DIP) rate (Harris et al., 2016) and growth rate (GR) inhibition (Hafner et al., 2016), both of which consider and correct for the variability in growth rate that is irrelevant to drug treatment, for example, basal differences in cellular growth rate across tumor lines, natural variation in experimental conditions, or differences in duration of experiments. This is accomplished via normalizing the net growth rate of the drug-treated cell population to that of the untreated control.
Among the key quantitative metrics extracted from dose-response analyses are potency and efficacy. Potency is a measure of drug activity expressed in terms of the dose required to produce an effect of given intensity (e.g. 50% inhibition in growth rate). The more potent a drug is, the lower dose of that drug is required to evoke a given response. Efficacy, on the other hand, is the maximum response achievable from an applied drug, for example, at its maximum tolerated dose (MTD). Conventionally, studies of cancer drug response have been primarily focused on variation in drug potency, assuming that it is the most important difference between effective and ineffective drugs or sensitive and resistant cells. Multiparametric studies, however, have demonstrated that different dose-response metrics encode distinct but complementary information regarding how cells respond to a drug (Fallahi-Sichani et al., 2013;Hafner et al., 2017); high potency is not necessarily associated with high efficacy and vice versa. Furthermore, some drugs (e.g. inhibitors of the Akt/mTOR pathway in breast cancer cell lines) exhibit unusually shallow dose-response curves (Fallahi-Sichani et al., 2013). Classical pharmacology has no ready explanation for this phenomenon, but single-cell analysis shows that it correlates with significant cell-to-cell variability in drug response. Therefore, it is proposed that multiple dose-response metrics and cell-to-cell variability in drug response should be investigated together in order to accurately evaluate the effectiveness of cancer drugs. The most effective cancer drugs, and consequently the most reliable pharmacogenomic predictions, may be those that are based on the highest efficacy at clinically relevant drug concentrations .

Statistical Models of Gene-Drug Associations
The development of pharmacogenomic models for predicting therapeutic response has been a key step toward delivering the promise of precision cancer medicine. Moreover, these models have helped researchers prioritize the preclinical investigation of candidate compounds as potential cancer therapeutics. Such computational models aim at associating quantitative metrics of drug response in tumor cells (as described above) with molecular features derived, most commonly, from genomic or transcriptomic data collected for the same (or similar) samples. There are, however, challenges associated with the development of such models, including questions regarding the types of datasets to be used for training and testing the models, the nature of computational approaches that would be most appropriate, and perhaps most importantly, the extent of generalizability of the developed models (Azuaje, 2016).
Among the computational approaches widely applied are supervised statistical and machine learning techniques (Azuaje, 2016). Common supervised models for finding correlations between continuous input and output data include linear regression and logistic regression. Classification models, such as support vector machines (SVM), random forests and k-nearest neighbors (KNN) models, are also widely used for mining patterns in categorical data (Amin et al., 2014;Chen et al., 2015;Cortés-Ciriano et al., 2015;Stetson et al., 2014). Predictions from these models have helped with the stratification of tumors into subtypes with specific genetic or transcriptional biomarkers that are predictive of their response to different classes of drugs.
It is worth noting that capturing gene-drug associations often requires analysis of a plethora of data acquired at a system level. Developing predictive and biologically meaningful models on such high-dimensional data requires particular consideration of model complexity and the trade-off between variance and bias. Therefore, many current algorithms for learning from large-scale biological data incorporate regularization to modulate model complexity. Regularization is a common approach to avoid overfitting by penalizing high-valued regression coefficients. This is achieved by adding penalties to more complex models and then sorting all possible models from least overfit to greatest, in order to identify the one with the best predictive power. Commonly used penalty terms include ℓ 1 and ℓ 2 (Friedman et al., 2010;Zou and Hastie, 2005). ℓ 1 regularization, used for example in lasso regression, imposes penalties equal to the absolute value of the magnitude of model coefficients, in a way that some coefficients can become zero and eliminated. ℓ 1 regularization therefore leads to sparse models, that is, models with fewer coefficients. ℓ 2 regularization, used for example in ridge regression, imposes penalties equal to the square of the magnitude of model coefficients. ℓ 2 regularization encourages coefficients to be small, but doesn't force them to be zero, and thus does not lead to sparse models. Elastic net models combine ℓ 1 and ℓ 2 methods. Among these methods, the elastic net regularization has been applied, in multiple large-scale studies, to develop predictive models of drug response across hundreds of cancer cell lines treated with hundreds of cancer drugs (Barretina et al., 2012;Basu et al., 2013;Garnett et al., 2012;Haverty et al., 2016). These penalized regression models have enabled the identification of cooperative interactions among genetic mutations, chromosomal copy number, and transcripts across the genome and signatures of response for each drug. They have also identified new genomic predictors of drug sensitivity among tumor cells of specific subtypes, and explained how oncogene dependency in tumor cells could be modified by additional tissue-specific or transcriptional biomarkers. Finally, they have revealed some of the frequently mutated genes that are associated with sensitivity to a broad range of cancer therapeutics.
Additional efforts to develop statistical models of gene-drug associations include algorithms that combine different standard techniques in order to enhance the overall performance of the models. For example, the combination of Bayesian inference modeling with ridge regression has been reported to have an improved predictivity in comparison with the cross-validated ridge regression, while requiring only a fraction of the computational time (Neto et al., 2014). The combination of the elastic net and principal component analysis (PCA) has been another example of such efforts (Park et al., 2014).

Network Models of Gene-Drug Associations
The broad spectrum of molecular alterations-which may include hundreds or thousands of genetic mutations-that define the tumorigenic states, make it unlikely that tumor stratification and prediction of therapeutic responses could be efficiently achieved only based on individual mutations or single gene expression biomarkers. As such, systematic analyses of tumors in a variety of drug-screening projects have shown that therapeutic vulnerabilities and oncogenic mutations do not necessarily overlap. In many cases, mutations or single-gene transcripts alone are poor predictors of drug sensitivity (Basu et al., 2013). The prediction power, however, could be improved when multiple variables, including genetic and transcriptomic profiles of tumors, are considered as interactive networks. For example, the integration of protein-protein interaction networks and prior knowledge of drug-target interactions with elastic net regression modeling has led to prioritization of genomic and transcriptomic features that are more likely to be associated with cancer drug responses (Lee et al., 2018;Wang et al., 2016).
In addition, tumor profiling studies have shown that many human tumors lack clinically actionable mutations, broadly defined as changes in DNA that, if detected in a patient's tumor, would be predicted to affect patient's response to treatment. In many cases, tumors carry alterations either in oncogenes that are not currently druggable (e.g. RAS oncogene) or in genes with poorly characterized therapeutic value. In these cases, integrating genomic and transcriptomic profiles via network analysis and identifying key tumorigenic signatures and their master regulators have enhanced the predictivity of pharmacological targets (Alvarez et al., 2018;Ergün et al., 2007;Gu et al., 2014). The advantage of network models (in comparison with gene-by-gene analysis approaches) is that they distinguish essential regulators of tumorigenic states from those that may be only a consequence of such states and thus representing less crucial targets (Califano and Alvarez, 2017).
A variety of network-based algorithms have been developed to computationally identify key regulators, drivers and pathways involved in tumorigenesis (Califano et al., 2012;Pe'er and Hacohen, 2011). Among the most commonly applied mechanistic approaches to infer master regulators of tumorigenic states are those that integrate enrichment analysis methods (built based on, for example, the Kolmogorov-Smirnov test (Daniel, 1990)) as in gene set enrichment analysis (GSEA) (Subramanian et al., 2005), with de novo identification of tissue-specific protein targets, which are generally referred to as regulons (Alvarez et al., 2016). These algorithms particularly account for target directionality, distinguishing among activated and repressed targets, regulon overlap between candidate master regulators, as well as target inference confidence. Such network-based analyses have led to the identification of key proteins that regulate specific tumor phenotypes (Castro et al., 2016), subtype-specific oncogenic drivers (Remo et al., 2015;Tovar et al., 2015), and cancer drug resistance (Piovan et al., 2013;Rodriguez-Barrueco et al., 2015).
When integrated with functional genetic or compound screens, network-based inference of master regulators has led to systematic identification of actionable targets in human tumors. Such analyses have defined compounds and drugs that could be prioritized on the basis of their ability to invert the activity of master regulator proteins that mechanistically regulate tumorigenic states (Alvarez et al., 2018). Such a paradigm for drug discovery-that is based on a broader tumor-checkpoint dependencycomplements conventional approaches to targeted therapies that are based on driver oncogene dependencies in tumor cells. The underlying hypothesis for these approaches is that tumors are more dependent on the coordinated activity of master regulator proteins that regulate their transcriptional state stability than on the individual oncogenes that initiate them (Califano and Alvarez, 2017).

Models of Adaptive Drug Response and Resistance
Genomic, transcriptomic and proteomic profiling of tumors provide multidimensional snapshots of the complexity and heterogeneity of human tumors. Although extremely helpful in guiding new actionable targets and biomarkers of initial drug sensitivity, such molecular profiles are not sufficient to predict the outcome of therapy for most cancer patients or even tumor cultures in laboratory settings. Despite the initial benefit, tumor responses to molecularly targeted therapies are often incomplete and followed by relapse and progression to a drug-resistant state, which is currently a major challenge in cancer treatment. There is, therefore, an urgent need to develop effective approaches to identify therapeutic targets that either prevent or overcome drug resistance.
Sub-clonal expansion of drug-resistant clones is a well-defined route of tumor relapse, where the emergence of drug resistance is attributed to the extensive amount of genomic diversity across individual cells within a tumor. Numerous studies, however, have shown that resistance could also be driven through nonmutational mechanisms and by a small subpopulation of drug-tolerant cells, which remain viable upon drug exposure, even when the majority of the cell population is rapidly killed. The emergence of these persister cells has been attributed, in some cases, to the enrichment of a small subpopulation of stem-like cells, which are intrinsically refractory to the effect of anticancer drugs (Rambow et al., 2018;Roesch et al., 2010;Sharma et al., 2010). Alternatively, drug-induced adaptive mechanisms can drive cell-to-cell heterogeneity in drug response. These mechanisms are associated with the action of compensatory signaling circuits and homeostatic processes that diminish the efficacy of drugs (Lito et al., 2012(Lito et al., , 2013. Adaptive drug resistance is often explained by feedback activation of pro-growth signaling pathways, through either cellautonomous mechanisms such as up-regulation or rewiring of mitogenic cascades, or noncell-autonomous changes in the microenvironment. Finally, epigenetic reprogramming following inhibition of oncogenic signaling may also induce drug-tolerant states, generating sub-populations of reversibly drug-resistant cells (Fallahi-Sichani et al., 2017;Shaffer et al., 2017), which constitute a reservoir from which stably resistant clones are ultimately selected and drive disease progression (Hata et al., 2016).
In all cases, drug resistance is associated with dynamic changes that drive molecularly distinct tumor cell states (Konieczkowski et al., 2018). Drug combinations that work cooperatively to prevent or suppress pro-growth or pro-survival signaling in these states can improve response rates and block the emergence of drug resistance . Predictive computational modeling could play an important role to guide rational design of effective drug combination strategies (Du and Elemento, 2015). Among the most commonly applied modeling approaches to predict effective drug combinations are: (i) kinetic modeling of protein-protein interactions and feedback regulatory mechanisms to reproduce signaling dynamics in silico, and thereby predicting the effect of various network alterations (Kirouac et al., 2013), (ii) integrating high-throughput, time-and dose-dependent measurements of drug-induced molecular changes with phenotypic outcomes of drug response through the application of statistical modeling and machine learning algorithms (e.g. partial least square regression analysis), and thereby predicting efficient ways to overcome adaptive mechanisms of drug resistance (Fallahi-Sichani et al., 2015;Lee et al., 2012), and (iii) network modeling of drug response in order to identify key resistance signaling pathways or master regulators that could be targeted by individual or synergistic drug combinations .

Single-Cell Models of Heterogeneous Drug Response
Recent advancements in single-cell technologies have substantially enhanced our understanding of cellular heterogeneity at the level of DNA, RNA, proteins, metabolites, and their consequences for heterogeneity in drug response in tumor cells. For example, single-cell gene expression profiles derived from cancer cell lines or cohorts of cancer patients have shown that all tumors include cells that express markers of drug resistance, whereas the frequencies of these cells vary from one tumor to another, predicting eventual tumor progression to resistance following drug treatments (Shaffer et al., 2017;Tirosh et al., 2016). Multiple studies have focused on evaluating single-cell gene expression profiles in tumors as they respond to anti-cancer therapies, during which populations of drug-resistant cells are selected or generated as a consequence of drug adaptation (Rambow et al., 2018). A common theme among these studies is the transition of cells into distinct transcriptional states, indicating epigenetic reprogramming, and ultimately the selection of pre-existing resistance genotypes.
Cells and the molecular processes underlying their behaviors, including their responses to drugs, are highly dynamic. Time-lapse live-cell microscopy experiments (using genetically encoded reporters), followed by single-cell tracking, have revealed how individual cells (even those generated from the same mother cells) exhibit time-dependent and reversible heterogeneity at the level of signaling and phenotypic responses to drug treatment (Fallahi-Sichani et al., 2017). Furthermore, highly multiplexed singlecell imaging technologies have provided information regarding the spatial distributions of different cell types as well as localization of individual genes and protein expressions in tumors (Lin et al., 2018). Altogether, these studies have demonstrated the new insights and advantages that single-cell measurements and analysis could offer as compared with bulk tumor profiling assays. Single-cell profiling of tumors at different levels has been rapidly leading to the generation of rich datasets. Although just recently collected, these datasets have already begun to provide new opportunities for computational model development, thereby improving our understanding of the tumor heterogeneity and the complexity of drug response (Azizi et al., 2018). We expect the next generation of drug discovery models to be focused on challenges of single-cell analysis.

Conclusion
In this article, we have discussed how recent advancements in high-throughput molecular profiling in combination with systems biology approaches have led to the discovery of novel targeted therapies, genomic and network-level predictors of drug response, and mechanisms of drug resistance in tumors. In parallel with technological developments and rapid generation of large-scale, highdimensional datasets, effective computational modeling methods are needed to capture different aspects of tumor complexity and heterogeneity. These methods will continue to improve our understanding of tumor behavior, thereby opening new avenues to rational design of more successful therapeutic strategies.