Understanding intratumor heterogeneity by combining genome analysis and mathematical modeling

Abstract

Cancer is composed of multiple cell populations with different genomes. This phenomenon called intratumor heterogeneity (ITH) is supposed to be a fundamental cause of therapeutic failure. Therefore, its principle‐level understanding is a clinically important issue. To achieve this goal, an interdisciplinary approach combining genome analysis and mathematical modeling is essential. For example, we have recently performed multiregion sequencing to unveil extensive ITH in colorectal cancer. Moreover, by employing mathematical modeling of cancer evolution, we demonstrated that it is possible that this ITH is generated by neutral evolution. In this review, we introduce recent advances in a research field related to ITH and also discuss strategies for exploiting novel findings on ITH in a clinical setting.

INTRODUCTION

Cancer is composed of multiple cell populations with different genomes. Each of the populations is called a clone (or subclone) and this phenomenon is called intratumor heterogeneity (ITH). ITH is observed in various types of cancers and is presumed to be a major cause leading to therapeutic resistance. If a tumor harbors a major clone sensitive to a specific anti‐cancer treatment, the tumor shrinks within a given period after the treatment. However, in most cases, a minor clone resistant to the chemotherapy exists in the tumor and predominantly regrows despite the intensive therapy. It is supposed that ITH can be generated by clonal branching during cancer evolution.

In 1976, Nowell proposed that clones acquiring somatic mutations were subject to natural selection in a stepwise manner, by which cancer could originate from a single normal cell. According to the hypothesis, cancer evolution can be regarded as Darwinian evolution of a unicellular organism that divides through the process of asexual reproduction. After this clonal evolution model was proposed, molecular biologists discovered the involvement of proto‐oncogenes and tumor suppressor genes in the process of carcinogenesis. In 1990, Fearon and Vogelstein integrated these discoveries with the clonal evolution model to propose the multistep tumorigenesis model. Specifically, in colorectal tumorigenesis, while accumulating multiple causal mutations of significant genes, including APC, KRAS, TP53 and SMAD4, a normal epithelial cell linearly transforms through a benign lesion into a malignant tumor. Since then, the view has been accepted that linear clonal evolution creates a uniform population of malignant cells, although ITH was shown to exist by single gene‐focused experiments. However, genomic studies employing next‐generation sequencers have recently demonstrated that branching evolution is more predominant and generates more extensive ITH than we have presumed hitherto. In this review, we focus only on genomic ITH in solid tumors, especially ITH of somatic single nucleotide mutations, although there exist non‐genomic ITH, including ITH in transcriptome and methylome. First, we introduce the recent genomic studies aiming to unveil ITH.

UNVEILING INTRATUMOR HETEROGENEITY BY GENOMIC ANALYSIS

