A boundary preserving numerical algorithm for the Wright-Fisher model with mutation

• Published in: BIT (Nordisk Tidskrift for Informationsbehandling) Vol. 52; no. 2; pp. 283 - 304
• Main Authors: Dangerfield, C. E., Kay, D., MacNamara, S., Burrage, K.
• Format: Journal Article
• Language: English
• Published: Dordrecht Springer Netherlands 01.06.2012; Springer
• Subjects: Computational Mathematics and Numerical Analysis; Exact sciences and technology; Global analysis, analysis on manifolds; Mathematical analysis; Mathematics; Mathematics and Statistics; Numeric Computing; Numerical analysis; Numerical analysis. Scientific computation; Numerical approximation; Numerical methods in probability and statistics; Sciences and techniques of general use; Sequences, series, summability; Topology. Manifolds and cell complexes. Global analysis and analysis on manifolds; 65C20; 92D99; Stochastic differential equations; Strong convergence; 65C30; Boundary preserving numerical algorithm; 65L20; Wright-Fisher model; Ion channels; Split step; Hölder condition; 60H35; Stochastic model; Differential equation; Probability distribution; Numerical method; Stochastic method; Algorithm; Euler scheme; Numerical approximation; Numerical analysis; Numerical solution; Integral equation; Stochastic equation
• ISSN: 0006-3835; 1572-9125
• DOI: 10.1007/s10543-011-0351-3

The Wright-Fisher model is an Itô stochastic differential equation that was originally introduced to model genetic drift within finite populations and has recently been used as an approximation to ion channel dynamics within cardiac and neuronal cells. While analytic solutions to this equation remain within the interval [0,1], current numerical methods are unable to preserve such boundaries in the approximation. We present a new numerical method that guarantees approximations to a form of Wright-Fisher model, which includes mutation, remain within [0,1] for all time with probability one. Strong convergence of the method is proved and numerical experiments suggest that this new scheme converges with strong order 1/2. Extending this method to a multidimensional case, numerical tests suggest that the algorithm still converges strongly with order 1/2. Finally, numerical solutions obtained using this new method are compared to those obtained using the Euler-Maruyama method where the Wiener increment is resampled to ensure solutions remain within [0,1].