Numerical method for multi-alleles genetic drift problem

Genetic drift describes random fluctuations in the number of genes variants in a population. One of the most popular models is the Wright-Fisher model. The diffusion limit of this model is a degenerate diffusion-convection equation. Due to the degeneration and convection, Dirac singularities will always develop at the boundaries as time evolves, i.e., the fixation phenomenon occurs. Theoretical analysis has proven that the weak solution of this equation, regarded as measure, conserves total probability and expectations. In the current work, we propose a scheme for 3-alleles model with absolute stability and generalize it to N-alleles case (N > 3). Our method can conserve not only total probability and expectations, but also positivity. We also prove that the discrete solution converges to a measure as the mesh size tends to zero, which is the exact measure solution of the original problem. The simulations illustrate that the probability density decays to zero first on the inner nodes, then also on the edge nodes except at the three vertex nodes, on which the density finally concentrates. The results correctly predict the fixation probability and are consistent with theoretical ones and with direct Monte Carlo simulations.