An Optimal Mass Transport Method for Random Genetic Drift

We propose and analyze an optimal mass transport method for a random genetic drift problem driven by a Moran process under weak-selection. The continuum limit, formulated as a reaction-advection-diffusion equation known as the Kimura equation, inherits degenerate diffusion from the discrete stochastic process that conveys to the blow-up into Dirac-delta singularities hence brings great challenges to both the analytical and numerical studies. The proposed numerical method can quantitatively capture to the fullest possible extent the development of Dirac-delta singularities for genetic segregation on one hand, and preserves several sets of biologically relevant and computationally favored properties of the random genetic drift on the other. Moreover, the numerical scheme exponentially converges to the unique numerical stationary state in time at a rate independent of the mesh size up to a mesh error. Numerical evidence is given to illustrate and support these properties, and to demonstrate the spatio-temporal dynamics of random generic drift.

1. Introduction. In population genetics, genetic drift describes random fluctuations in the numbers of gene variants (alleles) over time. Allele frequency, expressed as a percentage, measures the relative fraction of an allele at a particular locus in the population, and its change quantifies the intensity of the random genetic drift [17]. Typically, when genetic drift begins, it will continue until either i) the involved allele completely disappears from the population or ii) the allele establishes permanently at 100% frequency (called fixed). In either case, genetic drift causes gene variants to disappear because infrequently occurring alleles face a greater chance of being lost in a small population, or causes a new population genetically distinct from its original population such that initially rare alleles become much more frequent and even fixed [17,32,33,37]. Both events indicate that genetic drift can decrease the population's genetic diversity, and it plays a role in the evolution of new species.
Mathematical modeling of genetic drift dates back to the pioneering works of Ronald Fisher [15] and Sewall Wright [39,40,41]. The Wright-Fisher model employs a discrete stochastic process to model dynamics of finite populations at the individual level under the restrictions 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. It is later modified and extended by Patrick Moran [28,27] (allowing generations overlap) and Komoo Kimura [24,21,22,23] (allowing the genetic mutation to spread across the population). In particular, they show that, in the limit of a large population and weak selection, these processes can be approximated by the same diffusion approximation, namely the Kimura equation [22] which describes the probability of fixation of a mutant with frequency-independent fitness. The continuum framework makes a systematic qualitative and quantitative analysis of the new model possible thanks to the tools from modern analytical and numerical analysis.
In this paper, we consider the Moran process as a paradigm and introduce its large population limits with different drift-diffusion scalings assumption [11]. Consider the dynamics of a population with N individuals that are distinguished by two neutral alleles labeled A and B such that they do not affect the survival and reproduction ability of the individual. In light of the balance between selection and drift, Traulsen et al. [35] summarized the process into three simple steps: (a) selection-an individual is randomly selected for reproduction with a probability proportional to its fitness; (b) reproduction-the selected individual produces one identical offspring; (c) replacement-the offspring replaces a randomly selected individual in the population. The process is then repeated after each time step ∆t. The fitness or reproduction rate of an individual depends on environmental conditions as well as the type and frequency of its competitors. For each type, we associate a fitness function depending on the type frequency: In particular, when Ψ is constant, i.e., the fitness of an individual is genetically determined and not affected by interactions, one recovers the classical frequency-independent Moran process.
Under the assumption N −1 ∝ (∆t) µ and a suitable time rescaling, as N → ∞, one collects from [11,23,35] the following thermodynamical limit for the density f (x, t) of allele A when ν = µ = 1 as two special cases such that genetic drift is the only evolutionary force in (1.1a), and the evolutionary force includes genetic mutation, migration and selection in (1.1b). They correspond to the limits of discrete process with ν > µ = 1 2 and µ ∈ ( 1 2 , 1], ν = 1 − µ, respectively. The biological significance urges us to impose the non-flux boundary condition to each of (1.1) such that the following conservation holds where θ(x) is the fixation probability function that satisfies and it can be explicitly given by One notes that (1.4) recovers the conservation of total population for (1.1) with θ(x) ≡ 1, and conservation of mass center for (1.1a) with θ(x) = x.
In this work, we will take advantage of the fact that the following free energy admits (1.1) as a gradient flow with respect to a variation of optimal transport distances (see e.g. [36,6,4,10]) and that (1.1) can be rewritten as denotes the first variation of the energy E with a fixed mass constraint. We now recall several theoretically relevant results on the well-posedness of (1.1)-(1.2) by [12,11,34]. Let us denote by BM + ([0, 1]) the space of functions with positive Radon measure in [0,1].
According to Theorem 1.1, we always expect (1.1) to collapse into the Dirac-delta singularity regardless of the initial distribution. The long-time dynamics along with this spiky spatial profile well demonstrate and capture the formation of gene segregation, i.e., an asymptotic gene fixation with allele A (resp. allele B) occupying the whole population when L 0 f 0 (x)θ(x)dx = 0 (resp. 1).
More recently, the gradient flow structure (1.7) has been discussed in [10]. The analysis of the model is built upon the classical steepest descent variational schemes for nonlinear Fokker-Planck equations which are introduced in [20] and generalized in [1,2]. The transport distance between probability measures has to be adapted to the degenerate diffusion coefficient of (1.1). They introduced a suitable distance called the Shahshahani distance showing the convergence of the variational scheme for the specific case of replicator dynamics, i.e. (1.1) without the degenerate diffusion term. It is an open problem to show the convergence of this variational scheme including the degenerate diffusion term.
The spatial-temporal dynamics of the Kimura equation are well understood in the purely diffusive case. However, the singularity induced by the convergence towards Dirac-delta concentrations, which well model the biologically realistic gene segregation or fixation through the drift, imposes substantial challenges to their numerical approximation, in particular for long-time computations. Therefore, one of the crucial features of any numerical scheme to solve (1.7), while keeping its biological relevance, is to accurately capture the concentration phenomena at the discrete level. Lagrangian numerical schemes for these gradient flow interpretations of one-dimensional Fokker-Planck, thin-film and quantum drift equations have been devised by various authors, see e.g. [19,18,2,38,9,8,25,31,26,30]. We here build upon the strategy of solving an equivalent equation satisfied by the diffeomorphisms mapping the initial data to the solution at later times as performed in [2,7,8,3], see next section for details. We refer to [5] for a very recent review on the state of the art of this kind of numerical schemes. The main advantage of the present approach is being able to deal with Dirac-delta concentrations easily and rigorously now that we work in mass variables. This merit is also utilized by [2,8,3] to study blow-ups in Keller-Segel models for chemotactic cellular aggregation.
Other numerical schemes have recently been proposed in the literature to solve the Kimura equation (1.1). In [13] the authors introduced a numerical method based on Lagrangian maps that preserves the free energy decay of the system. They also analyze a convex-concave splitting approach which leads to an implicit method, and prove the unique solvability of this method. However, the equation is solved in the original variables therein, so they have to devise heuristic criteria to capture the concentration of mass towards the endpoints. Moreover, they do not take full advantage of the gradient flow structure (1.7) and the variational schemes at the discrete level. As an alternative, our approach is to construct numerical schemes directly on the optimal maps and can describe the Dirac-delta to the fullest extent. See details in Section 2. [42] performs an interesting "horse race" comparison of a serial of finite volume and finite element schemes for (1.1a). Their critical comparison of the long-time asymptotic performance urges carefulness in choosing a numerical method for this type of problem, especially when the main properties of the model are not kept by the scheme. We also want to mention that the genetic drift problem with multi-alleles, cast as a multi-dimensional PDE, is investigated in [43] through finite-difference methods, where the authors propose a numerical scheme with absolute stability and conserve several biologically/physically relevant quantities such as positivity, total probability, and conservation.
The rest of this paper is organized as follows. Section 2 introduces the evolution problem for the map in mass variables as in [2,8,5] together with the needed prerequisites in optimal transport theory. Since the free-energy functional shows again in the new variables a non-convex structure, we present our numerical method and then apply it for the genetic-drift problem by introducing the convex-splitting technique as the preprocessing step of the Euler implicit scheme. Section 3 is devoted to the analysis of the properties of the numerical scheme. We show that the discrete problem converges exponentially fast to a unique stationary state with a monotone drift V for any fixed discretization parameter modulo error terms in spatial discretization, and we characterize all possible stationary states and limit points otherwise. See Theorem 3.9 and Theorem 3.12 for our main results. Finally, Section 4 conducts several sets of numerical experiments to validate the shown properties and to accurately illustrate the long-time dynamics of random genetic drift.
2. Numerical Methods. In this section, we recast the random generic drift models into an evolution problem for map of the mass variable, and then propose a numerical scheme to solve the new equation. Since the associated new free-energy functional is non-convex, we introduce the convex-splitting technique as the preprocessing step of the Euler implicit scheme in its numerical solver.
2.1. Optimal transport and the Wasserstein distance. We first introduce a suitable Wasserstein distance in the probability space P([0, 1]) such that equation (1.1) can be interpreted as a gradient flow of the free energy (1.6). Due to the presence of the variable coefficient x(1 − x) in (1.1), the quadratic Wasserstein distance will not be based on the usual Euclidean distance but on the induced generalized Shahshahani distance dt, for x, y ∈ Ω = (0, 1). Chalub et al. [10,Lemma 10] prove that the infimum in (2.1) is achieved at a unique constant-speed geodesic, and d can be uniformly extended toΩ ×Ω as moreover, d defines a distance inΩ and the metric space (Ω, d) is Polish. This Shahshahani distance is locally equivalent to the Euclidean one in the interior, but behaves differently close to the boundary. This difference is reflected by the dynamics of replicator-diffusion equation (1.1), which is locally uniformly parabolic in the interior, but degenerate at the boundaries. We now present several concepts from transport theory essential for this paper. Let µ and ν be two absolutely continuous measures with respect to the Lebesgue in P([0, 1]), and T be a measurable map from [0, 1] → [0, 1]. We say that T transports µ onto ν and write ν = T µ if ν(B) = µ • T −1 (B) for any measurable set B ⊂ [0, 1]. We also say ν is the push-forward or the image measure of µ by T . To introduce the corresponding Wasserstein distance between µ and ν, one defines as soon as the source measure µ has no atoms. In fact, by Brenier's theorem, if µ is absolutely continuous with respect to the Lebesgue measure, then there exists a measurable nondecreasing map T such that ν = T µ. The proper definition of the associated Wasserstein distance needs a relaxed variational scheme given by where Π runs over the set of transference plans Γ between µ and ν. One important simplification in 1D is that the optimal plan can be characterized fully in terms of the inverse of cumulative distribution functions, see [36]. To be specific, let F and G be the cumulative of the 1D functions f and g, and define the pseudo-inverse By Brenier's theorem and the definition of the image measure, we have Then it is straightforward to obtain ϕ = G −1 • F , and the Wasserstein distance in (2.2) becomes The proof of this fact without relying to Brenier's theorem can be found in [36, Section 2.2]. By the definition of pseudo-inverse function, we can derive the following evolution equation satisfied by Φ(η, t) We would like to remark that, the Dirichlet boundary condition of Φ(1, t) = 1 applies here since one have from the strong maximum principle that f (x, t) in strictly positive in [0, 1] for all time t if f 0 (x) ≥, ≡ 0. This, on the other hand, indicates that Φ(η, t) must be strictly increasing in η, which is an important property to preserve for numerical schemes. Actually, according to the non-flux boundary condition of f (x, t), we have This implies the Dirichlet boundary condition holds Φ(0, t) = 0 and Φ(1, t) = 1 for all time t. However, a free boundary condition should be adopted if f (x, t) remains compactly supported when studying problems with diffusion degenerate inside the domain. Note that (2.4) has the following free energy Then we will connect this evolution problem for the map Φ(η, t) pushing forward the initial data f 0 to the solution f (x, t) at time t with the continuum limit of implicit Euler steps obtained as Euler-Lagrange conditions for suitable variational problems. Finally, let us restate Theorem 1.1 in terms of the map Φ(η, t).
In the sequel, we will design a numerical scheme capable of accurately capturing the long-time behavior described in Theorem 2.1.

