# Asymptotic Green’s function for the stochastic reproduction of competing
variants via Fisher’s angular transformation

Bhavin S. Khatri  The Francis Crick Institute, Mill Hill Laboratory,
The Ridgeway, London, NW7 1AA, U.K.

###### Abstract

The Wright-Fisher Fokker-Planck equation describes the stochastic dynamics of
self-reproducing, competing variants at fixed population size. We use Fisher’s
angular transformation, which defines a natural length for this stochastic
process, to remove the co-ordinate dependence of it’s diffusive dynamics,
resulting in simple Brownian motion in an unstable potential, driving variants
to extinction or fixation. This insight allows calculation of very accurate
asymptotic formula for the Green’s function under neutrality and selection,
using a novel heuristic Gaussian approximation.

Understanding the interplay between stochastic and deterministic forces in
systems with different reproducing variants is a theme that arises, and has
importance, in many different scientific fields Blythe2007  including language
evolution Baxter2006 , protein evolution Kimura1984 ; Akashi2012 , the
evolution of biodiversity HubbellBook2001 ; McKane2000 ; Volkov2003  and
population genetics Crow1970 ; Wright1945 . This article is concerned with a
fundamental question in population genetics, given the possibility of only two
reproducing variants, how does the probability distribution of gene frequency
$x(t)$ change over time, given it is known at a prior time point $x_{0}=x(0)$,
subject to small number fluctuations (genetic drift) and selection
(competition), though the analogous question may be posed in any of these
fields. We address this question in the context of the Wright-Fisher model,
which is the canonical model of stochastic dynamics incorporating both these
features.

The diffusion approximation Kimura1964 , of the Wright-Fisher model describes
the stochastic dynamics of gene frequency $x$ ($=n/N$, where $n$ is the number
of copies of the mutant allele and $N$ the total population):

$\displaystyle\frac{\partial G(x,x_{0};t)}{\partial
t}=-\frac{\partial}{\partial x}$
$\displaystyle\left(sx(1-x)G(x,x_{0};t)\right)$ (1)
$\displaystyle+\frac{1}{2N}\frac{\partial^{2}}{\partial
x^{2}}\left(x(1-x)G(x,x_{0};t)\right),$

where $G(x,x_{0};t)$ is the transition probability density, or Green’s
function, of gene frequency given an initial condition
$G(x,x_{0};0)=\delta(x-x_{0})$, and $s$ is the selection coefficient. This
equation is derived in the large $N$ limit from a Master equation of discrete
populations of each variant at a fixed $N$ Kimura1964 .

This Fokker-Planck equation has been studied extensively. In particular,
Kimura calculated a series solution for the neutral equation in terms of
Gegenbauer polynomials Kimura1954 ; Kimura1964 , which was later extended to
the multi-allele case by Baxter, et al Baxter2007 . For the case of selection
Kimura also calculate a series solution, however, the eigenvalues could not be
represented in closed-form in terms of the population size $N$ and selection
coefficient $s$ Kimura1964 . More recently, a number of methods have been
developed to calculate the Green’s function under selection, including a
numerical matrix approach Song2012  and perturbation theory based on a path-
integral formulation of the Wright-Fisher process Schraiber2014 . However, for
many practical applications, such as virus evolution, where population sizes
are large and generation times short, these approaches are not very practical
as a large number of terms is required for convergence at short times. The
solution of Voronka and Keller Voronka1975 , which uses an asymptotic ray
approximation of the Green’s function, is valid at short times and for models
of selection, neutrality and mutation, but their approach is not very
intuitive and unwieldy requiring switching between different solutions in a
time-dependent manner. We present a simple short-time asymptotic calculation
of the Green’s function in closed form for both neutrality and selection,
which has intuitive appeal as it exploits Fisher’s angular transformation,
which as we show is the natural co-ordinate for Wright-Fisher stochastic
dynamics Antonelli1977 .

Figure 1: Effective deterministic force in angular domain for Wright-Fisher
process, where $2NF(\theta)=-\cot\theta+Ns\sin\theta$.

