11institutetext: Independent researcher
11email:  22institutetext: ESILV - École Supérieure d’Ingénieurs Léonard de
Vinci, Paris, France 22email:
33institutetext: Sony Computer Science Laboratories Inc, Tokyo, Japan
33email:

# cCorrGAN: Conditional Correlation GAN for Learning Empirical Conditional
Distributions in the Elliptope

Gautier Marti 11 0000-0001-6497-5702 Victor Goubet 22 0000-0003-4078-4703
Frank Nielsen 33 0000-0001-5728-0726

###### Abstract

We propose a methodology to approximate conditional distributions in the
elliptope of correlation matrices based on conditional generative adversarial
networks. We illustrate the methodology with an application from quantitative
finance: Monte Carlo simulations of correlated returns to compare risk-based
portfolio construction methods. Finally, we discuss about current limitations
and advocate for further exploration of the elliptope geometry to improve
results.

###### Keywords:

Generative Adversarial Networks Correlation Matrices Elliptope Geometry
Empirical Distributions Quantitative Finance Monte Carlo Simulations.

##  1 Introduction

Since the seminal Generative Adversarial Networks (GANs) paper [10],
adversarial training has been successful in several areas which are typically
explored by machine learning researchers (e.g. image generation in computer
vision, voice and music generation in audio signal processing). To a smaller
extent, GANs were also applied for generating social and research citations
networks. We are not aware of existing generative adversarial networks for
sampling matrices from an empirical distribution defined in the elliptope, the
set of correlation matrices, except our previous work [18]. In this work, we
propose extensions of the original CorrGAN model keeping in mind matrix
information geometry [20].

This body of work can be motivated by applications in quantitative finance,
and possibly in other fields relying on correlation and covariance matrices.
Having access to a generative model of realistic correlation matrices can
enable large scale numerical experiments (e.g. Monte Carlo simulations) when
current mathematical guidance falls short. In the context of quantitative
finance, it can improve the falsifiability of statements and more objective
comparison of empirical methods. For example, paper A and paper B concurrently
claim their novel portfolio allocation methods are the state of the art, out-
performing the rest of the literature on well-chosen period and universe of
assets. Can we verify and obtain better insights than these claims? We will
briefly illustrate that the generative models presented in this paper can help
perform a more reliable comparison between portfolio allocation methods.
Moreover, by analyzing the simulations, we are able to extract insights on
which market conditions are more favorable to such or such portfolio
allocation methods.

##  2 Related Work

Research related to our work can be found in three distinct areas: generation
of random correlation matrices, generative adversarial networks (in finance),
information geometry of the elliptope in machine learning.

###  2.1 Sampling of Random Correlation Matrices

Seminal methods detailed in the literature are able to generate random
correlation matrices without any structure. These methods are sampling
uniformly (in a precise sense) in the set of correlation matrices, e.g. the
onion method [9]. Extended onion methods and vines can be used to generate
correlation matrices with large off-diagonal values [17], but they are not
providing any other particular structure by design. Other methods allow to
simulate correlation matrices with given eigenvalues [5], or with the
Perron–Frobenius property [14]. CorrGAN [18] is the first and so far only
method which attempts to learn a generative model from a given empirical
distribution of correlation matrices. Preliminary experiments have shown that
the synthetic correlation matrices sampled from the generative model have
similar properties than the original ones.

###  2.2 GANs in Finance

Generative Adversarial Networks (GANs) in Finance are just starting to be
explored. Most of the research papers (and many are still at the preprint
stage) were written after 2019. The main focus and motivation behind the use
of GANs in Finance is data anonymization [1], followed by univariate time
series (e.g. stock returns) modeling [8, 15, 24, 26]. Also under
investigation: modeling and generating tabular datasets which are often the
format under which the so-called ‘alternative data’ hedge funds are using is
structured; An example of such ‘alternative data’ can be credit card (and
other) transactions which can also be modeled by GANs [27]. To the best of our
knowledge, only our previous work [18] tackles the problem of modeling the
joint behaviour of a large number of assets.

###  2.3 Geometry of the Elliptope in Machine Learning

Except a recent doctoral thesis [4], the core focus of the research community
is on the Riemannian geometry of the SPD cone, e.g. [19, 2], SPDNet [13] and
its extension SPDNet with Riemannian batch normalization [3]. Concerning the
elliptope, a non-Riemannian geometry (Hilbert’s projective geometry) has also
been considered for clustering [21]. The Hilbert elliptope distance is an
example of non-separable distance which satisfies the information monotonicity
property [21].

##  3 Our Contributions

We suggest the use of GANs to learn empirical conditional distributions in the
correlation elliptope; We highlight that the correlation elliptope has
received less attention from the information geometry community in comparison
to the symmetric positive definite cone of covariance matrices; We highlight
where information geometric contributions could help improve on the proposed
models (the Euclidean Deep Convolutional Generative Adversarial Network
(DCGAN) is already competitive if not fully satisfying from a geometric point
of view); Finally, we illustrate the use of such generative models with Monte
Carlo simulations to understand the empirical properties of portfolio
allocation methods in quantitative finance.

##  4 The Set of Correlation Matrices

###  4.1 The Elliptope and its Geometrical Properties

Let $\mathcal{E}=\\{C\in\mathbb{R}^{n\times n}~{}|~{}C=C^{\top}$, $\forall
i\in\\{1,\ldots,n\\}$, $C_{ii}=1,\forall x\in\mathbb{R}^{n}$, $x^{\top}Cx\geq
0\\}$ be the set of correlation matrices, also known as the correlation
elliptope.

The convex compact set of correlation matrices $\mathcal{E}$ is a strict
subspace of the set of covariance matrices whose Riemannian geometry has been
well studied in the information geometry literature, e.g. Fisher-Rao distance
[25], Fréchet mean [19], PCA [12], clustering [23, 16], Riemannian-based
kernel in SVM classification [2]. However, $\mathcal{E}$ endowed with the
Riemannian Fisher-Rao metric is not totally geodesic (cf. Figures 2, 2), and
this has strong implications on the principled use of the tools developed for
the covariance matrices. For example, consider SPDNet [13] (with an eventual
Riemannian batch normalization layer [3]). Despite this neural network is fed
with correlation matrices as input, it will generate covariance matrices in
subsequent layers. How can we stay in the elliptope in a principled and
computationally efficient fashion?

In Figure 2, each $2\times 2$ covariance matrix is represented by a 3D point
$(x,y,z)$. The blue segment $(x=z=1)$ is the set of $2\times 2$ correlation
matrices. In green, the geodesic (using the Riemannian Fisher-Rao metric for
covariances) between correlation matrix $\gamma(0)=(1,-0.75,1)$ and
correlation matrix $\gamma(1)=(1,0.75,1)$. The geodesic
$t\in[0,1]\rightarrow\gamma(t)$ (and in particular the Riemannian mean
$\gamma(0.5)$) is not included in the blue segment representing the
correlation matrix space.

In Figure 2 we compare several means:

  * •

(M1) Euclidean mean $\frac{A+B}{2}$,

  * •

(M2) Riemannian barycenter (in general, a covariance matrix $\Sigma^{\star}$),

  * •

(M3)
$C^{\star}=\mathrm{diag}(\Sigma^{\star})^{-\frac{1}{2}}\Sigma^{\star}\mathrm{diag}(\Sigma^{\star})^{-\frac{1}{2}}$
($\Sigma^{\star}$ normalized by variance),

  * •

