Abstract
We study a situation that arises in the somatic evolution of cancer. Consider a finite population of replicating cells and a sequence of mutations: type 0 can mutate to type 1, which can mutate to type 2. There is no back mutation. We start with a homogeneous population of type 0. Mutants of type 1 emerge and either become extinct or reach fixation. In both cases, they can generate type 2, which also can become extinct or reach fixation. If mutation rates are small compared to the inverse of the population size, then the stochastic dynamics can be described by transitions between homogeneous populations. A “stochastic tunnel” arises, when the population moves from all 0 to all 2 without ever being all 1. We calculate the exact rate of stochastic tunneling for the case when type 1 is as fit as type 0 or less fit. Type 2 has the highest fitness. We discuss implications for the elimination of tumor suppressor genes and the activation of genetic instability. Although our theory is developed for cancer genetics, stochastic tunnels are general phenomena that could arise in many circumstances.
IN 1951, Hermann J. Muller proposed that cancer arises when a single cell receives multiple mutations (Muller 1951). This perspective, which has been confirmed by numerous studies of cancer genetics, is the primary motivation for the present article. We describe the dynamics of sequential mutations in the context of somatic evolution of cancer.
The majority of human cancers arise in epithelial tissues. These tissues are organized into small compartments of cells (Mintz 1971; Kovacs and Potten 1973; Bachet al. 2000). We analyze the population dynamics of one such compartment. Often compartments are subdivided into stem cells and differentiated cells (Cairns 2002). The former divide throughout the life of the individual and replenish the compartment, whereas the latter have a limited lifetime. This design gives rise to interesting population dynamics, which we consider elsewhere (Franket al. 2003; Michoret al. 2003a), but ignore here. In this article, we consider a population of equivalent cells. The analysis can be interpreted as applying to the subpopulation of stem cells within a compartment or to an effective population size defining the spatial scale of somatic selection within a tissue. We answer the following question: How long does it take to accumulate two mutations in a compartment of N cells? We consider three types of cells, labeled 0, 1, and 2, because they contain 0, 1, and 2 mutations, respectively. Type 1 and type 2 cells can have a reproductive rate (somatic fitness) that can differ from type 0 cells. In general, we are interested in the case where the fitness of type 1 is equal to or less than the fitness of type 0, while the fitness of type 2 is greater than the fitness of type 0. This scenario describes, for example, the inactivation of tumor suppressor genes or activation of genetic instability (Figure 1).
Tumor suppressor genes (TSGs) require inactivating mutations in both alleles (Knudson 1971, 2001; Moolgavkar and Knudson 1981). In a normal cell, according to our notation type 0, both alleles are functioning. In a type 1 cell, one allele has been inactivated, but the cell retains a normal phenotype. In a type 2 cell, both alleles are inactivated. Often cells with an inactivated TSG fail to undergo programmed cell death (apoptosis) and therefore have a tremendous selective advantage compared to normal cells. There are various mutational mechanisms that can inactivate TSGs, including point mutations, insertions, deletions, loss of whole chromosomes, or mitotic recombination. Often the inactivation of a TSG results in a local accumulation of cells (neoplasia), which provides a first step toward cancer. Therefore, we calculate how long it takes for a compartment of cells to inactivate a TSG and thereby initiate tumor growth. Many genetic pathways to cancer require the inactivation of several TSGs (Fearon and Vogelstein 1990; Kinzler and Vogelstein 1998; Luebeck and Moolgavkar 2002). In this context, our theory helps to calculate how long it takes to eliminate the next TSG in an already established population of precancer cells (sometimes called adenoma).
Another application of our theory is genetic instability (Lengaueret al. 1998). Since cancer cells accumulate many genetic changes, it was proposed that certain mutations increase the rate at which subsequent mutations occur (Loebet al. 1974; Jackson and Loeb 1998; Loeb 1998; Lengaueret al. 1998). This concept is called genetic instability. A main question of all of cancer biology is to what extent genetic instability is an early event and therefore a driving force of tumorigenesis or a late stage consequence (Kiberstis and Marx 2002; Sieberet al. 2002). Our theory addresses the following situation: the first mutation induces genetic instability and thereby increases the rate at which the second mutation might occur. How long does it take for the second mutation to occur, especially when the first mutation has a selective cost? This question is important when evaluating the possibility that genetic instability is a driving force of cancer progression (Nowaket al. 2002; Komarova et al. 2003a,b; Michoret al. 2003b).
Chromosomal instability (CIN) is a particular form of genetic instability that refers to increased rates of gaining or losing whole chromosomes or arms of chromosomes (Lengaueret al. 1998). In yeast, several hundred genes that contribute to maintaining the stability of chromosomes during cell division are known. Mutations in such genes can trigger the CIN phenotype. A number of human CIN genes have been discovered (Ragajapolanet al. 2003). We are interested in calculating the time until a population of CIN cells has inactivated a TSG. The first allele of the TSG is inactivated by a point mutation occurring with a probability of 10^{–7}/gene/cell division. The second allele of the TSG is inactivated by losing the chromosome (arm), which contains the locus of that gene. In CIN cells the probability for this event was measured to be 10^{–2}/chromosome/cell division (Lengaueret al. 1998). Hence, the mutation rates for the first and second event differ by many orders of magnitude.
Throughout the article we assume that the population size is less than the inverse of the mutation rate from type 0 to type 1. In this case, we can approximate the evolutionary dynamics as a Markov process describing the transition from a population that contains only type 0 cells to one that contains only type 1 cells to one that contains only type 2 cells. In addition there is the possibility that the population moves from all type 0 to all type 2, without ever visiting the all type 1 state. We refer to the latter transition as “stochastic tunneling.” In previous articles we have derived tunneling rates for various limits (Komarova et al. 2003a,b) and we have used tunnels to evaluate the role of chromosomal instability for cancer initiation (Nowaket al. 2002). Here we introduce a new stochastic description and provide a mathematical derivation of the tunneling rate that holds for a wide range of population sizes.
In the model section, we present the basic theory. In subsequent sections, we calculate the rate of tunneling, derive simple approximations for the rate of tunneling, compare our calculations with computer simulations, and then discuss simultaneous double mutations. In the discussion, we summarize our results.
MODEL
Consider a population of N asexually reproducing cells. There are three types of cells, denoted by 0, 1, and 2, respectively. Type 0 mutates to type 1 at rate u_{1}, while type 1 mutates to type 2 at rate u_{2}. Direct mutation from type 0 to type 2 is possible in principle, but can normally be neglected for reasons that we discuss later.
The relative reproductive rates (fitness) of cell types 0, 1, and 2 are given by 1, r, and a, respectively. Usually we are interested in the situation where r ≤ 1 and a > 1. Thus, the intermediate type, 1, is either neutral or selectively disadvantageous compared to type 0, while type 2 has a higher fitness than type 0.
We assume that cells reproduce asynchronously. Therefore we use a Moran process instead of the standard WrightFisher model. Each elementary step of the stochastic process consists of a birth and a death event. For birth, one of the N cells is chosen at random proportional to fitness. It will give rise to an offspring subject to possible mutation. For death, one of the N cells is chosen at random. The total population size, N, is strictly constant.
Time is measured in units of cell divisions. If cells divide on average once per day, then our timescale is given in days. At time t = 0, all N cells are of type 0. After some time, a single mutant cell of type 1 will be generated. This cell will lead to a lineage of type 1 cells that will either go to extinction or go to fixation. The probability of fixation starting from a single type 1 cell with fitness r is given by
If Nu_{1} ≪ 1, the mean waiting time until the appearance of a successful mutant is much longer than the average time required for the fixation of the mutant. Hence, we can describe the evolutionary dynamics as a simple Markovian jump from an all 0 population to an all 1 population. The transition rate between these two states is given by Nu_{1}ρ(r). The waiting time for the transition follows a negative exponential distribution with mean 1/Nu_{1}ρ(r).
After the fixation of type 1, a similar process describes the transition from an all 1 population to an all 2 population. If Nu_{2} ≪ 1, we can regard this transition as a Markovian jump occurring at a rate Nu_{2}ρ(a/r). The waiting time for the transition from all 1 to all 2 follows a negative exponential distribution with mean 1/Nu_{2}ρ(a/r). The probability of fixation of a type 2 cell in a population of N – 1 type 1 cells is ρ(a/r).
If the mutation giving rise to the second mutant is fast and if the fitness advantage of the second mutant is large, the second mutant might appear and spread before the fixation of the first mutant (Figure 2). This leads to a onestep transition from the all 0 population to the all 2 population without ever visiting an all 1 population. This phenomenon is called stochastic tunneling (Nowaket al. 2002; Komarova et al. 2003a,b). Tunneling is especially important when the intermediate mutant is deleterious and the population size is large, because then the expected waiting time for its fixation is very long.
Denote by x_{0}, x_{1}, and x_{2}, respectively, the probability that the system is in state all 0, all 1, and all 2. The Markovian jump process is described by the forward Kolmogorov differential equation
THE RATE OF TUNNELING
Whenever a new mutant cell of type 1 is produced in a population of type 0 cells, the lineage starting from it can temporally increase in abundance, but will eventually become extinct in the majority of cases. The probability of extinction is given by 1 –ρ(r). Tunneling occurs when the lineage of type 1 mutants, before extinction, produces a type 2 mutant that reaches fixation.
Suppose a new type 1 mutant is produced at time 0. Let Y(t) denote the number of cells of the lineage arising from this one mutant cell (Figure 2). For a given trajectory Y(t), the probability of appearance and fixation of a type 2 mutant is given by 1 – P, where
The calculation of the tunneling probability can be carried out by considering a doubly stochastic process. There are two sources of stochasticity: the first one is the stochasticity generating the trajectory of type 1 mutants, Y(t); the second one is the stochasticity generating the spread of type 2 mutants. If the average with respect to the former stochasticity is denoted by E[*], then the rate of tunneling can be expressed as
Let us define
EXPLICIT FORMULAS FOR THE RATE OF TUNNELING
We can derive explicit formulas for the tunneling rate by further examining Equation 7. According to the argument in appendix a, we can derive the following approximate formula:
When one of the terms within the square root symbol in Equation 10 is much larger than the other, we have simple expressions for two extreme cases,
COMPUTER SIMULATIONS
In Figures 3 and 4, we compare our analytical results with direct computer simulation of the Moran process. At each elementary step of the Moran process, one cell is chosen for reproduction, possibly with mutation, and one cell is chosen for elimination. Thus the total population size remains constant. For each parameter choice, we compute many independent runs of the stochastic process. Then we calculate the fraction of runs that have resulted in fixation of type 2 cells by a certain time. The results are compared with the solution of the Kolmogorov equation (Equation 1), using different formulas for the tunneling rate, R. There is perfect agreement if we use Equations 6 and 8 for calculating the tunneling rate. There is excellent agreement if we use the simplified Equation 11, except when N is very close to the critical value at which tunneling starts to become important. We also show that ignoring the tunnel totally fails once the population is larger than a critical size.
Figure 3 shows simulations of the neutral tunnel, where r = 1. Solid circles indicate the results of the computer simulation. The solid curve uses Equations 6 and 8 for the tunneling rate. The broken curve, which is in excellent agreement with the data, uses Equation 11. The other broken curve, the one that remains constant for larger values of N, neglects tunneling. This is accurate only if
Figure 4 shows simulations of tunneling through a deleterious intermediate mutant, r < 1. Again, the tunneling rate using Equations 6 and 8 (solid curve) is a perfect fit. Using Equation 11 also leads to an excellent agreement. Neglecting the tunnel is completely wrong if N is above a critical size, which is implicitly given by ru_{2}ρ(a)/(1 – r) >ρ(r).
The mathematical analysis assumes that the lineages starting from different mutants of type 1 behave independently, which is implicit in the calculation of the doubly stochastic process. This assumption holds if mutations from type 0 to type 1 occur infrequently, Nu_{1} ≪ 1. If this inequality does not hold, then a lineage starting from a single type 1 might not go extinct before the next mutation creating a type 1 mutant arises. This effect should cause a deviation from our results for very large N.
DIRECT MUTATION
In this article, we assumed that mutation from type 0 to type 1 and the mutation from type 1 to type 2 must occur in separate cell division events. We also must consider the possibility that these two mutations occur in a single event of cell division, causing a direct mutation from type 0 to type 2. If the mutation rates are u_{1} and u_{2}, then double mutants might arise at rate u_{1}u_{2}. Thus, direct mutation would result in the transition from an all 0 to an all 2 population at rate R_{direct} = Nu_{1}u_{2}ρ(a). In contrast, the tunneling rate is R = Nu_{1} ((r)/(1 – r))u_{2}ρ(a), when the population size is sufficiently large. These two rates have the same order of magnitude with respect to the mutation rate. Comparing these two rates, we conclude that direct mutation is faster than tunneling if r < 0.5. Otherwise tunneling is faster.
 Download figure
 Open in new tab
 Download powerpoint
 Download figure
 Download figure
 Open in new tab
 Download powerpoint
 Download figure