Fokker-Planck equations with co-ordinate dependent diffusion constants such as
Eqn.1 have the property that space is explored at different rates dependent on
the position in the domain; using this intuition Antonelli et al Antonelli1977
, suggested the natural definition of length for a stochastic process be
related to the differential
$\mathrm{d}\theta^{2}=\sum_{ij}g_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j}$, where
$g_{ij}$ is a metric tensor and taken to be the inverse of the covariance
matrix $g^{ij}$. In one-dimension, this is simply
$\mathrm{d}\theta^{2}=\mathrm{d}x^{2}/D(x)$, which represents the
(differential) mean square distance traversed in equal times. For the
diffusion constant of random drift, $D(x)=x(1-x)$, the stochastic distance is
simply

$\theta=\int^{x}\frac{\mathrm{d}x^{\prime}}{\sqrt{x^{\prime}(1-x^{\prime})}}=\cos^{-1}(1-2x).$
(2)

This is Fisher’s angular transformation Fisher1922 ; Fisher1930a . Applying
this transformation to the neutral Fokker-Planck equation ($s=0$), we get

$\frac{\partial p(\theta,t)}{\partial
t}=\frac{1}{2N}\left(\frac{\partial^{2}p(\theta,t)}{\partial\theta^{2}}+\frac{\partial}{\partial\theta}\left(\cot\theta
p(\theta,t)\right)\right),$ (3)

or the equivalent stochastic differential equation GardinerBook ; vanKampen

$\frac{\mathrm{d}\theta}{\mathrm{d}t}=-\frac{1}{2N}\cot(\theta)+\eta(t),$ (4)

where $\langle\eta(t)\rangle=0$ and
$\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t^{\prime})/N$. We see the
result of this transformation is a co-ordinate independent diffusion constant
in $\theta-$space, but now with an effective deterministic force
$f(\theta)=\cot(\theta)/2N$. This arises due to the co-ordinate dependent
diffusion constant $x(1-x)/N$ in the $x-$domain and drives the system towards
regions of decreasing diffusion constant. It is also exactly the spurious
drift term that arises in transforming between Ito and Stratonovich
descriptions of stochastic dynamics GardinerBook ; vanKampen . Note that in
the Langevin representation, Eqn.1 would have a multiplicative noise term,
which is transformed to additive noise in $\theta-$space in Eqn.4. Examining
the force in Fig.1, we see that it is unstable, on average driving a mutant
alelle to extinction if $\theta(0)<\pi/2$ and fixation if $\theta(0)>\pi/2$,
with a fixed point at $\theta=\pi/2$. To calculate an approximate solution of
the Green’s function, we make a Taylor expansion of the force about the fixed
point $\theta=\pi/2$ to linear order to give a linear stochastic differential
equation,
$\dot{\theta}=\frac{1}{2N}\left(\theta-\frac{\pi}{2}\right)+\eta(t)$. As we
will see this approximation works well even for initial frequencies of order
$10\%$ ($x(0)=0.1$), due to the non-linearity of the angular transformation,
which has the property of compressing the central range in $x-$space about
$x=1/2$ to a smaller central region in $\theta-$ space about $\theta=\pi/2$
(for example, $x=0.1\rightarrow\theta=0.64$ and $x=0.9\rightarrow\theta=2.5$).
Equivalently, this is an harmonic approximation of the effective potential
function in $\theta-$space, where
$\dot{\theta}=\partial_{\theta}U(\theta)+\eta(t)$ and
$U(\theta)=\frac{1}{2N}\ln\sin\theta$. As the resulting SDE is linear the
solution is straightforwardly computed as
$\theta(t)=\frac{\pi}{2}+\left(\theta_{0}-\frac{\pi}{2}\right)e^{t/2N}+\int_{0}^{t}\mathrm{d}t^{\prime}\eta(t^{\prime})e^{\frac{t-t^{\prime}}{2N}},$
where $\theta_{0}=\theta(0)$. As the solution is a sum of Gaussian random
variables $\eta$ the Green’s function for $\theta$ will be Gaussian with mean,
$\langle\theta(t)\rangle=\frac{\pi}{2}+\left(\theta_{0}-\frac{\pi}{2}\right)e^{t/2N}$,
since $\langle\eta(t)\rangle=0$ and variance,
$\langle\langle\theta^{2}(t)\rangle\rangle=e^{t/N}-1$, where the van Kampen
notation has been used,
$\langle\langle\theta^{2}\rangle\rangle=\langle\theta^{2}\rangle-\langle\theta\rangle^{2}$.
Note that the variance diverges for $t\gg N$ and the mean divergences for
$t\gg 2N$, to $-\infty$ when $\theta_{0}<\pi/2$ and to $+\infty$ for
$\theta_{0}>\pi/2$ and is fixed for all time at $\langle\theta\rangle=\pi/2$,
if $\theta_{0}=\pi/2$ the fixed point of the deterministic dynamics. The
Green’s function in $\theta-$space is then

