10
$\begingroup$

Kimura and Ohta (1968) showed that the expected time for a neutral allele to reach fixation (given that it will reach fixation) is

$$\bar t(p_0)=-4N\left(\frac{1-p_0}{p_0}\right)\ln(1-p_0),$$

where $p_0$ is the initial frequency and $N$ is the population size.

From their work, can we generalize this result to calculate the expected time to reach a frequency of $p_1$ (given that this frequency will be reached at some point), where $p_1$ is not necessarily equal to $1$?

$\endgroup$
  • 3
    $\begingroup$ In the above paper, a key quantity is u(p,t), the prob. that an allele fixes at time t given starting freq. p. They cite a previous paper showing that u(p,t) satisfies a certain PDE. In your case, you want to modify u(p,t) to be a sort of first-passage-time probability (i.e. the prob. that the allele freq. reaches p1 for the 1st time at time t given initial freq. p). In this case, u(p,t) may no longer satisfy the same PDE as before (check that other paper). In order to proceed, you need to find the PDE that your new u(p,t) satisfies. Look up first passage time problems, it may help. $\endgroup$ – A. Kennard Mar 29 '15 at 10:50
2
$\begingroup$

Plain English Answer:

I have written a computer function which simulates neutral evolution to solve this problem. It's not an exact mathematical answer, but it is basically the same approach as Kimura and Ohta took in (the second half of) their paper, except that my computer is more powerful than their one, so I could get far more precise estimates by simulating more populations than they did.

Figure One is a lattice plot of the relationship between P0 and the expected time to reach P1, for different values of P1 (by column) and different population sizes (by row). It is clear that the same relationships between P0, P1 and expected time to get from P0 to P1 are seen in every population size, but with larger populations, you should expect a longer time to get from P0 to P1.

When P0 and P1 are close together, you generally expect a smaller time to get from P0 to P1 than when P0 and P1 are far apart. You might think of this as a rule like 'it takes a longer time to travel further'. Expected time taken to get from P0 to P1 is smaller when P0 and P1 are on the same side of 0.5 as each other and P1 is further from 0.5 than P0 is, than when P0 is further from 0.5 than P1 is. You might think of this as a rule like 'it is quicker to travel when you're closer to fixation or extinction, than when you're in the middle frequencies'.

First-passage time lattice plot

Figure One: Expected time for a neutral allele to go from frequency P0 to frequency P1, given differences in P0, sub-plotted by P1 and population size. Modelled frequencies of P0 and P1 are 0.1, 0.3, 0.5, 0.7, and 0.9. Modelled population sizes are N = 10, N = 50, and N = 100. Each sub-plot represents one combination of P1 and population size, with P1 increasing from left to right, and population size increasing from bottom to top. All expectations are estimated from 10,000 simulated populations.

Some populations take longer than others to go from P0 to P1. Generally, most populations take a small number of generations to reach P1, and a smaller number take a long time – but a very few can take a very very long time.

I've attached my code, so if you know how to use R you can use it to estimate the expected time to get from P0 to P1 for any combination of P0, P1, and population size.

Additional details and assumptions for the technically inclined:

Relevant theory:

In a diploid, sexually-reproducing population of size $N$, there are $2N$ copies of a given locus. Each locus is occupied by some allele. All variation at our locus of interest is selectively neutral. For our purposes, we can consider that each locus is occupied either by our allele of interest ($A$), or another variant ($a$), ignoring variation in the 'other' alleles. The number of copies of $A$ in the population at the start of the process, $A | t=0$, is given by $2N*P0$.

Let us assume that population size is constant ($Nt+1 = Nt$ for all $t$), and that mating is random. Let us further assume that generations are distinct. Let $t = 1$ be the first generation born, and so on.

The expected number of copies of $A$ at time $t+1$ then follows a binomial distribution:

$At+1$ is distributed $Binom(2N, At/(2N))$

Since we have a transition rule for $At → At+1$, it is a fairly simple matter to simulate populations undergoing this form of evolution. I have written the R function neutralFPT to carry out this simulation and estimate the proportion of all populations which reach $P1$ at some point in their history, the expected time for a neutral allele to reach $P1$ from $P0$, given that it will reach $P1$, and the distribution of lengths of time taken to reach $P1$, given that a population will reach $P1$. The script is given in the final section of this answer.

Probability densities of first-passage times:

Probability densities of first-passage times over plausible values of $P0$, $P1$, and $N$ follow a similar structure – unimodal near the left margin, with a long right-tail of longer first-passage times (Figure 2).

First-passage time probability densities

Figure 2: probabiltiy-densities of first-passage times from P0 to P1 under different values of P0, P1, and 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.

Function: neutralFPT, with examples of use.

 ## 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.
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.