# Simple haploid selection model simulator
# p0 = initial allele frequency for "A" allele/genotype
# w1 = fitness for "A" allele/genotype
# w2 = fitness for "a" allele/genotype
# n = number of generations to follow
# This is an example of a deterministic model. It is deterministic in
# that no matter how many times we run it, if we start with the same set
# of starting conditions (p0, w1, w2) we will get exactly the same
# result. In the "real world" we also need to account for population
# size (and fluctations) and how that can influence allele frequencies
# via drift, and random variation (due to the environment for instance).
haploid.selection <- function(p0 = 0.01, w1 = 1, w2 = 0.9, n = 100) {
# Initialize some vectors to store allele frequencies (and its
# delta) and mean population fitness
p <- rep(NA,n) # a vector to store allele frequencies
delta_p <- rep(NA, n) # a vector to store changes in allele frequencies
w_bar <- rep(NA, n)
# starting conditions
p[1] <- p0 # starting allele frequencies
delta_p[1] <- 0 #change in allele frequency
w_bar[1] <- (p[1]*w1) + ((1-p[1])*w2)
# now we need to loop from generation to generation
for ( i in 2:n) {
w_bar[i - 1] <- (p[i - 1]*w1) + ((1-p[i - 1])*w2) # mean population fitness
p[i] <- (w1*p[i - 1])/w_bar[i - 1]
delta_p[i] <- p[i] - p[i-1]
}
if (any(p > 0.9999)) {
fixation <- min(which.max(p > 0.9999))
cat("fixation for A1 occurs approximately at generation:", fixation )
} else {
maxAlleleFreq <- max(p)
cat("fixation of A1 does not occur, maximum allele frequency is:", maxAlleleFreq )
}
# Plot the results
par(mfrow=c(2,2))
# 1. mean population fitness over time
plot(x = 1:n, y = w_bar,
xlab = "generations",
ylab = expression(bar(w)),
pch=20, col="red", cex = 2, cex.lab = 1.5,
main = paste("p0 = ", p0, "and s = ", (1 - (w2/w1))))
# 2. change in allele frequency over time
plot(x = 1:n, y = p,
xlab="generations",
ylab="Allele frequency (p)",
pch = 20, col = "red", cex.lab = 1.5)
# 3. plot of p[t+1] vs p[t]
p.1 <- p[-n]
p.2 <- p[-1]
plot(p.2 ~ p.1,
xlab = expression(p[t]),
ylab = expression(p[t+1]),
pch = 20, col = "red", cex = 2, cex.lab = 1.5)
# 4. plot of allele frequency change
plot(x = 2:n, y = delta_p[-1],
xlab = "generation",
ylab= expression(paste(Delta,"p")),
pch = 20, col = "red", cex = 2, cex.lab = 1.5)
}
haploid.selection(p0 = 0.001, w2 = 0.9, n = 1000)
haploid.selection(p0 = 0.000001, w2 = 0.9, n = 250)
haploid.selection(p0 = 0.000001, w2 = 0.65, n = 250)
haploid.selection(p0 = 1/(10^8), w1 = 1, w2 = 0.85 , n = 100)