The Cumulative Distribution Function Based Method for Random Drift Model

In the population genetics, the random genetic drift model describes that the number of gene variants (alleles) fluctuates randomly over time due to random sampling. The fraction of an allele in the population can be used to measure the intensity of the random genetic drift. In other words, the value of the fraction equals to zero or one, which means the allele disappears or is completely chosen in the system. Population genetics models, aiming at modeling genetic variability, had a natural start with discrete stochastic models at the individual level. The earlier mathematical model of random genetic drift, known as the Wright-Fisher model, was constructed by Ronald Fisher [13] and Sewall Wright [28, 29, 30]. The Wright-Fisher model is regarded as a discrete-time Markov chain under the assumption that the generations do not overlap and that each copy of the gene of the new generation is selected independently and randomly from the whole gene pool of the previous generation. Then Moran [21, 22] and Kimura [16, 17, 19, 20] derived the diffusion limit of random Wright-Fisher model. Chalub et. al spreaded the large population limit of the Moran process, and obtained a continuous model that may highlight the genetic-drift (neutral evolution) or natural selection [20, 4, 26].

We consider a population of constant size $N_{e}$ with a pair of two types A and B. Time $t$ is increased by the time step $\Delta t$, and the process is repeated. $x$ is the gene $A$ frequency at $t$ generation. Let $P_{N_{e},\Delta t}(t,x)$ denote the probability of finding a fraction $x$ of type $A$ individuals at time $t$ in population of fixed size $N_{e}$ evolving in discrete time, according to the Moran process. Then, in the limit of large population and small time steps, we postulate the existence of a probability density of states, that will depend on the precise way the limits are taken [4]. Namely:

In this paper, we study the following initial-boundary problem of $f(t,x)$:

where the function $M(x)$, which is typically a polynomial in $x$, incorporates the forces of migration, mutation, and selection, acting at time t[35, 20, 16].

Eq.(1.1) can be supplemented by the following conservation laws:

• (a)

  Mass conservation law:

  $$\int_{0}^{1}f(0,x)dt=\int_{0}^{1}f(t,x)dt=1.$$

  (1.2)

• (b)

  Moment conservation law:

  $$\frac{d}{dt}\int_{0}^{1}\theta(x)f(t,x)dx=0,$$

  (1.3)

  where $\theta(x)$ is the fixation probability function that satisfies

  $$x(1-x)\theta^{\prime\prime}+M(x)\theta^{\prime}=0,\ \theta(0)=0,\ \theta(1)=1.$$

  (1.4)

One notes that (1.3) recovers the conservation of expectation for pure drift ($M(x)=0$) with $\theta(x)=x$.

Based on different $M(x)$, the behavior of the solution has three cases:

The following theorem can be found in [5].

The appearance of the point measure $\delta_{0}$ ($\delta_{1}$) stands for that the fixation at gene $B$ ($A$) happens with a probability $a(t)$ ($b(t)$). The spatial temporal dynamics of the Kimura equation are well understood in the purely diffusive case and in only a relatively small number of population [18, 8, 7].

A well-known model is Muller’s ratchet, i.e, the fittest gene $B$ of individuals is eventually lost from the population and deleterious mutations ($B\rightarrow A$) slowly but irreversibly accumulate through rare stochastic fluctuations [24, 12]. In a finite asexual population, offspring inherit all the deleterious mutations their parents possess. Since these offspring also occasionally acquire new deleterious mutations, populations will tend to accumulate deleterious mutations over time. This effect is known as Muller’s ratchet.

There exists a unique steady solution $f_{\infty}$,

For the numerical simulation, the crucial features to solve (1.1) are that the numerical solution can keep the total probability conservation law and accurately capture the concentration phenomena at the discrete level. The total probability fails to keep the total probability by some classical numerical schemes [1, 27, 16]. A complete solution, whose total probability is unity, obtained by finite volume method (FVM) schemes in [35]. Xu et al. [33] compared a serial of finite volume and finite element schemes for the pure diffusion equation. Their critical comparison of the long-time asymptotic performance urges carefulness in choosing a numerical method for this type of problem, otherwise the main properties of the model might be not kept. In recent years, a variational particle method was proposed based on an energetic variational approach, by which a complete numerical solution can be obtained and the positivity of the solution can be kept [10]. However, some artificial criteria must be introduced to handle the concentrations near the boundary ends, even though it is designed based on the biological background. Recently, an optimal mass transport method based on pseudo-inverse of CDF is used to solve the model [3]. In this method, the feasible solution is strictly monotonous. However, for the cases of pure drift, natural selection and one-way mutation, the Dirac-$\delta$ concentrations must be developed, then the corresponding CDF must be discontinuous at the concentration points. This leads to a fact that its pseudo-inverse can not be strictly monotonous. However, numerical results were presented for the cases of pure drift and selection. So some manual intervention must be introduced in the numerical implements. The pure drift problem with multi-alleles, corresponding to a multidimensional PDE, is investigated in [34] by finite-difference methods, where the authors proposed a numerical scheme with absolute stability and several biologically/physically relevant quantities conserved, such as positivity, total probability, and expectation. So far, although quite a few numerical methods have been established for the pure drift and selection cases, efficient numerical works on the mutation case are not reported yet.

