Stochastic Tunnels in Evolutionary Dynamics

Corresponding author: Yoh Iwasa, Kyushu University, Fukuoka 812-8581, Japan., yiwasscb{at}mbox.nc.kyushu-u.ac.jp (E-mail)

Communicating editor: J. B. WALSH

	ABSTRACT

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

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 ; BACH et 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 (FRANK et al. 2003 ; MICHOR et 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 (Fig 1).

View larger version (46K): In this window In a new window Download PPT slide   Figure 1. Transition dynamics. (a) The basic architecture of our model contains three types, 0, 1, and 2, with fitness 1, r, and a, respectively. Mutation rates are given by u1 and u2. (b) As an example, we consider the inactivation of a tumor suppressor (TSP) gene. A wild-type cell, TSP+/+, mutates to a cell with one inactivated allele, TSP+/–, which mutates to a cell with two inactivated alleles, TSP–/–. The somatic fitnesses are given by 1, 1, and a > 1. (c) Another example is the emergence of chromosomal instability, CIN. We start with a cell where one allele of a TSP gene has been inactivated. This cell receives a mutation that triggers CIN. CIN might have a selective cost, r < 1, but greatly accelerates the loss of the second allele of the TSP gene. Such transitions are parts of general mutation-selection networks of cancer evolution.

Tumor suppressor genes (TSGs) require inactivating mutations in both alleles (KNUDSON 1971 , KNUDSON 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 (LENGAUER et al. 1998 ). Since cancer cells accumulate many genetic changes, it was proposed that certain mutations increase the rate at which subsequent mutations occur (LOEB et al. 1974 ; JACKSON and LOEB 1998 ; LOEB 1998 ; LENGAUER et 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 ; SIEBER et 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 (NOWAK et al. 2002 ; KOMAROVA et al. 2003A , KOMAROVA et al. 2003B ; MICHOR et 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 (LENGAUER et 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 (RAGAJAPOLAN et 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 (LENGAUER et 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 , KOMAROVA et al. 2003B ) and we have used tunnels to evaluate the role of chromosomal instability for cancer initiation (NOWAK et 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

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

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₁, while type 1 mutates to type 2 at rate u₂. 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 Wright-Fisher 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

This probability is 1/N if the mutant is neutral, r = 1, but is <1/N if the mutant is deleterious, r < 1 (DURRETT 2002 ; KOMAROVA et al. 2003A , KOMAROVA et al. 2003B ). When the population size N is finite, there is always a positive probability of fixation even for deleterious mutants.

If Nu₁ << 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₁(r). The waiting time for the transition follows a negative exponential distribution with mean 1/Nu₁(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₂ << 1, we can regard this transition as a Markovian jump occurring at a rate Nu₂(a/r). The waiting time for the transition from all 1 to all 2 follows a negative exponential distribution with mean 1/Nu₂(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 (Fig 2). This leads to a one-step transition from the all 0 population to the all 2 population without ever visiting an all 1 population. This phenomenon is called stochastic tunneling (NOWAK et al. 2002 ; KOMAROVA et al. 2003A , KOMAROVA et al. 2003B ). 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.

View larger version (29K): In this window In a new window Download PPT slide   Figure 2. Stochastic tunneling. In a homogeneous population of type 0 cells, type 1 mutants with fitness r 1 arise. They go to fixation or to extinction. In the latter case, type 1 cells might produce type 2 cells before becoming extinct. Type 2 cells have fitness a > 1 and are likely to reach fixation in the population. Hence, the system can go from all 0 to all 2 without having ever been in all 1.

Denote by x₀, x₁, and x₂, 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

(1)

In the following section, we calculate the rate of tunneling, R, from the all 0 state to the all 2 state. The solutions of system (1) with various expressions for R are compared with direct stochastic simulations of the Moran process.

	THE RATE OF TUNNELING

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

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 (Fig 2). For a given trajectory Y(t), the probability of appearance and fixation of a type 2 mutant is given by 1 – P, where

(2)

This expression is the zeroth term of a Poisson distribution. The exponent is given by the cumulative number of type 1 cells multiplied by the mutation rate from type 1 to type 2, ru₂, and the probability of fixation, (a), of one type 2 cell that arises in a population of type 0 cells. Note that using (a) is an approximation, because the type 2 cell arises in a population that contains a mixture of type 0 and type 1 cells. The probability that a particular trajectory, Y(t), does not generate a successful type 2 mutant is given by P.

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

(3)

Hence, the rate of tunneling is the product of three factors: the rate of producing type 1 mutants times the probability that the lineage of type 1 mutants will become extinct times the conditional expectation of initiating a successful type 2 lineage. The expectation is subject to the conditions that the initial number of type 1 cells is 1, Y(0) = 1, and that the lineage will eventually become extinct, Y() = 0. In this way, we exclude the contribution of the spread of type 2 mutants after fixation of type 1 mutants. Let us now calculate the conditional expectation. The type 1 lineage either goes to fixation, Y() = N, or goes to extinction, Y() = 0. There is no possibility for Y(t) to stay between 0 and N forever. Therefore, we have

If Y becomes fixed, Y() = N, and then the integral ₀ Y(t)dt diverges to infinity and P = 0. Since the second term in the right-hand side is zero, Equation 3 can be rewritten as

(4)

Note that the above expression for R is always positive, because P < 1 and P = 0 if the intermediate mutant becomes fixed, Y() = N.

Let us define

(5)

Note that V₀ = 1 and V_N = 0. The expectation in Equation 4 is equal to V₁. We can write

(6)

We have used Pr[Y() = 0] = 1 – (r) for the extinction probability. According to the Moran model, the transition from i to i + 1 occurs at rate (N – i)(ir/(ir + (N – i))). The transition from i to i – 1 occurs at rate i((N – i)/(ir + (N – i))). The system remains in i otherwise. To derive a recursive formula for V_i, we expand the function according to events occurring in a short time interval of length t.

By rearranging terms, dividing both sides by t, and taking the limit when t is very small, we obtain

(7)

This recursion holds for i = 1, 2, ... , N – 1. The boundary conditions are V₀ = 1 and V_N = 0. Using Equation 7, we can determine the solution of V₁. We can use an iterative approximation to obtain a numerical solution of Equation 7. For example, we start with the initial distribution, V_i = 1 – i/N for i = 0, ... , N, and calculate V_i in the next round by

(8)

When this iteration converges, the stationary distribution of V_i satisfies Equation 7. The value for V₁ together with Equation 6 gives the rate of tunneling.

	EXPLICIT FORMULAS FOR THE RATE OF TUNNELING

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

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:

(9)

As we pointed out earlier, 1 – V₁ – (r) in Equation 6 is always positive from the definition of a probability. If V₁ is replaced by the simpler formula shown above, then 1 – V₁ – (r) can become negative. Hence, we use approximation (9) only when Equation 6 is positive, and replace it by 0 otherwise. We have

(10)

The notation [x]₊ means max{0, x}.

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,

(11)

The first and second expressions provide excellent approximations for the rate of tunneling through a deleterious and a neutral mutant, respectively.

	COMPUTER SIMULATIONS

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

In Fig 3 and Fig 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 Equation 6 and Equation 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.

View larger version (16K): In this window In a new window Download PPT slide   Figure 3. Computer simulation and theory for the neutral tunnel, r = 1. The solid circles show the results of direct simulations of the Moran process. We performed 100,000 (top) and 10,000 (bottom) independent runs and show the fraction of runs that result in a homogeneous all 2 population after 1000 cell generations. One cell generation consists of N elementary steps of the Moran process. Cellular population sizes varied from N = 1 to 1000. The lines show the solution of the Kolmogorov forward equation. The solid line uses the tunneling rate as given by Equation 6 and Equation 8. The broken line that is in good agreement uses Equation 11. The other broken line indicates the solution neglecting the tunnel. Parameter values are u1 = 10–5, u1 = 2 x 10–3, r = 1, and a = 10.

View larger version (15K): In this window In a new window Download PPT slide   Figure 4. Computer simulation and theory for tunneling through a disadvantageous intermediate mutant, r < 1. The same as for Fig 3 is shown, but parameter values are u1 = 10–5, u1 = 10–2, r = 0.6, and a = 10.

Fig 3 shows simulations of the neutral tunnel, where r = 1. Solid circles indicate the results of the computer simulation. The solid curve uses Equation 6 and Equation 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 N < 1/.

Fig 4 shows simulations of tunneling through a deleterious intermediate mutant, r < 1. Again, the tunneling rate using Equation 6 and Equation 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₂(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. 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

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

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₁ and u₂, then double mutants might arise at rate u₁u₂. Thus, direct mutation would result in the transition from an all 0 to an all 2 population at rate R_direct = Nu₁u₂(a). In contrast, the tunneling rate is R = Nu₁((r)/(1 – r))u₂(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.

If the intermediate mutant is neutral, r = 1 (Fig 5), then tunneling, , is always faster than direct mutation, R_direct = Nu₁u₂(a), because the mutation rate u₂ is small.

View larger version (59K): In this window In a new window Download PPT slide   Figure 5. Transition rates in compartments. (a) The transition rates for r = 1. In a homogeneous compartment containing only type 0 cells, type 1 cells emerge and might take over the compartment. In this case, the compartment moves from state 0 to state 1 at a rate u1. Here u1 denotes the mutation rate producing a type 1 cell, N denotes the number of cells per compartment, and the probability of fixation of a neutral mutant is given by 1/N. Once the type 1 cell has taken over the compartment, a type 2 cell emerges and might be fixed. Then the compartment moves from state 1 to state 2 at a rate Nu2(a). Here, u2 denotes the mutation rate producing a type 2 cell, and (a) denotes the probability that a type 2 cell with fitness a > 1 reaches fixation in an all 1 compartment. Alternatively, a type 2 cell can be produced before the type 1 cell has taken over the compartment, such that the compartment moves from state 0 directly to state 2 without ever visiting state 1. This stochastic tunneling takes place at rate Nu1[u2(a) – 1/N]+. (b) The transition rates for r < 1. Here the compartment moves from state 0 to state 1 at a rate Nu1(r) and from state 1 to state 2 at a rate Nu2(a/r). The rate of the selected stochastic tunnel is given by Nu1[(ru2/(1 – r))(a) – (r)]+.

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

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

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 , ARMITAGE and DOLL 1957 ; FISHER 1958 ; COOK et al. 1969 ; BELL 1976 ; GOLDIE and COLDMAN 1979 , GOLDIE and COLDMAN 1985 ; MOOLGAVKAR and KNUDSON 1981 ; WHELDON 1988 ; CHAPLAIN 1995 ; BYRNE and CHAPLAIN 1996 ; OWEN and SHERRATT 1999 ; BAR-OR et al. 2000 ; HERRERO-JIMENEZ et 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₂(t) t²Nu₁u₂(a)/2 if N << 1/ and that x₂(t) tNu₁ if N >> 1/. Thus, the probability x₂(t) that the population consists of cells with two inactivated alleles of the tumor suppressor (TSP) gene at time t grows as a quadratic function of time for small population sizes and as a linear function of time for large population sizes. "Small" and "large" are defined by the rate of the second mutation. Chromosomal instability increases the rate of the second mutation (NOWAK et al. 2002 ) and therefore leads to the elimination of TSP genes in only one step also for rather small populations of cells. Note that the explicit equations for x₂(t) hold in the limit of low mutation rates or short timescales. Otherwise the solution of the Kolmogorov forward equation (Equation 1) must be used.

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₂(t) the probability that a population of N cells has made both mutational steps at time t. The probability x₂(t) is given by a convolution integral of two negative exponential distributions. For Nu₂t << 1, we obtain the approximation x₂(t) = t²Nu₁u₂(r)(a)/2. For Nu₂t >> 1, we obtain the approximation x₂(t) = tNu₁(r) as long as N is below a critical value. If instead N is greater than this critical value, we have x₂(t) = tNu₁u₂(r/(1 – r))(a). The critical value is implicitly given by u₂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 – r))ru₂(a), is a product of three factors. The first factor is the product of the population size N and the frequency of mutation-selection balance, u₁/(1 – r), indicating the expected number of type 1 mutants in the population. The second factor, ru₂, 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, . Why it is proportional to the square root of the second mutation rate is not easy to interpret.

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₂(a) when the type 1 mutant is deleterious and is equal to when the type 1 mutant is neutral (IWASA et al. 2003 , IWASA et al. 2004 ). These expressions multiplied by Nu₁ are the terms dominating in Equation 11 when the population size is large. Hence, tunneling is a phenomenon closely related to branching processes.

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 ; SAVILL et 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 rate-limiting steps in cancer progression and the effect of cellular population size and genetic instability.

	ACKNOWLEDGMENTS

We thank Natalia Komarova for many 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.

Manuscript received September 19, 2003; Accepted for publication December 11, 2003.

	APPENDIX 1

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

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

(A1)

The boundary conditions are (0) = 1 and (1) = 0. The explicit solution of Equation A1 can sometimes be found using special functions. For example, when the intermediate type is neutral (r = 1), we can prove

(A2)

where c = N²u₂(a) and is the modified Bessel function of the first kind (GRADSHTEYN and RYZHIK 2001 ). For the purpose of obtaining an approximate formula, we can use a heuristic argument as follows. Since we need to evaluate only (1/N) = V₁, we consider the solution of Equation A1 near y = 0, because we are concerned about the value of y = 1/N (<<1). If we neglect higher-order terms with respect to y, Equation A1 can be expressed as the sum of two exponential functions, (y) = C₊exp[₊y] + C_–exp[_–y], with exponents being the solution of

(A3)

This equation has two roots, one positive and one negative. Since (0) = 1 and (1) = 0, we expect that (y) decreases smoothly from (0) = 1. Then (y) decreases exponentially with the rate specified by the negative root. In effect, we assume (y) = exp[–fy] for small y. Replacing this in Equation A3 and setting y = 0, we have a quadratic equation for f, which can be solved as . This is very small because u₂ is small. Then V₁ in Equation 6 is approximated as

(A4)

	LITERATURE CITED

TOP ABSTRACT MODEL THE RATE OF TUNNELING EXPLICIT FORMULAS FOR THE... COMPUTER SIMULATIONS DIRECT MUTATION DISCUSSION APPENDIX 1 LITERATURE CITED

ARMITAGE, P. and R. DOLL, 1954  The age distribution of cancer and a multistage theory of carcinogenesis. Br. J. Cancer 8:1-12.[Medline]

ARMITAGE, P. and R. DOLL, 1957  A two-stage theory of carcinogenesis in relation to the age distribution of human cancer. Br. J. Cancer 11:161-169.[Medline]

BACH, S. P., A. G. RENEHAN, and C. S. POTTEN, 2000  Stem cells: the intestinal stem cell as a paradigm. Carcinogenesis 21:469-476.[Abstract/Free Full Text]

BAR-OR, R. L., R. MAYA, L. A. SEGEL, U. ALON, and A. J. LEVINE et al., 2000  Generation of oscillations by the p53-Mdm2 feedback loop: a theoretical and experimental study. Proc. Natl. Acad. Sci. USA 97:11250-11255.[Abstract/Free Full Text]

BELL, G. I., 1976  Models of carcinogenesis as an escape from mitotic inhibitors. Science 192:569-572.[Abstract/Free Full Text]

BYRNE, H. M. and M. A. J. CHAPLAIN, 1996  Growth of necrotic tumors in the presence and absence of inhibitors. Math. Biosci. 35:187-216.

CAIRNS, J., 2002  Somatic stem cells and the kinetics of mutagenesis and carcinogenesis. Proc. Natl. Acad. Sci. USA 99:10567-10570.[Abstract/Free Full Text]

CHAPLAIN, M. A. J., 1995  The mathematical modelling of tumour angiogenesis and invasion. Acta Biotheoretica 43:387-402.[CrossRef][Medline]

COOK, P. J., R. DOLL, and S. A. FELLINGHAM, 1969  A mathematical model for the age distribution of cancer in man. Int. J. Cancer 4:93-112.[Medline]

DURRETT, R., 2002 Probability models for DNA sequence evolution. Springer-Verlag, Berlin/Heidelberg, Germany/New York.

FEARON, E. R. and B. VOGELSTEIN, 1990  A genetic model for colorectal tumorigenesis. Cell 61:759-767.[CrossRef][Medline]

FISHER, J. C., 1958  Multiple mutation theory of carcinogenesis. Nature 181:651-652.[Medline]

FRANK, S. A., Y. IWASA, and M. A. NOWAK, 2003  Pattern of cell division and the risk of cancer. Genetics 163:1527-1532.[Abstract/Free Full Text]

GATENBY, R. A. and P. K. MAINI, 2003  Mathematical oncology: cancer summed up. Nature 421:321.[CrossRef][Medline]

GOLDIE, J. H. and A. J. COLDMAN, 1979  A mathematical model for relating the drug sensitivity of tumors to their spontaneous mutation rate. Cancer Treat. Rep. 63:1727-1733.[Medline]

GOLDIE, J. H. and A. J. COLDMAN, 1985  Genetic instability in the development of drug resistance. Semin. Oncol. 12:222-230.[Medline]

GRADSHTEYN, I. S., and I. M. RYZHIK, 2001 Table of Integrals, Series and Products, Ed. 5. Academic Press, New York.

HERRERO-JIMENEZ, P., A. TOMITA-MITCHELL, E. E. FURTH, S. MORGENTHALER, and W. G. THILLY, 2000  Population risk and physiological rate parameters for colon cancer. The union of an explicit model for carcinogenesis with the public health records of the United States. Mutat. Res. 447:73-116.[Medline]

HIGGS, P. G., 1998  Compensatory neutral mutations and the evolution of RNA. Genetica 102(103):91-101.

IIZUKA, M. and M. TAKEFU, 1996  Average time until fixation of mutants with compensatory fitness interaction. Genes Genet. Syst. 71:167-173.

INNAN, H. and W. STEPHAN, 2001  Selection intensity against deleterious mutations in RNA secondary structures and rate of compensatory nucleotide substitutions. Genetics 159:389-399.[Abstract/Free Full Text]

IWASA, Y., F. MICHOR, and M. A. NOWAK, 2003  Evolutionary dynamics of escape from biomedical intervention. Proc. R. Soc. Lond. Ser. B 270:2573-2578.[Medline]

IWASA, Y., F. MICHOR, and M. A. NOWAK, 2004  Evolutionary dynamics of invasion and escape. J. Theor. Biol. 226:205-214.[CrossRef][Medline]

JACKSON, A. L. and L. A. LOEB, 1998  The mutation rate and cancer. Genetics 148:1483-1490.[Abstract/Free Full Text]

KIBERSTIS, P. and J. MARX, 2002  The unstable path to cancer. Science 297:543.[CrossRef]

KIMURA, M., 1985  The role of compensatory neutral mutations in molecular evolution. J. Genet. 64:7-19.

KINZLER, K. W., and B. VOGELSTEIN, 1998 The Genetic Basis of Human Cancer. McGraw-Hill, Toronto.

KNUDSON, A. G., 1971  Mutation and cancer: statistical study of retinoblastoma. Proc. Natl. Acad. Sci. USA 68:820-823.[Abstract/Free Full Text]

KNUDSON, A. G., 2001  Two genetic hits to cancer. Nat. Rev. Cancer 1:157.[CrossRef][Medline]

KOMAROVA, N. L., A. SENGUPTA, and M. A. NOWAK, 2003a  Mutation-selection networks of cancer initiation: tumor suppressor genes and chromosomal instability. J. Theor. Biol. 223:433-450.[CrossRef][Medline]

KOMAROVA, N. L., C. LENGAUER, B. VOGELSTEIN, and M. A. NOWAK, 2003b  Dynamics of genetic instability in sporadic and familial colorectal cancer. Cancer Biol. Ther. 6:685-692.

KOVACS, L. and C. S. POTTEN, 1973  An estimation of proliferative population size in stomach, jejunum and colon of DBA-2 mice. Cell Tissue Kinet. 6:125-134.[Medline]

LENGAUER, C., K. W. KINZLER, and B. VOGELSTEIN, 1998  Genetic instabilities in human cancers. Nature 396:643-649.[CrossRef][Medline]

LITTLE, M. P. and E. G. WRIGHT, 2003  A stochastic carcinogenesis model incorporating genomic instability fitted to colon cancer data. Math. Biosci. 183:111-134.[CrossRef][Medline]

LOEB, L. A., 1998  Cancer cells exhibit a mutator phenotype. Adv. Cancer Res. 72:25-56.[Medline]

LOEB, L. A., C. F. SPRINGGATE, and N. BATTULA, 1974  Errors in DNA replication as a basis of malignant changes. Cancer Res. 34:2311-2321.[Abstract/Free Full Text]

LUEBECK, E. G. and S. H. MOOLGAVKAR, 2002  Multistage carcinogenesis and the incidence of colorectal cancer. Proc. Natl. Acad. Sci. USA 99:15095-15100.[Abstract/Free Full Text]

MICHALAKIS, Y. and M. SLATKIN, 1996  Interaction of selection and recombination in the fixation of negative-epistatic genes. Genet. Res. 67:257-269.[Medline]

MICHOR, F., M. A. NOWAK, S. A. FRANK, and Y. IWASA, 2003a  Stochastic elimination of cancer cells. Proc. R. Soc. Lond. 270:2017-2024.[Medline]

MICHOR, F., Y. IWASA, N. L. KOMAROVA, and M. A. NOWAK, 2003b  Local regulation of homeostasis favors chromosomal instability. Curr. Biol. 13:581-584.[CrossRef][Medline]

MINTZ, B., 1971  Clonal basis of mammalian differentiation. Symp. Soc. Exp. Biol. 25:45-70.

MOOLGAVKAR, S. H. and A. G. KNUDSON, 1981  Mutation and cancer: a model for human carcinogenesis. J. Natl. Cancer Inst. 66:1037-1052.

MULLER, H. J., 1951  Radiation damage to the genetic material. Sci. Progr. 7:93-165.

NOWAK, M. A., N. L. KOMAROVA, A. SENGUPTA, P. V. JALLEPALLI, and I. M. SHIN et al., 2002  The role of chromosomal instability in tumor initiation. Proc. Natl. Acad. Sci. USA 99:16226-16231.[Abstract/Free Full Text]

NORDLING, C. O., 1953  A new theory on the cancer inducing mechanism. Br. J. Cancer 7:68-72.[Medline]

OWEN, M. R. and J. A. SHERRATT, 1999  Mathematical modelling of macrophage dynamics in tumours. Math. Models Methods Appl. Sci. 9:513-539.[CrossRef]

PHILLIPS, P. C., 1996  Waiting for a compensatory mutation: phase zero of the shifting-balance process. Genet. Res. 67:271-284.[Medline]

PLOTKIN, J. B. and M. A. NOWAK, 2002  The different effects of apoptosis and DNA repair on tumorigenesis. J. Theor. Biol. 214:453-467.[CrossRef][Medline]

RAGAJAPOLAN, H., M. A. NOWAK, B. VOGELSTEIN, and C. LENGAUER, 2003  The significance of unstable chromosomes in colorectal cancer. Nat. Rev. Cancer 3:695-701.[CrossRef][Medline]

SAVILL, N. J., D. C. HOYLE, and P. G. HIGGS, 2001  RNA sequence evolution with secondary structure constraints: comparison of substitution rate models using maximum-likelihood methods. Genetics 157:399-411.[Abstract/Free Full Text]

SIEBER, O. M., K. HEINMANN, P. GORMANN, H. LAMLUM, and M. CRABTREE et al., 2002  Analysis of chromosomal instability in human colorectal adenomas with two mutational hits at APC. Proc. Natl. Acad. Sci. USA 99:16910-16915.[Abstract/Free Full Text]

STEPHAN, W., 1996  The rate of compensatory evolution. Genetics 144:419-426.[Abstract]

WHELDON, T. E., 1988 Mathematical Models in Cancer Research. IOP Publishing, Bristol, UK.

WODARZ, D. and D. KRAKAUER, 2001  Genetic instability and the evolution of angiogenic tumor cell lines. Oncol. Rep. 8:1195-1201.[Medline]

This article has been cited by other articles:

				R. Durrett and D. Schmidt Waiting for Two Mutations: With Applications to Regulatory Sequence Evolution and the Limits of Darwinian Evolution Genetics, November 1, 2008; 180(3): 1501 - 1509. [Abstract] [Full Text] [PDF]

				Y. Zhao and R. J. Epstein Programmed Genetic Instability: A Tumor-Permissive Mechanism for Maintaining the Evolvability of Higher Species through Methylation-Dependent Mutation of DNA Repair Genes in the Male Germ Line Mol. Biol. Evol., August 1, 2008; 25(8): 1737 - 1749. [Abstract] [Full Text] [PDF]

				M. M. Desai, D. Weissman, and M. W. Feldman Evolution Can Favor Antagonistic Epistasis Genetics, October 1, 2007; 177(2): 1001 - 1010. [Abstract] [Full Text] [PDF]

				K. Jain and J. Krug Deterministic and Stochastic Regimes of Asexual Evolution on Rugged Fitness Landscapes Genetics, March 1, 2007; 175(3): 1275 - 1288. [Abstract] [Full Text] [PDF]

				J. H. Bielas, K. R. Loeb, B. P. Rubin, L. D. True, and L. A. Loeb From the Cover: Human cancers express a mutator phenotype PNAS, November 28, 2006; 103(48): 18238 - 18242. [Abstract] [Full Text] [PDF]

				F. Michor, Y. Iwasa, and M. A. Nowak The age incidence of chronic myeloid leukemia can be explained by a one-mutation model PNAS, October 3, 2006; 103(40): 14931 - 14934. [Abstract] [Full Text] [PDF]