#####
# This is the Wright-Fisher model.
# Mutation or selection occur on a generational rate.
# Two alleles at one locus.
# See Neuhauser 2001.
# Note that this model is based on the binomial distribution (a discrete and non-negative distribution).
N <- 50 # Haploid population size.
ploid <- 1 # Ploidy level, usually haploid or diploid.
p <- 0.5 # Frequency of Allele A1 at generation i.
j <- round(p * (ploid*N)) # Number of gametes of type A1 at generation i.
gens <- 500 # Number of generations.
sims <- 8 # number of simulations to perform.
# Fitness parameters
# If w11 > w12 > w22 = directional selection.
# If w11 < w12 < w22 = directional selection.
# If w11, w22 < w12 = overdominance.
# If w11, w22 > w12 = underdominance.
w11 <- .33 # Relative fitness of homozygote 1
w12 <- .33 # Relative fitness of heterozygote
w22 <- .33 # Relative fitness of homozygote 2
wbar <- p^2 * w11 + 2*p*(1-p)*w12+ (1-p)^2*w22 # average fitness
# Mutation parameter.
# Note that in the classical bi-allelic model mutation does not
# result in novel alleles, it simple alters allelic frequencies.
u <- 0
#u <- 0.00001
label <- paste("X = Generation Number; ", "Simulations = ", sims, "; ", ploid,"N =", ploid*N, sep = "")
plot(1:(gens+1), seq(0,1, length.out=(gens+1)), type="n", ylim=c(0, 1), main="Wright-Fisher Model", xlab = label, ylab="Allele Frequency")
# theta <- 1.23 ; mtext(bquote(hat(theta) == .(theta)))
text(gens-50, y=0.9, bquote(mu == .(u)))
text(gens-50, y=0.85, bquote(omega[11] == .(w11)))
text(gens-50, y=0.8, bquote(omega[12] == .(w12)))
text(gens-50, y=0.75, bquote(omega[22] == .(w22)))
phi <- (p*(p*w11+(1-p)*w12))/wbar
psi <- phi*(1-u) +(1-phi)*u
#####
for (h in 1:sims)
{
for (i in 1:gens)
{
wbar[i] <- p[i]^2 * w11 + 2*p[i]*(1-p[i])*w12 + (1-p[i])^2*w22
phi[i] <- (p[i]*(p[i]*w11 + (1-p[i])*w12))/wbar[i]
psi[i] <- phi[i]*(1-u) +(1-phi[i])*u
j[i+1] <- rbinom(1, ploid*N, psi[i])
p[i+1] <- j[i+1]/(ploid*N)
}
lines(1:(gens+1), p, col=h)
}
#####