## Abstract

Results of Nowak and collaborators concerning the onset of cancer due to the inactivation of tumor suppressor genes give the distribution of the time until some individual in a population has experienced two prespecified mutations and the time until this mutant phenotype becomes fixed in the population. In this article we apply these results to obtain insights into regulatory sequence evolution in Drosophila and humans. In particular, we examine the waiting time for a pair of mutations, the first of which inactivates an existing transcription factor binding site and the second of which creates a new one. Consistent with recent experimental observations for Drosophila, we find that a few million years is sufficient, but for humans with a much smaller effective population size, this type of change would take >100 million years. In addition, we use these results to expose flaws in some of Michael Behe's arguments concerning mathematical limits to Darwinian evolution.

THERE is a growing body of experimental evidence that in Drosophila, significant changes in gene regulation can occur in a short amount of time, compared to divergence time between species. Ludwig *et al.* (1998, 2000, 2005) studied the evolution of the *even-skipped* stripe 2 enhancer in four Drosophila species (*Drosophila melanogaster*, *D. yakuba*, *D. erecta*, and *D. pseuodoobscura*). While expression is strongly conserved, they found many substitutions in the binding sites for bicoid, hunchback, Kruppel, and giant, as well as large differences in the overall size of the enhancer region. In addition, they uncovered several binding sites that have been gained and lost among these four species: a lineage-specific addition of the bicoid-3 binding site in *D. melanogaster* that is absent in the other species, a lineage-specific loss of the hunchback-1 site in *D. yakuba*, and the presence of an extra Kruppel site in *D. pseudoobscura* relative to *D. melanogaster*. These differences are nicely summarized in Figure 2B of Ludwig *et al.* (2005).

In a simulation study, Stone and Wray (2001) estimated the rate of *de novo* generation of regulatory sequences from a random genetic background. They found that for a given six-nucleotide sequence, the time until it arose in a 2-kb region in some individual was ∼5950 years for humans and 24 years for Drosophila. However, as MacArthur and Brookfield (2004) have already pointed out, there is a serious problem with Stone and Wray's computation. They assumed individuals in the population evolve independently, while in reality there are significant correlations due to their common ancestors.

Motivated by this simulation study, Durrett and Schmidt (2007) have recently given a mathematical analysis for regulatory sequence evolution in humans, correcting the calculation mentioned above. They assumed an effective population size of 10,000 and a per nucleotide mutation rate of μ = 10^{−8}. In this situation, the expected number of segregating sites in a 1-kb sequence is 1000(4*N*_{e}μ) = 0.4 so it makes sense to talk about a population consensus sequence. The authors defined this as the nucleotide at the site if there is no variability in the population and if the site is variable, the most frequent nucleotide at that site in the population. Using a generation time of 25 years, they found that in a 1-kb region, the average waiting time for words of length six was 100,000 years. For words of length eight, they found that the average waiting time was 375,000 years when there was a seven- of eight-letter match to the target word in the population consensus sequence (an event of probability ∼5/16) and 650 million years when there was not.

Fortunately, in biological reality, the match of a regulatory protein to the target sequence does not have to be exact for binding to occur. Biological reality is complicated, with the acceptable sequences for binding described by position weight matrices that indicate the flexibility at different points in the sequence. To simplify, we assume that binding will occur to any eight-letter word that has seven letters in common with the target word. If we do this, then the mean waiting time reduces to ∼60,000 years.

To explain the intuition behind the last result, consider the case of eight-letter words. If all nucleotides are equally likely and independent, then using the binomial distribution we see that a six of eight match to a given eight-letter target word has probability 63/16,384 ≈ 0.00385, so in a region of 1000 nucleotides, we expect to find 3.85 such approximate matches in the population consensus sequence. Simple calculations then show that the waiting time to improve one of these six of eight matches to seven of eight has a mean of 60,000 years. This shows that new regulatory sequences can come from small modifications of existing sequence.

Extending our previous work on the *de novo* generation of binding sites, this article considers the possibility that in a short amount of time, two changes will occur, the first of which inactivates an existing binding site, and the second of which creates a new one. This problem was studied earlier by Carter and Wagner (2002). In the next section, we present the model and then a simpler theoretical analysis based on work of Komarova *et al.* (2003) and Iwasa *et al.* (2004, 2005), who studied the onset of cancer due to the inactivation of tumor suppressor genes. Finally, we compare the theory with simulations and experimental results.

