The Cumulative Distribution Function-Based Method for Random Drift Model

We propose a numerical method to uniformly handle the random genetic drift model for pure drift, natural selection and mutation. For pure drift and natural selection cases, 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 the 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 revised finite difference method (rFDM), which can keep the total probability, is positivity-preserving and is unconditionally stable. The advantage of our algorithm is that it is straightforward and preserves the moment for pure drift and natural selection. For all cases involved, it can be proven that after a long time, the numerical solution converges to the numerical steady state. This is another advantage, and it means that the numerical scheme can catch the discontinuous jump of the CDF, and accurately predict the fixation probability for pure drift, natural selection and one-way mutation. Meanwhile, for the two-way mutation case, it can also catch the power law of the singularity accurately. The numerical results show the effectiveness of the scheme. Finally, by comparing rFDM with a standard finite difference method (sFDM), we find that for the pure drift problem, sFDM fails to keep the conservation of expectation and cannot predict the fixation probability accurately, because an artificial mutation is introduced near two boundary ends. For the two-way mutation case, sFDM fails to predict the correct power law of the singularity.