population genetics of fitness - Part I

f1 <- function(p = 0.5, t, s){
    q <- 1-p
    w_bar <- 1-t*p^2-s*q^2
    delta_p <- (p*q*2*(s-p*(s+t))) / (2*w_bar)
    
    return(list(delta_p, w_bar))
}

parameters <- data.frame(t = c(-0.5, 0  , 0.5, -0.5),
                         s = c( 0.5, 0.5, 0.5, -0.5))
models <- c("Additive", "Dominance", "Overdominance", "Underdominance")

par(mfrow = c(2,2))
ps <- seq(0.0, 1, 0.05)


for(i in 1:4){
    res <- sapply(ps, function(x){f1(p = x, t = parameters$t[i], s = parameters$s[i])})
    res <- data.frame(delta_p = unlist(res[1,]), w_bar = unlist(res[2,]))
    par(mar = c(5.1,4.1, 4.1, 4))
    plot(ps, res$delta_p, type = "o", col = "dodgerblue4", pch = 16, 
         ylab = NA, xlab = NA, axes = F, main = models[i],
         ylim = c(min(res$delta_p), max(res$delta_p)*1.2)) 
    limits <- round(seq( min(res$delta_p), max(res$delta_p), length.out = 10), 3)
    axis(side = 2, cex = 1.2, col = "dodgerblue4", col.axis = "dodgerblue4",
         at = limits)
    mtext(side = 2, "Variation in allele frequency", cex = 1, col = "dodgerblue4", line = 2.5)
    
    grid()
    abline(h = 0, col = "dodgerblue4", lty = 3)
    par(new = T)
    plot(ps, res$w_bar, lty = 2, lwd = 2, type = "b", col = "palegreen4", axes = F, xlab = NA,
         ylab = NA,  ylim = c(min(res$w_bar), max(res$w_bar)))
    
    mtext(side = 4, line = 2, "Average Fitness", cex = 1, col = "palegreen4",
          ylim = c(min(res$w_bar)*0.9, max(res$w_bar)))
    
    limits <- round(seq( min(res$w_bar), max(res$w_bar), length.out = 10), 3)
    axis(side=4, cex = 1.2, col = "palegreen4", col.axis = "palegreen4",
         at = limits)
    
    if( i == 3){
    abline(v = ps[which.max(res$w_bar)], lty = 3, col = "palegreen4", lwd = 2)
    } else if (i == 4) { 
        abline(v = ps[which.min(res$w_bar)], lty = 3, col = "palegreen4", lwd = 2)
    }
    
    axis(side = 1, cex = 1.2, col = "gray50", col.axis = "gray50")
    mtext(side = 1, "Allele frequency (p)", cex = 1, line = 2, col = "gray50")
}