The Cumulative Distribution Function Based Method for Random Drift Model

In this paper, we propose a numerical method to uniformly handle the random genetic drift model for pure drift with or without natural selection and mutation. For pure drift and natural selection case, the Dirac $\delta$ singularity will develop at two boundary ends and the mass lumped at the two ends stands for the fixation probability. For the one-way mutation case, known as Muller's ratchet, the accumulation of deleterious mutations leads to the loss of the fittest gene, the Dirac $\delta$ singularity will spike only at one boundary end, which stands for the fixation of the deleterious gene and loss of the fittest one. For two-way mutation case, the singularity with negative power law may emerge near boundary points. We first rewrite the original model on the probability density function (PDF) to one with respect to the cumulative distribution function (CDF). Dirac $\delta$ singularity of the PDF becomes the discontinuity of the CDF. Then we establish a upwind scheme, which keeps the total probability, is positivity preserving and unconditionally stable. For pure drift, the scheme also keeps the conservation of expectation. It can catch the discontinuous jump of the CDF, then predicts accurately the fixation probability for pure drift with or without natural selection and one-way mutation. For two-way mutation case, it can catch the power law of the singularity. %Moreover, some artificial algorithms or additional boundary criteria is not needed in the numerical simulation. The numerical results show the effectiveness of the scheme.