The advent of next‐generation sequencing technology enables us to obtain a comprehensive profile of somatic mutations in cancer; namely, cancer genome data. In various types of cancers, genomic studies involving a large number of different patients have created catalogs of driver genes, which play causal roles in tumorigenesis and progression, as shown in The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC). These studies have also demonstrated that remarkable heterogeneity of cancer genomes exists among patients, based on which the necessity of precision medicine is advocated. In addition to this intertumor heterogeneity, the presence of another kind of heterogeneity, ITH, has been revealed by genomic analysis. To analyze ITH, an approach called multiregion sequencing is commonly employed. Multiple samples obtained from physically separate regions within the tumor of a single patient are analyzed by next‐generation sequencing (Figure ). The samples from metastatic and recurrent lesions of the same patient are also sequenced along with a primary lesion. By multiregion sequencing, 2 categories of somatic single nucleotide mutations, “founder” and “progressor” mutations, are identified. These are present in all or some of the regions, respectively (they are also referred to by different terms in different studies: e.g. public and private mutations, respectively). Founder mutations are assumed to accumulate in the early phase of cancer evolution. The common ancestor clone that has acquired all the founder mutations then branches into subclones, which accumulate progressor mutations and contribute to the formation of ITH. Through these multiregion mutational profiles, we can infer an evolutionary history of the cancer by constructing a phylogenetic tree. For example, whole exome multiregion sequencing using multiple samples from primary and metastatic lesions in each of 10 renal cancer patients unveiled extensive ITH and divergent clonal evolution in renal cancer. Their study also revealed founder non‐silent mutations in some known driver genes such as VHL as well as progressor non‐silent mutations in other known driver genes such as SETD and BAP1. It is intriguing that, in some cases, independent mutations on different positions of the same driver gene were acquired by parallel evolution, strongly suggesting that a part of the ITH was generated by Darwinian evolution. We can also obtain information on ITH by performing deep sequencing in a sample from 1 region of a tumor (Figure ). An advantage of next generation sequencing technology is the ability to measure the allele frequency of a mutation by sequencing the genomic position of the mutation many times. Sequence reads of a higher coverage obtained by deep sequencing (typically from hundreds to thousands) enable us to calculate precise allele frequencies or to identify mutations of low allele frequencies. Based on mutant allele frequencies in a tumor, we can also evaluate the subclonal structure. If the tumor is clonal and diploid, the allele frequencies of all the mutations are 0.5. In realistic situations, because the tumor contains a small fraction of normal stromal cells, the mutation allele frequencies decrease to less than 0.5 and mutations whose genomic positions are subject to copy number alteration also show further deviated allele frequencies. Moreover, subclonal mutations are observed as mutations with low frequencies, which is not explained by the mixture of normal cells or copy number alteration. For example, in a study demonstrating whole genome deep sequencing of 21 primary breast tumors, the subclonal structures were explored by combining the information on mutant allele frequencies with allelic copy numbers. Furthermore, the order of mutations and copy number alterations was regarded as each phase in the evolutionary history of breast cancer. From the viewpoint of bioinformatics, exploring the clonal structure and inferring an evolutionary history from deep sequencing data is also a challenging task, for which multiple methods have been developed so far. Furthermore, in multiregion sequencing, local heterogeneity can be evaluated by obtaining a sufficient number of reads for each of the multiregion samples. Multiregion sequencing has been reported in various types of tumors, so far including brain tumors, breast cancers, colorectal cancers, esophageal cancers, head and neck cancers, hepatocellular carcinomas, lung cancers, melanomas, ovarian cancers, pancreatic cancers, prostate cancers and urothelial carcinomas, as summarized comprehensively in a previous review. We also performed whole exome multiregion sequencing in 9 cases of colorectal cancers to identify founder and progressor mutations in each case. The result obtained from 1 of the 9 cases is shown in Figure . Progressor mutations showed a mutational pattern which was geographically correlated with sampling locations. Moreover, we found that, in each region, founder mutations existed as clonal mutations while progressor mutations existed as subclonal mutations. This finding suggests that, even in each region, there existed ITH which was not captured by the resolution of multiregion sequencing. In addition, most of the mutations of known driver genes, such as APC and KRAS, were found to be founder mutations. Progressor mutations contained few driver mutations and parallel evolution was not confirmed, which is in contrast to the findings obtained in renal cancer. Moreover, frequent parallel evolutions were demonstrated in prostate cancer and low‐grade glioma but not in other kinds of cancer, suggesting that the origin of ITH cannot solely be explained by Darwinian evolution.

Multiregion sequencing. A, DNA samples from multiple regions of a single tumor were analyzed by next‐generation sequencing. B, Through multiregion mutation profiling, founder and progressor mutations were found as the common mutations in all the regions tested and only restricted regions, respectively. C, In a phylogenetic tree constructed from the multiregion mutation profile, the trunk and branches correspond to the founder and progressor mutations, respectively

Deep sequencing. A, A sample from a single region usually contains multiple types of cell populations with different genomes. B, Assume that the sample harbors normal cells, clone 1 with a clonal mutation, and clone 2 with the same clonal mutation and a subclonal mutation and all cell types are diploids without copy number alterations. C, By employing next‐generation sequencing, a sufficient number of sequence reads encompassing the positions of the target mutations are obtained to precisely estimate mutant allele frequencies which reflect the mixing proportion of each type in the whole cell populations

Multiregion sequencing of colorectal cancer. A, A schema of a multiregion sampling in a primary colorectal cancer and matched metastatic liver lesion. In this case, we obtained 20 samples from the primary lesion and 1 sample from the metastatic lesion. B, A multiregion mutation profile. The depth of red represents mutant allele frequency while the colors of sample labels were prepared so that the similarities of colors represent those of mutation patterns. C, A phylogenetic tree constructed from the multiregion mutation profile. The time when mutations in known driver genes of colorectal cancer is acquired is indicated along the tree

