In Part A, you're going to be running simulations based on the example we did in class: flower color in a field. Flower color is controlled by one locus that produces some amount of pigment. There are two kinds of pollinators in the field, bees and birds. You will control the preferences of the birds and the bees, you will control how many birds and how many bees, and you will control how much pigment an individual copy makes.
So there will be two alleles, a red allele and a white allele. Having two copies of the red allele will produce a lot of pigment, having two copies of the white allele will produce no pigment, and having one copy of the red allele will produce variable amounts of pigment.
By controlling the number of birds and bees, the preferences of the birds and bees, and the amount of pigment produced by a single copy of the red allele, you will attempt to produce scenarios that have a particular value of h, heritability, and s, fitness differences.
# 36 flowers: 12 each with 0, 1 and 2 red copiesg <- rep(0:2, each = 12)d <- 0.50 # redness of heterozygotecol <- c(0, d, 1)[g + 1] # 0 white ... 1 rednH <- 5; tH <- 0.00 # hummingbirds, and the palest flower they choosenB <- 5; tB <- 1.00 # bees, and the reddest flower they chooselikeH <- 0.1 + 0.9 * pmin(1, pmax(0, (col - tH + 0.2) / 0.2))likeB <- 0.1 + 0.9 * pmin(1, pmax(0, (tB - col + 0.2) / 0.2))# each visitor calls on a flower it chooses 8 times; a visit sets a seedseeds <- rpois(36, 8 * (nH * likeH + nB * likeB))w <- tapply(seeds, g, mean) # white, pink, reds <- 1 - w[1] / w[3] # red 1, pink 1 - h s, white 1 - sh <- (w[3] - w[2]) / (w[3] - w[1])
Selection and drift are both always operating. Differential reproduction can have a deterministic component, a part predictable from the heritable traits — the alleles — and a random component, driven by population fluctuations in breeding and the like.
Your goal in this exercise is to try and set up the simulation to produce a trajectory of the allele frequency over time that passes through the targets shown in red boxes. You can manipulate drift by controlling the population size and the level of inbreeding. You control the starting frequency, the heritability of the purple allele, and the fitness effects of the purple allele. You will try to get the starting populations to end up in some particular place.
A note that will be helpful: the strength of selection, s, controls how hard it selects against individuals with a particular trait. The heritability is going to influence the heterozygotes. Low-heritability traits are hard for selection to fully remove from a population, because the heterozygote has the same phenotype as the most fit genotype. Selection won't differentiate between the heterozygotes and the fit homozygotes. Drift will always be at play. And so to keep a frequency low, but not zero, requires low heritability and relatively weak drift. On the other hand, to keep a frequency in the middle requires balancing selection: high heritability, overdominance.
# individuals, one locus, fitness by genotypeN <- 400F <- 0.00 # held there every generationp0 <- 0.50 # how common purple is at the starth <- 0.5s <- 0.00w <- c(1, 1 - h * s, 1 - s) # yellow homozygote, heterozygote, purplesf <- 2 * F / (1 + F) # chance an offspring's two parents are one individualk <- round(2 * N * p0)g <- matrix(sample(rep(1:0, c(k, 2 * N - k))), 2) # two copies eachfor (t in 1:100) { fit <- w[colSums(g) + 1] p1 <- sample(N, N, TRUE, fit) p2 <- ifelse(runif(N) < sf, p1, sample(N, N, TRUE, fit)) g <- rbind(g[cbind(sample(2, N, TRUE), p1)], g[cbind(sample(2, N, TRUE), p2)])}mean(g) # how common purple is after 100 generations
For this part, you're going to set heritability and selection for four different alleles relative to one another. You're also going to set the population size and the inbreeding. There will be a plot of squares that shows the genotype when you combine any two of the alleles together to make a single genotype. You're going to watch the frequency of alleles in many populations, simulated using your set heritabilities and selective strengths, and drift, change over time. And they'll be color coded by which allele is still present at the end. A mixed or gray line indicates one in which no single allele is at fixation, that is, there remains variation. Your goal will be to get frequencies, not of one population, but of all the populations together, that average out to match the bars shown in the predict panel.
Remember, weaker heritability means less efficient selection. Stronger selection means more powerful selection. And stronger drift tends to destroy variation at random. And so you're essentially exploring, fundamentally, selection strength, selection efficiency, and the randomness of drift together to produce predictable outcomes. Weak or inefficient selection and strong drift will tend to destroy variation randomly. Strong or efficient selection and weak drift will tend to destroy variation predictably. And balancing selection, where the heterozygotes have an advantage, will tend to preserve variation in proportion to how much balancing there is relative to the strength of drift.
# four alleles, each with its own h and sh <- c(0.5, 0.5, 0.5, 0.5)s <- c(0, 0, 0, 0)N <- 100; F <- 0.00w <- outer(1:4, 1:4, function(i, j) ifelse(i == j, 1 - s[i], 1 - h[i] * s[i] - h[j] * s[j]))winner <- replicate(200, { g <- matrix(sample(rep_len(1:4, 2 * N)), 2) t <- 0 while (length(unique(c(g))) > 1 && t < 600) { fit <- w[cbind(g[1, ], g[2, ])] p1 <- sample(N, N, TRUE, fit) p2 <- ifelse(runif(N) < 2 * F / (1 + F), p1, sample(N, N, TRUE, fit)) g <- rbind(g[cbind(sample(2, N, TRUE), p1)], g[cbind(sample(2, N, TRUE), p2)]) t <- t + 1 } if (length(unique(c(g))) == 1) g[1] else 0 # 0 = still mixed})table(factor(winner, 0:4)) / 200
Mutation was covered in Biology 145. And while we have not gotten to expand on it yet in lecture, we have mentioned it repeatedly as where new alleles come from. For now, the mechanics of mutation and what exactly is happening aren't going to matter so much, as we are going to explore drift and selection for a variety of new alleles. This will also be our first time looking at the age of an allele, which is directly related to a lot of the things you've done prior, but it's sort of a new wrinkle.
For this activity, you're going to control drift here as the effective population size, just using a slider labeled individuals. You are going to control the heritability and the selection of new alleles as they arise. That said, the new alleles are going to be random. They'll have some expected heritability and some variation in that. They'll have some expected selective impact and some variation around that. So you will be predicting means and variances of new mutations, of new alleles entering the population. You will then look at the distribution of the alleles that persisted and alleles that newly arose, and compare where they show up.
This involves a number of rather complicated plots. That lets us do a few different things. For one, as I talked about in class, inbreeding causes homozygosity, and homozygosity can cause disorders. If you are homozygous for two alleles that are both negative, that's bad, but if you're heterozygous for a bad allele and a good allele, and that bad allele has low heritability, that's fine. As a result, the effect of inbreeding on fitness is going to be mediated through heritability. That is, it's only bad to be homozygous if what you're homozygous for is worse than the heterozygous or the other homozygous condition.
So you're going to create simulations and run through and explore what happens. Many things will be held constant for the targets you're trying to hit. Use the practice toggle generously, and try your best to watch the simulations and figure out deductively what is happening.
# new harmful mutations, each a new allele at its own locusN <- 100; U <- 0.3 # new mutations per individual, each generationhm <- 0.50; hsd <- 0.10 # h: average, spreadsm <- 0.10; ssd <- 0.05 # s: average, spreadarise <- function(k) data.frame(h = pmin(2, pmax(-1, rnorm(k, hm, hsd))), s = pmin(1, rgamma(k, shape = (sm / ssd)^2, scale = ssd^2 / sm)))muts <- arise(0); g1 <- g2 <- matrix(0, N, 0) # two copies per locuswbar <- homo <- c()for (t in 1:400) { g <- g1 + g2 w <- exp((g == 1) %*% log(1 - muts$h * muts$s) + (g == 2) %*% log(1 - muts$s)) wbar <- c(wbar, mean(w)); homo <- c(homo, mean(rowSums(g == 2))) gam <- function(p) { k <- runif(length(g1[p, ])) < 0.5; ifelse(k, g1[p, ], g2[p, ]) } a <- sample(N, N, TRUE, w); b <- sample(N, N, TRUE, w) g1 <- t(sapply(a, gam)); g2 <- t(sapply(b, gam)) k <- rpois(1, U * N); muts <- rbind(muts, arise(k)) # new alleles, one copy each new <- matrix(0, N, k); new[cbind(sample(N, k, TRUE), seq_len(k))] <- 1 g1 <- cbind(matrix(g1, N), new); g2 <- cbind(matrix(g2, N), matrix(0, N, k))}
Fitness is going to be in reference to a particular trait. And for a particular trait, different values of it may have different fitnesses, and that surface need not be a smooth line. We've already seen this in a simple case with heterozygosity, where it could be a flat line with a dip at the homozygote for one allele, or it could be a valley, with the homozygotes both more fit than the heterozygotes, or a balancing selection style peak, with the heterozygotes at higher fitness than either homozygote. However, for most traits, which are controlled by many loci, we'll instead get a continuous function, like the pupfish we discussed, or the crossbills. In that case, there may be multiple values of the trait that have high fitness, with intermediate values having low fitness. Further, where you start and what is best may not be the same.
For this activity, you're going to try and get populations to spread out along a fitness landscape by controlling the rate of drift, the number of individuals, and the relative strength of selection: how deep the difference between trait values is.
A thing to note is that selection is not a forward-thinking process. While there may be a best peak, a best value of the trait, if the space between them is deep, selection can't move you from a low peak to a high peak. Selection is only acting instantaneously, one generation on. Moving across a valley doesn't use selection. It uses drift. While the end result may be better for the population, crossing the valley is not something selection can do on its own readily.
# two loci; the trait is the number of + copies, 0 to 4N <- 40d <- 0.10 # how deep the valley isw <- c(1, 1 - d, 1 - d, 1 - d, 1.2) # fitness at 0, 1, 2, 3, 4 copiesmu <- 0.001crossed <- replicate(50, { g <- matrix(0, N, 4) # everyone on the low peak for (t in 1:1500) { if (mean(rowSums(g)) >= 3.5) break fit <- w[rowSums(g) + 1] a <- sample(N, N, TRUE, fit); b <- sample(N, N, TRUE, fit) for (l in 0:1) # one copy from each parent at each locus g[, 2*l + 1:2] <- cbind(g[cbind(a, 2*l + sample(2, N, TRUE))], g[cbind(b, 2*l + sample(2, N, TRUE))]) flip <- runif(4 * N) < mu; g[flip] <- 1 - g[flip] } mean(rowSums(g)) >= 3.5})mean(crossed) # the share that reached the high peak