Erwartete Zeit für ein neutrales Allel, um eine Frequenz von p1p1p_1 zu erreichen, wenn es bei der Frequenz p0p0p_0 beginnt

Einfache englische Antwort:

Ich habe eine Computerfunktion geschrieben, die neutrale Evolution simuliert, um dieses Problem zu lösen. Es ist keine exakte mathematische Antwort, aber es ist im Grunde der gleiche Ansatz, den Kimura und Ohta in (der zweiten Hälfte) ihrer Arbeit verfolgten, außer dass mein Computer leistungsfähiger ist als ihr, sodass ich durch Simulation viel genauere Schätzungen erhalten konnte mehr Bevölkerung als sie taten.

Abbildung 1 ist ein Gitterdiagramm der Beziehung zwischen P0 und der erwarteten Zeit bis zum Erreichen von P1 für unterschiedliche Werte von P1 (nach Spalte) und unterschiedliche Populationsgrößen (nach Zeile). Es ist klar, dass die gleichen Beziehungen zwischen P0, P1 und der erwarteten Zeit, um von P0 nach P1 zu gelangen, in jeder Populationsgröße zu sehen sind, aber bei größeren Populationen sollten Sie mit einer längeren Zeit rechnen, um von P0 nach P1 zu gelangen.

Wenn P0 und P1 nahe beieinander liegen, erwarten Sie im Allgemeinen eine kürzere Zeit, um von P0 zu P1 zu gelangen, als wenn P0 und P1 weit voneinander entfernt sind. Man könnte sich das in der Regel vorstellen wie: „Es dauert länger, bis man weiterfährt“. Die erwartete Zeit, die benötigt wird, um von P0 zu P1 zu gelangen, ist kleiner, wenn P0 und P1 auf derselben Seite von 0,5 liegen und P1 weiter von 0,5 entfernt ist als P0, als wenn P0 weiter von 0,5 entfernt ist als P1. Sie können sich dies als Regel vorstellen wie: „Es ist schneller zu reisen, wenn Sie näher an der Fixierung oder Auslöschung sind, als wenn Sie sich in den mittleren Frequenzen befinden“.

Abbildung 1: Erwartete Zeit für ein neutrales Allel, um von Frequenz P0 zu Frequenz P1 zu wechseln, bei gegebenen Unterschieden in P0, subgrafisch dargestellt nach P1 und Populationsgröße. Modellierte Frequenzen von P0 und P1 sind 0,1, 0,3, 0,5, 0,7 und 0,9. Die modellierten Populationsgrößen sind N = 10, N = 50 und N = 100. Jede Unterdarstellung stellt eine Kombination aus P1 und Populationsgröße dar, wobei P1 von links nach rechts zunimmt und die Populationsgröße von unten nach oben zunimmt. Alle Erwartungen werden anhand von 10.000 simulierten Populationen geschätzt.

Einige Populationen brauchen länger als andere, um von P0 zu P1 zu gelangen. Im Allgemeinen brauchen die meisten Populationen eine kleine Anzahl von Generationen, um P1 zu erreichen, und eine kleinere Anzahl braucht lange – aber sehr wenige können sehr, sehr lange dauern.

Ich habe meinen Code angehängt. Wenn Sie also wissen, wie man R verwendet, können Sie damit die erwartete Zeit schätzen, um von P0 nach P1 für eine beliebige Kombination aus P0, P1 und Bevölkerungsgröße zu gelangen.

Zusätzliche Details und Annahmen für technisch Interessierte:

Relevante Theorie:

In einer diploiden, sich sexuell fortpflanzenden Population von großer Größe$N$, es gibt$2N$Kopien eines bestimmten Locus. Jeder Locus ist von einem Allel besetzt. Alle Variationen an unserem interessierenden Ort sind selektiv neutral. Für unsere Zwecke können wir davon ausgehen, dass jeder Locus entweder von unserem interessierenden Allel ($A$) oder eine andere Variante ($a$), wobei Variationen in den „anderen“ Allelen ignoriert werden. Die Anzahl der Kopien von$A$in der Bevölkerung zu Beginn des Prozesses,$A | t=0$, ist gegeben durch$2N*P0$.

Nehmen wir an, dass die Populationsgröße konstant ist ($Nt+1 = Nt$für alle$t$), und diese Paarung ist zufällig. Nehmen wir weiter an, dass die Generationen verschieden sind. Lassen$t = 1$die erste Generation sein, die geboren wird, und so weiter.

Die erwartete Anzahl von Kopien von$A$zum Zeitpunkt$t+1$dann folgt eine Binomialverteilung:

$At+1$wird ausgeliefert$Binom(2N, At/(2N))$

Da haben wir eine Übergangsregelung für$At → At+1$, ist es ziemlich einfach, Populationen zu simulieren, die diese Form der Evolution durchlaufen. Ich habe die R-Funktion neutralFPT geschrieben, um diese Simulation durchzuführen und den Anteil aller Populationen abzuschätzen, die erreichen$P1$irgendwann in ihrer Geschichte die erwartete Zeit für das Erreichen eines neutralen Allels$P1$aus$P0$, da es erreichen wird$P1$, und die Verteilung der Zeitdauer bis zum Erreichen$P1$, da eine Bevölkerung erreichen wird$P1$. Das Skript finden Sie im letzten Abschnitt dieser Antwort.

Wahrscheinlichkeitsdichten der First-Passage-Zeiten:

