Wright-Fisher Simplicial Diffusion

Updated 19 December 2025

Wright-Fisher simplicial diffusion is a stochastic process on the probability simplex that generalizes classical diffusion to multi-allele settings.
It employs geometric embeddings, spectral expansions, and exact simulation methods to capture the combined effects of drift, mutation, and selection.
This framework underpins applications in population genetics and machine learning by providing insights into fixation probabilities and evolutionary dynamics.

Wright-Fisher simplicial diffusion is a class of continuous-time, continuous-state stochastic processes modeling the joint evolution of allele frequencies in a population at a single genetic locus, where the state space is the probability simplex. The process generalizes classical one-dimensional Wright-Fisher diffusion to the multi-allele, simplex-valued setting, capturing the combined effects of random genetic drift, mutation, and (optionally) selection in the diffusion limit of large population size. The simplex geometry, degenerate diffusion at the boundary, and links with spherical Brownian motion and orthogonal polynomials make Wright-Fisher simplicial diffusion a central object in mathematical population genetics, probability, and, more recently, machine learning.

1. Mathematical Formulation and Geometric Interpretation

Let $k$ denote the number of alleles. The state of the system at time $t$ is $x(t) = (x_1(t), \ldots, x_k(t))$ , with $x_i(t) \geq 0$ , $\sum_{i=1}^k x_i(t) = 1$ , i.e., $x(t) \in \Delta^{k-1}$ , the standard $(k-1)$ -simplex. The most general neutral Wright-Fisher diffusion in Itô form, with parent-independent mutation parameters, is given by

$dx_i(t) = \frac12\bigl(\varepsilon_i - \mu\,x_i(t)\bigr)\,dt + \sum_{j=1}^k \sqrt{x_i(t)\,x_j(t)}\,db_{ij}(t), \qquad i=1,\dots,k,$

where $b_{ij}(t)$ are standard Brownian motions with $b_{ji} = -b_{ij}$ , $\varepsilon_i > 0$ the per-allele mutation parameters, and $\mu = \sum_i \varepsilon_i$ . The diffusion matrix is degenerate at the boundary; $x_i(t)=0$ is an exit or reflecting boundary depending on $\varepsilon_i$ . The Fokker-Planck (forward Kolmogorov) equation is

$\frac{\partial p}{\partial t} = -\sum_{i=1}^k \frac{\partial}{\partial x_i}\Bigl[\,\tfrac12 (\varepsilon_i-\mu x_i) p\Bigr] + \frac12 \sum_{i,j=1}^k \frac{\partial^2}{\partial x_i\partial x_j} \Bigl[\,x_i (\delta_{ij}-x_j) p\,\Bigr].$

A geometric characterization is available: isotropic Brownian motion on the unit $(k-1)$ -sphere $S^{k-1}$ , projected via $x_i = y_i^2$ ( $y \in S^{k-1}$ ), yields Wright-Fisher simplicial diffusion with $\varepsilon_i = \tfrac12$ (Maruyama et al., 2015). This mapping underlies the connection between hyperspherical and simplicial diffusion and enables efficient simulation algorithms.

2. Spectral Theory and Eigenfunction Expansions

The spectral properties of Wright-Fisher simplicial diffusion are central to both theoretical and computational approaches. The infinitesimal generator is self-adjoint (in the reversible case) with respect to the Dirichlet stationary density

$\pi_{\mathrm{stat}}(x) = \frac{\Gamma(\mu)}{\prod_i \Gamma(\varepsilon_i)} \prod_{i=1}^k x_i^{\varepsilon_i-1}.$

The spectrum consists of (multivariate) polynomial eigenfunctions:

In the neutral, no-mutation case ( $\varepsilon_i=0$ ), the eigenfunctions are homogeneous polynomials $x^\alpha$ related to symmetric multivariate Jacobi polynomials, orthogonal with respect to the weight $w(x) = \prod_i x_i$ (Tran et al., 2012, Maruyama et al., 2015).
With positive $\varepsilon_i$ , total-degree- $n$ polynomials $Q_n(x,x')$ form reproducing kernels, as in the Griffiths expansion: $p(x,t|x',0) = \pi_{\mathrm{stat}}(x) \sum_{n=0}^\infty e^{-\lambda_n t} Q_n(x,x'),$ with eigenvalues $\lambda_n = [n(n-1) + \mu n]/2$ (Maruyama et al., 2015). Under the $x_i = y_i^2$ mapping, this expansion is termwise equivalent to a Gegenbauer polynomial expansion of the transition kernel for Brownian motion on $S^{k-1}$ .

The spectral expansions facilitate explicit calculation of transition densities, moments, fixation probabilities, and exact simulation strategies (Jenkins et al., 2015, Sant et al., 2023). In the $k=2$ (biallelic) case, the eigenfunctions reduce to Jacobi polynomials and the transition density to a mixture of Beta distributions (Tran et al., 2012, Jenkins et al., 2015).

3. Boundary Structure, Hierarchical Extensions, and Loss of Alleles

Boundary stratification is crucial for both the analytic and probabilistic understanding of Wright-Fisher simplicial diffusion. The diffusion operator is degenerate at the boundary $x_i=0$ , and classical boundary conditions (Dirichlet/Neumann) are typically ill-defined.

The hierarchical extension scheme constructs a global solution to the forward Fokker-Planck equation on the closed simplex $\overline{\Delta^{k-1}}$ by recursively gluing solutions on lower-dimensional faces (Hofrichter et al., 2014). When an allele frequency vanishes (allele loss), probability mass exits the interior and seeds a Wright-Fisher process on the corresponding face with one fewer allele. This process continues recursively, so that the global solution keeps track of all surviving alleles at each stratum, and ensures the correct moment evolution and total-mass conservation.

For the backward Kolmogorov equation, an analogous recursive interpolation stratifies the solution along faces, guaranteeing $C^2$ regularity up to codimension-one boundaries and uniqueness (Hofrichter et al., 2014). The spectral expansion (Green’s function) is accordingly hierarchically glued across the simplex strata.

This mechanism removes boundary singularities and accurately models natural phenomena such as fixation or extinction of alleles, consistent with the probabilistic structure of genetic drift and absorption (Hofrichter et al., 2014, Hofrichter et al., 2014).

4. Simulation and Numerics: Exact Path and Bridge Sampling

Exact sampling of Wright-Fisher simplicial diffusions is obviated by their spectral expansions. In the reversible, neutral, parent-independent mutation case, the transition density can be written as a mixture (in $k = 2$ as Beta, in $k > 2$ as Dirichlet mixtures) whose weights are determined by Kingman coalescent block-counting or similar series (Jenkins et al., 2015, Sant et al., 2023). The core algorithm uses efficient “alternating series” inversion to sample the latent block-count variable and then completes with standard Beta or Dirichlet sampling.

For diffusion bridges (paths conditioned on endpoints), the transition density becomes a higher-order mixture (over multiple latent variables), but is still amenable to exact sampling (Sant et al., 2023). Non-neutral processes, including selection, are handled by Girsanov-type change of measure and Poisson thinning (retrospective acceptance) to adjust the neutral proposals (Jenkins et al., 2015, Sant et al., 2023).

This simulation paradigm enables unbiased Monte Carlo estimation even in regimes inaccessible to direct SDE simulation, such as paths conditioned to avoid absorption or under frequency-dependent selection.

5. Path Integral Formulations and Perturbative Expansions

The Wright-Fisher simplex diffusion and its discrete-time predecessor admit path-integral (sum-over-paths) representations (Waxman, 2024, Schraiber, 2013). In the discrete model, the full transition probability from $\mathbf{a}$ to $\mathbf{z}$ in $t$ generations is a sum over all possible intermediate frequency sequences, weighted by multinomial sampling (mutation and selection via exponential weights). The diffusion limit leads to a functional integral over continuous trajectories $x:[0,t]\to\Delta^{k-1}$ , with an Onsager–Machlup action comprised of drift, diffusion, and selection costs.

For weak selection, the path integral admits a perturbative (Dyson/Feynman-diagram) expansion in the scaled selection coefficient (Schraiber, 2013). In this expansion, each “vertex” corresponds to a selective scattering event, and the neutral propagators are known. The series converges absolutely for all finite $t$ , and error bounds are explicit.

These path-integral and diagrammatic approaches are beneficial for theoretical analysis, as they provide tractable expressions for transition probabilities and suggest approximation schemes beyond standard PDE treatments.

6. Extensions, Generalizations, and Statistical Implications

Wrigh-Fisher simplicial diffusion encompasses a suite of extensions:

Negative mutation rates and time-reversal: Generalizations to processes with negative “mutation rates” describe non-ergodic, transient dynamics, with explicit exit distributions calculable via skew-products of Bessel-square processes (Pal, 2010).
Irreversible mutation and fluxes: When the mutation matrix is non-reversible, stationary densities concentrate on simplex edges, with asymmetry in the fluxes between alleles directly encodable (Burden et al., 2016).
Matrix-valued and bivariate extensions: Extra discrete components or phases introduce matrix-valued generators, which admit matrix-valued spectral decompositions and associated recurrence/transience theory (Iglesia, 2011).

Statistical inference is sensitive to the boundary classification of the diffusion: The mutual absolute continuity or singularity of path measures under different parameter values is controlled by whether the process can reach (or leave) the simplex boundary, strongly impacting maximum likelihood or Bayesian estimators (Jenkins, 2024). Separating times reflect when two models first diverge in likelihood, directly linked to moments of allele loss or fixation.

7. Applications: Modern Sequence Modeling and Diffusion Generative Models

Recent advances exploit Wright-Fisher simplicial diffusion as a mathematically principled foundation for generative modeling over discrete or simplex-valued data (e.g., DNA, protein, and language sequences) (Chandra et al., 17 Dec 2025). Wright-Fisher diffusion naturally interpolates between discrete, Gaussian, and simplex-valued settings, depending on population size and mutation parameter regimes. This connection enables unified loss functions and likelihoods across different data domains, with computational benefits in terms of both speed and numerical stability.

Practical implementations leverage the exact Dirichlet mixture sampling, sufficient-statistic parameterizations, and efficient score-matching algorithms to achieve improved stability and likelihoods for conditional sequence generation tasks, outperforming previous simplex-diffusion baselines in chromatin-accessibility–conditioned DNA design.

Key references: (Maruyama et al., 2015, Hofrichter et al., 2014, Hofrichter et al., 2014, Tran et al., 2012, Jenkins et al., 2015, Burden et al., 2016, Sant et al., 2023, Waxman, 2024, Pal, 2010, Iglesia, 2011, Chandra et al., 17 Dec 2025, Jenkins, 2024, Schraiber, 2013).