pop_gene_model

Enter the real-world data on polyandry and SR frequency

real_world_data <- data.frame(
  real_freq_SR   = c(10.2, 13.5, 11.5, 11.0, 7.5, 3.5, 0),
  freq_polyandry = c(0.58, 0.52, 0.73, 0.59, 0.73, 0.73, 0.92)
)

# The average of P1 and P2 in from GIRALDO-PEREZ et al, https://onlinelibrary.wiley.com/doi/pdf/10.1111/jeb.12792
average_paternity_SR_males <- 0.2105

# The estimated frequency of SR male progeny who inherit SR
real_world_k <- 0.96

Function to make the mating table

The “mating table” is list of all the possible mating types (e.g. ST/SR mother with a SR/Y father, etc etc), which lists the possible offspring genotypes and their associated frequencies. The precise values in the mating table are a function of \(k\) (i.e. the strength of segregation distortion) and expected paternity of SR males when the female is double-mated to an ST and a SR male. We assume that the meiosis is Mendelian except in SR males, and that all crosses except those with a SR father produce a 50:50 sex ratio.

make_mating_table <- function(k, paternity_of_SR_males, mating_types){

  offspring <- vector(mode = "list", nrow(mating_types))
  
  make_offspring <- function(row, type, zygote_freq){
    data.frame(mother = mating_types$mother[row],
               father1 = mating_types$father1[row],
               father2 = mating_types$father2[row],
               type, 
               zygote_freq, 
               stringsAsFactors = FALSE)
  }
  
  g0 <- c("STST_female", "STSR_female",      # possible offspring types
          "ST_male", "SR_male")
  g1 <- c("STSR_female", "SRSR_female",
          "ST_male", "SR_male")

  p1 <- c(0.5 *  (1 + k), 0.5 * (1 - k))
  p2 <- c(0.25 * (1 + k), 0.25 * (1 + k),
          0.25 * (1 - k), 0.25 * (1 - k))  
  
  p3 <- c(0.5 * (1 - paternity_of_SR_males), # ST male wins, he transmits his ST-X with 50% chance
          0.5 * paternity_of_SR_males * (1 + k), # SR male wins with probability sc, he transmits his SR-X 0.5*(1+k) chance
          1 - 0.5 * (1 - paternity_of_SR_males) - 0.5 * paternity_of_SR_males * (1 + k)) # sons created by the remaining two possibilities

  p4 <- c(0.25 * (1 - paternity_of_SR_males), # ST male wins, he transmits his ST-X with 50% chance, same X in female meiosis
          0.25 * (1 - paternity_of_SR_males) + 0.25 * paternity_of_SR_males * (1 + k), # Either male wins and transmits the other X than the female does
          0.25 * paternity_of_SR_males * (1 + k), # SR male wins, he transmits his SR-X 0.5*(1+k) chance, same X in female meiosis
          0, 0)
  p4[4:5] <- (1 - sum(p4[1:3])) / 2 # sons created by the remaining possibilities

  rep25 <- rep(0.25, 4)
  rep5  <- rep(0.5, 2)
  
  # Single mating types
  offspring[[1]]  <- make_offspring(1, c("STST_female", "ST_male"), rep5)
  offspring[[2]]  <- make_offspring(2, g0, rep25)
  offspring[[3]]  <- make_offspring(3, c("STSR_female", "SR_male"), rep5)
  offspring[[4]]  <- make_offspring(4, c("STSR_female", "ST_male"), p1)
  offspring[[5]]  <- make_offspring(5, g1, p2)
  offspring[[6]]  <- make_offspring(6, c("SRSR_female", "SR_male"), p1)
  
  # Double mating types - two ST males
  offspring[[7]]  <- make_offspring(7, c("STST_female", "ST_male"), rep5)
  offspring[[8]]  <- make_offspring(8, g0, rep25)
  offspring[[9]]  <- make_offspring(9, c("STSR_female", "SR_male"), rep5)
  
  # Double mating types - ST and SR males
  offspring[[10]]  <- make_offspring(10, c("STST_female", "STSR_female", "ST_male"), p3)
  offspring[[11]]  <- make_offspring(11, c("STST_female", 
                                           "STSR_female", 
                                           "SRSR_female", 
                                           "ST_male", 
                                           "SR_male"), p4)
  offspring[[12]]  <- make_offspring(12, c("STSR_female", "SRSR_female", "SR_male"), p3)


  # Double mating types - two SR males
  offspring[[13]]  <- make_offspring(13, c("STSR_female", "ST_male"), p1)
  offspring[[14]]  <- make_offspring(14, g1, p2)
  offspring[[15]]  <- make_offspring(15, c("SRSR_female", "SR_male"), p1)
  
  output <- do.call("rbind", offspring)
  names(offspring) <- paste(mating_types[,1], mating_types[,2], sep = " x ")
  output
}