(M4) Fréchet mean constrained to the space of correlation matrices,

  * •

(M5) projection of $\Sigma^{\star}$ onto the elliptope using the Riemannian
distance.

We observe that (M3) is close but not equal to (M4) and (M5), which seems to
be equivalent. We have not found yet how to compute (M4) or (M5) efficiently.
Alternatively, a recent doctoral dissertation [4] on Riemannian quotient
structure for correlation matrices might help progress in this direction.

Figure 1: $2\times 2$ correlation matrices (blue); Fisher-Rao geodesic (green)
between two matrices. The geodesic does not stay inside the set of correlation
matrices.

Figure 2: Visualizing various means.

With a better understanding on how to define and efficiently compute a
Riemannian mean for correlation matrices, we can think of adapting SPDNet to
correlation matrices, and eventually use it as a component of the generative
adversarial network architecture.

###  4.2 Financial Correlations: Stylized Facts

Since we aim at learning empirical (conditional) distributions, different
domains may have different type of typical correlation matrices whose set
describes a subspace of the full correlation elliptope. In the present work,
we will showcase the proposed model on correlation matrices estimated on
stocks returns. We briefly describe the known properties, also known as
stylized facts, of these correlation matrices. We leave for future work or
research collaborations to explore the empirical space of correlation matrices
in other domains.

The properties of financial correlation matrices are well-known and documented
in the literature. However, until CorrGAN [18], no single model was able to
generate correlation matrices verifying a handful of them at the same time:
(SF 1) Distribution of pairwise correlations is significantly shifted to the
positive; (SF 2) Eigenvalues follow the Marchenko–Pastur distribution, but for
a very large first eigenvalue; (SF 3) Eigenvalues follow the Marchenko–Pastur
distribution, but for $\approx 5\%$ of other large eigenvalues; (SF 4) Perron-
Frobenius property (first eigenvector has positive entries); (SF 5)
Hierarchical structure of correlation clusters; (SF 6) Scale-free property of
the corresponding Minimum Spanning Tree.

When markets are under stress, it is widely known among practitioners that the
average correlation is high. More precisely, the whole correlation structure
changes and becomes different than the ones typical of rallying or steady
markets. These characteristics tied to the market regime would be even harder
to capture with a fully specified mathematical model.

##  5 Learning Empirical Distributions in the Elliptope using Generative
Adversarial Networks

We suggest that conditional generative adversarial networks are promising
models to generate random correlation matrices verifying all these not-so-
well-specified properties.

###  5.1 CorrGAN: Sampling Realistic Financial Correlation Matrices

The original CorrGAN [18] model is essentially a DCGAN (architecture based on
CNNs). CNNs have 3 relevant properties which align well with stylized facts of
financial correlation matrices: shift-invariance, locality, and
compositionality (aka hierarchy). Based on our experiments, this architecture
yields the best off-the-shelf results most consistently. However, outputs
(synthetic matrices) need to be post-processed by a projection algorithm (e.g.
[11]) as they are not strictly valid (but close to) correlation matrices. This
goes against the end-to-end optimization of deep learning systems. Despite
CorrGAN yields correlation matrices which verify the stylized facts, we can
observe with PCA (or MDS, t-SNE) that the projected empirical and synthetic
distributions do not match perfectly. We also noticed some instability from
one trained model to another (well-known GAN defect). For example, in our
experiments, the average Wasserstein distance (using the POT package [7])
between any two (PCA-projected) training sets is
$\mu_{E}:=6.7\pm\sigma_{E}:=6.8$ (max distance := 15) whereas the average
distance between a training set and a generated synthetic set is
$\mu_{G}:=18.8\pm\sigma_{G}:=8$ (min distance := 8). Ideally, we would obtain
$\mu_{E}\approx\mu_{G}$ instead of $\mu_{E}\ll\mu_{G}$. We observed that
CorrGAN was not good at capturing the different modes of the distribution.
This motivated the design of a conditional extension of the model.

###  5.2 Conditional CorrGAN: Learning Empirical Conditional Distributions in
the Elliptope

We describe in Figure 4 the architecture of cCorrGAN, the conditional
extension of the CorrGAN model [18]. To train this model, we used the
following experimental protocol: 1. We built a training set of correlation
matrices labeled with one of the three classes $\\{$stressed, normal,
rally$\\}$ depending on the performance of an equally weighted portfolio over
the period of the correlation estimation; 2. We trained several cCorrGAN
models; 3. We evaluated them qualitatively (PCA projections, reproductions of
stylized facts), and quantitatively (confusion matrix of a SPDNet classifier
[3], Wasserstein distance between projected empirical and synthetic
distributions) to select the best model. Samples obtained from the best model
are shown in Figure 4.

Figure 3: Architecture of the conditional GAN (cCorrGAN) which generates
$80\times 80$ correlation matrices corresponding to a given market regime
$\\{$stressed, normal, rally$\\}$

Figure 4: Empirical (left) and cCorrGAN generated (right) correlation matrices
for three market regimes

##  6 Application: Monte Carlo simulations of risk-based portfolio allocation
methods

Monte Carlo simulations to test the robustness of investment strategies and
portfolio construction methods were the original motivation to develop
CorrGAN, and its extension conditional CorrGAN. The methodology is summarized
in Figure 5: 1. Generate a random correlation matrix from a given market
regime using cCorrGAN, 2. extract features from the correlation matrix, 3.
estimate the in-sample and out-of-sample risk on synthetic returns which
verify the generated correlation structure, 4. predict and explain (with
Shapley values) a) performance decay from in-sample to out-of-sample, b) out-
performance of a portfolio allocation method over another. The idea of using
Shapley values to explain portfolio construction methods can be found in [22].

Figure 5: Methodology for analysing the results of the Monte Carlo simulations

Using such methodology, enabled by the sampling of realistic random
correlation matrices conditioned on a market regime, we can find, for example,
that the Hierarchicaly Risk Parity (HRP) [6] outperforms the naive risk parity
when there is a strong hierarchy in the correlation structure, well-separated
clusters, and a high dispersion in the first eigenvector entries, typical of a
normal or rallying market regime. In the case of a stressed market where
stocks plunge altogether, HRP performs on-par (or slightly worse because of
the risk of identifying a spurious hierarchy) than the simpler naive risk
parity.

Acknowledgments: We thank ESILV students for helping us running experiments
for hundreds of GPU-days: Chloé Daems, Thomas Graff, Quentin Bourgue, Davide
Coldebella, Riccardo Rossetto, Mahdieh Aminian Shahrokhabadi, Marco Menotti,
Johan Bonifaci, Clément Roux, Sarah Tailee, Ekoue Edem Ekoue-Kouvahey, Max
Shamiri.

