Efficient Fully Implicit Time Integration Methods for Modeling Cardiac Dynamics

Implicit methods are well known to have greater stability than explicit methods for stiff systems, but they often are not used in practice due to perceived computational complexity. This paper applies the backward Euler (BE) method and a second-order one-step two-stage composite backward differentiation formula (C-BDF2) for the monodomain equations arising from mathematically modeling the electrical activity of the heart. The C-BDF2 scheme is an L-stable implicit time integration method and easily implementable. It uses the simplest forward Euler and BE methods as fundamental building blocks. The nonlinear system resulting from application of the BE method for the monodomain equations is solved for the first time by a nonlinear elimination method, which eliminates local and nonsymmetric components by using a Jacobian-free Newton solver, called Newton--Krylov solver. Unlike other fully implicit methods proposed for the monodomain equations in the literature, the Jacobian of the global system after the nonlinear elimination has much smaller size, is symmetric and possibly positive definite, which can be solved efficiently by standard optimal solvers. Numerical results are presented demonstrating that the C-BDF2 scheme can yield accurate results with less CPU times than explicit methods for both a single patch and spatially extended domains.

[1]  Andrew J. Pullan,et al.  Mathematically Modelling the Electrical Activity of the Heart: From Cell to Body Surface and Back Again , 2005 .

[2]  William L. Briggs,et al.  A multigrid tutorial , 1987 .

[3]  J. Butcher The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods , 1987 .

[4]  C. Kelley Solving Nonlinear Equations with Newton's Method , 1987 .

[5]  Joakim Sundnes,et al.  Computing the electrical activity in the heart , 2006 .

[6]  Linda R. Petzold,et al.  Using Krylov Methods in the Solution of Large-Scale Differential-Algebraic Systems , 1994, SIAM J. Sci. Comput..

[7]  Wenjun Ying,et al.  A multilevel adaptive approach for computational cardiology , 2005 .

[8]  N. Hooke Efficient simulation of action potential propagation in a bidomain , 1992 .

[9]  J. Keener,et al.  A numerical method for the solution of the bidomain equations in cardiac tissue. , 1998, Chaos.

[10]  Michael Mascagni,et al.  The backward Euler method for numerical solution of the Hodgkin-Huxley equations of nerve conduction , 1990 .

[11]  Begnaud Francis Hildebrand,et al.  Introduction to numerical analysis: 2nd edition , 1987 .

[12]  Yousef Saad,et al.  Convergence Theory of Nonlinear Newton-Krylov Algorithms , 1994, SIAM J. Optim..

[13]  NeuronsSteven J. Cox,et al.  A Fast Fully Implicit Backward Euler Solver for Dendritic Neurons , 2000 .

[14]  C. Kelley Iterative Methods for Linear and Nonlinear Equations , 1987 .

[15]  J. Verwer,et al.  Numerical solution of time-dependent advection-diffusion-reaction equations , 2003 .

[16]  William L. Briggs,et al.  A multigrid tutorial, Second Edition , 2000 .

[17]  D. Brandt,et al.  Multi-level adaptive solutions to boundary-value problems math comptr , 1977 .

[18]  R. C. Barr,et al.  Propagation model using the DiFrancesco-Noble equations , 1992, Medical and Biological Engineering and Computing.

[19]  C. Henriquez,et al.  Cardiac propagation simulation. , 1992, Critical reviews in biomedical engineering.

[20]  J. Butcher Numerical Methods for Ordinary Differential Equations: Butcher/Numerical Methods , 2005 .

[21]  M. Courtemanche,et al.  Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. , 1998, The American journal of physiology.

[22]  Henk A. van der Vorst,et al.  Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems , 1992, SIAM J. Sci. Comput..

[23]  Randolph E. Bank,et al.  Transient simulation of silicon devices and circuits , 1985 .

[24]  D DiFrancesco,et al.  A model of cardiac electrical activity incorporating ionic pumps and concentration changes. , 1985, Philosophical transactions of the Royal Society of London. Series B, Biological sciences.

[25]  Rodrigo Weber dos Santos,et al.  Parallel multigrid preconditioner for the cardiac bidomain model , 2004, IEEE Transactions on Biomedical Engineering.

[26]  James P. Keener,et al.  Mathematical physiology , 1998 .

[27]  Natalia A. Trayanova,et al.  Computational techniques for solving the bidomain equations in three dimensions , 2002, IEEE Transactions on Biomedical Engineering.

[28]  Q. Sheng Solving Linear Partial Differential Equations by Exponential Splitting , 1989 .

[29]  A. Garfinkel,et al.  An advanced algorithm for solving partial differential equation in cardiac conduction , 1999, IEEE Transactions on Biomedical Engineering.

[30]  Arun V. Holden,et al.  Computational biology of the heart , 1998, The Mathematical Gazette.

[31]  K. Stüben A review of algebraic multigrid , 2001 .

[32]  J. Cooley,et al.  Digital computer solutions for excitation and propagation of the nerve impulse. , 1966, Biophysical journal.

[33]  L. Shampine,et al.  Numerical Solution of Ordinary Differential Equations. , 1995 .

[34]  Xiao-Chuan Cai,et al.  A fully implicit parallel algorithm for simulating the non‐linear electrical activity of the heart , 2004, Numer. Linear Algebra Appl..

[35]  Paul J. Lanzkron,et al.  An Analysis of Approximate Nonlinear Elimination , 1996, SIAM J. Sci. Comput..

[36]  G. Bett,et al.  Computer model of action potential of mouse ventricular myocytes. , 2004, American journal of physiology. Heart and circulatory physiology.

[37]  D. Rose,et al.  Linear algebraic transformations of the bidomain equations: implications for numerical methods. , 1994, Mathematical biosciences.

[38]  A. Hodgkin,et al.  A quantitative description of membrane current and its application to conduction and excitation in nerve , 1990 .

[39]  D. Rose,et al.  Global approximate Newton methods , 1981 .

[40]  A. Tveito,et al.  An operator splitting method for solving the bidomain equations coupled to a volume conductor model for the torso. , 2005, Mathematical biosciences.

[41]  J. Trangenstein,et al.  Operator splitting and adaptive mesh refinement for the Luo-Rudy I model , 2004 .

[42]  Y. Saad,et al.  GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems , 1986 .

[43]  D. Rose,et al.  Analysis of a multilevel iterative method for nonlinear finite element equations , 1982 .

[44]  J. M. Watt Numerical Initial Value Problems in Ordinary Differential Equations , 1972 .

[45]  E. Hairer,et al.  Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems , 2010 .

[46]  Tasso J. Kaper,et al.  N th-order operator splitting schemes and nonreversible systems , 1996 .

[47]  PETER K. MOORE AN ADAPTIVE H-REFINEMENT FINITE ELEMENT METHOD FOR PARABOLIC DIFFERENTIAL SYSTEMS IN THREE SPACE DIMENSIONS , 2000 .

[48]  Kendall E. Atkinson An introduction to numerical analysis , 1978 .

[49]  J. Lambert Numerical Methods for Ordinary Differential Equations , 1991 .

[50]  Peter K. Moore,et al.  An Adaptive Finite Element Method for Parabolic Differential Systems: Some Algorithmic Considerations in Solving in Three Space Dimensions , 1999, SIAM J. Sci. Comput..

[51]  Yousef Saad,et al.  Iterative methods for sparse linear systems , 2003 .