make_mating_types <- function(){
  expand.grid(mother = c("STST_female",
                         "STSR_female",
                         "SRSR_female"),
              father1 = c("ST_male",
                          "SR_male"),
              father2 = c("none",
                          "ST_male",
                          "SR_male"),
              stringsAsFactors = FALSE) %>%
    mutate(sorted_paste = map_chr(1:n(), function(i) paste0(sort(c(mother[i], father1[i], father2[i])), collapse = ""))) %>%
    distinct(sorted_paste, .keep_all = TRUE) %>% select(-sorted_paste)
}

mating_types <- make_mating_types()

make_mating_table(k = 0.9, paternity_of_SR_males = 0.23, mating_types) %>%
  rename(zygote_type = type) %>%
  kable() %>%
  kable_styling() %>%
  scroll_box(height = "500px")

Function to determine the mating type frequencies for the population

Given a set of genotype frequencies, it is simple to calculate the frequencies of each mating type. First, we implement selection, such that each genotype is represented in the mating types according to the product of its frequency and its fitness. Second, we determine the frequency of cross between male genotype \(i\) and female genotype \(j\) as \(i\times j\).

find_mating_type_frequencies <- function(pop, mating_types, freq_polyandry){

  # Implement viability selection
  pop$prop <- pop$prop * pop$fitness
  pop$prop <- pop$prop / sum(pop$prop)
  
  single_maters <- 1:6
  double_maters <- 7:15
  
  # Mating type freqs for the single-mating females is freq_single_maters * freq_females * freq_males, post selection
  mating_types$mating_freq[single_maters] <- (1 - freq_polyandry) * 
    pop$prop[match(mating_types$mother[single_maters], pop$type)] * 
    pop$prop[match(mating_types$father1[single_maters], pop$type)]
  
  # Mating type freqs for the single-mating females is freq_double_maters * freq_females * freq_father1 * freq_father2, post selection
  # NB the ST + SR double matings are doubled, since they can occur with either male as father1
  mating_types$mating_freq[double_maters] <- freq_polyandry * 
    pop$prop[match(mating_types$mother[double_maters], pop$type)] * 
    pop$prop[match(mating_types$father1[double_maters], pop$type)] *
    pop$prop[match(mating_types$father2[double_maters], pop$type)] 
  
  mating_types$mating_freq[10:12] <- 2 * mating_types$mating_freq[10:12]
  mating_types$mating_freq <- mating_types$mating_freq / sum(mating_types$mating_freq) # ensure the mating freqs sum to one

  mating_types
}

Define the main simulation function

This function iterates over generations, implementing selection and reproduction each time, until either A) the SR alelle fixes, B) the SR alelle goes extinct, or C) the generation timer expires.

# Helper function to calculate the frequency of the SR allele
calc_prop_SR <- function(pop){
  sum(c(pop$prop[c(1, 4)], 2 * pop$prop[2])) / (sum(c(2 * pop$prop[c(2,4,5)], pop$prop[c(1,3)])))
}