## References

  * [1] Assefa, S.: Generating synthetic data in finance: opportunities, challenges and pitfalls. Challenges and Pitfalls (June 23, 2020) (2020)
  * [2] Barachant, A., Bonnet, S., Congedo, M., Jutten, C.: Classification of covariance matrices using a Riemannian-based kernel for BCI applications. Neurocomputing 112, 172–178 (2013)
  * [3] Brooks, D., Schwander, O., Barbaresco, F., Schneider, J.Y., Cord, M.: Riemannian batch normalization for SPD neural networks. In: Advances in Neural Information Processing Systems. vol. 32, pp. 15489–15500 (2019)
  * [4] David, P.: A Riemannian Quotient Structure for Correlation Matrices with Applications to Data Science. Ph.D. thesis, The Claremont Graduate University (2019)
  * [5] Davies, P.I., Higham, N.J.: Numerically stable generation of correlation matrices and their factors. BIT Numerical Mathematics 40(4), 640–651 (2000)
  * [6] De Prado, M.L.: Building diversified portfolios that outperform out of sample. The Journal of Portfolio Management 42(4), 59–69 (2016)
  * [7] Flamary, R., Courty, N.: POT Python Optimal Transport library (2017), https://pythonot.github.io/
  * [8] Fu, R., Chen, J., Zeng, S., Zhuang, Y., Sudjianto, A.: Time series simulation by conditional generative adversarial net. arXiv preprint arXiv:1904.11419 (2019)
  * [9] Ghosh, S., Henderson, S.G.: Behavior of the NORTA method for correlated random vector generation as the dimension increases. ACM Transactions on Modeling and Computer Simulation (TOMACS) 13(3), 276–294 (2003)
  * [10] Goodfellow, I.J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial networks. arXiv preprint arXiv:1406.2661 (2014)
  * [11] Higham, N.J.: Computing the nearest correlation matrix—a problem from finance. IMA journal of Numerical Analysis 22(3), 329–343 (2002)
  * [12] Horev, I., Yger, F., Sugiyama, M.: Geometry-aware principal component analysis for symmetric positive definite matrices. In: Asian Conference on Machine Learning. pp. 1–16. PMLR (2016)
  * [13] Huang, Z., Van Gool, L.: A Riemannian Network for SPD Matrix Learning. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 31 (2017)
  * [14] Hüttner, A., Mai, J.F.: Simulating realistic correlation matrices for financial applications: correlation matrices with the Perron–Frobenius property. Journal of Statistical Computation and Simulation 89(2), 315–336 (2019)
  * [15] Koshiyama, A., Firoozye, N., Treleaven, P.: Generative adversarial networks for financial trading strategies fine-tuning and combination. Quantitative Finance pp. 1–17 (2020)
  * [16] Lee, H., Ahn, H.J., Kim, K.R., Kim, P.T., Koo, J.Y.: Geodesic clustering for covariance matrices. Communications for Statistical Applications and Methods 22(4), 321–331 (2015)
  * [17] Lewandowski, D., Kurowicka, D., Joe, H.: Generating random correlation matrices based on vines and extended onion method. Journal of multivariate analysis 100(9), 1989–2001 (2009)
  * [18] Marti, G.: CorrGAN: Sampling realistic financial correlation matrices using generative adversarial networks. In: ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). pp. 8459–8463. IEEE (2020)
  * [19] Moakher, M.: A differential geometric approach to the geometric mean of symmetric positive-definite matrices. SIAM Journal on Matrix Analysis and Applications 26(3), 735–747 (2005)
  * [20] Nielsen, F., Bhatia, R.: Matrix information geometry. Springer (2013)
  * [21] Nielsen, F., Sun, K.: Clustering in Hilbert’s projective geometry: The case studies of the probability simplex and the elliptope of correlation matrices. In: Geometric Structures of Information, pp. 297–331. Springer (2019)
  * [22] Papenbrock, J., Schwendner, P., Jaeger, M., Krügel, S.: Matrix evolutions: Synthetic correlations and explainable machine learning for constructing robust investment portfolios. Available at SSRN (2020)
  * [23] Shinohara, Y., Masuko, T., Akamine, M.: Covariance clustering on Riemannian manifolds for acoustic model compression. In: 2010 IEEE International Conference on Acoustics, Speech and Signal Processing. pp. 4326–4329. IEEE (2010)
  * [24] Takahashi, S., Chen, Y., Tanaka-Ishii, K.: Modeling financial time-series with generative adversarial networks. Physica A: Statistical Mechanics and its Applications 527, 121261 (2019)
  * [25] Wells, J., Cook, M., Pine, K., Robinson, B.D.: Fisher-Rao distance on the covariance cone. arXiv preprint arXiv:2010.15861 (2020)
  * [26] Wiese, M., Knobloch, R., Korn, R., Kretschmer, P.: Quant GANs: Deep generation of financial time series. Quantitative Finance 20(9), 1419–1440 (2020)
  * [27] Zheng, Y.J., Zhou, X.H., Sheng, W.G., Xue, Y., Chen, S.Y.: Generative adversarial network based telecom fraud detection at the receiving bank. Neural Networks 102, 78–86 (2018)


# Diffusion approximations in population genetics and the rate of Muller’s
ratchet

Camila Bräutigam  Max Planck Institute for Mathematics in the Sciences,
Leipzig, Germany  Matteo Smerlak<EMAIL_ADDRESS>Max Planck Institute for
Mathematics in the Sciences, Leipzig, Germany

###### Abstract

Diffusion theory is a central tool of modern population genetics, yielding
simple expressions for fixation probabilities and other quantities that are
not easily derived from the underlying Wright-Fisher model. Unfortunately, the
textbook derivation of diffusion equations as scaling limits requires
evolutionary parameters (selection coefficients, mutation rates) to scale like
the inverse population size—a severe restriction that does not always reflect
biological reality. Here we note that the Wright-Fisher model can be
approximated by diffusion equations under more general conditions, including
in regimes where selection and/or mutation are strong compared to genetic
drift. As an illustration, we use a diffusion approximation of the Wright-
Fisher model to improve estimates for the expected time to fixation of a
strongly deleterious allele, i.e. the rate of Muller’s ratchet.

##  I Introduction

The Wright-Fisher (WF) model captures in simple mathematical form the combined
effects of natural selection, mutation, and genetic drift, and therefore plays
a foundational role in quantitative evolutionary theory and population
genetics. Unfortunately, the main quantities of interest such as fixation
probabilities cannot be computed in closed form within this model. An
established approximation strategy relates the discrete WF model to a
continuous diffusion process in allele frequency space obtained via a scaling
limit in which the population size and timescale go to infinity in a
correlated manner. This scaling transformation was first used by Fisher [1],
Wright [2] and Kimura [3], and more rigorous expositions were given by
Kolmogorov [4], Malécot [5] and Feller [6]. This approach plays a central role
in most textbooks expositions [7, 8, 9, 10], and it has been said that “the
standard diffusion approximation has permeated the field so thoroughly that it
shapes the way in which workers think about the genetics of populations” [11].

A key limitation of the WF diffusion equation, however, is its narrow range of
applicability: for the deterministic forces of selection and mutation to
survive the scaling limit, selection coefficients and mutation rates must
scale like the inverse effective population size $1/N$. This condition is
neither biologically well motivated (selection coefficients and mutation rates
do not covary with population size) nor universally applicable (viral
populations, for instance, can have large effective sizes
$\mathcal{O}(10^{7})$ [12]). Perhaps because the scaling limit of the WF model
is often presented as _the_ diffusion approximation of population genetics,
diffusion methods are sometimes believed to be restricted to the weak
selection-weak mutation regime of evolution: “In order to use diffusion
approximations it is necessary that the mutation rates be of no larger order
of magnitude than $1/N$” [13]. Outside this regime, the common practice is to
either assume neutral evolution and effectively set all selection coefficients
and mutation rates to zero, as in Kimura’s original work [3], or to use
deterministic equations [14], _i.e._ to neglect genetic drift altogether.