$\displaystyle G_{\theta}(\theta,\theta_{0};t)$
$\displaystyle=\frac{1}{\sqrt{2\pi(e^{t/N}-1)}}$
$\displaystyle\times\exp{-\frac{\left(\theta-\theta_{0}e^{t/2N}-\frac{\pi}{2}(1-e^{t/2N})\right)^{2}}{2(e^{t/N}-1)}}.$
(5)

Note the similarity of form to the Green’s function of an overdamped harmonic
oscillator, but with the difference that, as discussed, here the mean and
variance diverge DoiEdwards . This solution does not obey the boundary
conditions at $\theta=0$ and $\theta=\pi$, which are required to be absorbing
and specifically to go linearly to zero at these points; this is in order for
the solution in $x-$space to be finite at the boundaries, as required due to
the singularity of the diffusion constant at $x=0$ and $x=1$ Baxter2007 . The
method of images cannot be used in this case as the required image has it’s
forces reversed and so does not obey the original Fokker-Planck equation.
However, as we argue in the discussion, for many applications, including virus
evolution, the short time behaviour ($t\ll N$) is most relevant. Transforming
back to $x-$space, we have,

$\displaystyle G_{x}(x,x_{0};t)$
$\displaystyle=\left|\frac{\mathrm{d}\theta}{\mathrm{d}x}\right|G_{\theta}(\theta(x),\theta_{0}(x_{0});t)$
$\displaystyle=\frac{1}{\sqrt{2\pi
x(1-x)(e^{t/N}-1)}}\exp{-\frac{\left(\cos^{-1}(1-2x)-\cos^{-1}(1-2x_{0})e^{t/2N}-\frac{\pi}{2}(1-e^{t/2N})\right)^{2}}{2\left(e^{t/N}-1\right)}},$
(6)

where the Jacobian is
$\left|\frac{\mathrm{d}\theta}{\mathrm{d}x}\right|=2/\sin\theta=1/\sqrt{x(1-x)}$.

Figure 2: Comparison of approximate calculation of neutral Greens function
(solid lines – Eqn.Asymptotic Green’s function for the stochastic reproduction
of competing variants via Fisher’s angular transformation) and numerical
integration of stochastic differential equation that arises from diffusion
approximation (solid circles). a) initial frequency $x_{0}=0.01$, b)
$x_{0}=0.1$, c) $x_{0}=0.5$.

The results are plotted in Fig.2, at various times and initial conditions, as
solid lines and compared against numerical integration of the exact neutral
Wright-Fisher stochastic differential equation (Eqn.4). We see that, in
general, the approximation works very well for short times $t\ll N$ and when
the initial frequency $x_{0}$ is not too small. More precisely we would expect
the approximation to be good for sufficiently short times compared to the
expected time to fixation, which is $\langle
t^{*}\rangle=-2N(x_{0}\ln(x_{0})+(1-x_{0})\ln(1-x_{0}))$; for
$x_{0}=\\{0.01,0.1,0.5\\}$, $\langle t^{*}\rangle\approx\\{0.1N,0.7N,1.4N\\}$,
which is consistent with the results in Fig.2. The solution is more simple and
intuitive than that of Voronka et al, Voronka1975 , where it is clear the
behaviour of gene frequencies is essentially that of Brownian motion in an
unstable harmonic potential; the non-linearity in $x-$space arises purely from
working in the more natural co-ordinates of the angular transformation, where
in particular the argument of the exponential in the Gaussian solution is just
the square of the stochastic distance between $x_{0}$ and $x$. For example, it
is instructive that, as a consequence, even for a neutral process the mean of
allele frequencies moves towards the extinction or fixation boundary for any
initial frequency $x_{0}\neq 0.5$, as is seen clearly from the plots of the
Green’s functions in Fig.2 and from the solution of the mean
$\langle\theta(t)\rangle$; this is not obvious from the neutral Wright-Fisher
diffusion equation in $x-$space (Eqn.1 with $s=0$).