In this paper, we rewrite the system (1.1) to the following new one on the cumulative distribution function (CDF) $F(t,x)=\int_{0}^{x}f(t,s)ds$ as

Taking the CDF into account, the Dirac $\delta$ singularity at the boundary points of original PDF $f$ will change to the discontinuity of the CDF $F$ at the boundary points, and the fixation probability (lumped mass) will change to the height of the discontinuous jump. The singularity of negative power law will change to a bounded positive power law.

We will propose a upwind numerical scheme for (1.10) with a key revision near the boundary, which can handle the pure drift with or without selection and mutation, is unconditional stable, and keeps the total probability and positivity. It also keeps the conservation of expectation for pure drift. The numerical results show that the scheme can catch the height of discontinuity at the ends and predict accurately the fixation probability for the cases of pure drift, natural selection and one-way mutation. It also predict accurately the negative power of the power law for two-way mutation case.

The rest of this paper is organized as follows. In Section 2, we construct the numerical scheme for (1.10). Some numerical analysis is presented in Section 3. In Section 4, several numerical examples are presented to validate the theoretical results and to demonstrate the ability to trace the long-time dynamics of random genetic drift. Some discussions will be presented in Section 5 about the relations between the revised scheme and the standard finite difference method.

In this section, we introduce the revised finite difference method (rFDM) for (1.10) with the central difference method for diffusion term and the upwind scheme for convection term. Let $h=1/K$ be the spatial step size, and $x_{i}=ih$, $i=0,\cdots,{K}$, be the spatial grid points. Let $\tau$ be the temporal step size, and $t_{n}=n\tau$, $n=1,2,\cdots,N$, be the temporal grid points.

Define the first order difference as

and denote the diffusion coefficient by $a(x)=x(1-x)$. The upwind numerical scheme, referred as rFDM, reads as: Given $F^{n-1}$, $F^{n}$ solves the following linear algebra system,

with $F^{n}_{0}=0$ and $F^{n}_{K}=1$, for $n=1,\cdots,N$, and the key revision on the diffusion cofficient

Finally, the solution of the original equation (1.1) can be recovered from $F$ by, for $n=1,\cdots,N$,

In the above formula, central difference does not applied at near boundary points $i=1$ or $K-1$, thanks to the fact that the discontinuous jump may occur at the boundary points.

In this section, we will focus on the numerical analysis for rFDM (2.2), including the unconditional stability, positivity preserving, and conservation law of the total probability and expectation.

Let ${\bf B}=(b_{ij})$ be the matrix of the linear system. Then ${\bf B}$ is a tri-diagonal matrix.

Note that

and

Then ${\bf B}$ is a M-matrix, i.e., any entry of the inverse matrix ${\bf B}^{-1}$ is positive. So rFDM (2.2) is unconditionally stable and positivity preserving.

By the integral by parts, the expectation $\mathcal{E}(t)$ satisfies that

So we define a discrete expectation as

For pure drift case, $M(x)\equiv 0$, then we have, by (2.2) and (2.3), that

So the discrete system (2.2) keeps the conservation of the discrete expectation for pure drift case. Please note that, without the revised treatment (2.3), the conservation of the discrete expectation is invalid.

In this section, we will show the effectiveness of this algorithm by different numerical examples. In Example 1, we verify the local convergence. In Section 4.1, the pure drift $M=0$ and natural selection $M(x)=x(1-x)(\eta x+\beta)$ case are studied in Example 2 and Example 3, respectively. In Section 4.2, we consider the mutation case $M(x)=\gamma(1-x)-\mu x$, including one-way mutation model, such as the Muller’s ratchet model, in Example 4 and two-way mutation model in Example 5.

Example 1. Local convergence

Although the singularity may develop at boundary points, the solution in interior area is sufficiently smooth. To verify the correctness of the numerical scheme, we check the local convergence.

Let $\Omega^{\prime}\subset\Omega$ be the interior area and

Define the error $e$ of $F(t,x)$ in $\mathcal{L}^{2}(\Omega^{\prime})$ and $\mathcal{L}^{\infty}(\Omega^{\prime})$ mode as

and

where $\{F_{i}^{n}\}_{i=0}^{K}$ is the numerical solution of CDF model (1.10), and $\{F_{e,i}^{n}\}_{i=0}^{K}$ is the corresponding exact solution at time $t^{n}$, $n=1,\cdots,N$.

In this example, we take the initial probability density as

Table 1 shows the error and local convergence order of $F(t,x)$ for $M(x)=x(1-x)(-4x+2)$ and $M(x)=0.2(1-x)+0.4x$ in $\mathcal{L}^{2}([0.3,0.7])$ and $\mathcal{L}^{\infty}([0.3,0.7])$ mode with different space and time grid sizes at time $t=0.1$. The reference ”exact” solution is obtained numerically on a fine mesh with $h=1/100000$ and $\tau=1/100000$. The results show that the local convergence of the numerical scheme is 2nd order in space and 1st order in time for different $M(x)$ in the inner region $[0.3,0.7]$.