Nevertheless, there are important cases where genetic drift is small compared
to selection and mutation, and yet plays a key role in driving evolutionary
change. An example is _Muller’s ratchet_ [15, 16], _viz._ the irreversible
accumulation of slightly deleterious mutations in finite asexual populations.
This process has been invoked in the evolution of mitochondria [17], RNA
viruses [18], Y chromosomes [19], ageing [20], cancer [21], and sex [22].
Obtaining analytical expressions for the mean click time of the ratchet is
therefore an important goal for theory.

At first sight, diffusion theory seems ideally suited for this problem:
Muller’s ratchet is akin to a fixation problem for a two-type Wright-Fisher
process, where one type represents unloaded alleles, and the other type
comprises all other, deleterious mutations [23, 19, 24]. However, recent works
have reinforced the belief that diffusion theory is inapplicable under finite
selection strengths [25, 26]. An emerging consensus is that diffusion theory
just does not apply to Muller’s ratchet, at least not in the slow click
regime: “The fact that the extinction of the fittest class is due to […] a
rare, large fluctuation and not […] a typical fluctuation prohibits simple
diffusive treatments of the ratchet” [26].111Ironically, Eq. (23) of [26]
coincides with the large $N$ Laplace approximation of the mean click time for
the WF diffusion equation, Eq. (4) below.

This paper has two objectives. First, we stress that diffusion methods in
population genetics do _not_ in fact require $s,u,v\sim N^{-1}$ ($s$ selection
coefficient, $u,v$ mutation rates) to correctly capture the interplay between
mutation and selection in large populations; an alternative form of the WF
diffusion equation, obtained as an interpolation of the WF process rather than
as a scaling limit, provides an excellent approximation valid under the more
general conditions $s,u,v\ll 1\ll N$. The fact that diffusion theory can be
used to derive more general results than the usual classical Wright-Fisher
diffusion equation has been noted in the past (see e.g. [11] and references
therein), but has not been sufficiently appreciated in our view. Second, we
use the alternative approximation by diffusion to derive a new analytical
formula for the mean fixation time of a deleterious allele, a central quantity
in the context of Muller’s ratchet. We find that not only can diffusion theory
be applied to this problem, it actually provides a better approximation than
previous estimates in [25, 26] based on WKB expansions.

##  II Scaling limit vs. interpolation

We focus on the haploid WF model for two types $A$ and $B$. For a population
with effective size $N$, this is the Markov chain on the set
$\\{0,1/N,\cdots,(N-1)/N,1\\}$ defined by the binomial sampling formula

$NX_{n+1}\sim\textrm{Binom}(N,p(X_{n})).$ (1)

Here $X_{n}$ denotes the frequency of $A$ types at generation $n$, and $p(x)$
represents the probability that an offspring is of type $A$ given that the
frequency of $A$ in the parental population is $x$. The form of the function
$p(x)$ depends on the specifics of mutation-selection dynamics. For instance,
if $A$ has $(1+s)$ times more offspring than $B$, and the probability of a
mutation $A\to B$ (resp. $B\to A$) is $u$ (resp. $v$), we have

$p(x)=\frac{(1+s)(1-u)x+v(1-x)}{1+sx}.$ (2)

Other functional forms are more appropriate for viability selection, etc.

Figure 1: Stationary distributions as result of $10^{6}$ simulated generations
of the WF process (Eq. (1), yellow histograms) compared to analytical
predictions from the WF diffusion (Eq. (4), blue line) and the diffusion by
interpolation (Eq. (7), orange line). The parameters $s,u,v$ as given in the
subplots are kept constant when moving vertically downwards where population
size $N$ increases. Deviations of the WF diffusion from the underlying
discrete model become apparent for strong selection $N\gg 1/s^{2}$. Figure 2:
Relative errors in approximations of the establishment time of a rare mutant
subject to reversions ($s,u\neq 0$). The errors of the two analytical
approximations, WF diffusion (Eq. (4), in blue) and diffusive approximation by
interpolation (Eq. (7), in orange) are computed relative to results from
simulations of the discrete WF process (Eq. (1)).

In the limit of large population sizes, the standard diffusion approximation
considers the WF model on the slow timescale $T=nN$, and assumes that
selection and mutation are comparable in strength to genetic drift, i.e. $Ns$,
$Nu$, and $Nv$ have finite limits $\sigma,\mu,\nu$ as $N\to\infty$. Under
these conditions, it can be shown that the paths of the WF model on the slow
timescale converge to those of the diffusion process

$\displaystyle dX_{T}=[\sigma X_{T}(1-X_{T})-\mu
X_{T}+(1-\nu)X_{T}]dT+\sqrt{X_{T}(1-X_{T})}dW_{T}$ (3)

where $dW_{T}$ denotes a standard Wiener process and the Ito convention for
stochastic differential equations is used. When brought back to the fast
timescale $t=T/N$, this process becomes

$\displaystyle
dX_{t}=[sX_{t}(1-X_{t})-uX_{t}+(1-v)X_{t}]dt+\sqrt{\frac{X_{t}(1-X_{t})}{N}}dW_{t}.$
(4)

Eq. (4) above is the “WF diffusion equation” found in textbooks. By
construction, this approximation relies on the condition that all evolutionary
parameters scale as $1/N$ when $N$ is taken to infinity, corresponding to the
“weak selection-weak mutation” regime of evolution [7]. A rule of thumb due to
Nei states that selection can be considered small if $Ns^{2}\ll 1$ [27, 25].

There is, however, another way to approximate the WF model with a diffusion
process. By the central limit theorem, Eq. (1) may be written as

$X_{n+1}\sim p(X_{n})+\sqrt{\frac{p(X)(1-p(X_{n}))}{N}}\nu$ (5)

where $\nu$ approaches a standard normal variable if $N\gg 1$. The latter
equation is the Euler-Maruyama discretization of the Ito diffusion equation

$dX_{t}=[p(X_{t})-X_{t}]dt+\sqrt{\frac{p(X_{t})(1-p(X_{t}))}{N}}dW_{t}$ (6)

at times $t=0,1,\cdots$. If the fraction $X_{t}$ changes only slightly during
each time increment (i.e. if $s,u,v\ll 1\ll N$), then Eq. (6) is an accurate
_interpolation_ of the original process, Eq. (1). We simplify Eq. (6) further
as

$dX_{t}=[p(X_{t})-X_{t}]dt+\sqrt{\frac{X_{t}(1-X_{t})}{N}}dW_{t}.$ (7)

Versions of Eq. (7) appear in e.g. [24, 28], but, to our knowledge, have never
been used to derive analytical results valid beyond the weak selection-weak
mutation limit of evolutionary dynamics. A mathematical discussion of neutral
diffusions at large mutation rates can be found in [29, 30].

##  III Stationary distributions and establishment times

#### Stationary distributions

That the diffusive approximation in Eq. (7) has a broader domain of validity
than the standard WF diffusion in Eq. (4) can be seen by comparing their
stationary distributions. We recall that, given a stochastic differential
equation of the kind described above,

$dX_{t}=a(X_{t})\,dt+\sqrt{b(X_{t})}\,dW_{t},$ (8)

a stationary distribution $P_{\text{eq}}$ exists as long as both boundaries,
$x=0$ and $x=1$, are reflecting ($v\neq 0$, $u\neq 0$) and can be calculated
as [31]

$P_{\text{eq}}(x)\propto\frac{\exp{\left[2\int^{x}\frac{a(y)}{b(y)}\mathop{}\\!\mathrm{d}y\right]}}{b(x)}.$
(9)