For the case of selection, from Eqn.1 and using Fisher’s angular
transformation, the stochastic differential equation for $\theta$ is

$\frac{\mathrm{d}\theta}{\mathrm{d}t}=-\frac{1}{2N}\left(\cot(\theta)-Ns\sin(\theta)\right)+\eta(t),$
(7)

where $\eta(t)$ has the same moments as before. Note that the contribution of
selection to the effective force tends to zero as
$\theta\rightarrow\\{0,\pi\\}$, which agrees with the intuition that when an
allele is rare, the change in allele frequency is dominated by drift; in
particular, for $\theta\ll 1$, and $Ns\gg 1$,
$2Nf(\theta)\approx-1/\theta+Ns\theta$ and the forces of drift and selection
are roughly in balance when $Ns\sim 1/\theta^{2}=1/4x$, where Fisher’s angular
transformation is $\theta\approx\sqrt{4x}$ for small $x$ – in other words when
the allele frequency $x\ll(4Ns)^{-1}$ drift dominates.

When selection is weak ($Ns\ll 1$), the effective force in the angular domain
is only a weak perturbation on the neutral force (Fig.1) and the Green’s
functions differ little from neutrality, particularly at short times (not
shown). A similar linear expansion of the force can be carried out to
calculate an asymptotic expression for the Green’s function under weak
selection, as shown in the Supplementary Online Material; the resulting
expression has similar accuracy compared to numerical simulations as in the
neutral case.

In the regime of strong selection $Ns\gtrsim 1$, the above approach gives a
poor approximation, due to the non-linearity of the effective force in the
angular domain (Fig.1). Here we present a heuristic approach to solving Eq.7
approximately, for any value of $Ns$. The approach is to assume that the
Green’s function of the non-linear SDE can be approximated by a Gaussian
process, where: 1) the time-varying mean is calculated as a solution to the
deterministic dynamics of the SDE Eq.7, with initial condition $\theta_{0}$,
which we show below can be calculated exactly; and 2) the time-varying
variance $\langle\langle\theta^{2}(t)\rangle\rangle$ is dependent on the local
gradient of the force, which varies as a function of the deterministic
solution, $\lambda=f^{\prime}(\langle\theta(t)\rangle)$. In general, if an
exact solution is not available to the deterministic dynamics, an approximate
solution that makes a linear approximation of the effective force (Fig.1)
about the initial condition, also gives accurate results at short times (not
shown).

