mating <- function(NP, SR, Nov, EMs, Chy){ ### Initialize the population ### Nf <- as.integer(NP * (1-SR)) #number of females Nm <- NP - Nf #number of males males <- matrix(0, nrow=Nm, ncol=2) #each male is defined by a normally-distributed trait that regulates its mating success; the number of mates is stored in the second column of this matrix males[,1] <- rnorm(Nm) # the trait value for each male is stored in the first column of the matrix females <- matrix(NA, nrow=Nf, ncol=Nov) #each female is defined by the number of ovules she produces; the identity of the father of each seed is stored in the cells of the matrix ### Mating ### for(i in 1:Nf){ for(j in 1:Nov){ # randomly choose males (with replacement) to act as sires of each seed; probability of mating depends on the trait value for each male while(is.na(females[i,j])){ #every seed gets fertilized if(Chy == 0){ lek <- sample.int(Nm, EMs) #females encounter EMs males and choose the one with the highest trait value choice <- which.max(males[,1][lek]) #identifies which of the sampled males has the highest trait value females[i,j] <- lek[choice] # female mates with that male chosenmale <- lek[choice] males[chosenmale, 2] <- males[chosenmale,2] + 1 # that male gets an (additional) mating } else { # note that if Chy != 0 females choose males in direct proportion to their trait value, with the value of Chy defining the minimum trait value females will accept attract <- rnorm(1) # females have a threshold for acceptance for each mating attempt and will only accept males from the upper half of the trait distribution choice <- sample.int(Nm, 1) #males are sampled with replacement until an accepatable mate is found if(males[choice,1] >= attract + Chy){ #females mate with males that meet her threshold for acceptance; this procedure chooses males as linear function of their trait value; see the plot at the end of "main.R", although the variance increases with the mean... perhaps this is realistic however. females[i,j] <- choice males[choice, 2] <- males[choice,2] + 1 } } } } } #loops to count the number of full sibs per seed family, as in Fig. 1 b <- rep(0, length.out=length(females[,1])) c <- rep(0, length.out=length(females[,1])) for(k in 1:length(females[,1])){ for(i in 1:length(females[1,])){ for(j in 1:length(females[1,])){ if(j > i){ c[k] <- c[k] + 1 if(females[k,i] == females[k,j]) b[k] <- b[k] + 1 else b[k] <- b[k] } } } } rp <- b/c # calculate the average correlated paternity per population to return it to main rp <- mean(rp) # calculate the selection differential; equivalent to the standardized univariate selection gradient; calculated here using the Robertson-Price identity to account for differences in siring success among males with different phenotypes meansiring <- mean(males[,2]) relfit <- males[,2]/meansiring seldif <- as.numeric(cov(scale(males[,1]), relfit)) # calculate Is # from Klug et al. (2010): The opportunity for sexual selection (Is): A standardized measure of intra-sexual variation in mating success; measured as the square of the coefficient of variation in mating success for a given sex. Is <- (sd(males[,2])/mean(males[,2]))^2 ### return calculated values ### return(list(rp = rp, males = males, Is = Is, sg = seldif)) }