EXPLORING EVOLUTIONARY PRINCIPLES UNDERLYING INTRATUMOR HETEROGENEITY BY MATHEMATICAL MODELING

As described above, genome analysis has unveiled extensive ITH in various types of tumors. However, genome analysis is not sufficient for understanding how ITH is generated. To answer this question, mathematical modeling is a powerful tool when combined with genomic analysis. Next, we introduce some of the mathematical modeling studies that have sought to clarify the evolutionary principles underlying ITH.

Mathematical modeling can give us insights into system‐level principles underlying a phenomenon and is widely employed in the natural sciences and engineering disciplines as well as in the social sciences. In a mathematical modeling study, we start from expressing a target dynamic system as a mathematical model, using mathematical expressions such as differential equations. The mathematical model contains variables representing the system status and parameters specifying the system dynamics. After constructing the mathematical model, we usually try to express the variables as a function of the parameters and initial values of the variables. When a mathematical model is a simple differential equation, the model is analytically solvable. That is, it is possible to obtain the function by manipulating the mathematical expression. However, because most mathematical models are not analytically solvable, we use a computer simulation to numerically solve the models and analyze system dynamics. Therefore, computer simulation is currently an indispensable tool in mathematical modeling studies.

One of the mathematical models that is often used for modeling cancer evolution is an agent‐based model. An agent‐based model assumes there are system components called agents and defines rules that specify behaviors of the agents as well as interactions between the agents, and between the agents and environments. By employing an agent‐based model, we can construct a flexible model and easily analyze the system dynamics by simulation, if parameter and initial variable values are given. In mathematical modeling of cancer evolution, it is natural to assume each cell to be an agent and thereby ITH is easily expressed as the difference of internal states of agents. For example, in a pioneering model, agents were assumed to be cells which contained a few genes and proliferated while accumulating mutations. As a result, the computer simulation succeeded in reproducing ITH observed in single gene‐focused experiments. Since then, multiple mathematical modeling studies employing agent‐based models have been developed to shed light on the principles underlying the generation of ITH. For example, stem cell hierarchy may contribute to ITH and the interaction between cells as well as the turnover of cells in 3‐D space may affect the formation of ITH.

Because the existing models could not completely reproduce the extensive ITH revealed by our multiregion sequencing of colorectal cancer, we developed a new agent‐based model, the branching evolutionary process (BEP) model, to simulate heterogeneous cancer evolution. In a way similar to the other models, the BEP model assumes cells to be agents (Figure A). Each cell harbors n genes including d driver genes, while each cell divides and dies in a unit time with a probability p and q, respectively. When the cell divides, each gene is mutated with probability r, and if any driver genes are mutated, the division probability p increases 10^f‐fold per mutation. In the BEP model, f can be regarded as the strength of the driver genes. Given that a cell without mutations divides according to this rule, after the normal cell acquires the first driver mutation, which accelerates cell division, the proportion of the clone originating from the cell increases in a whole cell population. By repeating these steps, each cell gradually accumulates driver mutations as well as accompanying passenger mutations, which do not affect the cell division rate, and, finally, a tumor is formed with numerous mutations accumulated. Depending on parameter values in the course of the cancer evolution, each of the cancer cells can accumulate different combinations of mutations to generate different kinds of ITH. In Figures B,C, an example of snapshots of 2‐D tumor growth is shown, which was simulated using the BEP model with appropriate parameter values. In this example, driver mutations gradually accumulated in the cells and a clone with 4 mutations was selected during natural selection, which finally became dominant in the tumor.

Branching evolutionary process (BEP) model. A, Each cell has n genes (10 genes in this toy model) while each cell divides and dies in a unit time with a probability p and q, respectively. When the cell divides, each gene is mutated with probability r. If any of d driver genes (4 genes in this toy model) are mutated, the division probability p increases 10f‐fold per mutation. B, Evolutionary snapshots obtained by simulating 2‐D tumor growth based on the BEP model with an appropriate parameter setting. The region with the same color represents a clone with the same set of mutated genes. C, Single‐cell mutation profiles at 3 time points in the simulated tumor growth. Top colored bands represent clones, while the blue bands on the left represent driver genes

