Wright, Fisher, and the Weasel

[latexpage]
Richard Dawkins’s computer simulation algorithm explores how long it takes a 28-letter-long phrase to evolve to become the phrase “Methinks it is like a weasel”. The Weasel program has a single example of the phrase which produces a number of offspring, with each letter subject to mutation, where there are 27 possible letters, the 26 letters A-Z and a space. The offspring that is closest to that target replaces the single parent. The purpose of the program is to show that creationist orators who argue that evolutionary biology explains adaptations by “chance” are misleading their audiences. Pure random mutation without any selection would lead to a random sequence of 28-letter phrases. There are $27^{28}$ possible 28-letter phrases, so it should take about $10^{40}$ different phrases before we found the target. That is without arranging that the phrase that replaces the parent is the one closest to the target. Once that highly nonrandom condition is imposed, the number of generations to success drops dramatically, from $10^{40}$ to mere thousands.

Although Dawkins’s Weasel algorithm is a dramatic success at making clear the difference between pure “chance” and selection, it differs from standard evolutionary models. It has only one haploid adult in each generation, and since the offspring that is most fit is always chosen, the strength of selection is in effect infinite. How does this compare to the standard Wright-Fisher model of theoretical population genetics?

The Wright-Fisher model, developed 85 years ago, independently, by two founders of theoretical population genetics, is somewhat similar. But it has an infinite number of newborns, precisely $N$ of whom survive to adulthood.  In its simplest form, with asexual haploids, the adult survivors each produce an infinitely large, but equal, number of offspring. In that infinitely large offspring population, mutations occur exactly in the expected frequencies. Each genotype has a fitness, and the proportions of the genotypes change according to them, exactly as expected. Since there are an infinite number of offspring, there is no genetic drift at that stage. Then finally a random sample of $N$ of the survivors is chosen to become the adults of this new generation.

In one case we get fairly close to the Weasel. That is when $N = 1$, so that there is only one adult individual. Now for some theory, part of which I explained before, in a comment here. We first assume that a match of one letter to the target multiplies the fitness by a factor of $1+s$, and that this holds no matter how many other letters match. This is the fitness scheme in the genetic algorithm program by keiths. A phrase that has $k$ matches to the target has fitness $(1+s)^k$. We also assume that mutation occurs independently in each letter of each offspring with probability $u$, and produces one of the other 26 possibilities, chosen uniformly at random. This too is what happens in keiths’s program. If a letter matches the target, there is a probability $u$ that it will cease to match after mutation. If it does not match, there is a probability $u/26$ that it will match after mutation.

The interesting property of these assumptions is that, in the Wright-Fisher model with $N = 1$, each site evolves independently because neither its mutational process nor its selection coefficients depend on the the state of the other sites. So we can simply ask about the evolution of one site.

If the site matches, after mutation $1-u$ of the offspring will still match, and $u$ of them will not. Selection results in a ratio of $(1+s)(1-u)\,:\,u$ of matching to nonmatching letters. Dividing by the sum of these, the probability of ending up with a nonmatch is then $u$ divided by the sum, or
\[
q \ = \ u\big/\left((1+s)(1-u)+u\right) \ = \ u \big/ (1 + s(1-u))
\]

Similarly if a letter does not match, after mutation the ratio of nonmatches to matches is $\frac{u}{26}\,:\,1-\frac{u}{26}$. After selection the ratio is $(1+s)\frac{u}{26}\,:\,1-\frac{u}{26}$. Dividing by the sum of these, the probability of ending up with a match is
\[
p \ = \ \frac{u}{26}(1+s) \bigg/\left( \frac{u}{26}(1+s) + \left(1 – \frac{u}{26}\right)\right)
\ = \ u(1+s) \big/ (26 + us)
\]

Thus, what we see at one position in the 28-letter phrase is that if it is not matched, it has probability $p$ of changing to a match, and if it is matched, it has probability $q$ of changing to not being matched. This is a simple 2-state Markov process. I’m going to leave out the further analysis — I can explain any point in the comments. Suffice it to say that

1. The probability, in the long run, of a position matching is $p/(p+q)$.

2. This is independent in each of the 28 positions. So the number of matches will, in the long run, be a Binomial variate with 28 trials each with probability $p/(p+q)$.

3. The “relaxation time” of this system, the number of generations to approach the long-term equilibrium, will be roughly $1/(p+q)$.