Transforming the deterministic part of Eqn.7 back to $x-$space, we have
$\dot{x}=-\frac{1}{2N}(\frac{1}{2}(1-2x)-2Nsx(1-x)$, the solution to which is
of the form $x=C+A\tanh{(\gamma t/2+B)}$. Transforming back to $\theta-$space
and using the initial condition $\theta_{0}=\langle\theta(0)\rangle$, the
solution to the deterministic dynamics of Eqn.7 is:

$\langle\theta(t)\rangle=\cos^{-1}\left(-\frac{1+2N\gamma\tanh\left(\gamma
t/2-\alpha\right)}{2Ns}\right)$ (8)

where $\alpha=\tanh^{-1}\left(\frac{2Ns\cos{\theta_{0}}+1}{2N\gamma}\right)$
and the characteristic rate of change of the mean is
$\gamma=\sqrt{1+4N^{2}s^{2}}/2N$.

The next step is to calculate the variance, which we motivate by considering
the situation when the slope of the effective force is fixed to a constant
$\lambda$, which gives a Gaussian solution with variance
$\langle\langle\theta^{2}(t)\rangle\rangle=\frac{1}{2N\lambda}(e^{2\lambda
t}-1)$. The linearity of the force characterises the Gaussian distribution and
so if we assume that the effective deterministic force varies slowly over a
range of theta representing the width of the probability density, we can then
heuristically replace $\lambda$ with the local derivative of the effective
force $\lambda(\langle\theta(t)\rangle)$ in the variance. This approximates
the local spreading of the probability density being solely due to the local
derivative of the force giving a time varying variance:

$\langle\langle\theta^{2}(t)\rangle\rangle=\frac{1}{2N\lambda(\langle\theta(t)\rangle)}(e^{2\lambda(\langle\theta(t)\rangle)t}-1).$
(9)

Note that for strong selection, the derivative of the effective deterministic
force $\lambda(\langle\theta(t)\rangle)$ will be zero at certain times, as can
be seen from the plot of the deterministic force in Fig.1; at these time
points it is simple to see that the variance remains well behaved as
$\lim_{\lambda\rightarrow
0}\\{\langle\langle\theta^{2}\rangle\rangle\\}\rightarrow t/N$, as one would
expect if the deterministic force tends to a constant. Transforming back to
$x-$space, the Green’s function solution is:

$\displaystyle
G_{x}(x,x_{0};t)=\frac{\exp{-\frac{\left(\cos^{-1}(1-2x)-\langle\theta\rangle\right)^{2}}{2\langle\langle\theta^{2}\rangle\rangle}}}{\sqrt{2\pi
x(1-x)\langle\langle\theta^{2}\rangle\rangle}},$ (10)

where $\langle\theta\rangle$ and $\langle\langle\theta^{2}\rangle\rangle$ are
given by Eqns.8 and 9, respectively, where $\cos(\theta_{0})=1-2x_{0}$.

We plot the results for $Ns=10$ ($x_{0}=\\{0.1,0.5,0.9\\}$ in Fig.3; see
Supplementary Online Information for plots with $x_{0}=\\{0.01,0.99\\}$)
(Green’s functions for $Ns=1$ are plotted in the Supplementary Online
Information). We find that for both $Ns=1$ and $Ns=10$ the heuristic approach
and the integration of the Wright-Fisher SDE (Eqn.7) agree very well at short
times compared to the average expected time for fixation/extinction of a
mutant. This is true even when $x_{0}$ is very close to $0$ or $1$
(Supplementary Online Information) and is reasonably accurate to quite long
times ($t\sim\tau$) for an initial frequency of $x_{0}=0.1$ (Fig.3A). In
addition, we see that as well as capturing the broad behaviour of the time-
varying mean and variance, the insets of the figures show the Green’s
functions plotted on a log scale, demonstrating that the approximation is also
very accurate in the tails of the distribution at short times. Finally, for
very long times when $\langle\theta(t)\rangle$ tends to $0$ or $\pi$, the
variance of the heuristic solution Eqn.9 diverges, as $\lambda$ diverges at
the boundaries, and the approximation fails; this is indicated in those cases
where there is no heuristic solution plotted for a given time in each plot.

Figure 3: Comparison of approximate calculation of Greens function for $Ns=10$
(solid lines – Eqn.10) and numerical integration of stochastic differential
equation that arises from diffusion approximation (solid circles). a) initial
frequency $x_{0}=0.1$, b) $x_{0}=0.5$, c) $x_{0}=0.9$. Green’s functions are
plotted at times given by fractions of $\tau=\frac{1}{s}(1+\ln{(Ns)})$, which
is approximately the expected time to fixation of a mutant which survives
drift and then is driven to fixation by selection Desai2007 .

To conclude, we have calculated very accurate approximations of the 2-allele
Green’s function (or transition probability density) of population genetics
for arbitrary selection coefficient $s$ and population size $N$. A key
advantage and insight of the approach outlined in this paper, is that it
transforms a non-linear Fokker-Planck equation to a simple problem of Brownian
motion in an effective potential. Together with the heuristic Gaussian
approximation, this represents, to the author’s knowledge, a new general
approach for asymptotically solving Fokker-Planck equation’s with a co-
ordinate dependent diffusion constant in slowly-varying potentials (or
equivalently SDEs with multiplicative noise), where the solution to the mean
behaviour is known; indeed, in 1-dimension a PDE with co-ordinate dependent
diffusion can always be transformed to one with co-ordinate independent
diffusion Antonelli1977 ; Baxter2007 . For more than two variants the methods
detailed in Antonelli1977 , suggests via the metric tensor, a potential route
to finding solutions in higher dimensions.

These results have potential application to detecting selection in time-series
data of the composition of variants, in biological evolution, language
evolution and for species in ecosystems. In particular, as these results have
accuracy in the asymptotic short-time limit, they will be applicable to
studying selection from time-series of variants (haplotypes) in virus
evolution, since they have large effective population sizes and short
generation times, meaning even sampling virus populations infrequently (on the
time scale of many months or years) would be accurately modelled by the
results of this paper.

