Efficient cardiac simulations using the Runge--Kutta--Chebyshev method

Abstract

Heart disease is one of the leading causes of death in Canada, claiming thousands of lives each year. Cardiac electrophysiology that studies the electrical activity in the human heart has emerged as an active research field in response to the demand for providing reliable guidance for clinical diagnosis and treatment to heart arrhythmias. Computer simulation of electrophysiological phenomena provides a non-invasive way to study the electrical activity in the human heart and to provide quantitative guidance to clinical applications.

With the need to unravel underlying physiological details, mathematical models tend to be large and possess characteristics that are challenging to mitigate. In this thesis, we describe numerical methods for solving widely used mathematical models: the bidomain model and its simplified form, the monodomain model. The bidomain model is a multi-scale cardiac electrophysiology model that includes a set of reaction-diffusion partial differential equations (PDEs) with the reaction term representing cardiac cell models that describes the chemical reactions and flows of ions across the cell membrane of myocardial cells at the micro level and the diffusion term representing current propagation through the heart at the macro level. We use the method of lines (MOL) to obtain numerical solution of this model. The MOL first spatially discretizes the system of PDEs, resulting in a system of ordinary differential equations (ODEs) at each space point, and we obtain fully discrete solutions at each space-time point using time-integration methods for ODEs. In this thesis, we propose innovative numerical methods for the time integration of systems of ODEs based on the Runge--Kutta--Chebyshev (RKC) method. We implement and compare our methods with those used by current research on time integration of ODEs on three problems: time integration of individual cardiac cell models, time integration of the cell model of a monodomain problem, and time integration of spatially discretized tissue equation in a monodomain benchmark problem proposed by S. Niederer et al. in 2011.

Numerical methods in cardiac electrophysiology research for solving ODEs include the forward Euler (FE) method, the Rush--Larsen (RL) method, the backward Euler method, and the generalized RL method of first-order. We introduce multistage first-order RKC methods and multistage first-order RL methods that are constructed by replacing the FE method with multistage first-order RKC methods. We implement all the aforementioned methods and test their efficiencies in time integration of 37 cardiac cell models. We find introducing the multistage RKC and RL methods allows larger step sizes to meet prescribed numerical accuracy; the increased time steps sped up time integration of 19 cell models. We replace the FE method with two-stage RKC method in time integration of cell model in a monodomain model. We find the increased time step introduced by applying this method improved the entire solving process by up to a factor of 1.4. We also apply the RKC(2,1) method to time integration of the tissue equation from a monodomain benchmark problem. Results show we have decreased the execution time of this benchmark problem by a factor of two. We note the increase of time step is from stability improvement brought by the numerical method. We finally give a quantitative explanation of stability improvement from introducing multistage RKC and RL methods for solving systems of ODEs considered in this thesis.