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

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

Towards New High-Order Operator Splitting Time-Integration Methods

Operator splitting (OS) methods represent a powerful strategy to solve an extensive range of mathematical models in the form of differential equations. They have a long and celebrated history, having been successfully used for well over half a century to provide efficient numerical solutions to challenging problems. In fact, OS methods are often the only viable way to solve many problems in practice. The simplest, and perhaps, most well-known OS methods are Lie--Trotter--Godunov and the Strang--Marchuk methods. They compute a numerical solution that is first-, and second-order accurate in time, respectively. OS methods can be derived by imposing order conditions using the Campbell--Baker--Hausdorff formula. It follows that, by setting the appropriate order conditions, it is possible to derive OS methods of any desired order. An important observation regarding OS methods with order higher than two is that, according to the Sheng--Suzuki theorem, at least one of their defining coefficients must be negative. Therefore, the time integration with OS methods of order higher than two has not been considered suitable to solve deterministic parabolic problems, because the necessary backward time integration would cause instabilities. Throughout this thesis, we focus our attention on high-order (i.e., order higher than two) OS methods. We successfully assess the convergence properties of some higher-order OS methods on diffusion-reaction problems describing cardiac electrophysiology and on an advection-diffusion-reaction problem describing chemical combustion. Furthermore, we compare the efficiency performance of higher-order methods to second-order methods. For all the cases considered, we confirm an improved efficiency performance compared to methods of lower order. Next, we observe how, when using OS and Runge--Kutta type methods to advance the time integration, we can construct a unique extended Butcher tableau with a similar structure to the ones describing Generalized Additive Runge--Kutta (GARK) methods. We define a combination of methods to be OS-GARK methods. We apply linear stability analysis to OS-GARK methods; this allows us to conveniently analyze the stability properties of any combination of OS and Runge--Kutta methods. Doing so, we are able to perform an eigenvalue analysis to understand and improve numerically unstable solutions

A Comparison of the Bidomain and EMI Models in Refractory Cardiac Tissue

Computational cardiac modelling has made incredible strides over the past 40 years toward becoming an integral component of healthcare. The majority of cardiac modelling is accomplished using the bidomain or monodomain models, equations describing electrical conduction in cardiac tissue. These models use a volume averaging approach in which the structure of individual cells is disregarded; instead, cells are treated homogeneously as a continuum. Although this approach often provides an adequate view of cardiac activity at the macro level, there are situations where this approximation is insufficient, such as when discontinuities at the cellular level are implicated in a given disease or phenomenon. To address this, a more detailed tissue model has recently been developed: the extracellular-membrane-intracellular (EMI) model. The EMI model explicitly defines the extracellular, membrane, and intracellular compartments to form a highly detailed model of cardiac tissue. However, this additional level of detail also poses a high computational cost. This thesis investigates the trade-off that exists between the conventional bidomain model and the EMI model. To do this, we carry out a comparison study. This constitutes the first EMI comparison study that has been conducted outside of the research group that developed the model. Using both models, we find the currents required to trigger consecutive action potentials at varying time intervals. We then use these data points to construct refractory profiles for each model and compare these profiles against available experimental data. Our findings demonstrate that within the framework of this study, the behaviour of the EMI model is noticeably closer to experimental data than the behaviour of the bidomain model. These results have implications on the way we approach tissue model selection in the future, as well as for our general understanding of the refractory properties of cardiac tissue