#Compute the expected change in frequency of a risky versus certain reproductive strategy #assuming either fixed population size or variable population size #generate number of offspring produced by parents with "risky" genotype risky <- function(parents = 8){2*rbinom(1, parents, 0.5)} deltap <- fixmut <- function(s = 0.01, N = 10,...) { #compute and store all post-selection frequencies p.star <- rep(0,2*N) for(j in 1:(2*N-1)){ p.star[j] <- j*(1+s)/(j*s + 2*N) } #generate entries of the matrix Q #use the R function that computes the binomial probability density Q <- matrix(0, nrow = (2*N-1), ncol = (2*N-1)) for(i in 1:(2*N-1)){ for(j in 1:(2*N-1)) { Q[i, j] <- dbinom(j, 2*N, p.star[i]) } } #generate entries fo the vector r (probs. of fixation in one step) r <- rep(0,(2*N-1)) for(i in 1:(2*N-1)) { r[i] <- dbinom(2*N, 2*N, p.star[i]) } #compute u, vector of ultimate fixation probabilities u <- rep(0,(2*N-1)) I <- diag(2*N-1) #the identity matrix u <- solve(I - Q, r) u[1] #output the fixation probability of a new, unique mutation } ################################ ############## 2N = 4 ########## #neutral mutation neut <- fixmut(s = 0, N = 2) #should be 1/(2*N) neut #advantageous mutation with s = 0.01 fixmut(s = 0.01, N = 2) fixmut(s = 0.01, N = 2)/neut #fixation prob relative to that of a neutral mutation #deleterious mutation with s = -0.01 fixmut(s = -0.01, N = 2) fixmut(s = -0.01, N = 2)/neut #fixation prob relative to that of a neutral mutation ############## 2N = 40 ########## #neutral mutation neut <- fixmut(s = 0, N = 20) #should be 1/(2*N) neut #advantageous mutation with s = 0.01 fixmut(s = 0.01, N = 20) fixmut(s = 0.01, N = 20)/neut #fixation prob relative to that of a neutral mutation #deleterious mutation with s = -0.01 fixmut(s = -0.01, N = 20) fixmut(s = -0.01, N = 20)/neut #fixation prob relative to that of a neutral mutation ############## 2N = 400 ########## #neutral mutation neut <- fixmut(s = 0, N = 200) #should be 1/(2*N) neut #advantageous mutation with s = 0.01 fixmut(s = 0.01, N = 200) fixmut(s = 0.01, N = 200)/neut #fixation prob relative to that of a neutral mutation #deleterious mutation with s = -0.01 fixmut(s = -0.01, N = 200) fixmut(s = -0.01, N = 200)/neut #fixation prob relative to that of a neutral mutation