## THE MODEL

Consider a population of 2*N* haploid individuals. The reader should think of this as the chromosomes of *N* diploid individuals evolving under the assumptions of random union of gametes and additive fitness. However, since we use the continuous-time Moran model, it is simpler and clearer to state our results for haploid individuals.

We start with a homogeneous population of wild-type individuals. We have two sets of possible mutant genotypes *A* and *B*. Wild-type individuals mutate to type *A* at rate *u*_{1} and type *A* individuals mutate to type *B* at rate *u*_{2}. We assume there is no back mutation. We think of the *A* mutation as damaging an existing transcription factor binding site and the *B* mutation as creating a second new binding site at a different location within the regulatory region. We assign relative fitnesses 1, *r*, and *s* to wild-type, *A* mutant, and *B* mutant individuals, respectively. See Figure 1 for a diagram of our model. We used the word damage above to indicate that the mutation may only reduce the binding efficiency, not destroy the binding site. However, even if it does, the mutation need not be lethal. In most cases the *B* mutation will occur when the number of *A* mutants is a small fraction of the population, so most individuals with the *A* mutation will also carry a working copy of the binding site.

We could also assume that the mutations occur in the other order: *B* first and then *A*. This is also a two-stage process that falls into the general framework of our analysis below under the appropriate fitness assumptions. The problems in population genetics to be solved are as follows: How long do we have to wait (i) for a type *B* mutation to arise in some individual or (ii) for the *B* mutant to become fixed in the population? These problems were investigated by Komarova *et al.* (2003) and Iwasa *et al.* (2004, 2005) for tumor suppressor genes, whose inactivation can lead to cancer, with *A* being the inactivation of one copy of the gene and *B* the inactivation of the other. A nice account of these results can be found in Section 12.4 of Nowak's (2006) book on evolutionary dynamics. Here, we apply these results to estimate the waiting time for a switch between two transcription factor binding sites, as defined in the statement of our problem above.

First we need to describe the population genetics model we are considering. Rather than use the discrete-time Wright–Fisher model, we use the continuous-time Moran model. We prefer the Moran process because it is a birth-and-death chain, which means that the number of type *A* individuals increases or decreases by one on each event. Biologically the Moran model corresponds to a population with overlapping generations and in the case of tumor suppressor genes is appropriate for a collection of cells in an organ that is being maintained at a constant size.

As the reader will see from the definition, the Moran process as a genetic model treats *N* diploids as 2*N* haploids and replaces one chromosome at a time. In this context it is common to invoke random union of gametes and assume fitnesses are additive, but that is not necessary. Since homozygous mutants are rare, the fitness of an *A* mutant is its fitness in the heterozygous state. Supposing that the relative fitnesses have been normalized to have maximum value 1, the dynamics may be described as follows:

Each individual is subject to possible replacement at rate 1.

A copy is made of an individual chosen at random from the population.

Mutation changes the copy from wild type to

*A*with probability*u*_{1}and from*A*to*B*with probability*u*_{2}.If the relative fitness of the proposed new individual is 1 −

*q*after mutation (where*q*is the selection coefficient), then the replacement occurs with probability 1 −*q*. Otherwise, nothing happens.

For more on this model, see Section 3.4 of Ewens (2004).

## THEORETICAL RESULTS

#### Neutral case:

Returning to the problem, we first consider the case in which the fitness of the *A* mutant *r* = 1 and the population is of intermediate size as compared to the mutation rates: that is, one in which the population size and mutation rates satisfy(1)where for any numbers *a* and *b*, is read “*a* is much less than *b*” and means *a*/*b* is small.

Theorem 1. *If* *and* , *the probability P*(*t*) *that a B mutation has occurred in some member of the population by time t is*(2)*where* ≈ *is read* “*approximately.*” *If B mutants become fixed with probability* β, *then the result for the fixation time is obtained by replacing u*_{2} *by* β*u*_{2} *in* (2).

In words, the waiting time τ_{B} for the first *B* mutation is roughly exponential with mean while the waiting time *T _{B}* for

*B*to become fixed is roughly exponential with mean . Hence, if β < 1, this increases the waiting time by a factor of rather than the 1/β that one might naively expect. The last conclusion in the theorem should be intuitive since successful mutations (

*i.e.*, those that go to fixation) occur at rate β

*u*

_{2}. In each of the next three theorems, the results for the fixation time can be obtained by replacing

*u*

_{2}by β

*u*

_{2}in the waiting-time result.

*Sketch of proof*. The mathematical proof of this result involves some technical complications, but the underlying ideas are simple. Here and in what follows, readers not interested in the underlying theory can skip the proof sketches. A simple calculation, see Section 2 of Iwasa *et al.* (2005), shows that the probability a type *A* mutant will give rise to a type *B* mutant before its family dies out is asymptotically . Since type *A* mutations arise at rate 2*Nu*_{1}, the time σ_{B} until an *A* mutant arises that will have a descendant of type *B* is exponential with mean . The amount of time after σ_{B} it takes for the *B* mutant to appear, τ_{B} – σ_{B}, is of order . Since , the difference τ_{B} – σ_{B} is much smaller than σ_{B} and the result holds for τ_{B} as well.

To give some intuition about how the *B* mutation arises, we note that in the neutral case, *r* = 1, the number of copies of the *A* allele increases and decreases with equal probability, so if we ignore the transitions that do not change the number of mutant alleles, the result is an unbiased random walk. Since such a walk represents the winnings of a gambler playing a fair game, the probability that the number will rise to before hitting 0 is 1/*k*. When this occurs, the central limit theorem of probability theory implies that the number of steps required to return to 0 is of order *k*^{2} = 1/*u*_{2}, since this requires the random walk to move by *k*, and by the central limit theorem this will take time of order *k*^{2}. Since *B* mutants have probability *u*_{2}, there is a reasonable chance of having a *B* mutation before the number of *A* mutants returns to 0. ▪

Iwasa *et al.* (2004) call this *stochastic tunneling*, since the second mutant (type *B*) arises before the first one (type *A*) fixes. Carter and Wagner (2002) had earlier noted this possibility but they did not end up with a very nice formula for the average fixation time; see their (2.2) and the formulas for the constants given in their Appendix. The assumption implies that throughout the scenario we have just described, the number of type *A* mutants is a small fraction of the population, so we can ignore the probability that the *A* mutants become fixed in the population. This means that in an intermediate-sized population (1) with *r* = 1, *B* mutations arise primarily through stochastic tunneling.

In contrast to populations of intermediate size, populations that are small (compared to the mutation rates) have fixation of the type *A* mutation before the type *B* mutation arises.

Theorem 2. *If* *and* , *then the probability P*(*t*) *that a B mutation has occurred in some member of the population by time t is*(3)

*Sketch of proof*. To explain this, we note that the waiting times between *A* mutations are exponential with mean 1/(2*Nu*_{1}) and each one leads to fixation with probability 1/(2*N*), so the time we have to wait for the first *A* mutation that will go to fixation is exponential with mean 1/*u*_{1}. The condition implies that it is unlikely for the *B* mutation to appear before *A* reaches fixation. The average time required for an *A* mutation to reach fixation conditional on fixation, which is 2*N* by a result of Kimura and Ohta (1969), and the average time required for the *B* mutation to appear after fixation of the *A* mutation, which is 1/(2*Nu*_{2}), are each short compared to 1/*u*_{1} and can be ignored. ▪

When 2*N* and are about the same size, fixation of an *A* mutation and stochastic tunneling are both possible situations in which a type *B* mutation can arise, and the analysis becomes very complicated. Iwasa *et al.* (2005) obtained some partial results; see their Equation 15. Recent work of Durrett *et al.* (2008) addresses this borderline case and gives the following result.

Theorem 3. *If* *and* , *then the probability P*(*t*) *that a B mutation has occurred in some member of the population by time t is P*(*t*) ≈ 1 – exp(– α(γ)*u*_{1}*t*), *where*(4)Hence, the mean waiting time in this case is 1/(α(γ)*u*_{1}).

To summarize the first three theorems, if , then the waiting time for the first *B* mutant to appear in the population, τ_{B}, is approximately exponential under the following conditions:

#### Deleterious A mutants:

Suppose now that *A* mutants have fitness *r* < 1, and return to the case of intermediate population size defined by (1).

Theorem 4. *Suppose* *and that* , *where* ρ *is a constant that measures the strength of selection against type A mutants. The probability P*(*t*) *that the B mutation has occurred in some member of the population by time t is*(5)

The proof of this is somewhat involved so we refer the reader to Iwasa *et al.* (2005) for details. In words, the conclusion says that the waiting time in the nonneutral case is still exponential, but the mean has been multiplied by 1/*R*. Note that when ρ = 0, *R* = 1, which is the neutral case. Thus, *A* mutants are essentially neutral when ρ ≈ 0, which is true when . When ρ is large, [since (ρ + 2/ρ)^{2} = ρ^{2} + 4 + 4/ρ^{2}] and we have 1/*R* ≈ ρ. Therefore, as ρ increases the waiting time increases.

Kimura (1985) considered compensatory mutations that are related to the situation studied in Theorem 4. His model has four genotypes *AB*, *A*′*B*, *AB*′, and *A*′*B*′, where *A* and *B* are wild-type alleles with corresponding mutant alleles *A*′ and *B*′. The single mutant genotypes *A*′*B* and *AB*′ have fitness 1 – *s* while *AB* and *A*′*B*′ have fitness 1. Assuming , the mutation rate, he used diffusion theory to conclude that the average time for the fixation of the double mutant was, see his (16) and (17) and take *h* = 0,where *S* = 4*N*_{e}*s* and *V* = 4*N*_{e}*v*. Evaluating the expression above numerically, he concluded that the fixation time was surprisingly short. Note that his result covers a different range of parameters since Theorem 4 supposes . However, stochastic tunneling still occurs. Kimura shows that the frequency of single mutants remains small until the second mutation occurs.

## SIMULATION RESULTS

The results in the previous section are theorems about the limit as *N* → ∞, and their proofs are based on arguing that various complications can be ignored, so we now turn to simulations to show that the approximations are good for even relatively small values of *N*. We use a standard algorithm, described in the next paragraph, to simulate the continuous-time Markov chain *X*(*t*) that counts the number of *A* mutants in the population at time *t*. Readers not interested in the details of our simulation algorithm can skip the next paragraph.

Let *T*_{0} = 0 and for *m* ≥ 1 let *T _{m}* be the time of the

*m*th jump of

*X*(

*t*). If

*X*(

*T*) = 0, we let

_{m}*t*

_{m}_{+1}be exponential with mean 1/(2

*Nu*

_{1}) and set

*T*

_{m}_{+1}=

*T*+

_{m}*t*

_{m}_{+1}and

*X*(

*T*

_{m}_{+1}) = 1. If

*X*(

*T*) =

_{m}*k*with 1 ≤

*k*< 2

*N*, then we let

*t*

_{m}_{+1}be exponential with mean 1/(

*p*+

_{k}*q*+

_{k}*r*), where

_{k}*p*is the rate of jumps to

_{k}*k*+ 1,

*q*is the rate of jumps to

_{k}*k*– 1, and

*r*is the rate an

_{k}*A*mutant replaces an

*A*mutant, as defined in the following equations. Note that the second term in

*p*accounts for new

_{k}*A*mutants that enter the population:We set

*X*(

*T*

_{m}_{+1}) =

*k*+ 1 with probability

*p*/(

_{k}*p*+

_{k}*q*+

_{k}*r*),