run_simulation <- function(generations,
                           k,
                           paternity_of_SR_males,
                           freq_polyandry,
                           w_STSR_female,
                           w_SRSR_female,
                           w_SR_male,
                           initial_freq_SR,
                           mating_types){
  
  # Make the initial population
  pop <- data.frame(
    type = c("STST_female",
             "STSR_female",
             "SRSR_female",
             "ST_male",
             "SR_male"),
    prop = c(0.5 * (1 - initial_freq_SR) ^ 2,
             (1 - initial_freq_SR) * initial_freq_SR,
             0.5 * initial_freq_SR ^ 2,
             0.5 * (1 - initial_freq_SR),
             0.5 * initial_freq_SR), 
    fitness = c(1, w_STSR_female, w_SRSR_female, 1, w_SR_male),
    stringsAsFactors = FALSE) %>% 
    arrange(type) 

  # Make the mating table for this specific parameter space
  mating_table <- make_mating_table(k, paternity_of_SR_males, mating_types)
  
  # Iterate over generations
  for(i in 1:generations){
    
    # Find the mating type frequencies, find the offspring frequencies, and renormalise the frequencies to sum to one
    prop_col <- mating_table %>%
      left_join(find_mating_type_frequencies(pop, mating_types, freq_polyandry), 
                by = c("mother", "father1", "father2")) %>%
      mutate(offspring_freq = mating_freq * zygote_freq) %>%
      group_by(type) %>%                         # sorts by type?
      summarise(prop = sum(offspring_freq)) %>%
      mutate(prop = prop / sum(prop)) %>%
      pull(prop)
    
    # Calculate the frequency of SR. 
    # Quit early if SR fixed/extinct-ish
    prop_SR <- calc_prop_SR(pop %>% mutate(prop = prop_col))
    # if(is.na(prop_SR)) print(pop)   # FOR DEBUGGING
    pop <- pop %>% mutate(prop = prop_col)
        
    if(prop_SR > 0.99) return(pop) 
    if(prop_SR < 0.0001) return(pop)
  }
  pop
}

Define a helper function to run the main simulation on a data frame of parameters

A helper function that runs run_simulation() on a data frame of parameter values.

run_many_simulations <- function(parameters){
  
  mating_types <- make_mating_types() %>% mutate(mating_freq = NA)
  
  # Carefully pass all the parameters to run_simulation()
  run_one_simulation <- function(row, parameters, mating_types){
    
    genotypes <- run_simulation(
      generations = parameters$generations[row],
      k = parameters$k[row],
      paternity_of_SR_males = parameters$paternity_of_SR_males[row],
      freq_polyandry = parameters$freq_polyandry[row],
      w_STSR_female = parameters$w_STSR_female[row],
      w_SRSR_female = parameters$w_SRSR_female[row],
      w_SR_male = parameters$w_SR_male[row],
      initial_freq_SR = parameters$initial_freq_SR[row],
      mating_types = mating_types
    )
    
    output <- data.frame(parameters[row, ]) %>% as_tibble()
    output$genotypes <- list(genotypes %>% select(-fitness))
    output
  }
  
  # Loop over all the different parameter spaces
  lapply(1:nrow(parameters), run_one_simulation, parameters = parameters, mating_types = mating_types) %>%
    do.call("rbind", .) %>% 
    mutate(prop_SR = map_dbl(genotypes, calc_prop_SR))
}

Define the parameter space to investigate

Used to define the parameter spaces that are plotted in Figure X.

resolution <- 101
gen <- 100
initial_freq <- 0.1

parameters <- rbind(
  
   expand.grid(
    generations = gen,
    k = real_world_k,
    paternity_of_SR_males = seq(0, 0.5, length = resolution),
    freq_polyandry = seq(0, 1, length = resolution),
    w_STSR_female = 0.88, # Fitness from our Table of results
    w_SRSR_female = 0.42, # Fitness from our Table of results
    w_SR_male = 1,
    initial_freq_SR = initial_freq,
    facet = "Recessive cost in females")
  
) %>% as_tibble()

