branching process - senecio7の日記

Uecker H & Hermisson J 2011 On the fixation process of a beneficial mutation in a variable environment. Genetics 188:915-930.

N_e,t = N_t / (2ξ(t, N_t) + b(t, N_t) + d(t, N_t) + s(t, N_t)) ... (5)

p₀(n₀, t) = P_n₀(0, t) = (1 − 1 / (A + B)(t))^n₀ = [∫₀^tμ(τ)e^−ρ(τ)dτ / (1 + ∫₀^tμ(τ)e^−ρ(τ)dτ)]^n₀ ... (13)

A + B = lim_t→∞(A + B)(t) = 1 + ∫₀^∞μ(t)e^−ρ(t))dt = ∫₀^∞λ(t)e^−ρ(t))dt ... (15)
the last equality valid for lim_t→∞exp(−ρ(t)) = 0
this condition is met in all examples considered below

for a single mutant (n₀ = 1)
p_fix = 1 / (A + B) = 2 / [1 + ∫₀^∞ (λ + μ)(t) exp(− ∫₀^t(s + r)(τ)dτ)dt] ... (16a)
= 2 / [1 + ∫₀^∞ (N(0)/N_e(t)) exp(− ∫₀^ts(τ)dτ)dt] ... (16b)
∫₀^tr(τ)dτ = ∫₀^t(N(⋅)(τ)/N(τ))dτ = ln(N(t)/N(0))
a similar expression was also derived by Ohta and Kojima (1968)
restricted to a constant population size and a Poisson offspring distribution, their result is, however, less general

time to reach an intermediate frequency x_c
a mutation that has survived the initial phase will on average have grown faster during this phase than predicted by the deterministic model
the frequency will thus be larger than if it had always grown following the deterministic path
it is well described by the deterministic path that is not started from a single individual, but from an (on average) larger "effective initial population size"
it consists of subsuming all the stochasticity of the path under this effective initial population size as a single random variable and then modeling the path deterministically
the procedure consists of two steps
in a first step the distribution of the effective initial population size is estimated via a branching process
in a second step the deterministic approximation of the full birth-death model is used to describe the allele frequency path starting from the effective initial population size

ν_t := n_t / E[n_t] = n_t exp(−ρ(t)) ... (33)
ν_t describes all the stochasticity that has accumulated in the branching process until time t
for t → ∞, ν_t converges (almost surely) to a positive random variable ν that summarizes the entire stochasticity of the process
n_t = ν exp(ρ(t)) in the limit t → ∞
we can interprete ν as the random initial population size of an ensemble of deterministically growing paths that approximate the original process n_t for large t

lim_t→∞exp(ρ(t))ln(B(t)/(A + B)(t)) = −1/(A + B)
we obtain in the limit t → ∞ the stationary distribution
P(ν ≤ ν₀ | not extinct) = 1 − exp(−p_fixν₀) ... (40)
this effective initial mutant population size may now be used as the starting value for the deterministic solution of the full model

we can express the random variable T_{x_c} defined via x(T_{x_c}) = x_c in terms of ν
T_{x_c} = β⁻¹(ln(x_c(N₀ − ν) / (1 − x_c)ν)), ν < x_cN₀ ... (43)
[error fixed]
P(T_{x_c} ≤ T) = exp[− p_fixN₀ / (((1 − x_c) / x_c)exp(β(T))] ... (44)
Equation 44 constitutes the main result of this section

for the fixation probability, we have seen in Equation 16a that there is a second way to summarize the dynamics in terms of the absolute rate of increase (s + r)(t) and the total rate of birth and death events
this is not possible for the distribution of T_{x_c}, linked to the fact that variable selection and demographic change are not equivalent in this case
while the selection function s(t) has a strong influence on the deterministic part of the frequency path x(t), changes in the total population size (with equal effect on mutants and residents) have no influence on the latter

constant selection and constant population size
P(T_{x_c} ≤ T) = exp[−(s / (ξ + s))(N / (((1 − x_c) / x_c)exp(sT) + 1)] ≈ exp[−2sN / (((1 − x_c) / x_c)exp(sT) + 1)] ... (49)
T_{x_c} can be interpreted as the age of a derived allele that is currently found at frequency x_c in the population
this is a consequence of the time-homogeneity of the model
the distribution of the time to reach fixation at x = 1 starting from a frequency (1 − x_c) equals the distribution of T_{x_c}
the times needed from 0 to 0.5 and from 0.5 to 1 follow the same distribution
both random variables are independent
p^~(t_fix) = ∫₀^t_fixp(T_1/2)p(t_fix − T_1/2)dT_1/2 ... (50)
an alternative expression for the distribution of t_fix has been derived before by Wang and Rannala (2004), using diffusion techniques
our approximation (Equation 50) is simpler than the series expansion in terms of Gegenbauer polynomials, using the eigenvalues of the coefficients of the spheroidal harmonics obtained by these authors
it is nevertheless highly accurate if selection is not very weak or the population size small

simulations
we performed individual based computer simulations for which we use a Gillespie algorithm
events happen at rate (2 ξ(t, N_t) + s(t, N_t)) n_t (N_t − n_t) / N_t + b(t, N_t) N_t + d(t, N_t) N_t
ξ(t, N_t) s(t, N_t), b(t, N_t) and d(t, N_t) are assumed to be constant between events
once the time of an event is fixed, which kind of event takes place is determined using the respective probabilities
subsequently, the values of N_t and n_t are updated and the rates are set to their new values
additional update steps between events did not change the results

the impact of the ecological dynamics on the genetics of adaptation has been studied for the so-called moving optimum model

even if the external environment is constant, individual alleles in a population can experience variable selection pressures if multiple selected alleles segregate in the population and if these alleles interfere due to either epistasis or linkage

the limitations of the branching approach are the usual ones
the fate of the mutation must be decided while the mutant frequency is small and the independence assumption of the branching process applies
deviations from the simulation results are found in the case of mutations that are almost neutral on average, with fixation probabilities p_fix ≲ 10 / N