_{k}*X*(

*T*

_{m}_{+1}) =

*k*with probability

*r*/(

_{k}*p*+

_{k}*q*+

_{k}*r*), and

_{k}*X*(

*T*

_{m}_{+1}) =

*k*– 1 with probability

*q*/(

_{k}*p*+

_{k}*q*+

_{k}*r*). In the first two cases there is a probability

_{k}*u*

_{2}of a

*B*mutation. We stop the simulation the first time a

*B*mutant appears or

*X*(

*T*) = 2

_{m}*N*. If an

*A*mutant goes to fixation, we add an exponential with mean 1/(2

*Nu*

_{2}) to the final time to simulate waiting for the

*B*mutation to appear.

Let *n* = 2*N*. Since our aim is to show that the theoretical predictions work well even for small values of *n*, we will, in most cases, consider the values *n* = 1000 and *n* = 10,000. Table 1 gives the seven simulation scenarios we study. Table 2 compares the predicted mean time from Theorem 1 with the average time found in 10,000 replications of each simulation. In making predictions for the examples, we consider that any numbers *a* and *b* satisfy if *a*/*b* ≤ 1/10. In case 3 our assumption (Equation 1) about intermediate population size holds, and we can see that the simulated mean is very close to the predicted mean. In cases 1 and 5, we replace the upper and lower inequalities in (1) by equality, respectively, so in each case one of the two assumptions is not valid. Cases 2 and 4 are intermediate, meaning that the upper and lower inequalities in (1), respectively, do not quite hold since the ratios are 1/4. Yet, cases 2 and 4 show good agreement with the predicted mean.

The last two cases are specific examples related to regulatory sequence evolution in Drosophila and humans, which we consider in more detail in the next section. The Drosophila effective population size is too large to use the true value in the simulations, but this is feasible for humans. In addition, we multiply *u*_{2} by for these special cases so that the ratio *u*_{1}/*u*_{2} = 30. This comes from our assumption that a mutation at any position in the binding site will damage it, but to create a new binding site we require one position to mutate to the correct letter. Given a 10-letter target binding site, then *u*_{1} is 3 × 10 times bigger than *u*_{2}.

In Figure 2, we plot the observed waiting time/predicted mean for case 3 and see a good fit to the exponential distribution, which agrees with our theoretical prediction. Figure 3 corresponds to case 1 and shows that the tail of the distribution looks exponential, but the simulated mean time is ∼1.5 times larger than the predicted mean. This is caused by the fact that since *u*_{1} = 1/(2*N*), the time τ_{B} – σ_{B} we have to wait for the *B* mutant to be produced is of the same order of magnitude as σ_{B}. Hence, the total waiting time τ_{B} is significantly larger than σ_{B}.

To explain the observed shape of the distribution, recall from the sketch of the proof of Theorem 1 that σ_{B} has exactly an exponential distribution. Adding the independent random variable τ_{B} – σ_{B}, which we assume has density *g*(*s*), yields the following distribution for τ_{B},since the integrand is close to *g*(0) for all *s* ∈ [0, *t*], and consequently the exponential fit is not good for small *t*. Wodarz and Komarova (2005) have done an exact calculation of the waiting time in the branching process approximation of the Moran model, which applies to case 1. As Figure 3 shows, the computation matches the Moran model simulation exceptionally well.

Figure 4 corresponds to case 5 and yields a simulated mean of ∼78% of its predicted value. The curve looks exponential, but it has the incorrect mean. In this case, Theorem 3 shows that fixation of an *A* mutation and stochastic tunneling are both possible scenarios in which a *B* mutation can arise, producing a shorter waiting time. More specifically, the assumptions of Theorem 3 hold with γ = 1, which gives α = 1.433, and Figure 4 shows that the exponential distribution with this rate fits the simulated data reasonably well.

## EXPERIMENTAL RESULTS

In the following two examples, we consider a 10-nucleotide binding site and suppose that transcription factor binding requires an exact match to its target. We assume that any mutation within the binding site will damage it (*A* mutation) and that at least one 10-nucleotide sequence exists within the regulatory region that can be promoted to a new binding site by one mutation (*B* mutation). Our previous results show, as MacArthur and Brookfield (2004) observed earlier, that the existence of these so-called “presites” is necessary for the evolution of new binding sites on a reasonable timescale (Durrett and Schmidt 2007).

#### Drosophila:

We assume a per nucleotide mutation rate of 10^{−8} per generation, a simplification of the values that can be found in the classic article of Drake *et al.* (1998) and the recent direct measurements of Haag-Liautard *et al.* (2007). If transcription factor binding involves an exact match to a 10-nucleotide target, then inactivating mutations have probability *u*_{1} = 10^{−7} and those that create a new binding site from a 10-letter word that does not match the target in one position have probability *u*_{2} = × 10^{−8}. If the target word is 6–9 nucleotides long or inexact matches are possible, then these numbers may change by a factor of 2 or 3. Such details are not very important here, since our aim is to identify the order of magnitude of the waiting time.

We set the effective population size *N* = 2.5 × 10^{6}, which agrees with the value given on p. 1612 of Thornton and Andolfatto (2006). To apply Theorem 1 we needThe ratio of the left number to the middle number ≈ 1/300, but the ratio of the middle number to the one on the right is 1/2, which says that is not a valid assumption. Ignoring this for a moment, Theorem 1 predicts a mean waiting time ofwhich translates into 3460 years if we assume 10 generations per year.

Since the assumption is not valid, we use our simulation result for case 6, which has a small population size with parameter values of 2*Nu*_{1} = 0.5 and similar to the Drosophila example, to see what sort of error we expect (see Tables 1 and 2). We see that in the simulation, the observed mean is ∼25% higher than the theoretical mean, so adding 25% to the prediction gives a mean waiting time of 4325 years.

A second and more important correction to our prediction is that Theorem 1 assumes that the *A* mutation is neutral and the *B* mutation is strongly advantageous. If we make the conservative assumption that the *B* mutation is neutral, then the fixation probability β = 1/2*N* = 2 × 10^{−7}, and by Theorem 1 the waiting time increases by a factor of to ∼9 million years. If the *B* mutation is mildly advantageous, *i.e.*, *s* – 1 = 10^{−4}, then β ≈ 10^{−4} and the waiting time increases by a factor of only 100 to 400,000 years.

If we assume that *A* mutants have fitness *r* < 1 where , then Theorem 4 implies that the waiting time is not changed, but if , then the waiting time is increased by a factor if ρ is large. If we use the value of 1 – *r* = 10^{−4}, the increase is roughly a factor of 2. From this we see that if both mutations are almost neutral (*i.e.*, relative fitnesses *r* ≈ 1 – 10^{−4} and *s* ≈ 1 + 10^{−4}), then the switch between two transcription factor binding sites can be done in <1 million years. This is consistent with the results for the *even-skipped* stripe 2 enhancer mentioned earlier.

#### Humans:

We now show that two coordinated changes that turn off one regulatory sequence and turn on another without either mutant becoming fixed are unlikely to occur in the human population. We assume a mutation rate of 10^{−8}, again see Drake *et al.* (1998), and an effective population size of *N* = 10^{4} because this makes the nucleotide diversity 4*N*_{e}μ close to the observed value of 0.1%. If we again assume that transcription factor binding involves an exact match to a 10-nucleotide target, then inactivating mutations have probability *u*_{1} = 10^{−7}, and those that create a new binding site from a 10-letter word that does not match the target in one position have probability *u*_{2} = 3.3 × 10^{−9}. For the assumptions of Theorem 1 to be valid we needThe ratio of the middle number to the one on the right is 1/500, but the ratio of the left number to the middle one ≈ 1.

Ignoring for the moment that one of the assumptions is not satisfied, Theorem 1 predicts a mean waiting time ofMultiplying by 25 years per generation gives 216 million years.

As shown in Tables 1 and 2, we have simulation results for humans using the exact parameters above. In 10,000 replications, the simulation mean is 6.46 million generations, which is only ∼75% of the predicted value. Multiplying by 0.75 reduces the mean waiting time to 162 million years, still a very long time. Our previous work has shown that, in humans, a new transcription factor binding site can be created by a single mutation in an average of 60,000 years, but, as our new results show, a coordinated pair of mutations that first inactivates a binding site and then creates a new one is very unlikely to occur on a reasonable timescale.

To be precise, the last argument shows that it takes a long time to wait for two prespecified mutations with the indicated probabilities. The probability of a seven of eight match to a specified eight-letter word is 8(3/4)(1/4)^{7} ≈ 3.7 × 10^{−4}, so in a 1-kb stretch of DNA there is likely to be only one such match. However, Lynch (2007, see p. 805) notes that transcription factor binding sites can be found within a larger regulatory region (10^{4} – 10^{6} bp) in humans. If one can search for the new target sequence in 10^{4} – 10^{6} bp, then there are many more chances. Indeed since (1/4)^{8} ≈ 1.6 × 10^{−5}, then in 10^{6} bp we expect to find 16 copies of the eight-letter word.

#### The edge of evolution?

Our final example of waiting for two mutations concerns the emergence of chloroquine resistance in *P. falciparum*. Genetic studies have shown, see Wooton *et al.* (2002), that this is due to changes in a protein PfCRT and that in the mutant strains two amino acid changes are almost always present—one switch at position 76 and another at position 220. This example plays a key role in the chapter titled “The mathematical limits of Darwinism” in Michael Behe's book, *The Edge of Evolution* (Behe 2007).

Arguing that (i) there are 1 trillion parasitic cells in an infected person, (ii) there are 1 billion infected persons on the planet, and (ii) chloroquine resistance has arisen only 10 times in the past 50 years, he concludes that the odds of one parasite developing resistance to chloroquine, an event he calls a *chloroquine complexity cluster* (CCC), are ∼1 in 10^{20}. Ignoring the fact that humans and *P. falciparum* have different mutation rates, he then concludes that “On the average, for humans to achieve a mutation like this by chance, we would have to wait a hundred million times ten million years” (Behe 2007, p. 61), which is 5 million times larger than the calculation we have just given.

Indeed his error is much worse. To further sensationalize his conclusion, he argues that “There are 5000 species of modern mammals. If each species had an average of a million members, and if a new generation appeared each year, and if this went on for two hundred million years, the likelihood of a single CCC appearing in the whole bunch over that entire time would only be about 1 in 100” (Behe 2007, p. 61). Taking 2*N* = 10^{6} and μ_{1} = μ_{2} = 10^{−9}, Theorem 1 predicts a waiting time of 31.6 million generations for one prespecified pair of mutations in one species, with having reduced the answer by a factor of 31,600.

We are certainly not the first to have criticized Behe's work. Lynch (2005) has written a rebuttal to Behe and Snoke (2004), which is widely cited by proponents of intelligent design (see the Wikipedia entry on Michael Behe). Behe and Snoke (2004) consider evolutionary steps that require changes in two amino acids and argue that to become fixed in 10^{8} generations would require a population size of 10^{9}. One obvious problem with their analysis is that they do their calculations for *N* = 1 individual, ignoring the population genetic effects that produce the factor of . Lynch (2005) also raises other objections.

## CONCLUSIONS

For population sizes and mutation rates appropriate for Drosophila, a pair of mutations can switch off one transcription factor binding site and activate another on a timescale of several million years, even when we make the conservative assumption that the second mutation is neutral. This theoretical result is consistent with the observation of rapid turnover of transcription factor binding sites in Drosophila and gives some insight into how these changes might have happened. Our results show that when two mutations with rates *u*_{1} and *u*_{2} have occurred andthen the first one will not have gone to fixation before the second mutation occurs, and indeed *A* mutants will never be more than a small fraction of the overall population. In this scenario, the *A* mutants with fitness *r* are significantly deleterious if is large, a much less stringent condition than the usual condition that 2*N*(1 – *r*) is large. Also, the success probability of the *B* mutant is dictated by its fitness relative to the wild type rather than relative to the *A* mutant. This follows because the fraction of *A* mutants in the population is small when the *B* mutant arises, and hence most individuals are wild type at that time.

The very simple assumptions we have made about the nature of transcription factor binding and mutation processes are not crucial to our conclusions. Our results can be applied to more accurate models of binding site structure and mutation processes whenever one can estimate the probabilities *u*_{1} and *u*_{2}. However, the assumption of a homogeneously mixing population of constant size is very important for our analysis. One obvious problem is that Drosophila populations undergo large seasonal fluctuations, providing more opportunities for mutation when the population size is large and a greater probability of fixation of an *A* mutation during the recurring bottlenecks. Thus, it is not clear that one can reduce to a constant size population or that the effective population size computed from the nucleotide diversity is the correct number to use for the constant population size. A second problem is that in a subdivided population, *A* mutants may become fixed in one subpopulation, giving more opportunities for the production of *B* mutants or perhaps leading to a speciation event. It is difficult to analyze these situations mathematically, but it seems that each of them would increase the rate at which changes occur. In any case one would need to find a mechanism that changes the answer by a significant factor to alter our qualitative conclusions.

## Acknowledgments

We thank Eric Siggia for introducing us to these problems and for many helpful discussions and Jason Schweinsberg for his collaboration on related theoretical results. We are also grateful to Nadia Singh, Jeff Jensen, Yoav Gilad, and Sean Carroll who commented on previous drafts and to two referees (one anonymous and Michael Lynch) who helped to improve the presentation. Both authors were partially supported by a National Science Foundation/National Institutes of General Medical Sciences grant (DMS 0201037). R.D. is also partially supported by a grant from the probability program (DMS 0202935) at the National Science Foundation. D.S. was partially supported by a National Science Foundation graduate fellowship at Cornell and after graduation was a postdoc at the Institute for Mathematics and its Applications in Minnesota, 2007–2008.

## Footnotes

Communicating editor: J. B. Walsh

- Received September 30, 2007.
- Accepted August 19, 2008.

- Copyright © 2008 by the Genetics Society of America