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 with or without natural selection and mutation. For pure drift and natural selection case, the Dirac δ 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 δ 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 δ singularity of the PDF becomes the discontinuity of the CDF. Then we establish a revised finite difference method(rFDM), 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. The numerical results show the effectiveness of the scheme. We also compare rFDM with a standard finite difference method(sFDM). We find that, for pure drift problem, sFDM fails to keeps the conservation of expectation and can not predict the fixation probability, because there is artificial mutation introduced near two boundary ends.