###### Acknowledgements.

I thank Richard A. Goldstein for initially suggesting the problem and for
useful discussions. I also thank Richard Blythe for useful comments on the
manuscript. This work was supported by The Francis Crick Institute which
receives its core funding from Cancer Research UK, the UK Medical Research
Council and the Wellcome Trust.

## References

  * (1) R. A. Blythe, A. J. McKane, Journal of Statistical Mechanics: Theory and Experiment 2007, P07018 (2007).
  * (2) G. J. Baxter, R. A. Blythe, W. Croft, A. J. McKane, Physical Review E 73, 046118 (2006).
  * (3) M. Kimura, The neutral theory of molecular evolution (Cambridge University Press, 1984).
  * (4) H. Akashi, N. Osada, T. Ohta, Genetics 192, 15 (2012).
  * (5) S. P. Hubbell, The Unified Neutral Theory of Biodiversity and Biogeography (Princeton University Press, 2001).
  * (6) A. McKane, D. Alonso, R. V. Solé, Physical Review E 62 (2000).
  * (7) I. Volkov, J. R. Banavar, S. P. Hubbell, A. Maritan, Nature 424, 1035 (2003).
  * (8) J. F. Crow, M. Kimura, et al., An introduction to population genetics theory. (1970).
  * (9) S. Wright, Proc Natl Acad Sci U S A 31, 382 (1945).
  * (10) M. Kimura, Journal of Applied Probability 1, 177 (1964).
  * (11) M. Kimura, Genetics 39, 280 (1954).
  * (12) G. J. Baxter, R. A. Blythe, A. J. McKane, Math Biosci 209, 124 (2007).
  * (13) Y. S. Song, M. Steinrücken, Genetics 190, 1117 (2012).
  * (14) J. G. Schraiber, Theoretical population biology 92, 30 (2014).
  * (15) R. Voronka, J. B. Keller, Mathematical Biosciences 25, 331 (1975).
  * (16) P. L. Antonelli, C. Strobeck, Advances in Applied Probability pp. 238–249 (1977).
  * (17) R. A. Fisher, Proceedings of the Royal Society of Edinburgh 42, 321 (1922).
  * (18) R. A. Fisher, Proceedings of the Royal Society of Edinburgh 50, 205 (1930).
  * (19) C. Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences (Springer, 2009).
  * (20) N. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, 1981).
  * (21) M. Doi, S. Edwards, The Theory of Polymer Dynamics (Oxford University Press, 1986).
  * (22) M. M. Desai, D. S. Fisher, Genetics 176, 1759 (2007).

Supplemental Materials

## Green’s function under weak selection

For the case of selection, using Eqn.1 in the main text Fisher’s angular
transformation $\theta=\cos^{-1}(1-2x)$, the stochastic differential equation
for $\theta$ is

$\frac{\mathrm{d}\theta}{\mathrm{d}t}=-\frac{1}{2N}\left(\cot(\theta)-Ns\sin(\theta)\right)+\eta(t),$
(S1)

where $\langle\eta(t)\rangle=0$ and
$\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t^{\prime})/N$. We first try
to solve this by expanding the effective deterministic force,

$f(\theta)=\frac{1}{2N}(\cot(\theta)-Ns\sin(\theta)),$ (S2)

about its zero $\theta^{\ast}$ and solve the resulting linear SDE as for the
neutral case. This only works well for the case of weak selection, as the
effective deterministic force becomes increasingly non-linear for all values
of $\theta$ when selection is strong (Fig.1 main text). The solution to
$f(\theta)=0$ is

$\theta^{\ast}=\cos^{-1}\left(\frac{-1+\sqrt{1+4N^{2}s^{2}}}{2Ns}\right),$
(S3)

which for weak selection, $4Ns\ll 1$, is simply,

$\theta^{\ast}\approx\cos^{-1}\left(Ns\right).$ (S4)

Intuitively, this makes sense, as for $s>0$, this gives $\theta^{\ast}<\pi/2$,
so the dividing point (separatrix) between initial conditions that result in
deterministic dynamics giving fixation of the mutant allele is shifted to
smaller values of $\theta_{0}$ compared to neutral ($\theta_{0}=\pi/2$); the
converse is true, for $s<0$, where $\theta^{\ast}>\pi/2$. The force expanded
to linear order is $f(\theta)=\lambda(\theta-\theta^{\ast})$, where

