# Load necessary packages
library(ggplot2); library(compiler); library(reshape2)
# Function to find the frequency of the SGE, polyandry and males. Just finds the SGE frequency if type = SGEonly (for speed)
freq = function(pop,type = "all"){
if(type == "SGEonly") return(sum(c(pop[c(3:4,7:8),"N"], pop[9:10,"N"]*2)) / sum(c(pop[1:4,"N"], pop[5:10,"N"]*2)))
freqA = sum(pop[seq(2,10,by=2),"N"]); freqSGE = sum(c(pop[c(3:4,7:8),"N"], pop[9:10,"N"]*2)) / sum(c(pop[1:4,"N"], pop[5:10,"N"]*2))
SR = sum(pop[1:4,"N"]) / sum(pop[,"N"]) # Sex ratio
return(c(freqA,freqSGE,SR))}; freq = cmpfun(freq)
# Initialise the population.
# freqP is the initial frequency of the polyandry allele
# freqSGE is the initial frequency of the X* chromosome
# Alleles are coded as follows: Sex chromosomes: #0=Y, 1=X, 2=X*; polyandry allele = 1, monandry allele = 0
init.SGE = function(freqP, freqSGE){
pop = cbind(N=1, as.matrix(expand.grid(P=0:1,sex.chrom=c(10,20,11,12,22))))
pop[,"N"] = 0.5* c((1-freqP)*(1-freqSGE), freqP*(1-freqSGE),
(1-freqP)*freqSGE, freqP*freqSGE,
(1-freqP)*(1-freqSGE)*(1-freqSGE), (freqP)*(1-freqSGE)*(1-freqSGE),
(1-freqP)*freqSGE*(1-freqSGE)*2, (freqP)*freqSGE*(1-freqSGE)*2,
(1-freqP)*freqSGE*freqSGE, (freqP)*freqSGE*freqSGE )
return(pop)}; init.SGE = cmpfun(init.SGE)
# Little function needed to save space in brood.comp function below.
seqer = function(x) seq(from = x, to = x + 18, by = 6); seqer = cmpfun(seqer)
# Calculate the composition of monandrous and polyandrous broods for all possible types of crosses
# C is the parameter C in the paper (i.e. competitiveness of X*Y males' sperm)
# D is the parameter C in the paper (i.e. the strength of meiotic drive) (d=1 total drive, 0.5 moderate drive, 0 no drive)
# Returns a list containing two dataframes: One containing the offspring genotype frequencies for all crosses involving a monandrous female,
# And another containing the offspring genotype frequencies for all crosses involving a polyandrous female
brood.comp = function(C,D){
D2 = c(0.5-D/2,0.5+D/2); D4 = c(rep(0.25-D/4,2),rep(0.25+D/4,2)); D8 = c(rep(0.125-D/8,4),rep(0.125+D/8,4))
# first matings
broods = matrix(0,nrow=24,ncol=10) # 24 cross types, 10 offspring types
broods[1,c(1,5)] = 0.5 # Cross 1: 0 10 male x 0 11 female = 0 10 (1), 0 11 (5)
broods[2,c(1,2,5,6)] = 0.25 # Cross 2: 0 10 male x 1 11 female = 0/1 10 (1:2), 0/1 11 (5:6)
broods[3,c(1,3,5,7)] = 0.25 # Cross 3: 0 10 male x 0 12 female = 0 10 (1), 0 20 (3), 0 11 (5), 0 12 (7)
broods[4,c(1:8)] = 0.125 # Cross 4: 0 10 male x 1 12 female = 0/1 10 (1:2), 0/1 20 (3:4), 0/1 11 (5:6), 0/1 12 (7:8)
broods[5,c(3, 7)] = 0.5 # Cross 5: 0 10 male x 0 22 female = 0 20 (3), 0 12 (7)
broods[6,c(3,4,7,8)] = 0.25 # Cross 6: 0 10 male x 1 22 female = 0/1 20 (3:4), 0/1 12 (7:8)
broods[7,c(1,2,5,6)] = 0.25 # Cross 7: 1 10 male x 0 11 female = 0/1 10 (1:2), 0/1 (5:6)
broods[8,c(2,6)] = 0.5 # Cross 8: 1 10 male x 1 11 female = 1 10 (2), 1 11 (6)
broods[9,c(1:8)] = 0.125 # Cross 9: 1 10 male x 0 12 female = 0/1 10 (1:2), 0/1 20 (3:4), 0/1 11 (5:6), 0/1 12 (7:8)
broods[10,c(2,4,6,8)] = 0.25 # Cross 10: 1 10 male x 1 12 female = 1 10 (2), 1 20 (4), 1 11 (6), 1 12 (8)
broods[11,c(3,4,7,8)] = 0.25 # Cross 11: 1 10 male x 0 22 female = 0/1 20 (3:4), 0/1 12 (7:8)
broods[12,c(4,8)] = 0.5 # Cross 12: 1 10 male x 1 22 female = 1 20 (4), 1 12 (8)
broods[13,c(1,7)] = D2 # Cross 13: 0 20 male x 0 11 female = 0 10 (1), 0 12 (7)
broods[14,c(1,2,7,8)] = D4 # Cross 14: 0 20 male x 1 11 female = 0/1 10 (1:2), 0/1 12 (7:8)
broods[15,c(1,3,7,9)] = D4 # Cross 15: 0 20 male x 0 12 female = 0 10 (1), 0 20 (3), 0 12 (7), 0 22 (9)
broods[16,c(1:4,7:10)] = D8 # Cross 16: 0 20 male x 1 12 female = 0/1 10 (1:2), 0/1 20 (3:4), 0/1 12 (7:8), 0/1 22 (9:10)
broods[17,c(3,9)] = D2 # Cross 17: 0 20 male x 0 22 female = 0 20 (3), 0 22 (9)
broods[18,c(3,4,9,10)] = D4 # Cross 18: 0 20 male x 1 22 female = 0/1 20 (3:4), 0/1 22 (9:10)
broods[19,c(1,2,7,8)] = D4 # Cross 19: 1 20 male x 0 11 female = 0/1 10 (1:2), 0/1 12 (7:8)
broods[20,c(2,8)] = D2 # Cross 20: 1 20 male x 1 11 female = 1 10 (2), 1 12 (8)
broods[21,c(1:4,7:10)] = D8 # Cross 21: 1 20 male x 0 12 female = 0/1 10 (1:2), 0/1 20 (3:4), 0/1 12 (7:8), 0/1 22 (9:10)
broods[22,c(2,4,8,10)] = D4 # Cross 22: 1 20 male x 1 12 female = 1 10 (2), 1 20 (4), 1 12 (8), 1 22 (10)
broods[23,c(3,4,9,10)] = D4 # Cross 23: 1 20 male x 0 22 female = 0/1 20 (3:4), 0/1 22 (9:10)
broods[24,c(4,10)] = D2 # Cross 24: 1 20 male x 1 22 female = 1 20 (4), 1 22 (10)
monandrous = broods[seq(1,23,by=2),]
# Second mating
polyandrous = matrix(0,nrow=48,ncol=10) # X* males lose out against X males
paternity = c(0.5,0.5, 1/(1+C), 1/(1+C)); paternity2 = c(C/(1+C), C/(1+C), 0.5,0.5) # First mating, then 4x possible second male
polyandrous[1:4,] = paternity*(broods[rep(2,4),]) + (1-paternity)*(broods[seqer(2),]) # 0 10 male x 1 11 female
polyandrous[5:8,] = paternity*(broods[rep(4,4),]) + (1-paternity)*(broods[seqer(4),]) # 0 10 male x 1 12 female
polyandrous[9:12,] = paternity*(broods[rep(6,4),]) + (1-paternity)*(broods[seqer(6),]) # 0 10 male x 1 22 female
polyandrous[13:16,] = paternity*(broods[rep(8,4),]) + (1-paternity)*(broods[seqer(2),]) # 1 10 male x 1 11 female
polyandrous[17:20,] = paternity*(broods[rep(10,4),]) + (1-paternity)*(broods[seqer(4),]) # 1 10 male x 1 12 female
polyandrous[21:24,] = paternity*(broods[rep(12,4),]) + (1-paternity)*(broods[seqer(6),]) # 1 10 male x 1 22 female
polyandrous[25:28,] = paternity2*(broods[rep(14,4),]) + (1-paternity2)*(broods[seqer(2),]) # 0 20 male x 1 11 female
polyandrous[29:32,] = paternity2*(broods[rep(16,4),]) + (1-paternity2)*(broods[seqer(4),]) # 0 20 male x 1 12 female
polyandrous[33:36,] = paternity2*(broods[rep(18,4),]) + (1-paternity2)*(broods[seqer(6),]) # 0 20 male x 1 22 female
polyandrous[37:40,] = paternity2*(broods[rep(20,4),]) + (1-paternity2)*(broods[seqer(2),]) # 1 20 male x 1 11 female
polyandrous[41:44,] = paternity2*(broods[rep(22,4),]) + (1-paternity2)*(broods[seqer(4),]) # 1 20 male x 1 12 female
polyandrous[45:48,] = paternity2*(broods[rep(24,4),]) + (1-paternity2)*(broods[seqer(6),]) # 1 20 male x 1 22 female
return(list(monandrous,polyandrous))}; brood.comp = cmpfun(brood.comp)
# Mate the population
# Arguments: a dataframe containing the genotype frequencies of the current population (pop)
# The monandrous and polyandrous crosses genotype frequencies (made with the brood.comp function)
# The cost of polyandry 'cost.poly' (called p in the paper)
# The cost to X*X* homozygote females 'cost.X.fem' (called h in the paper)
# The effect of X* on female fecundity in heterozygote females (confusingly called h here - never actually implemented or discussed in the paper, and set to 0 in all simulations, meaning that heterozygous (XX*) females have the same fitness as XX females throughout. Note that this h is not the same as the h in the paper!!)
mating.sge = function (pop, monandrous, polyandrous,cost.poly,cost.X.fem, h){
male.freqs <- pop[1:4, "N"] / sum(pop[1:4, "N"]) # Assumes that all females get fertilised, irrespective of male frequency
freq.mono.crosses <- expand.grid(f=pop[c(5,7,9),"N"],m=male.freqs)
freq.mono.crosses <- freq.mono.crosses[,1]*freq.mono.crosses[,2]*rep(c(1, 1-cost.X.fem*h, 1-cost.X.fem),4) # females with one or two X* might pay a fecundity cost
freq.poly.crosses <- expand.grid(m2=male.freqs, f=pop[c(6,8,10),"N"], m1=male.freqs)
freq.poly.crosses <- freq.poly.crosses[,1]*freq.poly.crosses[,2]*freq.poly.crosses[,3]*(1-cost.poly)*rep(rep(c(1,1-cost.X.fem*h,1-cost.X.fem),each=4),4) # fecundity selection against polyandrous females
pop[,"N"] <- colSums(rbind(monandrous * freq.mono.crosses, polyandrous * freq.poly.crosses)) # new pop REPLACES old one
pop[,"N"] <- pop[,"N"] /sum(pop[,"N"]) # convert back to frequencies
return(pop) }; mating.sge = cmpfun(mating.sge)
# Test function - this simply makes a human-readable table of all the genotype frequencies produced by all the possible mating combinations for a given combination of C and D
# Not used in the simulation, but useful for checking for bugs in coding the rather complex mating functions
brood.calculator = function(C, D){
broods = brood.comp(C=C,D=D)
broods = rbind(broods[[1]], broods[[2]])
firstM = c(rep("0 10", 3), rep("1 10", 3), rep("0 20", 3), rep("1 20", 3), rep("0 10", 12), rep("1 10", 12), rep("0 20", 12), rep("1 20", 12))
secondM = c(rep("None", 12), rep(c("0 10","1 10","0 20","1 20"), 12))
female = c(rep(c("0 11","0 12","0 22"), 4), rep(rep(c("1 11","1 12","1 22"), each=4), 4))
broods<-data.frame(broods,firstM,secondM,female)
names(broods)[1:10] <- c("0 10","1 10","0 20","1 20", "0 11", "1 11", "0 12", "1 12", "0 22", "1 22")
return(broods)}
# Function that initialises the population, and then mates it for 'gen' generations
# If ploty = TRUE, the function makes a plot of genotype (and male) frequencies over time
sge.sim = function(gen, freqP, freqSGE, C,D, cost.poly,cost.X.fem, h, ploty=F){
pop = init.SGE(freqP=freqP, freqSGE=freqSGE); broods = brood.comp(C=C,D=D)
if(ploty==T) { freqs = matrix(0,ncol=3,nrow=gen+1); freqs[1,] = freq(pop)
for(ii in 1:gen) {pop = mating.sge(pop=pop, monandrous=broods[[1]], polyandrous=broods[[2]],cost.poly=cost.poly,cost.X.fem=cost.X.fem, h=h)
freqs[ii+1,] = freq(pop)}
freqs = data.frame(0:gen, freqs); names(freqs) = c("Generation","Polyandry","X*","Males"); freqs <- melt(freqs, id="Generation")
p=ggplot(data=freqs,aes(x=Generation, y=value, colour=variable)) + geom_line(aes(linetype=variable)) + theme_bw(base_size=16) +labs(y="Frequency") +theme(legend.position = "none", panel.grid.major=element_blank(), panel.grid.minor=element_blank())
print("Final frequency of polyandry, X* and males respectively: ")
print(freq(pop))
return(p)}
else for(ii in 1:gen) pop = mating.sge(pop=pop, monandrous=broods[[1]], polyandrous=broods[[2]],cost.poly=cost.poly,cost.X.fem=cost.X.fem, h=h)
return(freq(pop))};sge.sim = cmpfun(sge.sim)
# PLAY WITH THE SIMULATION HERE by adjusting the model's parameters
# Red line = frequency of polyandry allele, blue line = frequency of males, green line = frequency of X*
sge.sim(gen=1000, freqP=.3, freqSGE=.001, C=0.4, D=1, cost.poly=0.02,cost.X.fem=1, h=0, ploty=T)