Run the simulation

cores <- 7

run_all <- function(){
  if(!file.exists("data/simulation_output.rds")){ # if there is no data yet....
    
    # Divide the job into 100 equal-ish chunks so that we can see the % completion as it proceeds
    split_parameters <- split(parameters, sample(1:100, nrow(parameters), replace = TRUE))
    n_chunks <- length(split_parameters)
    for(i in 1:n_chunks){
      print(paste(i, "out of", n_chunks))
      plan("multicore")
      chunk_size <- floor(nrow(split_parameters[[i]]) / cores)
      simulation_output <- split(split_parameters[[i]], 
                                 c(rep(1:(cores-1), each = chunk_size), 
                                   rep(cores, nrow(split_parameters[[i]]) - chunk_size*(cores-1)))) %>%
        future_lapply(run_many_simulations) %>% do.call("rbind", .)
      saveRDS(simulation_output, paste("data/simulation_output_", i, ".rds", sep = ""))
    }
    
  } else { # if there is already some data...
    
    simulation_output <- readRDS("data/simulation_output.rds")
    
    if(nrow(simulation_output) != nrow(parameters)){
      pasted_done <- apply(simulation_output[, names(simulation_output) %in% names(parameters)], 1, paste0, collapse = "_")
      pasted_to_do <- apply(parameters, 1, paste0, collapse = "_")
      trimmed_parameters <- parameters[!(pasted_to_do %in% pasted_done), ]
      
      # Divide the job into 100 equal-ish chunks so that we can see the % completion as it proceeds
      split_parameters <- split(trimmed_parameters, sample(1:10, nrow(parameters), replace = TRUE))
      n_chunks <- length(split_parameters)
      for(i in 1:length(split_parameters)){
        print(paste(i, "out of", n_chunks, sep = ""))
        plan("multicore")
        chunk_size <- floor(nrow(split_parameters[[i]]) / cores)
        simulation_output <- split(split_parameters[[i]], 
                                   c(rep(1:(cores-1), each = chunk_size), 
                                     rep(cores, nrow(split_parameters[[i]]) - chunk_size*(cores-1)))) %>%
          future_lapply(run_many_simulations) %>% do.call("rbind", .)
        saveRDS(simulation_output, paste("data/simulation_output_surplus_", i, ".rds", sep = ""))
      }
    }
  }
  
  simulation_output <- list.files(path = "data", pattern = "simulation_output_", full.names = TRUE) %>%
    lapply(readRDS) %>% do.call("rbind", .)
  saveRDS(simulation_output, "data/simulation_output.rds")
}

run_all()
simulation_output <- readRDS("data/simulation_output.rds")
  
# use to delete and re-run:    unlink("data/simulation_output.rds")
# use to delete individual files:     unlink(list.files(path = "data", pattern = "simulation_output_", full.names = TRUE))

Make Figure 3

figure_3 <- simulation_output %>% 
  filter(w_SR_male == 1) %>%
  mutate(k = paste("Strength of drive:", k)) %>%
  ggplot(aes(100 * freq_polyandry, 100 * paternity_of_SR_males, fill = 100 * prop_SR)) +
  geom_blank() +
  geom_raster() + 
  stat_contour(aes(z = 100 * prop_SR), colour = "grey10", binwidth = 10) +
  geom_hline(yintercept = 50, colour = "black", linetype = 3) +
  labs(x = "% females that mate twice", 
       y = "% paternity of SR males in a double mating") +
  annotate("pointrange", x = 100 * real_world_data$freq_polyandry, y = 21.05, 
           pch = 21, alpha = 0.7,
           fill = "tomato", size = .3, ymin = 20, ymax = 20, colour = "grey10") +
  scale_fill_distiller(palette = "YlGnBu", direction =  1, name = "% SR at equilibrium") + #, limits = c(0,100)) +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) +   # , labels = c(0, 0.25, 0.5, 0.75, 1)
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        text = element_text(family = "Lato", size = 12),
        strip.background = element_rect(colour = "black", fill = "grey90", size = .8),
        legend.position = "top")