$\lambda=f^{\prime}(\theta^{*})=\frac{1}{2N}\left(\frac{1}{1-\cos^{2}(\theta^{*})}+Ns\cos(\theta^{*})\right),$
(S5)

which again for weak selection is approximately

$\lambda=f^{\prime}(\theta^{\ast})\approx\frac{1}{2N}\left(1+2N^{2}s^{2}\right).$
(S6)

The stochastic equation of motion is then

$\frac{\mathrm{d}\theta}{\mathrm{d}t}=\lambda\left(\theta-\theta^{\ast}\right)+\eta(t).$
(S7)

The analysis then proceeds in the same way as for the neutral case giving the
Green’s function in $x-$space as

$\displaystyle G_{x}(x,x_{0};t)=\sqrt{\frac{N\lambda}{\pi x(1-x)(e^{2\lambda
t}-1)}}\exp{-\frac{\left(\cos^{-1}(1-2x)-\cos^{-1}(1-2x_{0})e^{\lambda
t}-\theta^{\ast}(1-e^{\lambda t})\right)^{2}}{\left(e^{2\lambda
t}-1\right)/N\lambda}},$ (S8)

In the regime where $Ns\ll 1$, stochastic simulations and Eqn.S8 agree well at
short times with a similar accuracy (not shown) as shown in Fig.2 in the main
text; the results show the Green’s function under weak selection are only a
small perturbation on the neutral Green’s function Eqn.6 in main text, and
only diverge significantly for long times where this approximation, in any
case fails.

## Supplementary Figures

Figure S1: Comparison of approximate calculation of Greens function for $Ns=1$
(solid lines – Eqn.10 in main text) and numerical integration of stochastic
differential equation that arises from diffusion approximation (solid
circles). a) initial frequency $x_{0}=0.01$, b) $x_{0}=0.99$.

Figure S2: Comparison of approximate calculation of Greens function for $Ns=1$
(solid lines – Eqn.10 in main text) and numerical integration of stochastic
differential equation that arises from diffusion approximation (solid
circles). a) initial frequency $x_{0}=0.1$, b) $x_{0}=0.5$, c) $x_{0}=0.9$.

We plot the results for $Ns=1$ ($x_{0}=\\{0.01,0.99\\}$ in Fig.S1 and
$x_{0}=\\{0.1,0.5,0.9\\}$ in Fig.S2). As discussed in the main text, we find
that for both $Ns=1$ and $Ns=10$ the heuristic approach and the integration of
the Wright-Fisher SDE (Eqn.7 in main text) agree very well at short times
compared to the average expected time for fixation/extinction of a mutant.

For $Ns=1$, we expect the mean time to fixation/extinction to be of order
$\sim N=1/s$ and so the Green’s functions are plotted at different times $t$
which are fractions of $N$. We see that the time range from zero that the
approximation is accurate decreases as the initial frequency $x_{0}$ is nearer
to either of the boundaries, but for sufficiently short times, even when
$x_{0}=0.01$ or $x_{0}=0.99$ (Fig.S1), the heuristic solution is very
accurate. The main difference in the Green’s functions for $Ns=1$ and $Ns=10$
are that the distributions are more narrow about their peak for $Ns=10$, which
is as expected as under stronger selection as the dynamics will be more
deterministic; we see that the heuristic approximation captures this behaviour
very accurately.

In Fig.S3, we have plotted the approximate heuristic Green’s function for
strong selection ($Ns=10$), which initial frequencies of
$x_{0}=\\{0.01,0.99\\}$. We see that the approximation is again very good for
sufficiently short times compared to the expect time to fixation or
extinction.

Figure S3: Comparison of approximate calculation of Greens function for
$Ns=10$ (solid lines – Eqn.10 in main text) and numerical integration of
stochastic differential equation that arises from diffusion approximation
(solid circles). a) initial frequency $x_{0}=0.01$, b) $x_{0}=0.99$. Green’s
functions are plotted at times given by fractions of
$\tau=\frac{1}{s}(1+\ln{(Ns)})$, which is approximately the expected time to
fixation of a mutant which survives drift and then is driven to fixation by
selection.