Tiempo esperado para que un alelo neutral alcance una frecuencia de p1p1p_1 cuando comienza en la frecuencia p0p0p_0

Respuesta en inglés sencillo:

He escrito una función de computadora que simula la evolución neutral para resolver este problema. No es una respuesta matemática exacta, pero es básicamente el mismo enfoque que Kimura y Ohta tomaron en (la segunda mitad de) su artículo, excepto que mi computadora es más poderosa que la de ellos, por lo que podría obtener estimaciones mucho más precisas al simular más poblaciones que ellos.

La figura uno es un diagrama de celosía de la relación entre P0 y el tiempo esperado para llegar a P1, para diferentes valores de P1 (por columna) y diferentes tamaños de población (por fila). Está claro que las mismas relaciones entre P0, P1 y el tiempo esperado para llegar de P0 a P1 se observan en todos los tamaños de población, pero con poblaciones más grandes, debe esperar un tiempo más largo para llegar de P0 a P1.

Cuando P0 y P1 están muy juntos, por lo general se espera un tiempo menor para llegar de P0 a P1 que cuando P0 y P1 están muy separados. Puede pensar en esto como una regla como 'se necesita más tiempo para viajar más lejos'. El tiempo esperado para ir de P0 a P1 es menor cuando P0 y P1 están en el mismo lado de 0,5 y P1 está más lejos de 0,5 que P0, que cuando P0 está más lejos de 0,5 que P1. Puede pensar en esto como una regla como 'es más rápido viajar cuando está más cerca de la fijación o extinción que cuando está en las frecuencias medias'.

Figura uno: tiempo esperado para que un alelo neutral pase de la frecuencia P0 a la frecuencia P1, dadas las diferencias en P0, subgraficado por P1 y el tamaño de la población. Las frecuencias modeladas de P0 y P1 son 0,1, 0,3, 0,5, 0,7 y 0,9. Los tamaños de población modelados son N = 10, N = 50 y N = 100. Cada subparcela representa una combinación de P1 y el tamaño de la población, con P1 aumentando de izquierda a derecha y el tamaño de la población aumentando de abajo hacia arriba. Todas las expectativas se estiman a partir de 10.000 poblaciones simuladas.

Algunas poblaciones tardan más que otras en pasar de P0 a P1. En general, la mayoría de las poblaciones tardan un pequeño número de generaciones en llegar a P1, y un número menor tarda mucho tiempo, pero muy pocas pueden tardar mucho tiempo.

Adjunté mi código, por lo que si sabe cómo usar R, puede usarlo para estimar el tiempo esperado para llegar de P0 a P1 para cualquier combinación de P0, P1 y tamaño de la población.

Detalles adicionales y suposiciones para los técnicamente inclinados:

Teoría relevante:

En una población diploide con reproducción sexual de tamaño$N$, existen$2N$copias de un locus dado. Cada locus está ocupado por algún alelo. Toda variación en nuestro lugar de interés es selectivamente neutral. Para nuestros propósitos, podemos considerar que cada locus está ocupado por nuestro alelo de interés ($A$), u otra variante ($a$), ignorando la variación en los 'otros' alelos. El número de copias de$A$en la población al inicio del proceso,$A | t=0$, es dado por$2N*P0$.

Supongamos que el tamaño de la población es constante ($Nt+1 = Nt$para todos$t$), y que el apareamiento es aleatorio. Supongamos además que las generaciones son distintas. Dejar$t = 1$ser la primera generación nacida, y así sucesivamente.

El número esperado de copias de$A$en el momento$t+1$luego sigue una distribución binomial:

$At+1$esta distribuido$Binom(2N, At/(2N))$

Como tenemos una regla de transición para$At → At+1$, es una cuestión bastante simple simular poblaciones que experimentan esta forma de evolución. He escrito la función R neutralFPT para realizar esta simulación y estimar la proporción de todas las poblaciones que alcanzan$P1$en algún momento de su historia, el tiempo esperado para que un alelo neutral alcance$P1$de$P0$, dado que alcanzará$P1$, y la distribución de los períodos de tiempo necesarios para alcanzar$P1$, dado que una población alcanzará$P1$. El guión se da en la sección final de esta respuesta.

Densidades de probabilidad de los tiempos del primer paso:

Densidades de probabilidad de tiempos de primer paso sobre valores plausibles de$P0$,$P1$, y$N$siguen una estructura similar: unimodal cerca del margen izquierdo, con una cola derecha larga de tiempos de primer paso más largos (Figura 2).

Figura 2: densidades de probabilidad de los tiempos del primer paso de P0 a P1 bajo diferentes valores de P0, P1 y 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.

Función: neutralFPT, con ejemplos de uso.

 ## 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.