To explore the principles of ITH generation, we performed a large number of BEP simulations using a supercomputer with various parameter settings to find conditions leading to the extensive ITH observed in our genomic analyses. As a result, when cancer evolution was simulated with an assumption of a high mutation rate (i.e. with a large r value), followed by computer simulation of multiregion sequencing, we could reproduce mutation profiles similar to those obtained by our multiregion sequencing of colorectal cancers (Figure A,B). That is, irrespective of the presence of founder mutations, progressor mutations contributed to the formation of a heterogeneous mutation profile, which was geographically correlated with sampling locations. Moreover, we could also reconstruct local heterogeneity, as illustrated by the finding that progressor mutations existed as subclonal mutations in each region. Intriguingly, while driver mutations were acquired as founder mutations, progressor mutations contained few driver mutations and most of them consisted of neutral mutations that did not affect the cell division rate. This suggests that, after the appearance of the common ancestor clone with accumulated driver mutations, extensive ITH was generated by neutral evolution. In neutral evolution, which was introduced by Motoo Kimura, neutral mutations that are neither advantageous nor disadvantageous for survival and growth prevail in a population by chance. This contrasts with classical Dawinian evolution, which is caused by natural selection; therefore, when introduced, the neutral evolution theory was strongly opposed by supporters of Darwinian evolution theory. However, it is currently well‐accepted that most of the genetic diversity among biological species is generated by neutral evolution. Moreover, single‐cell mutation profiles of the simulated tumor suggest that the tumor consists of a large number of minute clones with numerous neutral mutations accumulated (Figure C). This extensive ITH generated by neutral mutation could be a fundamental cause of therapeutic failure. Because whether a mutation is neutral or not depends on the surrounding environment, an environmental change induced by a specific therapy can convert a neutral mutation which has no effect before the therapy to a driver mutation leading to therapeutic refractoriness (resistant mutation). This means that any type of therapy can potentially generate a resistant clone with a resistant mutation, which leads to the tumor relapse even if the therapy is temporarily effective (Figure ).

A computer‐simulated tumor with extensive intratumor heterogeneity (ITH) generated by neutral evolution. A, A tumor depicted by branching evolutionary process (BEP) simulation with an assumption of a high mutation rate. B, A simulated multiregion mutation profile of the simulated tumor. Cell populations in the regions labeled with A‐H (A) were extracted and their averaged mutation profiles were obtained. Note that the simulated profile is similar to the real one obtained from the colorectal cancer (Figure B), and that driver mutations consisted only of founder mutations. C, A simulated single‐cell mutation profile of the simulated tumor, suggesting the existence of numerous clones that cannot be detected by multiregion sequencing

Acquisition of therapeutic resistance by neutral evolution. Among numerous subclones generated by neutral evolution, a resistant clone emerges and expands at an accelerated pace, thereby leading to tumor relapse

Similar to our study, there is a recent report demonstrating the presence of uniformly high ITH in colorectal tumors unveiled by genome analysis through multiregion and deep sequencing. Moreover, the authors found that progressor mutations identified in some regions were often scattered in geographically separated regions as subclonal mutations of low allele frequencies. This phenomenon could be reproduced by computer simulation, and the Big Bang model was proposed. In the Big Bang model, a large number of subclones are generated in the early phase of cancer evolution, and then these subclones expand without natural selection, while partially mixing, to eventually show uniformly high ITH in every region of the tumor. Another study involving liver cancers demonstrated that a statistic value indicative of clonal diversity in multiregion sequencing data was consistent with a theoretical value derived analytically from the neutral evolution model. The claim that neutral evolution accumulated numerous mutations in the tumor is consistent with ours. An approach to test the neutrality using deep sequencing data has also been proposed. Mathematical modeling proved that allele frequencies of new mutations decreased inversely with tumor size in neutral evolution, based on the fact that the new mutations in the genome only represent labels of individual cell lineages. In the case where Darwinian evolution works, both subclonal driver and passenger mutations are inherited at a higher frequency than expected under neutral evolution, thereby generating signatures of clonal outgrowth. That is to say, there exist “too many” mutations at a high frequency. Based on this logic, we can determine whether the distribution of allele frequencies in deep sequencing data originates from neutral or Darwinian evolution. By applying this approach to whole exam sequencing data of multiple cancer types, the degree of contribution of natural selection was examined in each cancer type. In another recent study, this approach was developed further to analyze multiregion sequence data, where the distributions of allele frequencies in 2 separate regions of a single tumor were examined together to discriminate between Darwinian and neutral evolution. The researchers analyzed multiregion sequencing data to find evolutionary modes specific to different tumor types while simulation of cancer evolution and multiregion sequencing were employed to benchmark the proposed approach. In addition to the cancer type‐specific observation of parallel evolution, these results showed that the neutral evolution model was evident in some types of tumors while the contribution of Darwinian evolution was prominent in other types. However, the reason why the evolutionary principles underlying ITH differ in different types of cancers remains unclear. We expect that mathematical modeling will be a promising and powerful tool for gaining an insight into this issue.

