############### Individual based model to simulate sex specific arrival times ############### ########## ########## ########## This model is based on a code written by Hannah Kokko ########## ########## and described in: ########## ########## Kokko et al. (2006) Why do female migratory ########## ########## birds arrive later than males? Journal of Animal Ecology. 75 ########## ########## ########## ############################################################################################### ########### Authors: Simeon Lisovski & Markus Ritz ############################################ ############################################################################################### ## Parameters: ## NM: Number of South Polar Skua males ## NF: Number of South Polar Skua females ## alpha: Male mate choice parameter (1 = no preference, > 1 prevalence towards Brown skua females) ## EPCprob: Extra pair paternity rate ## sigma: ratio of allelic values in hybrids (1 = all hybrids inherit allelic value from South Polar Skua father) ## mutprob: Allel mutation probability ## mortality: mortality rate ## NrT: numbers of Territories ## QT: maximal quality of each terrirory ## plot: logical, if TRUE a plot with the arrival times will be produced skua.IBM <- function(NM, NF, BS, alpha, EPCprob, sogma, mutprob, mortality, NrT, QT, plot=TRUE) { # indices # index=1; a=2; b=3; arrtime=4; aEPC=5; bEPC=6 # Females and Males contain: # index; genes a and b; arrtime Fe <- matrix(c(1:NF, 40*runif(NF*2)), ncol=3, nrow=NF) Ma <- matrix(c(1:NM, 40*runif(NM*2)), ncol=3, nrow=NM) BS <- matrix(ncol=2, nrow = BS) if(nrow(BS)>0) BS[,1] <- 1:nrow(BS) Tr <- matrix(c(sort(round(runif(NrT, 1, QT)), decreasing=T), rep(NA, 2*NrT)), nrow=NrT, ncol=3) Fe <- cbind(Fe, matrix(ncol=3, nrow=NF)) Ma <- cbind(Ma, matrix(ncol=1, nrow=NM)) data <- cbind(NA, NA, NA) generation <- 0 repeat{ Fe[ , 4] <- NA Ma[ , 4] <- NA Fe[ , c(5,6)] <- NA if(!is.null(BS)) BS[, 2] <- NA if(nrow(Fe)>NF) Fe <- Fe[1:NF, ] if(nrow(Ma)>NM) Ma <- Ma[1:NM, ] Tr[ , 2] <- NA Tr[ , 3] <- NA mcurve <- rep(NA, 40) fcurve <- rep(NA, 40) for( t in 0:40) { # migration: everyone in state 0 moves up if they so wish move <- which(is.na(Ma[ , 4]) & runif(nrow(Ma)) < (1/ (1 + exp( -2 * (t - Ma[ , 2]))))) if(length(move)>0) Ma[move, 4] <- t move <- which(is.na(Fe[ , 4]) & runif(nrow(Fe)) < (1/ (1 + exp( -2 * (t - Fe[ , 3]))))) if(length(move)>0) Fe[move, 4] <- t # females who have arrived but not yet chosen an EPC mate do so epc <- which(!is.na(Fe[ , 4]) & is.na(Fe[ , 5])) if(length(epc)>0) { f <- which(!is.na(Ma[ , 4])) if(length(f)>0) { f <- f[sample(length(f),length(epc), replace=T)] Fe[epc, c(5,6)] <- Ma[f, c(2,3)] } } mcurve[t+1] <- mean(Ma[, 4], na.rm=T) fcurve[t+1] <- mean(Fe[, 4], na.rm=T) # some mortality malemort <- exp(-mortality*t) dead <- which(!is.na(Ma[, 4]) & runif(nrow(Ma))0) { Ma <- Ma[-dead, ] ind <- setdiff(Tr[, 2], Ma[, 1]) if(length(ind)>0) Tr[which(Tr[,2]%in%ind), 2] <- NA ind <- setdiff(BS[,2], Ma[, 1]) if(length(ind)>0) BS[which(BS[,2]%in%ind), 2] <- NA } femalemort <- exp(-mortality*t) dead <- which(!is.na(Fe[, 4]) & runif(nrow(Fe))0) { Fe <- Fe[-dead, ] ind <- setdiff(Tr[, 3], Fe[, 1]) if(length(ind)>0) Tr[which(Tr[,3]%in%ind), 3] <- NA } # males without a territory look for single BS females f1 <- setdiff(Ma[!is.na(Ma[, 4]), 1], c(Tr[, 2], BS[, 2])) f2BS <- which(is.na(BS[, 2])) f2SPS <- which(is.na(Tr[, 2]) & !is.na(Tr[, 3])) if(length(f1)>0 & (length(f2BS)>0 | length(f2SPS)>0)) f1 <- f1[runif(length(f1))<((alpha/2)^(log(nrow(BS)/(nrow(BS)+nrow(Fe)), 0.5)))] # if(length(f1)>0 & (length(f2BS)>0 | length(f2SPS)>0)) f1 <- f1[runif(length(f1))0 & length(f2BS)>0){ if(length(f1)==length(f2BS)) BS[f2BS, 2] <- f1 if(length(f1) > length(f2BS)) BS[f2BS, 2] <- f1[1:length(f2BS)] if(length(f1) < length(f2BS)) BS[f2BS[1:length(f1)], 2] <- f1 } # males without a territory look for single territorial females f1 <- setdiff(Ma[!is.na(Ma[, 4]), 1], c(Tr[, 2], BS[,2])) if(any(is.na(f1))) f1 <- f1[-which(is.na(f1))] f2 <- which(!is.na(Tr[, 3]) & is.na(Tr[, 2])) if(length(f1)>0 & length(f2)>0) { if(length(f1)==length(f2)) Tr[f2, 2] <- f1 if(length(f1) > length(f2)) Tr[f2, 2] <- f1[1:length(f2)] if(length(f1) < length(f2)) Tr[f2[1:length(f1)], 2] <- f1 } # those males without a territory look for territories f1 <- setdiff(Ma[!is.na(Ma[, 4]), 1], c(Tr[, 2], BS[,2])) if(any(is.na(f1))) f1 <- f1[-which(is.na(f1))] f2 <- which(is.na(Tr[, 2])) if(length(f1)>0 & length(f2)>0) { if(length(f1)==length(f2)) Tr[f2, 2] <- f1 if(length(f1) > length(f2)) Tr[f2, 2] <- f1[1:length(f2)] if(length(f1) < length(f2)) Tr[f2[1:length(f1)], 2] <- f1 } # females without a territory look for single territorial males f1 <- setdiff(Fe[!is.na(Fe[, 4]), 1], Tr[, 3]) if(any(is.na(f1))) f1 <- f1[-which(is.na(f1))] f2 <- which(!is.na(Tr[, 2]) & is.na(Tr[, 3])) if(length(f1)>0 & length(f2)>0) { if(length(f1)==length(f2)) Tr[f2, 3] <- f1 if(length(f1) > length(f2)) Tr[f2, 3] <- f1[1:length(f2)] if(length(f1) < length(f2)) Tr[f2[1:length(f1)], 3] <- f1 } # those females without a territory look for territories f1 <- setdiff(Fe[!is.na(Fe[, 4]), 1], Tr[, 3]) if(any(is.na(f1))) f1 <- f1[-which(is.na(f1))] f2 <- which(is.na(Tr[, 3])) if(length(f1)>0 & length(f2)>0) { if(length(f1)==length(f2)) Tr[f2, 3] <- f1 if(length(f1) > length(f2)) Tr[f2, 3] <- f1[1:length(f2)] if(length(f1) < length(f2)) Tr[f2[1:length(f1)], 3] <- f1 } } # end of migration data <- rbind(data, cbind(mean(Ma[, 2]), mean(Fe[, 3]), sum(!is.na(BS[,2])))) if(plot==TRUE){ if(generation/2 == floor(generation/2)) { plot(0:40, mcurve, ylim=c(0, max(c(mcurve, fcurve), na.rm=T)), type="l", col="blue", xlab="days", ylab="mean arrival") lines(0:40, fcurve, col="red") } } newmale <- matrix(ncol=4) newfemale <- matrix(ncol=4) for(i in 1:nrow(Tr)) { if(!is.na(Tr[i, 2]) & !is.na(Tr[i, 3]) & Tr[i, 1]>0) { f1 <- which(Ma[, 1]==Tr[i, 2]) f2 <- which(Fe[, 1]==Tr[i, 3]) # with certain probability some of male offspring are extra-pair... newest <- matrix(rep(Ma[f1,], each=Tr[i, 1]), nrow=Tr[i, 1]) # with 0.5 probability male genes are inherited from mum allelfrommom <- matrix(runif(nrow(newest)*2)<0.5, ncol=2) newest[allelfrommom[, 1], 2] <- Fe[f2, 2] newest[allelfrommom[, 2], 3] <- Fe[f2, 3] EPC <- which(runif(nrow(newest))0) newest[EPC, c(2,3)] <- matrix(rep(Fe[f2, c(5, 6)], each=length(EPC)), ncol=2) newmale <- rbind(newmale, newest) # with certain probability some of female offspring are extra-pair... newest <- matrix(rep(Ma[f1,], each=Tr[i, 1]), nrow=Tr[i, 1]) # EPC <- which(runif(nrow(newest))0) newest[EPC, c(2,3)] <- matrix(rep(Fe[f2, c(5, 6)], each=length(EPC)), ncol=2) # with 0.5 probability male genes are inherited from mum allelfrommom <- matrix(runif(nrow(newest)*2)<0.5, ncol=2) newest[allelfrommom[, 1], 2] <- Fe[f2, 2] newest[allelfrommom[, 2], 3] <- Fe[f2, 3] EPC <- which(runif(nrow(newest))0) newest[EPC, c(2,3)] <- matrix(rep(Fe[f2, c(5, 6)], each=length(EPC)), ncol=2) newfemale <- rbind(newfemale, newest) } } # hybrids if(nrow(BS)>0) { for(i in 1:nrow(BS)) { if(!is.na(BS[i, 2])) { ratio <- round(QT*BSxSPS_ratio, 0) if(ratio>0) { f1 <- which(Ma[, 1]==BS[i, 2]) newest <- matrix(rep(Ma[f1, ], each=ratio), nrow=ratio) newmale <- rbind(newmale, newest) newest <- matrix(rep(Ma[f1,], each=ratio), nrow=ratio) newfemale <- rbind(newfemale, newest) } } } newmale <- newmale[-1, ] newfemale <- cbind(newfemale[-1, ], matrix(ncol=2, nrow=nrow(newfemale)-1)) if(nrow(newmale)>0){ newmale[ , 1] <- seq(max(Ma[,1])+1, max(Ma[,1])+nrow(newmale)) mut <- matrix(runif(nrow(newmale)*2)0){ newfemale[ , 1] <- seq(max(Fe[,1])+1, max(Fe[,1])+nrow(newfemale)) mut <- matrix(runif(nrow(newfemale)*2)=50){ mod <- lm(c(data[(nrow(data)-9):(nrow(data)),2] - data[(nrow(data)-9):(nrow(data)),1])~c(1:10)) if(coef(mod)[2]>=(-0.005) & coef(mod)[2]<=0.005) STOP <- TRUE } if(STOP) break } return(data[,2]-data[,1]) }