Exploring the interaction of drift, mutation, and migration

Background

Although it is very simple, Sewall Wright’s infinite island model of migration is a useful way to explore how drift, migration, and mutation interact to influence the genetic structure of populations that are divided into subpopulations. In this model we assume that

1. Every allele introduced into a population differs from any of the resident alleles.

2. Every population receives the same fraction, \(m\), of migrants.¹

Given those assumptions and assuming an infinite alleles model of mutation, we showed that²

\[ F_{ST} = \frac{1}{4N_e(m + \mu) + 1} \quad . \]

A more realistic model of migration is the finite island model of migration.³ It’s still very simple, but it adds one thing that makes it a bit more realistic:

It assumes that there are a finite number of populations exchanging migrants with one another. Every population still receives the same fraction of migrants, \(m\), but we now also specify that the fraction of migrants coming from any one population into this one is \(m/(k - 1)\), where \(k\) is the number of popuations.

The formula for \(F_{ST}\) under the finite island model is a lot more complicated than the one above. If you’re interested, you can look in the R code below.⁴ Let’s compare the predictions for a few values of \(m\) and \(k\) with \(N_e = 100\) and \(\mu = 10^{-3}\).

Comparing the predictions

rm(list = ls())

infinite_island <- function(n_e, m, mu) {
  return(1/(4*n_e*(m + mu) + 1))
}

lambda <- function(mu, m, k) {
  m_prime <- m*k/(k - 1)
  l <- (1/2)*((1 - mu)^(-2)*(1- m_prime)^(-2) - 1)
  return(l)
}

gamma <- function(mu, m, k) {
  m_prime <- m*k/(k - 1)
  g <- (1 + m_prime*(2 - m_prime)/(k*mu*(2 - mu)*(1 - m_prime)^2))^(-1)
  return(g)
}

finite_island <- function(n_e, m, mu, n_islands) {
  l <- lambda(mu, m, n_islands)
  g <- gamma(mu, m, n_islands)
  h_s <- 4*n_e*l*g/(4*n_e*l*g + 1)
  h_b <- h_s*(1 + 4*n_e*l)/(4*n_e*l)
  h_t <- ((n_islands - 1)/n_islands)*h_b + (1/n_islands)*h_s
  f_st <- (h_t - h_s)/h_t
  return(f_st)
}

n_e <- 100
mu <- 10^(-4)
k_vals <- rep(c(5, 25, 250, 2500), 2)
m_vals <- c(rep(0.01, 4), rep(0.05, 4))
dat <- data.frame(k = k_vals,
                  m = m_vals,
                  Infinite = rep(NA, length(k_vals)),
                  Finite = rep(NA, length(m_vals)),
                  Ratio = rep(NA, length(k_vals)))
for (i in 1:nrow(dat)) {
  dat$Infinite[i] <- infinite_island(n_e, dat$m[i], mu)
  dat$Finite[i] <- finite_island(n_e, dat$m[i], mu, dat$k[i])
  dat$Ratio[i] <- dat$Infinite[i]/dat$Finite[i]
}
round(dat, 3)

Comparing predictions to simulations

From the table above, it’s apparent that predictions of the infinite island and finite island models are rather different when the number of islands is small. Let’s run a simulation of the finite island model and see how well the simulation results match up to each prediction.

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.4     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   2.0.1     ✓ forcats 0.5.1
── Conflicts ──────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(hierfstat)

initialize <- function(k, n_loci) {
  p <- array(dim = c(n_loci, k, 2))
  for (j in 1:n_loci) {
    for (i in 1:k) {
      p[j, i, 1] <- runif(n = 1, min = 0, max = 1)
      p[j, i, 2] <- 1 - p[j, i, 1]
    }
  }
  return(p)
}

make_symmetric_matrix <- function(m, k) {
  M <- diag(x = 1 - m, nrow = k)
  for (i in 1:(k - 1)) {
    for (j in (i+1):k) {
      M[i, j] <- m/(k-1)
      M[j, i] <- M[i, j]
    }
  }
  return(M)
}

sample_two_alleles <- function(p) {
  allele_1 <- ((runif(1) < p) + 1)*10
  allele_2 <- (runif(1) < p) + 1
  return(allele_1 + allele_2)
}