If the intermediate mutant is neutral, r = 1 (Figure 5), then tunneling,
We note that double mutation could be mechanistically impossible when we consider tunneling via a genetic instability mutation. Suppose the first mutation occurs in a gene that triggers genetic instability, which then increases the rate of the second mutation. The rate increase, however, may come into effect only during the next cell division. In this case, tunneling with a fast second step is much faster than a double mutation with two slow steps.
DISCUSSION
We have provided a theory that allows us to calculate the rates at which mutations happen in the somatic evolution of cancer. This is the continuation of an ongoing effort to understand dynamics of cancer progression (Nordling 1953; Armitage and Doll 1954, 1957; Fisher 1958; Cooket al. 1969; Bell 1976; Goldie and Coldman 1979, 1985; Moolgavkar and Knudson 1981; Wheldon 1988; Chaplain 1995; Byrne and Chaplain 1996; Owen and Sherratt 1999; BarOret al. 2000; HerreroJimenezet al. 2000; Wodarz and Krakauer 2001; Plotkin and Nowak 2002; Gatenby and Maini 2003; Little and Wright 2003).
The inactivation of a tumor suppressor gene requires, first, a neutral mutation and, second, an advantageous mutation (“advantageous” for the cell, not for the organism). How long does it take to eliminate a tumor suppressor gene in a population of N cells? We find that x_{2}(t) ≈ t^{2}Nu_{1}u_{2}ρ(a)/2 if
Genetic instability requires a mutation that might be deleterious for the cell. This mutation, however, increases the mutation rate and therefore the chance that this cell undergoes an advantageous mutation. Thus the essence of the emergence of genetic instability is described by a pathway from type 0 to 1 and 2, where type 1 is less fit than type 0, but type 2 is more fit than type 0. Denote by x_{2}(t) the probability that a population of N cells has made both mutational steps at time t. The probability x_{2}(t) is given by a convolution integral of two negative exponential distributions. For Nu_{2}t ≪ 1, we obtain the approximation x_{2}(t) = t^{2}Nu_{1}u_{2}ρ(r)ρ(a)/2. For Nu_{2}t ≫ 1, we obtain the approximation x_{2}(t) = tNu_{1}ρ(r) as long as N is below a critical value. If instead N is greater than this critical value, we have x_{2}(t) = tNu_{1}u_{2}(r/(1 – r))ρ(a). The critical value is implicitly given by u_{2}r/(1 – r) =ρ(r)/ρ(a). If N is near the critical value, then the more accurate tunneling equations should be used.
The explicit formulas for the rate of tunneling are useful for demonstrating the effect of various parameters and for providing an intuitive understanding. For example, the approximate formula for the deleterious tunnel, R = N(u_{1}/(1 – r))ru_{2}ρ(a), is a product of three factors. The first factor is the product of the population size N and the frequency of mutationselection balance, u_{1}/(1 – r), indicating the expected number of type 1 mutants in the population. The second factor, ru_{2}, denotes the mutation rate from type 1 to type 2. The third factor, ρ(a), is the probability of successful fixation of the second mutant once it appears in the population. Such a simple and intuitive understanding is not available for the corresponding formula when the type 1 mutant is neutral,
When the population size is very large, the fixation of either deleterious or neutral mutants is difficult. As a consequence, the relative importance of tunneling increases. In the limit of a large total population size, we can formulate a branching process in continuous time. For example, suppose that cells are of two types: type 1 and type 2 have reproductive rates of r and a, and both types have death rate of 1. There is a small mutation rate from type 1 to type 2, but no back mutation. Using the theory of multitype branching processes, we can calculate the probability for a single cell of type 1 to produce a lineage of type 2 cells that proliferates exponentially in time. This probability is equal to (r/(1 –r))u_{2}ρ(a) when the type 1 mutant is deleterious and is equal to
Although we were motivated by the somatic evolution of cancer, stochastic tunneling may be useful also in other evolutionary processes. An important simplification is that genetic recombination can be neglected in the somatic evolution of cancer, while it is important for sexually reproducing organisms.
Tunneling is important when the intermediate mutants have a lower fitness than the wild type, but the second mutant is advantageous. There is some similarity to compensatory neutral evolution (Kimura 1985; Iizuka and Takefu 1996; Michalakis and Slatkin 1996; Phillips 1996; Stephan 1996). For example, in helical regions of RNA molecules, a change in one site destabilizes the secondary structure. This mutation can be compensated by a mutation in the opposite site of the base pair (Higgs 1998; Innan and Stephan 2001; Savillet al. 2001). In this case, the second mutant has the same fitness as the original type, while in the situation discussed in this article, the second mutant has a greater fitness than the wild type. A more important difference is that the studies of compensatory neutral mutations focus on the effect of recombination. The evolution of a fitter genotype via intermediate deleterious genotype(s) is more difficult in the presence of recombination than in its absence. For somatic evolution of cancer, we can neglect recombination.
The results presented here provide a step toward constructing a comprehensive mathematical theory for the somatic evolution of cancer. We need to understand the nature of ratelimiting steps in cancer progression and the effect of cellular population size and genetic instability.
Acknowledgments
We thank Natalia Komarova formany helpful discussions. This work was completed during Y.I.'s visit to the Institute for Advanced Study, Princeton, New Jersey. Support from the David and Lucile Packard Foundation and Jeffrey Epstein is gratefully acknowledged.
APPENDIX: DERIVING EXPLICIT FORMULAS FOR THE RATE OF TUNNELING
Let φ(i/N) = V_{i} and y = i/N. Calculating the Taylor expansion of φ(y), Equation 7 is rewritten as
Footnotes

Communicating editor: J. B. Walsh
 Received September 19, 2003.
 Accepted December 11, 2003.
 Copyright © 2004 by the Genetics Society of America