In Fig. 1 we compare the analytical expressions from both diffusive
approximations for $P_{\text{eq}}$ with simulations of the WF process with
$p(x)$ from Eq. (2). Results are shown for a range of parameters $s,u,v$, both
within and outside the weak selection-weak mutation regime defined by
$s,u,v=\mathcal{O}(1/N)$.

Figure 3: Accuracy of different analytical expressions for the mean time
between clicks of Muller’s ratchet in the slow click regime, as a function of
$N$ for $S=U=0.01$ (left) and as a function of $S=2U$ for $N=100$ (right). Our
result (Eq. (10), red line), derived from the diffusion interpolation (Eq.
(7)) through Laplace’s method correctly approximates Wright-Fisher simulations
(mean click time over $10^{3}$ replicates, represented as black dots) of the
two-type reduced model of Muller’s ratchet, including when selection and
mutation are strong ($S=U\gg 1/N$).

#### Establishment times with one-way mutations

The classical fixation problem in population genetics considers the fate of a
single selectively advantageous mutant ($s>0$) subject to finite population
size fluctuations (genetic drift) only. Mutation rates are zero ($u=v=0$).
Diffusive approximations can be used in that case to derive an analytical
estimate for the probability of fixation at the boundaries $x=0$ or $x=1$ and
corresponding expected times to fixation. For completeness, we report the
results of the classical fixation problem in Appendix A, but note that the
differences between predictions from the two diffusive schemes are only minor.

Instead, we focus here on a modified fixation problem in which the
advantageous mutant is additionally subjected to one-way mutations ($u\neq 0$,
$v=0$). As in the classical setting, when a single mutant $X_{0}=1/N$ appears
in the population, the probability of immediate extinction is high (due to
genetic drift). Still, there is a finite probability of reproduction and
establishment. Establishment means here reaching a metastable state $x_{c}>0$,
defined by $p(x_{c})=x_{c}$. The boundary $x=1$ is no longer absorbing due to
reversions at a finite rate $u$ and the mutant will almost certainly go
extinct in finite time, since $x=0$ represents the only absorbing boundary.
The intermittent dynamics may still be of interest and can be characterised by
the probability and time of establishment.

In Fig. 2 we compare the relative errors of mean establishment times
$\overline{t}_{\textrm{est}}$ between WF simulations and predictions derived
from the two diffusive approximations (Eq. (4) and Eq. (7)). Details for the
calculations are found in the Appendix. Here, again, the relative errors for
the prediction of the WF diffusion increase outside the weak selection regime
$s^{2}>1/N$, unlike the estimates derived from the diffusive approximation by
interpolation (Eq. (7)). Details of the calculations can be found in the
Appendix A.

##  IV Rate of Muller’s ratchet

Haigh’s model of Muller’s ratchet [32] considers a population of fixed size
$N$ where deleterious mutations accumulate at a constant rate $U$, such that
an individual carrying $k$ mutations has reduced fitness
$\left(1-S\right)^{k}$. When $Ne^{-U/S}>1$, fixations are rare and clicks of
the ratchet are large deviations (rare fluctuations exponentially suppressed
in $N$) this is commonly referred to as the slow click regime. Although
significant progress was made in the last decade [33], Haigh’s model is
difficult to study analytically, and a common approach is to approximate it by
a two-type WF model [24, 25, 26]. In this simpler model, one class corresponds
to the fittest (unmutated) individuals, with unit fitness, and all other
mutants are collectively assigned a reduced fitness $(1-s)$ with
$s=(1-e^{-U})/(1-e^{-U/S})$; the effective mutation rate away from the fittest
class is then $u=1-e^{-U}$. Assuming fecundity selection following [25, 26],
we arrive at the WF model in Eq. (1) with sampling probability
$p(x)=x(1-u)/(1-s+sx)$. In contrast to Eq. (2) the selection coefficient is
redefined as a coefficient of negative selection, hence the normalisation term
of mean fitness is adjusted accordingly.

We used classical expressions for hitting times of one-dimensional diffusions
[7] combined with Laplace’s approximation of integrals to derive from Eq. (7)
an analytical expression for the expected time between clicks of the ratchet;
technical details are provided in the Appendix. On a logarithmic scale, this
expression reads

$\log\mathbb{E}(T_{\textrm{click}})\underset{N\to\infty}{\sim}2N\log\left[\left(\frac{u}{s}\right)^{u}\left(\frac{1-u}{1-s}\right)^{1-u}\right],$
(10)

compared to the approximation that is found by [26]

$\log\mathbb{E}(T_{\textrm{click}})\sim 2N[s-u+u\log(u/s)],$ (11)

and corresponds to the expected hitting time of the WF diffusion (Eq. (4)
instead of Eq. (7)). Fig. 3 compares the accuracy of our result with earlier
heuristic diffusion approximations [19, 34] as well as more recent WKB
expansions of Moran models [26]. The near-perfect agreement in the large $N$
limit—including when $S,U\gg 1/N$—confirms that diffusion methods are more
general than scaling limits and can be successfully applied outside the weak
selection-weak mutation regime of evolution.

##  V Conclusion

Because it does not depend on the mode (fecundity or viability) of selection,
on generations being overlapping or non-overlapping, etc., the classical WF
diffusion equation (4) is sometimes described as “universal”. What this means
is simply that, in the special regime where $s,u,v$ are all very small
relative to the finite size fluctuations ($s,u,v=\mathcal{O}(1/N)$), these
differences are immaterial for the long-term dynamics of the population. But
that is no longer true in regimes where selection and/or mutation are much
stronger than genetic drift. In these regimes, diffusion theory does not break
down; instead, other, non-universal diffusion approximations of the WF model
can be used to obtain accurate results for fixation times and other quantities
of interest.

Importantly, diffusion equations such as Eq. (7) can be used to estimate the
timescales of large deviations, as observed e.g. in the slow clicking regime
of Muller’s ratchet. This conclusion is in contrast with previously expressed
concerns regarding the use of diffusion theory in the context of large
deviations [35, 36, 37, 38]: “The fact that the extinction of the fittest
class is due to such a rare, large fluctuation and not the cause of a typical
fluctuation prohibits simple diffusive treatments of the ratchet” [26]. Our
results show, instead, that the analytical approximation of Muller’s ratchet
within a diffusive framework improves previous results in the strong selection
regime $s>1/\sqrt{N}$. We expect that applications of diffusive approximations
beyond the ones here presented might allow further progress in the analysis of
evolutionary processes that are well captured by the discrete-time Wright-
Fisher model, but not by its textbook diffusion approximation.

## Acknowledgments

We thank Maseim Kenmoe for help with Laplace integrals, Aleksander Klimek and
Anton Zadorin for comments on the manuscript, and Sophie Pénisson for pointing
us to the mathematical references [29, 30]. Funding for this work was provided
by the Alexander von Humboldt Foundation in the framework of the Sofja
Kovalevskaja Award endowed by the German Federal Ministry of Education and
Research.

##  Appendix A Appendix

###  A.1 Derivation of Muller’s ratchet rate formula

Eq. (7) from the main text describes the dynamics of the stochastic variable
$X_{t}$ in terms of a SDE. Moving to an alternative dynamical description in
terms of the probability density function $\rho(x,t)$ of the stochastic
variable $x$ at time $t$, the corresponding drift-diffusion equation (also
Fokker-Planck equation or Kolmogorov forward equation) takes the following
form

$\frac{\partial\rho(x,t)}{\partial t}=-\frac{\partial}{\partial
x}\left[a(x)\,\rho(x,t)\right]+\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\left[b(x)\,\rho(x,t)\right],$ (12)

