On a Modified Random Genetic Drift Model: Derivation and a Structure-Preserving Operator-Splitting Discretization

TL;DR

The paper addresses the degeneracy-induced ill-posedness of the Kimura diffusion by replacing the domain with $(\delta,1-\delta)$ and introducing Robin-type boundaries, supplemented by boundary variables $a(t)$ and $b(t)$ to preserve mass and capture fixation dynamics. It derives the modified Kimura system via Energetic Variational Approaches, resulting in a bulk equation $\partial_t\rho = \partial_{xx}^2 (x(1-x)\rho)$ with coupled boundary dynamics, and selects $\alpha=2(1-\delta)$ to approximately conserve the first moment. A structure-preserving, Lagrangian–Eulerian operator-splitting discretization is developed, combining a Lagrangian interior step with an Eulerian boundary step to ensure mass conservation, positivity, and first-moment preservation. Numerical results illustrate the approach's ability to reproduce key Kimura features such as fixation and equilibrium behavior while achieving higher efficiency through a particle-based interior scheme. The work provides a robust computational framework and lays groundwork for extending to more complex, higher-dimensional genetic drift models.

Abstract

One of the fundamental mathematical models for studying random genetic drift is the Kimura equation, derived as the large-population limit of the discrete Wright-Fisher model. However, due to the degeneracy of the diffusion coefficient, it is impossible to impose a suitable boundary condition that ensures the Kimura equation admits a classical solution while preserving biological significance. In this work, we propose a modified model for random genetic drift that admits classical solutions by modifying the domain of the Kimura equation from $(0, 1)$ to $(δ, 1 - δ)$ with $δ$ being a small parameter, which allows us to impose a Robin-type boundary condition. By introducing two additional variables for the probabilities in the boundary region, we effectively capture the conservation of mass and the fixation dynamics in the original model. To numerically investigate the modified model, we develop a hybrid Eulerian-Lagrangian operator splitting scheme. The scheme first solves the flow map equation in the bulk region using a Lagrangian approach with a no-flux boundary condition, followed by handling the boundary dynamics in Eulerian coordinates. This hybrid scheme ensures mass conservation, maintains positivity, and preserves the first moment. Various numerical tests demonstrate the efficiency, accuracy, and structure-preserving properties of the proposed scheme. Numerical results demonstrate the key qualitative features of the original Kimura equation, including the fixation behavior and the correct stationary distribution in the small-$δ$ limit.

Figures (5)

• Figure 1: Density evolution for $\rho^{1}_{0}(X)$ at (a) $t = 0.1$, (b)$t = 0.2$, and (c) $t = 1.5$; and for $\rho^{2}_{0}(X)$ at (d) $t = 0.1$, (e) $t = 0.2$, and (f) $t= 1.5$ with $h = \frac{1 - 2\delta}{150}$, $\tau = \frac{1}{10000}$, $\alpha = 2(1-\delta)$, and $N = 150$
• Figure 2: Comparison of the Lagrangian mass functions with the Eulerian solutions: (a) $a(t)$ (upper left) and (b) $b(t)$ (upper right) with initial density $\rho_{0}^{1}$; (c) $a(t)$ (lower left) and (d) $b(t)$ with initial density $\rho^{2}_{0}$. $h = \frac{1 - 2\delta}{1200}$, $\tau = \frac{1}{10000}$, $\alpha = 2(1 - \delta)$, and $N = 1200$.
• Figure 3: Mass and first moment deviation $\mathcal{M}_{1}(t) - \mathcal{M}_{1}(0)$ evolution for $\rho_{0}^{1}$ in (a) and (b), and for $\rho_{0}^{2}$ in (c) and (d), with $\alpha = 2(1 -\delta)$.
• Figure 4: Energy and Particle Number evolution for $\rho_{0}^{2}$ with $\alpha = 2(1 -\delta)$. The initial decay of the energy is due to the diffusion of the particles within the domain. When the particles start moving towards the boundary, the energy increases and approaches the steady state, which suggests that all the particles have been absorbed into the boundary.
• Figure 5: Evolution of mass functions with $a(t)$ and $b(t)$ with the initial density $\rho_{0}^{1}$, different values of $\delta$, and $\alpha = 2(1-\delta)$. Here, $a_{ref}$ and $b_{ref}$ represent the mass functions with $\delta = 10^{-4}$.

Theorems & Definitions (12)

• Remark 3.1
• Proposition 3.1
• proof
• Remark 3.2
• Proposition 3.2
• proof
• Remark 3.3
• Remark 4.1
• Remark 4.2
• Proposition 4.1