Discretization for Euler-Lagrange of Replicator-Diffusion (1.1).
We now consider the spatio-temporal discretization for the Euler-Lagrange problem (2.4) of the full Replicator-Diffusion equation (1.1). Throughout this paper, we assume that the time step τ and space step h are constant in the discretization. In terms of the Wasserstein distance (2.3), the Jordan-Kinderlehrer-Otto (JKO) steepest descent scheme implies that finding the inverse distribution function in (2.4) corresponds to solve the following for a fixed time step τ > 0 In light of (2.2), we find that Hence an approximated Euler-Lagrange equation (2.4) associated to this minimization problem is If we denote Φ k i = Φ(ih, kτ ) for i = 0, 1, · · · , N and N h = 1, k ∈ N, our full finite difference discretization of (2.6) is the following implicit scheme: with the Dirichlet boundary condition Φ k+1 The solution at each time step is computed by an iterative Newton's procedure with initial diffeomorphism obtained by a preprocessing step, and one does not need a CFL condition for this implicit in time discretization.

Preprocessing
Step via a Convex Splitting Technique. System (2.4) can be viewed as a gradient flow associated with the energy functional (2.5). This energy has a non-convex structure, and it prohibits the direct application of the proposed implicit scheme due to the singularity of the numerical scheme on the boundary and local convergence of Newton's method. However, the merits brought by the convex-splitting technique enable us to construct a numerical scheme in the preprocessing step and then calculate the initial diffeomorphism for Newton's method in the Euler implicit scheme.
To this end, we write V = V c − V e for V c and V e being smooth convex functions and then apply the convex splitting method of Eyre [14] to obtain the following semi-discrete scheme Here, the map from Φ k to Φ k+1 is an optimal transport in the sense that Φ k+1 minimizes the functional where W is a convex functional explicitly given by Let us introduce the discrete space domain and its closureQ := Q ∂Q with boundary The full finite difference discretization of (2.8) is formulated as follows: with the boundary condition Φ k 0 = 0 and Φ k N = 1 for each k. At each time step, we solve a system of nonlinear equations by a damped Newton's iteration. For notational simplicity, let us denote the nonlinear functional F in (2.11) as then we will calculate the Jacobian matrix DF of (2.12) and determine the Newton update γ k+1,n+1 by where index k corresponds to the temporal discretization and n to the Newton iteration.
3. Numerical Analysis. This section analyzes the implicit Euler scheme (2.7) and the convex splitting scheme (2.11). We first collect several important properties of the splitting technique applied in [13] that includes the unique solvability, convergence of the Damped Newton method (2.13) and dissipation of energy. Then we prove that there exists a lower bound for the discrete energy, which implies the existence of numerical solutions. Finally, in light of the fact that any steady state must be a Heaviside-type step function, we show the convergence of of the implicit scheme (2.7).

