In a previous lab we looked at how diploid populations with two alleles for a gene, with frequencies p and q, evolved based on the relative fitness of three phenotypes, AA, Aa and aa. We called the relative fitnesses of these phenotypes \(w_{11}\), \(w_{12}\) and \(w_{22}\) respectively. To simulate the evolution of a diploid population with three alleles for a gene we’ll need some new nomenclature.
We’ll call our three alleles A, B and C the frequencies of these alleles p, q and r.
We now have six possible phenotypes: AA, AB, AC, BB, BC and CC. We’ll call the relative fitnesses of these phenotypes \(w_{11}\), \(w_{12}\), \(w_{13}\), \(w_{22}\), \(w_{23}\) and \(w_{33}\).
The relative frequencies of these phenotypes are \(p^2\), \(2pq\), \(2pr\), \(q^2\), \(2qr\) and \(r^2\) and thus we can redefine \(\bar{w}\) as:
\[\bar{w} = w_{11} p^2 + w_{12} 2pq + w_{13} 2pr + w_{22} q^2 + w_{23} 2qr + w_{33} r^2\]
We’ll also need a new procedure for calculating the next generation of allele frequencies:
new.pqr <- function(p, q, r=1-p-q, w11=1, w12=1, w13=1, w22=1, w23=1, w33=1){
wbar <- w11*p^2 + w12*2*p*q + w13*2*p*r + w22*q^2 + w23*2*q*r + w33*r^2
new.p <- (w11*p^2 + w12*p*q + w13*p*r)/wbar
new.q <- (w22*q^2 + w12*p*q + w23*q*r)/wbar
new.r <- (w33*r^2 + w13*p*r + w23*q*r)/wbar
return(c(new.p, new.q, new.r))
}Let’s try it out:
new.pqr(p = 0.4, q=0.2)## [1] 0.4 0.2 0.4new.pqr(p = 0.4, q=0.2, w11=2)## [1] 0.4827586 0.1724138 0.3448276new.pqr(p = 0.4, q=0.2, w11=0.5, w22=2)## [1] 0.3333333 0.2500000 0.4166667new.pqr(p = 0.4, q=0.2, w11=1, w12=0.8, w13=0.6, w22=0.4, w23=0.2, w33=0)## [1] 0.6060606 0.1818182 0.2121212Let’s look 10 generations in the future.
p = 0.4; q=0.2; r=1-p-q; w11=1; w12=0.8; w13=0.6; w22=0.4; w23=0.2; w33=0
for (i in 1:10) {
p <- new.pqr(p, q, 1-p-q, w11, w12, w23, w22, w23, w33)[1]
q <- new.pqr(p, q, 1-p-q, w11, w12, w23, w22, w23, w33)[2]
r <- 1-p-q;
}
p; q; rWe’ll also need a function that returns the frequencies of phenotypes given p, q and r and assuming random mating:
pheno.freqs <- function(p, q, r=1-p-q){
return(c(p^2, 2*p*q, 2*p*r, q^2, 2*r*q, r^2))
}pheno.freqs(0.5, 0.3)## [1] 0.25 0.30 0.20 0.09 0.12 0.04pheno.freqs(0.1, 0.1)## [1] 0.01 0.02 0.16 0.01 0.16 0.64pheno.freqs(0.1, 0.9)## [1] 0.01 0.18 0.00 0.81 0.00 0.00Finally, let’s write a function that plots the phenotype and allele frequencies over time
library(ggplot2)
plot.allele.history <- function(p, q, r=1-p-q, w11=1, w12=1, w13=1, w22=1, w23=1, w33=1, gens=100){
history <- data.frame(generation=numeric(), cat=character(), Type=character(), frequency=numeric())
for (i in 1:gens){
phenos <- data.frame(generation=rep(i,6), cat=rep("phenotypes",6), Type=c("AA", "AB", "AC", "BB", "BC", "CC"), freq=pheno.freqs(p, q))
alleles <- data.frame(generation=rep(i,3), cat=rep("alleles",3), Type=c("p", "q", "r"), freq=c(p, q, r))
history <- rbind(history, phenos,alleles)
new.freqs <- new.pqr(p, q, 1-p-q, w11, w12, w13, w22, w23, w33)
p <- new.freqs[1]; q <- new.freqs[2]; r <- new.freqs[3]
}
g <- ggplot(history, aes(generation, freq, color=Type))+geom_line()+ylim(c(0,1))+
facet_grid(cat~.)
return(g)
}Let’s try this out!
plot.allele.history(p=0.35, q=0.25, w11=0.8, w22=0.8)
plot.allele.history(p=0.35, q=0.25, w11=0.8, w22=0.8, w33=0.8)
plot.allele.history(p=0.35, q=0.25, w12=0.8, w23=0.8)Can you create…
A scenario where all three alleles survive in different frequencies?
A scenario where two of the three alleles survive in the equal frequencies?
A scenario where two of the three alleles survive in different frequencies?
Challenge:
Write equations for the equilibrium allele frequencies. Can you simplify this set of equations to describe what must be true of the allele frequencies when equilibrium is reached?