where in our setting $a(x)=p(x)-x$ and $b(x)=x(1-x)/N$.
We are interested in the expected time for the stochastic variable $x$ to hit
the boundary at $x=0$, given an initial position $x=x_{0}$ (to be specified
below). In the scope of diffusion theory this is referred to as the mean first
passage time $\overline{t}(x_{0})$, which is described by a second order
differential equation [31]

$a(x_{0})\frac{\mathop{}\\!\mathrm{d}\overline{t}(x_{0})}{\mathop{}\\!\mathrm{d}x_{0}}+\frac{1}{2}b(x_{0})\frac{\mathop{}\\!\mathrm{d}^{2}\overline{t}(x_{0})}{\mathop{}\\!\mathrm{d}x_{0}^{2}}=-1.$
(13)

In the setting of Muller’s ratchet, the diffusion of Eq. (12) is bounded by
definition to the interval $[0,1]$ (since $x$ describes a relative frequency).
Solving Eq. (13), subject to the appropriate boundary conditions for
reflection at $x=1$ and absorption at $x=0$, namely

$\frac{\mathop{}\\!\mathrm{d}\overline{t}(x)}{\mathop{}\\!\mathrm{d}x}\rvert_{x=1}=0\qquad\text{and}\qquad\overline{t}(0)=0,$

the expected click time of Muller’s ratchet can be expressed in terms of the
following equation

$\begin{split}\overline{t}(x_{0})&=\int_{0}^{1}\frac{2}{b(x)}\exp{\left(\Psi(x)\right)}\,\times\\\
&\int_{0}^{\min[x,x_{0}]}\exp{\left(-\Psi(y)\right)}\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}x,\end{split}$
(14)

where $\Psi(x)=\int^{x}\frac{2a(z)}{b(z)}\mathop{}\\!\mathrm{d}z$. Eq. (14)
cannot be solved analytically, but instead the proportionality of the argument
in the exponential $\Psi(x)\propto N$ to the population parameter $N$, renders
the integral equation suitable for a Laplace approximation. We therefore
expand $\Psi(x)$ in a Taylor series around a critical point $x_{c}$ and
truncate the expansion by assuming
$(x-x_{c})=\mathcal{O}\left(\frac{1}{\sqrt{N}}\right)$, so that

$\begin{split}\Psi(x)&=\Psi(x_{c})+\Psi^{\prime}(x_{c})\left(x-x_{c}\right)\\\
&+\frac{1}{2}\Psi^{\prime\prime}(x_{c})\left(x-x_{c}\right)^{2}+\mathcal{O}\left(\frac{1}{\sqrt{N}}\right),\end{split}$
(15)

since $N(x-x_{c})^{3}=\mathcal{O}\left(\frac{1}{\sqrt{N}}\right)$. Hence, the
exponential integral is approximately reduced to a Gaussian integral on a
bounded interval, which can be evaluated in terms of error functions.
Generally speaking,

$\begin{gathered}\int_{a}^{b}\exp\left(\Phi(x)\right)\approx\int_{a}^{b}\exp\left[\Phi(x_{c})+\Phi^{\prime}(x_{c})\,(x-x_{c})+\frac{1}{2}\Phi^{\prime\prime}(x_{c})\,(x-x_{c})^{2}\right]\mathop{}\\!\mathrm{d}x\\\
=C\times\,\int_{a}^{b}\exp\left[\left(\sqrt{\frac{\Phi^{\prime\prime}(x_{c})}{2}}(x-x_{c})+\frac{\Phi^{\prime}(x_{c})}{\sqrt{2\Phi^{\prime\prime}(x_{c})}}\right)^{2}\right]\mathop{}\\!\mathrm{d}x\\\
=C\times\sqrt{\frac{2}{\Phi^{\prime\prime}(x_{c})}}{i\mkern
1.0mu}\int_{{i\mkern
1.0mu}\frac{(a-x_{c})\Phi^{\prime\prime}(x_{c})+\Phi^{\prime}(x_{c})}{\sqrt{2\Phi^{\prime\prime}(x_{c})}}}^{{i\mkern
1.0mu}\frac{(b-x_{c})\Phi^{\prime\prime}(x_{c})+\Phi^{\prime}(x_{c})}{\sqrt{2\Phi^{\prime\prime}(x_{c})}}}\exp(-\chi^{2})\,\mathop{}\\!\mathrm{d}\chi\\\
=C\sqrt{\frac{\pi}{2\Phi^{\prime\prime}(x_{c})}}\left(\operatorname{erfi}\left[\frac{(b-x_{c})\Phi^{\prime\prime}(x_{c})+\Phi^{\prime}(x_{c})}{\sqrt{2\Phi^{\prime\prime}(x_{c})}}\right]-\operatorname{erfi}\left[\frac{(a-x_{c})\Phi^{\prime\prime}(x_{c})+\Phi^{\prime}(x_{c})}{\sqrt{2\Phi^{\prime\prime}(x_{c})}}\right]\right),\end{gathered}$
(16)

where
$C=\exp\left[\Phi(x_{c})-\frac{\Phi^{\prime}(x_{c})^{2}}{2\Phi^{\prime\prime}(x_{c})}\right]$
and $\operatorname{erfi}(x)$ refers to the imaginary error function.
Considering in Eq. (14) an additional expansion of the non-exponential
function around the critical point
$g(x)=\frac{2}{b(x)}=g(x_{c})+\mathcal{O}\left(\frac{1}{\sqrt{N}}\right)$, as
well as an approximation of the limit $\min[x,x_{0}]\approx\min[x_{c},x_{0}]$
we notice that the previously nested integrals decouple and we can deduce the
full approximation by applying the Laplace approximation (Eq. (16)) to both
integrals separately

$\begin{gathered}\overline{t}(x_{0})=\int_{0}^{1}\frac{2}{b(x)}\exp{\left(\Psi(x)\right)}\int_{0}^{\min[x,x_{0}]}\exp{\left(-\Psi(y)\right)}\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}x\\\
\approx
g(x_{c})\frac{\pi}{2\sqrt{-\Psi^{\prime\prime}(y_{c})\Psi^{\prime\prime}(x_{c})}}\exp\left[\Psi(x_{c})-\frac{\Psi^{\prime}(x_{c})^{2}}{2\Psi^{\prime\prime}(x_{c})}-\Psi(y_{c})+\frac{\Psi^{\prime}(y_{c})^{2}}{2\Psi^{\prime\prime}(y_{c})}\right]\times\\\
\left(\operatorname{erfi}\left[\frac{-(\min[x_{c},x_{0}]-y_{c})\Psi^{\prime\prime}(y_{c})-\Psi^{\prime}(y_{c})}{\sqrt{-2\Psi^{\prime\prime}(y_{c})}}\right]-\operatorname{erfi}\left[\frac{-(0-y_{c})\Psi^{\prime\prime}(y_{c})-\Psi^{\prime}(y_{c})}{\sqrt{-2\Psi^{\prime\prime}(y_{c})}}\right]\right)\times\\\
\left(\operatorname{erfi}\left[\frac{(1-x_{c})\Psi^{\prime\prime}(x_{c})+\Psi^{\prime}(x_{c})}{\sqrt{2\Psi^{\prime\prime}(x_{c})}}\right]-\operatorname{erfi}\left[\frac{(0-x_{c})\Psi^{\prime\prime}(x_{c})+\Psi^{\prime}(x_{c})}{\sqrt{2\Phi^{\prime\prime}(x_{c})}}\right]\right).\end{gathered}$
(17)