Several Properties of the Convex Splitting Scheme.
Lemma 3.1. The numerical scheme (2.11) is uniquely solvable in space Q given by (2.10). Proof. We first introduce the discretization of (2.9) as is a convex function on the closed convex setQ and J N (y) = +∞ on boundary ∂Q, it is straightforward to obtain that there exists a unique minimizer x ∈ Q.
To show the unique solvability of the scheme, it suffices to prove that x ∈ Q is the minimizer of J N (y) if and only if it is a solution of (2.11). Suppose that x ∈ Q minimizes J N (y). Define j 1 ( ) := J N (x + (y − x)). Then there exists 0 > 0 small enough such that x + (y − x) ∈ Q for any ( , y) ∈ (− 0 , 0 ) ×Q. Now that j 1 ( ) achieves its minimizer at zero, one has j 1 (0) = 0 hence h Therefore, x is a solution of (2.11).
To prove the "only if" part, we assume that x ∈ Q solves the scheme (2.11) for any y ∈ Q. Then where the last inequality holds since m − ln m − 1 > 0 for any positive m, and V c is convex function with positive second derivative. This completes the proof.
We now introduce the following concept to prove the convergence of the Newton's method of the scheme (2.11).
Definition 3.2. ( [29]) Let G be a finite-dimensional real vector space and Q be an open nonempty convex subset of G. Then a convex function Λ ∈ C 3 : Q → R is called self-concordant on Q if there exists a constant a 0 > 0 such that the following inequality holds for all x ∈ Q and all u ∈ G: where (D k Λ(x)[u 1 , · · · , u k ] is its k-th differential taken at x alone the collection of direction (u 1 , · · · , u k ).
Then we have the following theorem.
Since both the linear and quadratic functions have zero third-order derivatives, one can easily find that J 1 N (y) is selfconcordant for any a 0 in Q.
We proceed to prove that j 2 (ξ) := J 2 N (y + ξu) is a self-concordant function of ξ along every line u in Q. To this end, we have from direct calculations that Since V c (x) is convex and smooth, for any y i ∈ Q there exists a constant M v > 0 such that Let us define then we proceed to find that where the fourth inequality follows from the following Now, let us choose with M v satisfies (3.2), then J 2 N (y) is self-concordant for y ∈ Q with parameter a 0 . This implies that J N (y) is selfconcordant and the Newton's iteration is convergent thanks to Theorem 2.2.3 in [29].
Let {Φ i } i=0,1,...,N be any strictly increasing sequence with Φ 0 = 0 and Φ N = 1. Then the discrete free energy functional is defined as Then the following theorem states that the energy dissipation is preserved through the discretization.
Lemma 3.4. The discrete energy dissipation for the evolution of the discrete energy (3.4) holds for the scheme (2.11).
Proof. We rewrite According to the inequality ln( x y ) ≤ (x − y) 1 y and convexity of V c , we can obtain Then it follows from the summation-by-parts that Similarly, in light of the convexity of E 2 N we find that On the other hand, let us rewrite scheme (2.11) into the following then one infers from (3.6)-(3.7) that which is claimed.

3.2.
Convergence of the Euler Implicit Scheme. We proceed to study convergence of the discrete solution of (2.7) to a unique steady state in the long time. For this purpose we first prove the existence of this steady state. and our argument begins with the following observation. Proof. The trick of the proof is to find the "median" of the sequence. Let k 0 be the integer such that Therefore, one finds that These inequalities, together with the uniform boundedness |V (x)| ≤ M in [0, 1], enable us to estimate where we apply the fact 0 < Φ k 0 +1 − Φ k 0 < 1 for the second inequality and the Jensen inequality for the third one.
As a consequence we establish the existence of a stationary state realized as a minimizer of the free energy functional. Notice that lower bound of the discrete energy functional (3.5) diverges to −∞ as h → 0, and this is consistent with the fact that the continuous energy functional (2.5) is not bounded from below. Now, we are ready to show that the steady state of (2.4) must be a Heaviside step function.
Lemma 3.6. Denote {Φ i } i=0,1,...,N as the discrete solution of (2.7). Define {F k i } k∈N i=0,1,··· ,N as Then for any fix time step size τ > 0, the following statements hold: (i) the conservation law holds for the scheme (2.7) as follows If V (·) does not change sign in [0, 1], the steady state of the numerical scheme (2.7) must be {Φ * i } i=0,1,...,N . Proof. Summing over i, one can find from (2.7) that which implies that the conservation law in (i) holds for the scheme. Suppose that V (·) is of one sign in [0, 1]. Let {Φ * i } i=0,1,...,N be an equilibrium of scheme (2.7) such that Let m 1 be the largest spatial index for Φ * m1 = 0 and m 2 be the smallest for Φ * m2+2 = 1. Then (3.12) implies ) ≤ 0 leads a contradiction to V (·) > 0 since we then must have Φ * i ≡ 1 for i = m 1 + 1, ..., N . One also gets a contradiction if V (·) < 0. Therefore, the definitions of m 1 and m 2 imply that in either case the steady state is a nondecreasing step function from zero to one, whereas the jump location is determined by the conservation law (3.9). Now that {F k i } k∈N i=0,1,··· ,N is bounded in time thanks to (3.8), the Bolzano-Weierstrass theorem finds a subsequence convergence to {F * i } i=0,1,...,N and then the steady state can be written as (3.10) with a = 1 and m determined by (3.11). Furthermore, for any > 0, we can choose τ small enough such that for all n 2 > n 1 > 1 which implies that {F k i } k∈N i=0,1,··· ,N is Cauchy sequence with respect to k. Therefore, as τ goes to zero, we have the unique existence of {F * i } i=0,1,...,N , which implies the existence and uniqueness of the steady state as expected. Remark 3.7. For diffusion equation (1.1a), (3.13) and (3.14) imply that m 1 = m 2 as V ≡ 0, hence there exists a ∈ (0, 1) such that the steady state is given by (3.10). Moreover, the conservation of mass center within the discrete scheme implies that it the integer-valued function. The assumption V (·) being one sign is technical and our numerics suggest that the conclusion still holds if otherwise.
..,N be the discrete solution of (2.7). Define {F k i } k∈N i=0,1,··· ,N as follows Then for any smooth function V (·) the following statements are true: (i) the conservation law (3.9) holds for the scheme (2.7); Proof. The proof of the first statement is quite similar to that for Lemma 3.6 and it is omitted. To show that (3.10) under (3.11) is a steady state of (2.7), we first recall that the discrete function {Φ * i } i=0,1,...,N is an equilibrium if and In order to show that (3.10) supplemented by (3.11) is a steady state of the scheme (2.7), one only needs to verify the condition of Φ * m+1 . If V (·) is strictly positive or negative, (3.11) implies that a = 0 or a = 1. If V (·) = 0, the steady state is uniquely determined by (3.15). If V (·) changes sign in [0, 1], then (2.7) implies Substituting Φ * m = 0, Φ * m+1 = a and Φ * m+2 = 1 into (3.16) gives that V (a) = 0, which is expected. If {F k i } k∈N i=0,1,··· ,N has the limit {F * i } i=0,1,...,N , the parameter m can be determined by (3.11). In either case, we prove the statement in (ii). We next analyze the long time properties of the numerical scheme concerning its convergence to the steady state. We will first show the exponential decay for the diffusion equation (1.1a) for any fixed time step size modulo error terms. As for the replicator-diffusion equation (1.1) we will be able to obtain a similar result in the case of monotone drift fitness potentials.
for h small enough and C 0 both depending only on the distance of the initial data to the steady state Φ 0 − Φ * , that is, the fully discrete numerical scheme converges exponentially to the steady state in time modulo O(h) terms.
Proof. To show the convergence to the steady state, we recall from Lemma 3.6 that the discrete solution takes the form 0, for i = 0, 1, · · · , m, a, for i = m + 1, 1, for i = m + 2, · · · , N, with a and m to be determined. The time evolution of the L 2 -distance between Φ k and the stationary state Φ * reads To estimate A n , B n and C n , we apply (2.7) with V (·) = 0 and κ = 2 and deduce from the summation-by-part that where we skip the index k + 1 for simplicity. On the other hand, the conservation of mass center implies that Therefore, for any fixed τ > 0, there exists some k 0 ∈ N such that b k+1 ≤ 2 for any k > k 0 , and Finally, let us choose the space step small such that h ≤ Φ 0 − Φ * 2 /5. Then one concludes that 1 1+2τ k ≈ exp(−2kτ ) ≈ exp(−2t) for τ small, and this completes the proof.
Remark 3.10. When h goes to zero, we can deduce a posteriori uniqueness of the equilibrium from the convergence property. Indeed, letΦ * be another equilibrium state and set Φ k = Φ k+1 =Φ * in (3.18) with h = 0. Then we have from the evolution that Φ * − Φ * 2 ≤ 0 which proves the uniqueness of the stationary solution.
Remark 3.11. In the continuum case one can linearize the equation around the steady states (i.e., the Heaviside step function) and collects the following eigenvalue problem (3.19) −ϕ = λ x(1−x) ϕ, x ∈ (0, 1), ϕ(0) = ϕ(1) = 0, the principle eigenvalue of which determines the exponential convergence rate of (1.1) to the unique steady state (e.g., Theorem 3 in [11]). One notes that λ in (3.19) has a Rayleigh's variational quotient as and Theorem 262 in [16] implies that λ ≥ 2, with the principal eigenvalue λ = 2 achieved at ϕ(x) = Cx(1 − x) for some C ∈ R. This indicates that the exponential convergence rate in time of our discrete scheme in Theorem 3.9 is optimal since the decay rate cannot be better for generic initial data. Our numerical simulations in the coming section support this conclusion.
The second main result of this paper goes as follows.
Theorem 3.12. The solution of the numerical scheme (2.7) for the Replicator-Diffusion equation (1.1) with a monotone fitness potential (V either non-positive or non-negative) converges exponentially to the steady state of (1.1) with exponential rate κ (the diffusion rate) as the time step size τ goes to zero modulo O(h) terms. More precisely, the numerical scheme (2.7) for the general replicator-dynamics equation (1.1) with a monotone fitness potential satisfies for h small enough and C 0 both depending only on the distance of the initial data to the steady state Φ 0 − Φ * and the bound of V .
Proof. First of all, Lemma 3.8 implies that the sequence given by (3.10)-(3.11) is an equilibrium of (2.7). With |V (·)| ≤ M 1 in [0, 1] for some constant M 1 > 0, the same calculations in (3.17) lead to Therefore, one finds In light of the concavity of x(1 − x), we can obtain from Taylor's expansion that We next show the convergences of , then summing equation (2.7) with respect to i yields The fact that 0 ≤ c 1 k ≤ m + 1 for each k ∈ N implies that {c 1 k } is Cauchy with the limit a. Similarly, since 0 ≤ d k ≤ N − 1 and d k+1 − d k < 0, one can show that is Cauchy and its limit is (N + a − m − 2). Case ii). If V (·) < 0. Then we denote c 2 k+1 := , then summing equation (2.7) with respect to i yields The fact that 0 ≤ c 2 k ≤ N − m − 2 for each k ∈ N implies that {c 2 k } is Cauchy with the limit (N − m − 2). Similarly, since 0 ≤ d k ≤ N − 1 and d k+1 − d k > 0, one can show that is Cauchy and its limit is (N + a − m − 2), then {c 1 k } is Cauchy with the limit a.
Since c 1 k → a and d k → (N + a − m − 2) with 0 ≤ a ≤ 1, there exists some k 1 ∈ N such that for any k > k 1 Therefore from (3.21) we deduce that 1