Wahrscheinlichkeitsdichten von First-Passage-Zeiten über plausiblen Werten von$P0$,$P1$, und$N$folgen einer ähnlichen Struktur – unimodal in der Nähe des linken Randes, mit einem langen rechten Schwanz mit längeren First-Passage-Zeiten (Abbildung 2).

Abbildung 2: Wahrscheinlichkeitsdichten der First-Passage-Zeiten von P0 nach P1 bei unterschiedlichen Werten von P0, P1 und N. A: P0 = 0,5, P1 = 0,9, N = 50; B: P0 = 0,9, P1 = 0,5, N = 50; C: P0 = 0,5, P1 = 0,9, N = 500.

Funktion: neutralFPT, mit Anwendungsbeispielen.

 ## this function simulates the number of generations taken for a
 # neutral allele to go from a starting proportion, P0, to a final
 # proportion, P1, in a random population, given that it will at some
 # point reach P1. Allele frequencies are not modelled after the point
 # when P1 is reached.
## neutralFPT was generated in response to Stack Exchange user Remi.b, at
 # http://biology.stackexchange.com/questions/30812/expected-time-for-a-neutral-allele-to-reach-a-frequency-of-p-1-when-starting-a
## neutralFPT was written by Shane Baylis, 2015
 # for R version 3.2.2

neutralFPT <- function(P0, P1, N, niter) {

tOut <- c(rep(NaN, niter))
 # create a vector of blank t-values
statOut <- c(rep(NaN, niter))
 # create a vector of population status values (reached P1, or didn't)

if(P0 == P1) stop("P0 and P1 are set to the same value!")
if(P0 == 0 | P0 == 1) stop("P0 is set to zero or one, so its frequency
                           cannot change!")
if(P0 < 0 | P0 > 1) stop("P0 must be between zero and one")
if(P1 < 0 | P1 > 1) stop("P1 must be between zero and one")
## work out whether you're heading upwards or downwards
if(P1 > P0) { # i.e., our target is above us
    for (i in 1:niter) {
        NAllele <- round(2*(P0*N))
        Target <- round(2*(P1*N))
        t <- 0
        while (NAllele < Target && NAllele != 0 && NAllele != 2*N) {
            t <- t+1
            NAllele <- rbinom(1, 2*N, (NAllele/(2*N)))
        }
        if(NAllele >= Target) {
            statOut[i] <- 1 ## 1 indicates that P1 occurred
            tOut[i] <- t
        }else{
            statOut[i] <- 0
            tOut[i] <- Inf
        }
    }
}else{ ## i.e., our target is below us
    for (i in 1:niter) {
       NAllele <- round(2*(P0*N))
       Target <- round(2*(P1*N))
       t <- 0
       while (NAllele > Target && NAllele != 0 && NAllele != 2*N) {
           t <- t+1
           NAllele <- rbinom(1, 2*N, (NAllele/(2*N)))
       }
       if(NAllele <= Target) {
           statOut[i] <- 1 ## 1 indicates that P1 occurred
           tOut[i] <- t
       }else{
           statOut[i] <- 0
           tOut[i] <- Inf
       }
   }
}
successes <- sum(statOut) # the number of populations in which
                          # P1 was reached
propSuccesses <- successes / niter # the proportion of populations
                                   # in which P1 was reached
successTimes <- subset(tOut, statOut == 1)
expectedFPT <- mean(successTimes)
medianFPT <- median(successTimes)
outs <- list(successes=successes, propSuccesses=propSuccesses,
             successTimes=successTimes,expectedFPT=expectedFPT,
             medianFPT=medianFPT, trials=niter)
return(outs)
} # function close       

## neutralFPT examples #################################################

sim <- neutralFPT(0.5, 0.9, 500, 10000)
sim$expectedFPT # Numeric. shows the expected (i.e., mean) first-passage
                # time from P0 to P1, in generations, given that a
                # population reached P1.
sim$medianFPT # Integer. Shows the median first-passage time from P0 to P1,
              # in generations, given that a population reached P1.
sim$propSuccesses # Integer. The proportion of simulated populations which
                  # reached P1. Other populations reached fixation
                  # or extinction of A without reaching P1.
sim$successTimes # Vector. The number of generations taken to reach P1,
                 # for all populations which reached P1.
hist(sim$successTimes, xlab="First passage time (generations)",
     main=paste(sim$successes, "successes from", sim$trials,
         "populations"))  ## histogram of first-passage times

PZero <- c(rep(c(0.1, 0.3, 0.5, 0.7, 0.9), 15))
POne <- c(rep(c(rep(0.1,5),rep(0.3,5),rep(0.5,5),rep(0.7,5),rep(0.9,5)),3))
PopSize <- c(rep(10,25),rep(50,25),rep(100,25))
FPTs <- c(rep(NaN, length(starts)))
testFrame <- data.frame(PZero, POne, PopSize, FPTs)
testFrame <- subset(testFrame, starts != ends)
for(s in 1:nrow(testFrame)) {
    sim <- with(testFrame, neutralFPT(PZero[s],POne[s],PopSize[s],10000))
    testFrame$FPTs[s] <- sim$expectedFPT
}  ## estimates expected first passage times from P0 to P1 for a variety
    # of values of P0, P1, and population size. Outputs the estimates to
    # a table called testFrame.

require(lattice)
with(testFrame, dotplot(FPTs~PZero|POne*PopSize,
   main="Expected First-Passage Time by P0, P1, and Population Size",
   xlab="starting frequency (P0)")) ## generates the lattice plot used as
                                     # Figure One.