Stochastic Tunnels in Evolutionary Dynamics
Yoh Iwasa, Franziska Michor, Martin A. Nowak

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).

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.

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 u1, while type 1 mutates to type 2 at rate u2. 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 ρ(r)=11r11rN. 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,b). When the population size N is finite, there is always a positive probability of fixation even for deleterious mutants.

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.

If Nu1 ≪ 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 Nu1ρ(r). The waiting time for the transition follows a negative exponential distribution with mean 1/Nu1ρ(r).

After the fixation of type 1, a similar process describes the transition from an all 1 population to an all 2 population. If Nu2 ≪ 1, we can regard this transition as a Markovian jump occurring at a rate Nu2ρ(a/r). The waiting time for the transition from all 1 to all 2 follows a negative exponential distribution with mean 1/Nu2ρ(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 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 (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 x0, x1, and x2, 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 x.0=Nu1ρ(r)x0Rx0x.1=Nu1ρ(r)x0Nu2ρ(a)x1x.2=Nu2ρ(a)x1+Rx0. (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

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 P=exp[ru2ρ(a)0Y(t)dt]. (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, ru2, 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 R=Nu1Pr[Y()=0]E[1PY(0)=1,Y()=0]. (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 E[PY(0)=1]=E[PY(0)=1,Y()=0]Pr[Y()=0]+E[PY(0)=1,Y()=N]Pr[Y()=N]. If Y becomes fixed, Y(∞) = N, and then the integral 0Y(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 R=Nu1{Pr[Y()=0]E[PY(0)=1]}. (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 Vi=E[PY(0)=i]fori=0,1,2,3,,N. (5) Note that V0 = 1 and VN = 0. The expectation in Equation 4 is equal to V1. We can write R=Nu1{1V1ρ(r)}. (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 (Ni)(ir/(ir + (Ni))). The transition from i to i – 1 occurs at rate i((Ni)/(ir + (Ni))). The system remains in i otherwise. To derive a recursive formula for Vi, we expand the function according to events occurring in a short time interval of length Δt. Vi=exp[ru2ρ(a)iΔt]E[exp[ru2ρ(a)ΔtY(t)dt]Y(0)=i]=(1ru2ρ(a)iΔt)jE[exp[ru2ρ(a)ΔtY(t)dt]Y(Δt)=j]×Pr(ji,inΔt)+O((Δt)2)=(1ru2ρ(a)iΔt){rVi+1+Vi1}i(Ni)ir+(Ni)Δt+(1ru2ρ(a)iΔt)Vi(1(r+1)i(Ni)ir+(Ni)Δt)+O((Δt)2). By rearranging terms, dividing both sides by Δt, and taking the limit when Δt is very small, we obtain 0=ru2ρ(a)Vi+{r(Vi+1Vi)+(Vi1Vi)}Niir+(Ni). (7) This recursion holds for i = 1, 2,..., N – 1. The boundary conditions are V0 = 1 and VN = 0. Using Equation 7, we can determine the solution of V1. We can use an iterative approximation to obtain a numerical solution of Equation 7. For example, we start with the initial distribution, Vi = 1 – i/N for i = 0,..., N, and calculate Vi in the next round by Vi=(rVi+1+Vi1)/(ru2ρ(a)ir+(Ni)Ni+r+1),i=1,2,3,,N1. (8) When this iteration converges, the stationary distribution of Vi satisfies Equation 7. The value for V1 together with Equation 6 gives the rate of tunneling.

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: V1=1(1r)+(1r)2+2(1+r)ru2ρ(a)1+r+[higher-order terms]. (9) As we pointed out earlier, 1 – V1 –ρ(r) in Equation 6 is always positive from the definition of a probability. If V1 is replaced by the simpler formula shown above, then 1 – V1 –ρ(r) can become negative. Hence, we use approximation (9) only when Equation 6 is positive, and replace it by 0 otherwise. We have RNu1[(1r)+(1r)2+2(1+r)ru2ρ(a)1+rρ(r)]+. (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, R={Nu1[r1ru2ρ(a)ρ(r)]+if1r2u2ρ(a)Nu1[u2ρ(a)1N]+if1r2u2ρ(a).} (11) The first and second expressions provide excellent approximations for the rate of tunneling through a deleterious and a neutral mutant, respectively.

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 Equations 6 and 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 × 10–3, r = 1, and a = 10.

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 N<1u2ρ(a) .

Figure 4.

—Computer simulation and theory for tunneling through a disadvantageous intermediate mutant, r < 1. The same as for Figure 3 is shown, but parameter values are u1 = 10–5, u1 = 10–2, r = 0.6, and a = 10.

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 ru2ρ(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, Nu1 ≪ 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 u1 and u2, then double mutants might arise at rate u1u2. Thus, direct mutation would result in the transition from an all 0 to an all 2 population at rate Rdirect = Nu1u2ρ(a). In contrast, the tunneling rate is R = Nu1 ((r)/(1 – r))u2ρ(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.

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)]+.

If the intermediate mutant is neutral, r = 1 (Figure 5), then tunneling, R=Nu1u2ρ(a) , is always faster than direct mutation, Rdirect = Nu1u2ρ(a), because the mutation rate u2 is small.

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; Bar-Oret al. 2000; Herrero-Jimenezet 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 x2(t) ≈ t2Nu1u2ρ(a)/2 if N1u2ρ(a) and that x2(t)tNu1u2ρ(a) if N1u2ρ(a) . Thus, the probability x2(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 (Nowaket 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 x2(t) hold in the limit of low mutation rates or short time-scales. 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 x2(t) the probability that a population of N cells has made both mutational steps at time t. The probability x2(t) is given by a convolution integral of two negative exponential distributions. For Nu2t ≪ 1, we obtain the approximation x2(t) = t2Nu1u2ρ(r)ρ(a)/2. For Nu2t ≫ 1, we obtain the approximation x2(t) = tNu1ρ(r) as long as N is below a critical value. If instead N is greater than this critical value, we have x2(t) = tNu1u2(r/(1 – r))ρ(a). The critical value is implicitly given by u2r/(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(u1/(1 – r))ru2ρ(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, u1/(1 – r), indicating the expected number of type 1 mutants in the population. The second factor, ru2, 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, R=Nu1u2ρ(a) . 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))u2ρ(a) when the type 1 mutant is deleterious and is equal to u2ρ(a) when the type 1 mutant is neutral (Iwasa et al. 2003, 2004). These expressions multiplied by Nu1 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; 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 rate-limiting 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) = Vi and y = i/N. Calculating the Taylor expansion of φ(y), Equation 7 is rewritten as 1+(r1)y1yru2ρ(a)ϕ(y)=1rNdϕdy+1+r2N2d2ϕdy2. (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 ϕ(y)k=0[c(1y)]k+1(k+1)!k!=c(1y)k=0(c(1y))2k+1(k+1)!k!=c(1y)I1(2c(1y)), (A2) where c = N2u2ρ(a) and I1(z)=k=0((z2)2k+1k!(k+1)!) 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) = V1, 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] + Cexp[λy], with exponents being the solution of 1+r2N2λ21rNλru2ρ(a)=0. (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 fN=((1r)+(1r)2+2(1+r)ru2ρ(a))(1+r) . This is very small because u2 is small. Then V1 in Equation 6 is approximated as V1=exp[fN]1(1r)+(1r)2+2(1+r)ru2ρ(a)1+r+[higher-order terms]. (A4)

Footnotes

  • Communicating editor: J. B. Walsh

  • Received September 19, 2003.
  • Accepted December 11, 2003.

LITERATURE CITED

View Abstract