figure_3 %>% ggsave(filename = "figures/figure_3.pdf", height = 4.9, width = 4.5)
figure_3

Make Figure 4

parameters <- rbind(
  
   expand.grid(
    generations = gen,
    k = real_world_k, 
    paternity_of_SR_males = average_paternity_SR_males, 
    freq_polyandry = real_world_data$freq_polyandry,
    w_STSR_female = 0.88, # Fitness from our Table of results
    w_SRSR_female = 0.42, # Fitness from our Table of results
    w_SR_male = c(0.8, 0.9, 1),
    initial_freq_SR = initial_freq,
    facet = "Recessive cost in females")
  
) %>% as_tibble()

observed_vs_expected_data <- 
  run_many_simulations(parameters) %>%
  left_join(real_world_data, by = "freq_polyandry") %>%
  rename(predicted_p_SR = prop_SR)

figure_4 <- observed_vs_expected_data %>%
  mutate(w_SR_male = paste("Fitness of SR males:", w_SR_male)) %>%
  ggplot(aes(real_freq_SR, 100 * predicted_p_SR, fill = 100*freq_polyandry)) +
  geom_abline(slope = 1, intercept=0, linetype=2) + 
  geom_point(pch=21, size = 3) + 
  scale_x_continuous(limits =c(0,35)) +
  scale_y_continuous(limits =c(0,35)) +
  scale_fill_distiller(palette = "Reds", direction = 1, name = "% polyandry") + 
  labs(x="Observed frequency of SR (%)",
       y="Predicted frequency of SR (%)") + 
  facet_wrap(~w_SR_male) + 
  theme_bw() + 
  theme(strip.background = element_blank(),
        text = element_text(family = "Lato", size = 12)
)

figure_4 %>% ggsave(filename = "figures/figure_4.pdf", height = 3.3, width = 8)
figure_4

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] showtext_0.5-1     showtextdb_2.0     sysfonts_0.7.2    
[4] ggplot2_3.1.0      kableExtra_0.9.0   future.apply_1.0.1
[7] future_1.11.1.1    dplyr_0.8.0.1      purrr_0.3.2       

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5   xfun_0.6           listenv_0.7.0     
 [4] colorspace_1.3-2   htmltools_0.3.6    viridisLite_0.3.0 
 [7] yaml_2.2.0         rlang_0.3.4        pillar_1.3.1.9000 
[10] glue_1.3.1.9000    withr_2.1.2        RColorBrewer_1.1-2
[13] plyr_1.8.4         stringr_1.4.0      munsell_0.5.0     
[16] gtable_0.2.0       workflowr_1.3.0    rvest_0.3.2       
[19] codetools_0.2-15   evaluate_0.13      labeling_0.3      
[22] knitr_1.22         parallel_3.5.1     curl_3.3          
[25] highr_0.8          Rcpp_1.0.1         readr_1.1.1       
[28] backports_1.1.2    scales_1.0.0       jsonlite_1.6      
[31] fs_1.3.1           hms_0.4.2          digest_0.6.18     
[34] stringi_1.4.3      rprojroot_1.3-2    grid_3.5.1        
[37] tools_3.5.1        magrittr_1.5       lazyeval_0.2.2    
[40] tibble_2.0.99.9000 crayon_1.3.4       pkgconfig_2.0.2   
[43] xml2_1.2.0         assertthat_0.2.1   rmarkdown_1.10    
[46] httr_1.4.0         rstudioapi_0.10    R6_2.4.0          
[49] globals_0.12.4     git2r_0.23.0       compiler_3.5.1