Simulating substitution rate of neutral mutations

Question

I am trying to computationally simulate a population based on the Wright-Fisher model I would like to get to the classic result of the neutral theory of molecular evolution that the rate of neutral substitution equals the rate of mutation. However, the results of my simulations never show that. I will briefly explain my simulation process and how I calculate the mutation and substitution rates. Please let me know if there is something wrong in my procedure.

I have a non-structured haploid population that grows in discrete time steps. I start with a population of N individuals in an environment with carrying capacity K. I use a K-allele model of mutation. At each generation, every individual will reproduce n offsprings and it will immediately die afterwards. A u fraction of offsprings will be chosen randomly to have a single mutation in each generation (u is mutation rate per gene per generation). At each generation, I record how many mutations have occurred and how many mutations have reached fixation in the population (reached frequency of more than 99% of individuals). Thus, mutation rate and substitution rate per generation is calculated as the number of mutations and number of substitutions divided by population size, respectively. With this set up, even with hundreds of repetitions, the rate of fixation will never be as big as mutation rate u.

Where do you think is the problem? I have struggled with this a lot and I appreciate any help.

Edit: My code is written in Python. This my pseudo-code:

Parameters: population size, genome_size, mutation_rate, generation_number, fixation_threshold.

'''
Population_size: Here it is 1500
Genome_size: the number of loci that can be mutated. Each locus can be mutated to any of ACTG. Here 10^4
Mutation_rate: number of mutations per genome per generation. Here 0.05
Generation_number: The number of generations that the population grows. Here 1000
Fixation_threshold: The fraction of individuals having an allele in the population so that the allele is considered fixed. Here 0.99
'''

The initial population is created. All individuals are the same without any mutations
A vector to hold the number of mutations and a vector to hold the number of fixations is created

For each generation:
    I determine how many individuals in each generation should be mutated. Here 8.
    I randomly choose  (uniform with replacement) the individuals to be mutated 
    For each individual in the population:
        It will produce an offspring. The offspring inherits the genome of its parent, i.e. the mutations of the parent.
        If the individual is among the ones to be mutated, the offspring will receive a mutation that is randomly picked.
            The mutation can be a back mutation
        The individual will die
    I record the number of mutations that occurred (it will always be 8 in this case)
    I record any mutation that has a frequency in the population of more than 0.99 of individuals. This is the number of fixations

I run this code for at least 100 times. I expect that the average number of fixations to be equal to the average number of mutations in each generation, or they converge over generations. Because the neutral theory of Evolution by Kimura predicts that for neutral sites in the genome, substitution rate is equal to the mutation rate.

Answer 1

This just got too big for a comment and I think might actually qualify as an answer.

Can you give me a source for what you're reading? According to this, which matches the definitions you are using, K=(Nu)/N=u where K is the substitution rate, u is the mutation rate in the population, Nu is the individual mutation rate, and N is the population size.

I went back to Kimura's original paper in Nature in 1968 and he derives that u(p), the substitution rate in the population, is (almost exactly) the same as the mutation rate in an individual(p). But this follows from a different paper(prolific guy) where he shows that the fixation probability for a given allele is roughly 1/N, where N is the population size. Your individual mutation rate in your model is not u, but Nu.

You are confusing the mutation/substitution rate in the population, which is used for molecular clocks and other things, and the individual mutation rate. If the molecular mutation rate, is say, 1/20 (as you assume in your model) then you are expecting a fixation rate per mutation of 1, essentially. If the substitution rate is going to be the same as the molecular mutation rate, every single mutation must reach fixation. If the same locus mutates twice in the same generation(quite possible with enough individuals), then both must somehow reach fixation, even though the first to reach fixation will necessarily destroy the other.

I hope that helps you sort out why Nu will never equal K, and helps the results of your model make sense.

(sidebar: Nu isn't the raw molecular mutation rate but the individual mutation rate, which is actually higher because while growing to adulthood the average human gamete is involved in about 50 cell divisions, which leads to a few more mutations per person.)

Answer 2

I don't fully understand how your model works

I don't fully understand how you model your population. Is it a standard $oop$ (object-oriented programming) where you simulate each individual? Or is it a simulation where you already use some mathematical model? You could eventually copy-paste your code. it is probably not very long, right? What language did you use? Do you have only one locus? Do you simulate several positions within these locus? Do you allow for back mutations?

It seems that you allow the population size to change over time and I don't quite understand how your parameter $K$ does limit the growth rate though. You should write down some equations about that. In any case, change in population size will modify the probability of mutation fixation.

Here is probably some important modifications to do on your code

It is a bit hard to comment on your simulations without actually seeing some equations or seeing the code but I really think that you'll need two important modifications in your code.

1. fix population size $N$.
2. To simulate reproduction, randomly choose (uniform distribution) individuals (with replacement) from the population. The chosen individuals form the new generation. For example:

In R:

new_pop = sample(old_pop, N, replace=T)

or directly

pop = sample(pop, N, replace=T)

In Python:

from random import randrange
new_pop = []
for i in xrange(N):
    new_pop.append(old_pop[randrange(N)])

When you get your results, don't forget that you're talking about a stochastic process and you may want to run enough simulation or enough independent loci in order to be able to compare with the predictions from mathematical models.