# 3.12.2009 # Simulation des 2-Typ Wright-Fisher-Modells # mit (gen-weiser) Selektion und Mutation WFstep <- function(N, s, nu1, nu2, p) { # simuliere eine neue Generation: # 2N Ind, aktueller Anteil a=p, # Sel.vorteil von Typ a=s/2N, Mut.rate a->A=nu1, Mut.rate A->a=nu2 ptilde<-(p*(1+s/(2*N))*(1-nu1/(2*N))+(1-p)*nu2/(2*N))/(1+p*s/(2*N)) rbinom(1, size=2*N, prob=ptilde)/(2*N) } # # Simulation und Darstellung von Pfaden: # # Diese Parameter gerne variieren: N <- 100 # (diploide) Populationsgroesse Tmax <- 1000 # Anzahl zu simulierender Generationen x0 <- 0.5 # Anteil Typ a-Individuen am Anfang nu1 <- 0.2 nu2 <- 0.3 s <- 2 anteil.a <- numeric(Tmax) anteil.a[1] <- x0 for (t in 2:Tmax) { anteil.a[t]<- WFstep(N, s, nu1, nu2, anteil.a[t-1]) } plot(1:Tmax, anteil.a, type="l", ylim=c(0,1), xlab="Zeit [Generationen]", ylab="Anteil Typ a", main=c("Simulation des WF-Modells", paste("2N=",2*N,", s=",s,", nu1=",nu1,", nu2=",nu2,sep=""))) # Zum Gleichgewicht: # # Empirische Verteilung der besuchten Zustaende zwischen T0 und T1 # darstellen # N <- 100 T0 <- 3000 T1 <- 30000 x0 <- 0.5 nu1 <- 0.2 nu2 <- 0.1 s <- 0 anteil.a <- numeric(T1); anteil.a[1] <- x0 for (t in 2:T1) { anteil.a[t]<- WFstep(N, s, nu1, nu2, anteil.a[t-1]) } hist(anteil.a[T0:T1], prob=TRUE,xlim=c(0,1), xlab="Anteil a", main=c(paste("Empir. Vert. des Anteils a zw. Gen.",T0,"und",T1), paste("2N=",2*N,", s=",s,", nu1=",nu1,", nu2=",nu2,sep=""))) f <- function(x) exp(2*s*x)*x^(2*nu2-1)*(1-x)^(2*nu1-1) fnormierung <- integrate(f, lower=1/N, upper=1-1/N) dichte <- function(x) f(x)/fnormierung$value curve(dichte, add=TRUE) # # Fixationsw'keit im Fall mit Selektion # per Simulation schaetzen # wofixiert <- function(x, s, N) { while (1) { if (x==0 || x==1) return(x) x<-WFstep(N, s, 0, 0, x) } } N <- 100 s <- 1 x0 <- 0.3 versuche <- 1000 empir.fix.wk <- mean(replicate(versuche, wofixiert(x0, s, N))) empir.fix.wk # Zum Vergleich das theoretische Limesergebnis der WF-Diffusion (1-exp(-2*s*x0))/(1-exp(-2*s))