For the ratchet click time in the slow-clicking regime, we assume the initial
position to be precisely at the deterministic equilibrium, i.e. where
$a(x_{0})=x_{0}\left(\frac{1-u}{1-s+sx_{0}}-1\right)\overset{!}{=}0$, which
means the initial point $x_{0}$ and the critical point $x_{c}$ of the
potential $\Psi(x)$ coincide, $x_{0}=x_{c}=1-u/s$. It is from this, that the
condition for the slow clicking regime can be deduced, namely $u<s$, which
implies the existence of a deterministic equilibrium.
We arrive finally at the following expression for the expected click time in
terms of our evolutionary parameters $N,s,u$

$\displaystyle\begin{gathered}\mathbb{E}(T_{\textrm{click}})\approx\frac{\pi
s^{2}\exp\left(2N\left(\frac{(s-u)^{2}}{-2s^{2}+4su-2u}+u\log\left(\frac{u}{s}\right)+(1-u)\log\frac{1-u}{1-s}\right)\right)}{2u\sqrt{\frac{s^{4}}{(s-1)^{2}u}-\frac{s^{2}}{u-1}}(s-u)}\times\\\
\left(\operatorname{erf}\left(\frac{(s-u)\sqrt{-\frac{N(s-1)^{2}}{s^{2}-2su+u}}}{s-1}\right)-\operatorname{erf}\left(\frac{(2s-1)(s-u)^{2}\sqrt{-\frac{N(s-1)^{2}}{s^{2}-2su+u}}}{(s-1)^{2}s}\right)\right)\times\\\
\left(\operatorname{erfi}\left(\frac{N\left(1-\frac{u}{s}\right)}{\sqrt{\frac{N(u-1)u}{s^{2}}}}\right)+\operatorname{erfi}\left(\frac{Nu}{s\sqrt{\frac{N(u-1)u}{s^{2}}}}\right)\right).\end{gathered}$
(21)

On a logarithmic scale, Eq. (21) yields in the limit of large populations

$\log\mathbb{E}(T_{\textrm{click}})\underset{N\to\infty}{\sim}2N\log\left[\left(\frac{u}{s}\right)^{u}\left(\frac{1-u}{1-s}\right)^{1-u}\right],$
(22)

in accordance with Eq. (10) from the main text.

Figure 4: Errors in predicted fixation probabilities to the classical fixation
problem over varying selection coefficients $s$ and population sizes
$N=100,1000$. Approximations by [39] (Eq. (28), in grey), [40] in the WF
diffusion scheme (Eq. (27), in blue), and from the new diffusion scheme (Eq.
(29), in orange) are compared to relative number of fixations in $10^{8}$
($N=100$) and $10^{7}$ ($N=1000$) WF simulations.

An analogous procedure, in which we replace $a(x)=sx(1-x)-ux$ from the WF
diffusion, Eq. 4, can be used to derive an estimate for Muller’s ratchet rate
from the WF diffusion

$\displaystyle\begin{gathered}\mathbb{E}(T_{\textrm{click}})\approx\frac{\pi
s}{2(s-u)u}\exp{\frac{-N(s-u)^{2}}{u}}\times\exp{\left(2N[s-u+u\log\Big{[}\frac{u}{s}\Big{]}]\right)}\times\\\
\quad\left[\mathrm{erf}\left(\sqrt{\frac{-N}{u}}\frac{(s-u)^{2}}{s}\right)-\mathrm{erf}\left(\sqrt{\frac{-N}{u}}(s-u)\right)\right]\times\quad\left[\mathrm{erfi}\left(\sqrt{-Nu}\left(\frac{s}{u}-1\right)\right)+\mathrm{erfi}\left(\sqrt{-Nu}\right)\right]\,.\end{gathered}$
(25)

For $\sigma<1-\frac{u}{s}$, where $\sigma^{2}=\frac{u}{2Ns^{2}}$ denotes the
variance of the Gaussian, the approximation can be reduced to

$\mathbb{E}(T_{\textrm{click}})\approx\sqrt{\frac{\pi}{Nu}}\frac{s}{(s-u)^{2}}\exp{\big{(}2N(s-u+\log\Big{[}\frac{u}{s}\Big{]})\big{)}},$
(26)

which is the same equation as derived by [26] through WKB approximations of
the associated Moran process (after a rescaling of time and population size
$t^{\prime}\rightarrow t/2N$ and $N^{\prime}\rightarrow 2N$). The result in
Eq. (26) coincides on a logarithmic scale with Eq. (11), as described in the
main text.

###  A.2 Fixation Probability and Fixation Times

Additionally, we here report the comparison of the estimated fixation
probabilities and fixation times for the Wright-Fisher model without mutations
$u=v=0$. In Fig. 4 we compare both diffusive approximations, Eq. (4) and Eq.
(7), as well as the formula reported in [39]. We notice that the corrections
in the strong selection regime are minor. Analytical formulae are compared to
mean fixation probabilities from a total of $10^{8}$ (for $N=100$) and
$10^{7}$ (for $N=1000$) Wright-Fisher simulations (Eq. (1)). The estimate from
the WF diffusion, Eq. (4) is the classical formula for the fixation problem as
presented by Crow and Kimura [40]

$\pi_{\textrm{CK}}=\frac{1-e^{-2s}}{1-e^{-2Ns}}.$ (27)

More recently, it was shown in [39] that

$\pi_{\textrm{SH}}=\frac{1-(1+s)^{-2}}{1-(1+s)^{-2N}}$ (28)

is slightly more accurate that the Crow-Kimura formula, among other good
properties.

Making use of the diffusive scheme by interpolation, Eq. (7), given an initial
point $x=x_{0}$, the probability of fixation at $x=1$ is approximated by

$\pi(x_{0})=\frac{(1+s)^{2N}(1+x_{0}s)^{-2N}((1+x_{0}s)^{2N}-1-x_{0}s)}{(1+s)^{2N}-1-s},$
(29)

which reduces to the classical fixation problem for $x_{0}=1/N$.

Using Eq. (13) with the appropriate boundary conditions
$\overline{t}(0)=\overline{t}(1)=0$ (two absorbing boundaries), the mean first
hitting time at $x=1$, conditional on fixation at $x=1$, with initial
condition $x=x_{0}$ can be computed through the following integral [31]

$\begin{split}\overline{t}(x_{0})=\frac{\pi^{0}(x_{0})}{\pi^{1}(x_{0})}\int_{0}^{x_{0}}\frac{2\pi^{1}(x)}{b(x)\Psi(x)}\int_{0}^{x}\Psi(y)\,\mathop{}\\!\mathrm{d}y\,\mathop{}\\!\mathrm{d}x\\\
+\int_{x_{0}}^{1}\frac{2\pi^{1}(x)}{b(x)\Psi(x)}\,\int_{x}^{1}\Psi(y)\,\mathop{}\\!\mathrm{d}y\,\mathop{}\\!\mathrm{d}x,\end{split}$
(30)

with $\Phi(x)=\int^{x}\frac{-2a(z)}{b(z)}\mathop{}\\!\mathrm{d}z$. Herein,
$\pi^{1}(x)$ and $\pi^{0}(x)$ denote the fixation probabilities at $x=1$ and
$x=0$, respectively, as a function of an initial point $x_{0}$. The relative
errors in estimates deriving from the differing $a(x)$ in the WF diffusion
(Eq. (4)) and the diffusive approximation by interpolation (Eq. (7)) are
compared in Fig. 5 relative to mean times of fixation from WF simulations for
varying parameters $s,N$.