PERSPECTIVE

As described so far, substantial progress has been made in the understanding of ITH. As a next step, we should apply this understanding to a clinical setting to tackle the problems related to therapeutic resistance. Although a large number of molecular target drugs are currently available or under development, drug‐resistance frequently appears during chemotherapy in most of the drugs, which leads to therapeutic failure. Mathematical modeling assuming simple ITH has shown that the regrowth of tumors can be delayed or prevented by adjusting a therapeutic regimen. The normal clinical practice is that an anti‐cancer drug is continuously administered to a cancer patient in the maximum tolerated dose, if possible. Given that a tumor is composed of major and minor clones which are sensitive and resistant to the chemotherapy, respectively, the tumor temporarily shrinks because the major chemosensitive clone is eradicated. However, as the major chemosensitive clone disappears, the minor resistant clone can grow freely due to the release from growth competition; that is, the competitive relationship between the 2 clones is dissolved. On the other hand, if we can keep the 2 clones in a competitive state while controlling the total tumor volume within an acceptable threshold, the survival of the patient can be prolonged as compared with the routine continuous administration. Thus, the usefulness of several alternative treatment regimens has been proposed by mathematical modeling. For example, in “adaptive therapy”, the initial dose of an anti‐cancer drug is high and then the dosage is decreased as the tumor shrinks to eventually maintain the sensitive clone at a level sufficient to suppress the growth of the resistant clone. Similarly, in “metronomic therapy” an anti‐cancer drug is continuously administered in a low dose. Treatment regimens using a combination of drugs are also proposed as exemplified by “double bind therapy.” In this therapy, 2 different kinds of drugs are administered alternately to expose each of the 2 clones to their respective effective drugs, which can theoretically keep the tumor volume at a constant level. When combined with novel findings obtained by recent advances in genome analysis, these treatment strategies can be further refined. In particular, liquid biopsy based on circulating tumor DNA (ctDNA) appears to be an essential tool for monitoring the efficacy of the treatment and for detecting recurrences in the early phases. ctDNA is also referred to as a tumor‐derived portion of cell‐free DNA (cfDNA), which is all non‐encapsulated DNA circulating in the bloodstream. By applying digital PCR or deep sequencing to cfDNA extracted from patients' plasma, we can non‐invasively detect mutations in ctDNA. The allele frequencies of the mutations in ctDNA are supposed to reflect the real‐time clonal proportions in the whole tumor, including primary and metastatic lesions, which can provide an opportunity for tracking the temporal dynamics of ITH during a therapeutic course. For example, the time‐series data acquired by digital PCR of ctDNA showed the emergence of a resistant KRAS mutation during anti‐EGFR therapy in colorectal cancer patients. Moreover, mathematical analysis of the time‐series data suggested that the KRAS mutation already existed in the tumor before the initiation of chemotherapy, which is consistent with the view derived from our neutral evolution model (Figure ). In a recent study designed to combine multiregion sequencing of the tumor and deep sequencing of ctDNA, first, surgical samples of early‐stage lung cancer were analyzed by multiregion sequencing to identify target mutations which can trace existent clones in the primary tumors. Then, deep sequencing of the target mutations was performed using time‐series ctDNA samples after the surgery. The researchers succeeded in unveiling the clonal evolutionary dynamics during the treatment period to the eventual tumor relapse. These time‐series analyses utilizing ctDNA‐based liquid biopsy are expected to give us deeper insight into the heterogeneous cancer evolution responsible for therapeutic failure and are indispensable for the realization of anti‐cancer therapies based on mathematical modeling.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.