Ожидаемое время для нейтрального аллеля, чтобы достичь частоты p1p1p_1 при запуске с частоты p0p0p_0

Простой английский ответ:

Я написал компьютерную функцию, которая имитирует нейтральную эволюцию, чтобы решить эту проблему. Это не точный математический ответ, но в основном это тот же подход, который Кимура и Охта использовали во (второй половине) своей статьи, за исключением того, что мой компьютер мощнее, чем их, поэтому я мог получить гораздо более точные оценки, моделируя больше населения, чем они.

Рисунок 1 представляет собой решетчатый график взаимосвязи между P0 и ожидаемым временем достижения P1 для разных значений P1 (по столбцам) и разных размеров популяции (по строкам). Понятно, что одни и те же отношения между P0, P1 и ожидаемым временем перехода от P0 к P1 наблюдаются для каждого размера популяции, но с большими популяциями вы должны ожидать более длительное время перехода от P0 к P1.

Когда точки P0 и P1 расположены близко друг к другу, вы, как правило, ожидаете, что путь от P0 до точки P1 займет меньше времени, чем когда точки P0 и P1 находятся далеко друг от друга. Вы можете думать об этом как о правиле, похожем на «дальнее путешествие занимает больше времени». Ожидаемое время, необходимое для перехода от P0 к P1, меньше, когда P0 и P1 находятся на одной стороне 0,5 друг от друга, а P1 дальше от 0,5, чем P0, чем когда P0 дальше от 0,5, чем P1. Вы можете думать об этом как о правиле вроде «Путешествовать быстрее, когда вы ближе к фиксации или исчезновению, чем когда вы находитесь на средних частотах».

Рисунок 1. Ожидаемое время перехода нейтрального аллеля от частоты P0 к частоте P1 с учетом различий в P0, нанесенных на график по P1 и размеру популяции. Смоделированные частоты P0 и P1 составляют 0,1, 0,3, 0,5, 0,7 и 0,9. Смоделированные размеры популяции: N = 10, N = 50 и N = 100. Каждый подучасток представляет собой одну комбинацию P1 и размера популяции, при этом P1 увеличивается слева направо, а размер популяции увеличивается снизу вверх. Все ожидания оцениваются по 10 000 смоделированных популяций.

Некоторым популяциям требуется больше времени, чем другим, чтобы перейти от P0 к P1. Как правило, большинству популяций требуется небольшое количество поколений, чтобы достичь P1, и меньшему количеству требуется много времени, но очень немногим может потребоваться очень и очень много времени.

Я приложил свой код, поэтому, если вы знаете, как использовать R, вы можете использовать его для оценки ожидаемого времени перехода от P0 к P1 для любой комбинации P0, P1 и численности населения.

Дополнительные детали и предположения для технически подкованных:

Соответствующая теория:

В диплоидной, размножающейся половым путем популяции размером$N$, Существуют$2N$копии данного локуса. Каждый локус занят некоторым аллелем. Все вариации в нашем локусе интереса избирательно нейтральны. Для наших целей мы можем считать, что каждый локус занят либо интересующим нас аллелем ($A$), или другой вариант ($a$), игнорируя вариации «других» аллелей. Количество копий$A$в популяции в начале процесса,$A | t=0$, дан кем-то$2N*P0$.

Предположим, что численность популяции постоянна ($Nt+1 = Nt$для всех$t$), и это спаривание случайно. Предположим далее, что поколения различны. Позволять$t = 1$родиться в первом поколении и так далее.

Ожидаемое количество копий$A$вовремя$t+1$затем следует биномиальное распределение:

$At+1$распространяется$Binom(2N, At/(2N))$

Поскольку у нас есть правило перехода для$At → At+1$, довольно просто смоделировать популяции, подвергающиеся этой форме эволюции. Я написал функцию нейтрального FPT R, чтобы выполнить это моделирование и оценить долю всех популяций, которые достигают$P1$в какой-то момент их истории ожидаемое время для достижения нейтрального аллеля$P1$из$P0$, учитывая, что он достигнет$P1$, а также распределение времени, необходимого для достижения$P1$, учитывая, что население достигнет$P1$. Сценарий приведен в последнем разделе этого ответа.

Плотности вероятностей времени первого прохождения:

Плотности вероятностей времени первого прохождения над правдоподобными значениями$P0$,$P1$, а также$N$следуют аналогичной структуре - унимодальной около левого края, с длинным правым хвостом более длительного времени первого прохода (рис. 2).

Рис. 2. Плотности вероятностей времени первого прохождения от P0 до P1 при различных значениях P0, P1 и N. A: P0 = 0,5, P1 = 0,9, N = 50; Б: Р0 = 0,9, Р1 = 0,5, N = 50; В: Р0 = 0,5, Р1 = 0,9, N = 500.

Функция: нейтральный FPT, с примерами использования.

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