Numerical Studies.
We now present several sets of numerical simulations to illustrate and verify the analytical properties of the proposed scheme and to demonstrate further spatio-temporal dynamics within the generic drift which are not captured by the analysis. Note that the time-dependent solution converges to the Dirac-delta function in the long-time limit. In our map variables, this corresponds to the convergence to a stationary map Φ ∞ (η) given by a step function. To deal with the issue that discrete values near the boundary are too close to be distinguished from each other under machine precision, we set the following additional boundary criteria throughout our simulations

Diffusion Equation.
We first numerically study the purely diffusive model (1.1a). The initial data are chosen to be either the symmetric f 0 (x) = π π−2 (1 − sin(πx)) or asymmetric f 0 (x) = 2x so we test the robustness of the proposed scheme with respect to initial data regarding the long-time dynamics. We remind the reader that the conserved quantity (1.4) determines asymptotically the proportion of the mass allocated to each end of the interval, that is the jump location in mass variables, see Theorem 2.1. To study the convergence of the numerical scheme, we consider the discrete L 1 error for Φ k,τ i,h , i = 0, 1, ..., N as follows is obtained with mesh size h = 1/1600 and time increment τ = 0.0005. Notice that this quantity is equivalent to an approximation of the Monge-Kantorovich or 1-Wasserstein transport distance between the associated densities. Moreover, the experimental order of convergence (EOC) in spatial-discretization h is defined as and the experimental order of convergence (EOC) in time-discretization τ is defined as EOC τ := log 2 (E τ ) − log 2 (E τ /2 ). Table 4.1 presents the errors and EOCs in spatial-discretization computed at a final time T = 0.992, a multiple of τ = 0.016, with constant time increment τ = 0.0005. We can see that difference errors decrease as the mesh is refined and a finer spatial-discretization results in higher accuracy. Besides, the orders of convergence in spatial-discretization indicate a convergence order of greater than one in both cases, and EOCs increase for a finer mesh. The results from both symmetric and asymmetric initial data indicate that the proposed scheme is robust to initial data. Table 4.2 shows the errors and EOCs in time-discretization with a fixed space mesh h = 1/1600. Similar robust convergence witnesses that the method converges as the time mesh shrinks.     1a) subject to symmetric initial data on the left and asymmetric initial data on the right. Both experiments capture the dominance of dynamics by the Dirac-delta as stated in Theorem 2.1. In particular, one finds that the degeneracy of diffusion at the endpoints makes the spatio-temporal dynamics of (1.1a) relatively simple in light of the proposed mass transport method now that the singularities can be quantified.
Next we present the evolution of (1.1a) out of the initial data f 0 (x) = π−2 π (1 − sin(πx)) and f 0 (x) = 2x with the step size h = 1/999 and time step τ = 0.001. Figure 4.1 captures the convergence to the steady-state of the numerical scheme given by a Heaviside step function, with at most one intermediate value, as stated in Lemma 3.6 or Lemma 2.1. We note that the jump locations η 0 = 1 2 and 1 3 for these two initial data, respectively, and this, together with the conservation of center of mass, are well preserved by our method. One finds that the (exponential) decay rate of the discrete scheme approaches 1 as space or time mesh shrinks. These findings agree well with Theorem 3.9.
transport distance between the associated densities. We can find that the convergence is faster than the rate obtained in Theorem 3.9 -exponential convergence in time modulo O(h) terms. This confirms the derivation of Theorem 3.9 and the existence of time modulo O(h) terms. We can also see that the convergence rates depend on τ and h, and as either one shrinks, the convergence of the L 2 -norm fits better e −t initially, which is consistent with the theoretical analysis, however, the decay rate is dramatically influenced by the O(h)-error in Theorem 3.9 as time evolves.