Figure 5: Relative errors in approximations of fixation times for varying
selection coefficients $s$ and population sizes $N=100,1000$ as calculated
from Eq. (30), for the WF diffusion (in blue) and the new diffusion (in
orange). Errors are computed relative to simumlated mean fixation times over
$10^{8}$ ($N=100$) and $10^{7}$ ($N=1000$) runs of of WF simulations.

###  A.3 Mean establishment time with reversions

As presented in Fig. 2 of the main text, analytical estimates between the two
diffusive approximations for the Wright-Fisher model with selection and one-
way mutations $v=0$ are compared. These formulae can be derived equivalently
to the above described fixation times. The probability of establishment, given
an inital point $x_{0}$, is given by

$\pi^{x_{c}}(x_{0})=\frac{\int_{0}^{x_{0}}\Phi(y)\mathop{}\\!\mathrm{d}y}{\int_{0}^{x_{c}}\Phi(y)\mathop{}\\!\mathrm{d}y},$
(31)

where as before $\Phi(x)=\int^{x}\frac{-2a(z)}{b(z)}\mathop{}\\!\mathrm{d}z$
and $x_{c}$ denotes the metastable state defined by $p(x_{c})=x_{c}$.

From this the expected time to establishment , assuming the initial frequency
$x_{0}=\frac{1}{N}$, is computed by

$\begin{split}\overline{t}_{\text{est}}=\frac{2\left(1-\pi^{x_{c}}\left(\frac{1}{N}\right)\right)}{\pi^{x_{c}}\left(\frac{1}{N}\right)}\int_{0}^{\frac{1}{N}}\int_{0}^{x}\frac{\Phi(y)\pi^{x_{c}}(x)}{b(x)\Phi(x)}\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}x\\\
+\int_{\frac{1}{N}}^{x_{c}}\int_{x}^{x_{c}}\frac{2\Phi(y)\pi^{x_{c}}(x)}{b(x)\Phi(x)}\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}x\,.\end{split}$
(32)

###  A.4 WF simulations

Stochastic simulations were run on Wolfram Mathematica according to the
dynamics from Eq. (1), while varying the evolutionary parameters $N,s,u,v$ as
given in the main text. Due to limited computation times only a constrained
set of parameters are tested that do not exceed computational feasibility.

## References

  * Fisher [1923] R. A. Fisher, Proc. R. Soc. Edinb. 42, 321 (1923).
  * Wright [1931] S. Wright, Genetics 16, 97 (1931).
  * Kimura [1955] M. Kimura, Proc. Natl. Acad. Sci. U.S.A. 41, 144 (1955).
  * Kolmogoroff [1931] A. Kolmogoroff, Math. Ann. 104, 415 (1931).
  * Malécot [1945] G. Malécot, C. R. Acad. Sci. Paris 221, 340 (1945).
  * Feller [2015] W. Feller, in  _Selected Papers II_ (Springer International Publishing, 2015) pp. 121–140.
  * Ewens [2004] W. J. Ewens,  _Mathematical Population Genetics_ (Springer New York, 2004).
  * Rice [2004] S. Rice,  _Evolutionary theory: mathematical and conceptual foundations_ (Sinauer Associates, 2004).
  * Nagylaki [2013] T. Nagylaki,  _Introduction to theoretical population genetics_ , Vol. 21 (Springer Science & Business Media, 2013).
  * Crow [2017] J. F. Crow,  _An introduction to population genetics theory_ (Scientific Publishers, 2017).
  * Wakeley [2005] J. Wakeley, Genetics 169(1): 1-7 (2005).
  * Lumby  _et al._ [2020] C. K. Lumby, L. Zhao, J. Breuer, and C. J. Illingworth, eLife 9, e56915 (2020).
  * Ewens [1965] W. J. Ewens, Biometrics 21, 386 (1965).
  * Bürger [2000] R. Bürger,  _The Mathematical Theory of Selection, Recombination, and Mutation_ (John Wiley & Sons, 2000).
  * Muller [1932] H. J. Muller, Am. Nat. 66, 118 (1932).
  * Felsenstein [1974] J. Felsenstein, Genetics 78, 737–756 (1974).
  * Loewe [2006] L. Loewe, Genet. Res. 87, 133 (2006).
  * Chao [1990] L. Chao, Nature 348, 454 (1990).
  * Gordo and Charlesworth [2001] I. Gordo and B. Charlesworth, Genet. Res. 78, 149 (2001).
  * Govindaraju  _et al._ [2020] D. R. Govindaraju, H. Innan, and R. A. Veitia, Trends Genet. 36, 395 (2020).
  * López  _et al._ [2020] S. López, E. L. Lim, S. Horswell, K. Haase, A. Huebner, M. Dietzen, T. P. Mourikis, T. B. Watkins, A. Rowan, S. M. Dewhurst, _et al._ , Nat. Genet. 52, 283 (2020).
  * Hartfield and Keightley [2012] M. Hartfield and P. D. Keightley, Integr. Zool. 7, 192 (2012).
  * Stephan  _et al._ [1993] W. Stephan, L. Chao, and J. G. Smale, Genet. Res. 61, 225 (1993).
  * Waxman and Loewe [2010] D. Waxman and L. Loewe, J. Theor. Biol. 264, 1120 (2010).
  * Assaf and Mobilia [2011] M. Assaf and M. Mobilia, J. Theor. Biol. 275, 93 (2011).
  * Metzger and Eule [2013] J. J. Metzger and S. Eule, PLoS Comput. Biol. 9, e1003303 (2013).
  * Nei [2005] M. Nei, Mol. Biol. Evol. 22, 2318 (2005).
  * Manhart  _et al._ [2012] M. Manhart, A. Haldane, and A. V. Morozov, Theor. Popul. Biol. 82, 66 (2012).
  * Norman [1975] M. F. Norman, SIAM J. Appl. Math. 29, 225 (1975).
  * Ethier and Norman [1977] S. N. Ethier and M. F. Norman, Proc. Natl. Acad. Sci. U.S.A. 74, 5096 (1977).
  * Van Kampen [1992] N. G. Van Kampen,  _Stochastic processes in physics and chemistry_ , Vol. 1 (Elsevier, 1992).
  * Haigh [1978] J. Haigh, Theor. Popul. Biol. 14, 251 (1978).
  * Neher and Shraiman [2012] R. A. Neher and B. I. Shraiman, Genetics 191, 1283 (2012).
  * Jain [2008] K. Jain, Genetics 179, 2125 (2008).
  * Doering  _et al._ [2005] C. R. Doering, K. V. Sargsyan, and L. M. Sander, Multiscale Model. Simul. 3, 283 (2005).
  * Ovaskainen and Meerson [2010] O. Ovaskainen and B. Meerson, Trends Ecol. Evol. 25, 643 (2010).
  * Hanggi  _et al._ [1984] P. Hanggi, H. Grabert, P. Talkner, and H. Thomas, Phys. Rev. A 29, 371 (1984).
  * Kessler and Shnerb [2007] D. A. Kessler and N. M. Shnerb, J. Stat. Phys. 127, 861 (2007).
  * Sella and Hirsh [2005] G. Sella and A. E. Hirsh, Proc. Natl. Acad. Sci. U.S.A. 102, 9541 (2005).
  * Kimura [1962] M. Kimura, Genetics 47, 713 (1962).