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.