4. At what value of $s$ do we expect selection to overcome genetic drift? There is no value where drift does not occur, but we can roughly say that when the equilibrium frequency $p/(p+q)$ is twice as great as the no-selection value of 1/27, that selection is then having a noticeable effect. The equation for this value is a quadratic equation in $s$ whose coefficients involve $u$ as well. If $u$ is much smaller than $s$, then this value of $s$ is close to 0.44.

5. The probability of a match at a particular position, one which started out not matching, will be after $t$ generations
\[
P(t) \ = \ \frac{p}{p+q} (1\:-\:(1-p-q)^t)
\]
and the probability of a nonmatch given that one started out matching is
\[
Q(t) \ = \ \frac{q}{p+q} (1\:-\:(1-p-q)^t)
\]

I don’t have the energy to see whether the runs of keiths’s program verify this, but I hope others will.

291 thoughts on “Wright, Fisher, and the Weasel

  1. keiths: Yes. I see a better match between my N=10 and N=20 runs:

    Generation: 1318300
    Total population size: 5000 Number of survivors per generation: 10
    Selection: ONCoefficient: 0.30 Mutation rate: 0.010
    Mean Hamming distance:13.42 Standard deviation:2.20

    I get mean 13.40 and s.d. 2.19.

    keiths:

    Generation: 1599300
    Total population size: 5000 Number of survivors per generation: 20
    Selection: ONCoefficient: 0.15Mutation rate: 0.005
    Mean Hamming distance:12.76 Standard deviation:2.17

    Mean 12.84 and s.d. 2.25. The entire plot is not very interesting — I’m just dumping in all of the data at present, meaning that there are 10 + 20 = 30 curves — so I’m going with just the early generations. If you look closely, you can distinguish the “heavier” group of 20 parents in the second run from the group of 10 parents in the first run.

    Number of parents: 10
    ************************************************************************
    Number of sites L : 28
    Number of possible alleles A : 27
    Number of parents : 10
    Number of offspring per parent : 500
    Mutation rate u : 0.01
    Selection parameter s : 0.3
    p = u * (1 + s) / (A-1 + u * s) : 0.000499942314348
    q = u / (1 + s * (1 – u)) : 0.0077101002313
    Assume equilibrium at generation: 10000
    W-F probability of fit allele at site P = p/(p+q) : 0.0608939979992
    Expected num. of fit alleles in parent L*P : 1.70503194398
    Mean observed num. of fit alleles in parent : 14.6011827261
    Expected std. dev. of observed nums. of fit alleles: 1.26538758181
    Std. dev. of observed numbers of fit alleles : 2.19382371378
    ************************************************************************
    Number of parents: 20
    ************************************************************************
    Number of sites L : 28
    Number of possible alleles A : 27
    Number of parents : 20
    Number of offspring per parent : 250
    Mutation rate u : 0.005
    Selection parameter s : 0.15
    p = u * (1 + s) / (A-1 + u * s) : 0.0002211474669
    q = u / (1 + s * (1 – u)) : 0.00435066347618
    Assume equilibrium at generation: 10000
    W-F probability of fit allele at site P = p/(p+q) : 0.0483719623697
    Expected num. of fit alleles in parent L*P : 1.35441494635
    Mean observed num. of fit alleles in parent : 15.1601330201
    Expected std. dev. of observed nums. of fit alleles: 1.13529698209
    Std. dev. of observed numbers of fit alleles : 2.24921122137
    ************************************************************************
    CPU times: user 23min 37s, sys: 7.93 s, total: 23min 45s
    Wall time: 23min 50s

  2. Tom:

    Would you please check this?

    Sorry, Tom, but with those parameters the fitnesses would range from 1 to 10,001^10,000, which is far beyond the range of the double-precision fitness variable I use.

    Augmenting the precision doesn’t seem worth it to me.

  3. keiths:

    Mean Hamming distance: 13.42 Standard deviation: 2.20

    Tom:

    I get mean 13.40 and s.d. 2.19.

    keiths:

    Mean Hamming distance: 12.76 Standard deviation: 2.17

    Tom:

    Mean 12.84 and s.d. 2.25.

    Pretty good agreement.

  4. Separate figures now for min-mean-max of two runs. My quick-and-dirty code “likes” to process offspring in big batches, not split into 200 little fragments.

    keiths: And as expected, an even better match between N=100 and N=200:

    Generation: 525400
    Total population size: 5000 Number of survivors per generation: 100
    Selection: ON Coefficient: 0.10 Mutation rate: 0.010
    Mean Hamming distance:11.20 Standard deviation:2.05

    Mean 11.13 and s.d. 2.04 for the last 900,000 of 1,000,000 generations.

    keiths: Generation: 526000
    Total population size: 5000 Number of survivors per generation: 200
    Selection: ON Coefficient: 0.05 Mutation rate: 0.005
    Mean Hamming distance:11.27 Standard deviation:2.03

    Mean and s.d. 2.04.

    Number of parents: 100
    ************************************************************************
    Number of sites L : 28
    Number of possible alleles A : 27
    Number of parents : 100
    Number of offspring per parent : 50
    Mutation rate u : 0.01
    Selection parameter s : 0.1
    Assume equilibrium at generation: 10000
    Mean observed num. of fit alleles in parent : 16.86679268
    Std. dev. of observed numbers of fit alleles : 2.04897128063
    ************************************************************************
    Number of parents: 200
    ************************************************************************
    Number of sites L : 28
    Number of possible alleles A : 27
    Number of parents : 200
    Number of offspring per parent : 25
    Mutation rate u : 0.005
    Selection parameter s : 0.05
    Assume equilibrium at generation: 10000
    Mean observed num. of fit alleles in parent : 16.9037339205
    Std. dev. of observed numbers of fit alleles : 2.03701360181
    ************************************************************************
    CPU times: user 2h 24min 32s, sys: 1min 11s, total: 2h 25min 44s
    Wall time: 2h 26min 11s

  5. keiths: Sorry, Tom, but with those parameters the fitnesses would range from 1 to 10,001^10,000, which is far beyond the range of the double-precision fitness variable I use.

    Augmenting the precision doesn’t seem worth it to me.

    I’ve wondered all along just what it is that you value. I’d genuinely like to know.

  6. Tom,

    I’ve wondered all along just what it is that you value. I’d genuinely like to know.

    I answered that already:

    The modified Weasel achieves its purpose. It shows that even with a tiny effective population size, drift needn’t dominate selection. It also provides a nice visual demonstration of the population being pushed toward equilibrium, where drift and selection cancel each other out on average, then wandering stochastically above and below that equilibrium point.

    It also validates, and is validated by, the predictions of Joe’s Wright-Fisher model.

  7. In one case we get fairly close to the Weasel. That is when N = 1, so that there is only one adult individual.

    Why look any further?

  8. Joe Felsenstein,

    Would you please, now that the Weaseling has died down, recommend one of the many simulations, available on the Web, of the Wright-Fisher model? Or perhaps you would like to share one of your own implementations. (I suspect that you’ve done various of them over the decades, and also that you’ve made assignments of them in your teaching.)

    In 2016, one needs a compelling reason not to run computational experiments within an interactive environment. Have you, in planning future development of PHYLIP, had a look at Jupyter (the offshoot of IPython)? The project recently obtained over $6 million in funding.

  9. Mung:

    [Me, in the Original Post] In one case we get fairly close to the Weasel. That is when N = 1, so that there is only one adult individual.

    Why look any further?

    Mung may have forgotten what the purpose of this thread was: it was to use the Wright-Fisher model to get a good sense of what the effect of moderate amounts of selection is, and what the effect of different numbers of surviving adults is. So we of course need to look further than $N = 1$, and further than having infinitely strong selection.

  10. Tom English,

    I presume you mean multi-locus selection simulations. For a single-locus program using the Wright-Fisher model to simulate natural selection, migration, and mutation, and genetic drift, the interactive program PopG, distributed by our lab here should be good.

    Don’t have time to write using a new language. I have a two-page multilocus selection program, one I used to show a creationist troll that his assertions were wrong. I could post that here and then you all could put a better interface on it.

  11. Let me qualify what I just said. If you had some number of loci, say 10 of them, and they were not linked, so that there was free recombination between them, and if the fitness scheme was multiplicative, so that the fitness of a genotype that had $k$ favorable alleles was $(1+s)^k$, then each at each locus the gene frequency would change essentially indepndently.

    So a single-locus simulation program like PopG will show what happens. 10 loci would act essentialy the same as 10 independent runs of PopG, or one run of 10 independent populations. There doesn’t seem to be much point in posting a 10-locus program.

  12. Joe Felsenstein: I presume you mean multi-locus selection simulations.

    Yes, I do. Some days back, I Googled Wright-Fisher model simulation, and got about 456,000 results. I decided not to wade through what was available, and instead to demonstrate just how little is required to do the the particular job at hand. (OK, I did just sample what’s out there, and found that I can’t tell for sure what would have worked for us.)

    Joe Felsenstein: I could post that here and then you all could put a better interface on it.

    I’ve always worked on guts, not skin. I’m interested now in learning some leading-edge techniques for presentation online. That’s bound to take some time.

  13. I looked at the first page of the Google search. Some of the papers were about teaching programs — one was an assignment to a class asking them to write the program. There was an interesting (if flawed) paper on simulating conditioned WF processes.

    All were one-locus programs.

    I am not sure whether you are asking me to post my multilocus program, or to let you write your own. I will post it (in a comment a few minutes from now) so as to let anyone who wants to play with it do so. As I indicated, the results will be little different from what you would get running PopG once for each locus.

  14. Delayes in posting a multilocus program. I found the one that I had and discovered that it was both (a) a modification of a program someone else had sent me, and (b) didn’t really do selection for favorable alleles at all loci.

    I will think about doing a completely new version. But for now folks should use PopG for one locus at a time.

  15. Tom:

    Some days back, I Googled Wright-Fisher model simulation, and got about 456,000 results. I decided not to wade through what was available, and instead to demonstrate just how little is required to do the the particular job at hand.

    Ah, now I’m beginning to see why Tom is so grumpy regarding the modified Weasel, and what’s behind comments like this:

    But I’m infinitely more interested in general theoretical results than in whacking my Weasel in public.

    He’s got a bad case of Weasel envy. 🙂

    P.S. No one said it was difficult, Tom. I’m sure you could do a fine full-featured Weasel, long and furry.

  16. Joe Felsenstein: I am not sure whether you are asking me to post my multilocus program, or to let you write your own.

    I’ve got something else I’d rather address right now. When you first contacted me, way back when, I was working on a response to a Weasely paper by a Templeton Foundation grant recipient:

    David H. Glass, University of Ulster, “Parameter Dependence in Cumulative Selection,” 2014.

    His Templeton project was “Explaining and Explaining Away.”

    You’ve set things up beautifully with this post. Glass claims that there’s biological relevance in his study. I think you’ll want to say a thing or two about that. It’s a very brief paper. I’m about to start an OP.

    I’d love to have keiths join in, but he’s pretty busy with real work, keeping the likes of kairosfocus and Mung in check.

  17. Hmm. Northern Ireland does seem to be an epicenter of creationism.

    Glass seems surprised that if you have half as much mutation, it can take twice as long to reach the target.

    Will await your OP.

  18. Tom English: I’d love to have keiths join in, but he’s pretty busy with real work, keeping the likes of kairosfocus and Mung in check.

    It’s that insidious plot to force his children to hear about intelligent design in science classes. It’s becoming more and more likely to become a reality as each day passes.

    Can you just imagine the damage if they taught Weasel Testing in CS courses?

  19. Tom,

    Responding to Glass might require us to “whack our Weasels in public”, as you put it. Are you up for that?

    I’d love to have keiths join in, but he’s pretty busy with real work, keeping the likes of kairosfocus and Mung in check.

    Glass’s misconceptions have a lot in common with those of KF, Mung, and Sal. His concern with a separate “CCM” (correct character mutation rate) mirrors KF’s “latching” obsession. And like Mung and Sal, he wants to overstate the problem that drift presents to selection.

  20. keiths: Glass’s misconceptions have a lot in common with those of KF, Mung, and Sal. His concern with a separate “CCM” (correct character mutation rate) mirrors KF’s “latching” obsession.

    The equations that I gave, above, allow there to be a separate Correct Character Mutation rate. Just replace the mutation rate from Match to Not, $u/26$ by that rate. The alphabet size only comes into my math in that term.

  21. Mung: It’s that insidious plot to force his children to hear about intelligent design in science classes. It’s becoming more and more likely to become a reality as each day passes.

    Can you just imagine the damage if they taught Weasel Testing in CS courses?

    It is evident from Glass’s paper that there is a lack of teaching of population genetics in the science courses he took.

  22. Joe Felsenstein: The equations that I gave, above, allow there to be a separate Correct Character Mutation rate.Just replace the mutation rate from Match to Not, $u/26$, by

    Woops. Reverse that. Make it: Replace the mutation rate from Nonmatch to Match, $u/26$, …

  23. keiths: Responding to Glass might require us to “whack our Weasels in public”, as you put it. Are you up for that?

    Not really. I’m very lazy in documenting code that I plan to use only once. And I don’t keep notes. I find now that I wrote, 16 months ago, an alternative main program that tests the routines I’m using. So I have to figure out what the testing actually was, before I can decide whether I’m actually satisfied with it.

    I did leave myself one interesting note, and it’s one that you might want to ponder. The speedup is relative to a “slow_mutate” routine that is much faster than what you’re doing. By the way, you need to get rid of your sort, if you haven’t done so already.

  24. Tom,

    I did leave myself one interesting note, and it’s one that you might want to ponder. The speedup is relative to a “slow_mutate” routine that is much faster than what you’re doing.

    There are lots of things I would have changed if my goal were to minimize execution time, but given that my program typically reaches equilibrium in a matter of seconds, I haven’t seen the need.

    By the way, you need to get rid of your sort, if you haven’t done so already.

    I left that code there in case I later wanted to implement a pure Weasel option where the fittest offspring was always selected. I’ve since decided that it’s better to maintain separate versions — a pure Weasel and a modified drift/selection Weasel.

  25. keiths: I left that code there in case I later wanted to implement a pure Weasel option where the fittest offspring was always selected.

    The heap comes early in sophomore-level courses on data structures. Should you take the radical leap to C++, you’ll have access to standard algorithms for heap manipulation. It’s really no harder to do selection right than to make excuses for doing it wrong — unless you’ve invested a lot of ego in what you’ve done wrong.

    keiths: I haven’t seen the need.

    I’ve always seen the need, because I see the Weasel as a special (very limited) case of evolutionary computation, not evolutionary computation as a generalization of the Weasel. I see pitfalls galore in seeing things the latter way. I’ll of course let Joe speak to analogous concerns, if any, about theoretical evolutionary genetics.

    I’ve previously had the need for arbitrary-precision calculations in evolutionary computation, and I saw immediately the need for them in a simulation with multiplicative fitness. An “obvious” (only after reading Joe’s post) response to Glass is to redo his study with multiplicative fitness contributions. So what do you think now about calculating $(1 + s)^{100}$ for largish $s$?

  26. Wow. You actually are experiencing Weasel envy!

    Don’t worry, Tom — we still think you’re smart, even if your program hasn’t been getting the attention you would have liked.

  27. keiths: Don’t worry, Tom — we still think you’re smart, even if your program hasn’t been getting the attention you would have liked.

    No one who could read my program would make much of it. It dawns on me only now that you’re someone who can’t read my program. I’m actually a bit on the slow side, don’t you think?

    You’ve derided Mung for not understanding C, and now you play games with me, when you actually don’t understand my code?

    Perhaps you should tell us what languages you can handle.

  28. An interesting open question is the extent to which more traditional programming languages are evolution friendly. As one might expect, languages like C and Java are much better suited for humans than evolution.

    – Kenneth De Jong

  29. That’s pitiful, Tom.

    You aren’t the only competent programmer around here, and you’ll need to come to grips with that fact, like it or not.

  30. Mung,

    keiths admits Tom is a competent programmer.

    Of course I do. I’m not threatened by Tom’s competence.

    It’s sad that he seems to be threatened by mine.

  31. keiths: That’s pitiful, Tom.

    Can I just remind folks that Noyau is available for those wishing to step outside the rules perimeter.

  32. Alan,

    “That’s pitiful, Tom” is not a violation of the rules, and in any case your intervention is neither needed nor required.

  33. keiths:
    Alan,

    “That’s pitiful, Tom” is not a violation of the rules, and in any case your intervention is neither needed nor required.

    Didn’t say it was, Keiths. Stay within the rules and there’ll be no need for me to intervene.

  34. (replying on a more relevant thread, from the Noyau(2) thread)

    Mung:
    OMagain has nothing, as usual. Would you be able to produce if I put up some money, OMagain?

    A quantitative measure of cumulative selection and software tests demonstrating how much cumulative selection is provided by the Dawkins Weasel program.

    How much would your time and effort be worth to do that?

    One might think that Mung is proposing a bet. Does Mung envisage betting on how much cumulative selection occurs in a run of Dawkins’s Weasel simulation? Perhaps this time there will be a clear bet proposed.

    My own view is that it is worth no one’s money to ask whether the Weasel simulates cumulative selection. It just does. Does someone here say it doesn’t? maybe Mung?

  35. Mung,

    Here’s a “quantitative measure of cumulative selection” that is similar to the one that Joe proposed:

    1. Run my basic Weasel n times using a set of parameters chosen to give a reasonable convergence time. Record the number of generations it takes to reach the target phrase each time, and take the average.

    2. Repeat the process using the same parameters but with selection turned off.

    Compare the results.

    Having been spanked again and again by the Weasel, are you willing to bet against it, Mung?

  36. Mung,
    I don’t read your posts so if you are serious then someone else will have to pass on your message.

    How much would your time and effort be worth to do that?

    I’m happy to apply my usual rate card. When you specify the project in sufficient detail so a proper estimate can be made be sure to find a way to let me know.

    However, there will be no mates rates for you!

Leave a Reply