make_sample <- function(p, n_sample) {
  n_pops <- dim(p)[2]
  n_loci <- dim(p)[1]
  names <- paste("Locus", seq(1:n_loci), sep = "")
  pops <- numeric(n_pops*n_sample)
  ct <- 1
  for (k in 1:n_pops) {
    population <- paste("Population", k, sep = "")
    for (i in 1:n_sample) {
      pops[ct] <- population
      ct <- ct + 1
    }
  }
  locus_data <- matrix(nrow = n_pops*n_sample, ncol = n_loci)
  for (i in 1:n_loci) {
   ct <- 1
    for (k in 1:n_pops) {
      for (n in 1:n_sample) {
        locus_data[ct, i] <- sample_two_alleles(p[i, k, 1])
        ct <- ct + 1
      }
    }
  }
  df <- data.frame(locus_data)
  colnames(df) <- paste("Locus", seq(1:n_loci), sep = "")
  df$Population <- pops
  df <- relocate(df, Population)
  return(df)
}

simulate <- function(n_e, m, mu, n_islands, n_loci, n_gen, n_sample) {
  p <- initialize(n_islands, n_loci)
  M <- make_symmetric_matrix(m, n_islands)
  V <- make_symmetric_matrix(mu, 2)
  g_st <- numeric(n_loci)
  p_star <- array(dim = c(n_loci, n_islands, 2))
  for (i in 1:n_gen) {
    for (k in 1:n_loci) {
      p_star[k, , ] <- M %*% p[k, , ] %*% V
      for (j in 1:n_islands) {
        p[k, j, ] <- rmultinom(1, 2*n_e, p_star[k, j, ])/(2*n_e)
      }
      mu_p <- mean(p[k, , 1])
      g_st[k] <- ((n_islands - 1)/n_islands)*var(p[k, , 1])/(mu_p*(1 - mu_p))
    }
  }
  df <- make_sample(p, n_sample)
  return(wc(df)$FST)
}

run_simulation <- function(n_e, m, mu, n_islands, n_loci, n_gen, 
                           n_sample, n_repetitions)
{
  f_st <- numeric(n_repetitions)
  for (i in 1:n_repetitions) {
    if (n_repetitions > 50) {
      cat(".", sep="")
      if ((i %% 50) == 0) {
        cat(i, "\n", sep ="")
      }
    }
    f_st[i] <-simulate(n_e, m, mu, n_islands, n_loci, n_gen, n_sample)
  }
  return(f_st)
}

plot_simulation <- function(f_st, n_e, m, mu, n_islands) {
  df <- data.frame(F_st = f_st)
  p <- ggplot(df, aes(x = F_st)) +
    geom_histogram(bins = 20, alpha = 0.4) +
    geom_vline(xintercept = infinite_island(n_e, m, mu),
               linetype = "dashed",
               color = "red") +
    geom_vline(xintercept = finite_island(n_e, m, mu, n_islands),
               linetype = "dashed",
               color = "blue") +
    ggtitle(paste("N_e = ", n_e, ", m = ", m, ", mu = ", mu, 
                  ", n_islands = ", n_islands, sep = "")) +
    theme_bw()
  return(p)
}

f_st <- run_simulation(100, 0.05, 1e-4, 5, 10, 1000, 25, 100)

..................................................50
..................................................100

plot_simulation(f_st, 100, 0.05, 1e-4, 5)

The dashed red line is the prediction from the infinite island model. The dashed blue line is the prediction from the finite island model.

Exercise 6

The last plot shows that neither the infinite island nor the finite island prediction captures the results of this particular simulation very well, it only does so on average and there’s a lot of variation from one replicate to the next. This outcome reflects the “evolutionary” or “genetic” sampling we discussed when comparing Nei’s \(G_{ST}\) with Weir and Cockerham’s \(\theta\). This weeks lab exercise will explore that variability a bit further and explore condidtions that affect the amount of variability and how useful the infinite island approximation is. To do so, you’ll want to use the following code.⁵ You’ll be changing some or all of

• n_loci - The number of loci included in the genetic sample,
• n_sample - The number of (diploid) individuals sampled in each population,
• n_gen - The number of generations for which to run the simulation, and
• n_repetitions - The number of times to repeat the simulation for any set of the first three choices. This item corresponds to the number of evolutionary or genetic sampling events included in your simulation.

For each choice of these four parameters, the simulation will run through all combinations of N_e = (25, 200, 250, 500), m = (0.05, 0.01, 0.005, 0.001), mu = (0.001, 0.0001), and n_islands = (10, 25, 50).

n_loci <- 10
n_sample <- 25
n_gen <- 1000
n_repetitions <- 10