Replicator-Diffusion Equation.
We proceed to study the full Replicator-Diffusion equation (1.1) that includes the forward Kimura equation with frequency selection. In particular, we restrict our attention to the choice that V (x) = αx + β for some constant α and β. Throughout the simulations we fix f 0 (x) = 6x(1 − x) with the space mesh size h = 1/999 and time step size τ = 10 −3 , and then study the variation of the diffusion rate and fitness potential, starting with κ = 2, V (x) = 4x + 2 and κ = 4, V (x) = −4x + 2 as examples.
One observes from Figure 4.3-(a) that the numerical solution converges to the step function with the discrete conservation laws holding. As shown in Figure 4.3-(b), there exist some jumps in the dissipation of the entropy, which can be attributed to the additional boundary criteria (4.1). One can observe the smearing of these jumps by decreasing the tolerance in (4.1). However, this tolerance value cannot be chosen too small either for Newton's iterations to converge. Finally, Figure 4.3-(c) shows the decay in the L 2 -distance or 2-Wassertein transport distance towards equilibration. One again finds that the (exponential) decay rate of the discrete scheme approaches 1 as the space or time mesh shrinks for short times. This finding agrees well with Theorem 3.12. We also observe in Figure 4.3-(Top) that the asymmetric initial data leads to an interesting phenomena of different equilibration time scales, a first transient slower time scale in which all the mass from the left is pushed towards 0 followed by a faster decay towards the equilibrium for large times. We remark that the proofs of our main results rely on the assumption that V (x) is of one sign for technique purpose, and the numerics in the bottom of Figure 4.3 favor the application of our scheme otherwise. However, we have to exclude the relaxation of this assumption in this work, unfortunately.
4.3. Spatial-Temporal Dynamics. We conclude this section by presenting several additional sets of numerical experiments on the full model (1.1) to illustrate and verify our main results, and to capture the spatio-temporal dynamics of the random genetic drift when analytical tools lack. For these purposes, we choose featured initial data to showcase the dynamics of generic shrift to the variation of the potential function. For simplicity, we assume κ = 2 and consider two different fitness potential V (x) = 2 and V (x) = −3x + 1. The step size h = 1/999 and time step τ = 0.001 are fixed for all the numerical experiments. Each row in Figure 4.4 shows the evolution of the inverse cumulative distribution function Φ of (1.1) corresponding to the initial data depicted in each column. The dynamics are soon dominated by the Dirac-delta concentrations at the ends of the interval in each case, though one observes relatively different evolutionary transient phenomena such as merging of bumps and different equilibration speeds. All these findings point to the computational/mathematical evidence that the evolutionary changes at the molecular level are caused by random genetic drift when described by a Moran process. Finally, we report in Table 4.3 the error of the computed jump location against the theoretical one given in Theorem 2.1. They add another layer in the robustness of our proposed scheme as the location in Φ determines the portion of the alleles in the genetic aggregation.