High order time-stepping methods for cardiac electrophysiology models

Cardiac simulations require to solve as efficiently as possible the monodomain or bidomain equations coupled to ionic models. Ionic models can be very large systems of nonlinear, stiff differential equations. They are solved on spatial discretization grids of the heart, that for usual simulations count millions of degrees of freedom. These models are stiff, which prompts for using implicit solvers. But, since part of the equations is semi-linear, several efficient methods have been proposed, following the first idea of [Rush, Larsen, IEEE TBME, 1978], that consists in solving implicitly only the semi-linear equations. Anyway, all current methods, based on a 1st order linearization, still require time-steps of the order of 20 µs: the equations are solved 50 times to compute 1 ms of real time, 500 000 times for 10 s. Following a work of [Perego, Venziani, ETNA, 2009], we developped a new whole family of very high order methods dedicated to cardiac solvers, and to more general stiff problems. It is based upon the emerging techniques of exponential integrator as developped in [Hochbruck, Ostermann, Acta Numerica, 2010]. In this work , we propose a family of accurate explicit methods of the Exponential Adams Bashforth type. Beside in-depth analysis, convergence proofs, and error estimates, we give the expression of such time-stepping schemes for an approximation of order 2, 3, and 4, and the general expression of arbitrary order k. We will present synthetic numerical results for cardiac ionic models. In conclusion, the 3rd and 4th order method could run with a time-step of 100µs, with still a very good accuracy, and no additional computational cost per iteration. In future work, we plan to combine these methods with high-order spatial discretization, and provide numerical techniques to improve the relevance and overall resolution of cardiac models.