df <- data.frame(N_e = NA, m = NA, mu = NA, n_islands = NA, Observed = NA,
                 Infinite = NA, Finite = NA)
for (n_e in c(25, 100, 250, 500)) {
  for (m in c(0.05, 0.01, 0.005, 0.001)) {
    for (mu in c(0.001, 1e-4)) {
      for (n_islands in c(10, 25, 50)) {
        cat("N_e = ", n_e, ", m = ", m, ", mu = ", mu, 
            "n_islands = ", n_islands, "\n" )
        f_st <- run_simulation(n_e, m, mu, n_islands, n_loci, n_gen, 
                               n_sample, n_repetitions)
        df <- add_row(df, 
                      N_e = rep(n_e, n_repetitions),
                      m = rep(m, n_repetitions), 
                      mu = rep(mu, n_repetitions), 
                      n_islands = rep(n_islands, n_repetitions),
                      Observed = f_st,
                      Infinite = infinite_island(n_e, m, mu),
                      Finite = finite_island(n_e, m, mu, n_islands))
      }
    }
  }
}
df <- filter(df, !is.na(N_e))

Running one set of simulations for a particular choice of n_loci, n_sample, n_gen, and n_repetitions takes close to half an hour on my MacBook. You’ll get a printout periodically that lets you know what the current set of parameters being simulated is, and you’ll see the set or parameters that’s already been completed. Because running the simulation for one set of parameters takes so long, you can’t possibly explore all of the combinations of those parameters you might be interested in, but I’d like you to explore at least five and address the following questions:⁶

1. Are there conditions under which the predictions of the infinite island and finite island model are similar enough that we don’t need to worry about the added mathematical complexity of the finite island model?

2. How close is the relationship between observed estimates of \(F_{ST}\) from the simulations and those predicted from the infinite island model? from the finite island model?

3. Under what conditions is the observed variation in \(F_{ST}\) relatively small?

A few hints

The easiest way to run different combinations of parameters is to copy the code block above and replace n_loci, n_sample, n_gen, and n_repetitions with new choices. To keep track of what you’re doing, I recommend doing that for one combination of parameters, examining the results as described below, and moving on to a new set of parameter (with a new code block) after that.

You can plot the observed \(F_{ST}\) versus the \(F_{ST}\) predicted from either model using the code below.⁷ In addition to showing how well the simulations match the predictions, you can get a feel for how much variabiity there is among the simulations by seeing how scatterd the points are.

## for the infinite island model
##
p <- ggplot(df, aes(x = Infinite, y = Observed)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", 
              color = "red") + theme_bw()
p
## for the finite island model
##
p <- ggplot(df, aes(x = Finite, y = Observed)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", 
              color = "red") + theme_bw()
p

If you’d like a more quantitative way to compare the models,⁸ a simple way is to run linear regressions of Observed on Infinite and Expected like this:

summary(lm(Observed ~ Infinite, data = df))
summary(lm(Observed ~ Finite, data = df))

If predictions were perfect the (Intercept) term in the regression would be 0, and the other term (Infinite in the first regression and Finite in the second) would be 1. Given the variability of the simulations, you won’t see numbers that are exactly 0 or 1. You can judge whether they are “close” by noticing whether the number in the Estimate column is more than two standard errors away from what you predict.

The residual standard error below the table measures how much error there is around the regression, i.e., it gives you a way to compare quantitatively how much variation there is in \(F_{ST}\).

---

1. Notice that \(m\) is the backwards migration rate, i.e., the fraction of a population composed of migrants that come from another population, not the fraction of residents that migrate to another population.↩︎

2. Actually, we didn’t show what happens with both migration and mutation, but from what we saw you can probably guess where this came from.↩︎

3. Latter, B.D.H. 1973. The island model of population differentiation: a general solution. Genetics 73:147-157 doi: 10.1093/genetics/73.1.147↩︎

4. You won’t find the formula for \(F_{ST}\) in Latter, but the formulas in the R code follow from the results there.↩︎

5. Notice that in order to run this code, you’ll also need to run the code above so that the relevant functions are loaded in R’s memory.↩︎

6. Keep in mind that you’ll only be able to answer these questions for the set of simulations that you complete. I’ll be grading your responses on the simulation results that you show, not on what an exhaustive examination of all possibilities might be.↩︎

7. The diagonal red dashed line is the 1:1 line. Every point should fall on this line if observed exactly matches the prediction.↩︎

8. And I’d encourage you to take a more quantitative approach↩︎