Numerical methods for simulation of electrical activity in the myocardial tissue

Mathematical models of electric activity in cardiac tissue are becoming increasingly powerful tools in the study of cardiac arrhythmias. Considered here are mathematical models based on ordinary differential equations (ODEs) and partial differential equations (PDEs) that describe the behaviour of this electrical activity. Generating an efficient numerical solution of these models is a challenging task, and in fact the physiological accuracy of tissue-scale models is often limited by the efficiency of the numerical solution process. In this thesis, we discuss two sets of experiments that test ideas for making the numerical solution process more efficient. In the first set of experiments, we examine the numerical solution of four single cell cardiac electrophysiological models, which consist solely of ODEs. We study the efficiency of using implicit-explicit Runge-Kutta (IMEX-RK) splitting methods to solve these models. We find that variable step-size implementations of IMEX-RK methods (ARK3 and ARK5) that take advantage of Jacobian structure clearly outperform most methods commonly used in practice for two of the models, and they outperform all methods commonly used in practice for the remaining models. In the second set of experiments, we examine the solution of the bidomain model, a model consisting of both ODEs and PDEs that are typically solved separately. We focus these experiments on numerical methods for the solution of the two PDEs in the bidomain model. The most popular method for this task, the Crank-Nicolson method, produces unphysical oscillations; we propose a method based on a second-order L-stable singly diagonally implicit Runge-Kutta (SDIRK) method to eliminate these oscillations. We find that although the SDIRK method is able to eliminate these unphysical oscillations, it is only more efficient for crude error tolerances

A multiresolution space-time adaptive scheme for the bidomain model in electrocardiology

This work deals with the numerical solution of the monodomain and bidomain models of electrical activity of myocardial tissue. The bidomain model is a system consisting of a possibly degenerate parabolic PDE coupled with an elliptic PDE for the transmembrane and extracellular potentials, respectively. This system of two scalar PDEs is supplemented by a time-dependent ODE modeling the evolution of the so-called gating variable. In the simpler sub-case of the monodomain model, the elliptic PDE reduces to an algebraic equation. Two simple models for the membrane and ionic currents are considered, the Mitchell-Schaeffer model and the simpler FitzHugh-Nagumo model. Since typical solutions of the bidomain and monodomain models exhibit wavefronts with steep gradients, we propose a finite volume scheme enriched by a fully adaptive multiresolution method, whose basic purpose is to concentrate computational effort on zones of strong variation of the solution. Time adaptivity is achieved by two alternative devices, namely locally varying time stepping and a Runge-Kutta-Fehlberg-type adaptive time integration. A series of numerical examples demonstrates thatthese methods are efficient and sufficiently accurate to simulate the electrical activity in myocardial tissue with affordable effort. In addition, an optimalthreshold for discarding non-significant information in the multiresolution representation of the solution is derived, and the numerical efficiency and accuracy of the method is measured in terms of CPU time speed-up, memory compression, and errors in different norms.Comment: 25 pages, 41 figure

Adaptive step ODE algorithms for the 3D simulation of electric heart activity with graphics processing units

In this paper we studied the implementation and performance of adaptive step methods for large systems of ordinary differential equations systems in graphics processing units, focusing on the simulation of three-dimensional electric cardiac activity. The Rush-Larsen method was applied in all the implemented solvers to improve efficiency. We compared the adaptive methods with the fixed step methods, and we found that the fixed step methods can be faster while the adaptive step methods are better in terms of accuracy and robustness. (c) 2013 Elsevier Ltd. All rights reserved.This work has been partially funded by Universitat Politecnica de Valencia through Programa de Apoyo a la InvestigaciOn y Desarrollo (PAID-06-11) and (PAID-05-12), by Generalitat Valenciana through projects PROMETEO/2009/013 and Ayudas para la realizacion de proyectos de I+D para grupos de investigacion emergentes GV/2012/039, and by Ministerio Espafiol de Economia y Competitividad through project TEC2012-38142-004.García Mollá, VM.; Liberos Mascarell, A.; Vidal Maciá, AM.; Guillem Sánchez, MS.; Millet Roig, J.; González Salvador, A.; Martínez Zaldívar, FJ.... (2014). Adaptive step ODE algorithms for the 3D simulation of electric heart activity with graphics processing units. Computers in Biology and Medicine. 44:15-26. https://doi.org/10.1016/j.compbiomed.2013.10.023S15264

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

A TIME-AND-SPACE PARALLELIZED ALGORITHM FOR THE CABLE EQUATION

Electrical propagation in excitable tissue, such as nerve fibers and heart muscle, is described by a nonlinear diffusion-reaction parabolic partial differential equation for the transmembrane voltage $V(x,t)$, known as the cable equation. This equation involves a highly nonlinear source term, representing the total ionic current across the membrane, governed by a Hodgkin-Huxley type ionic model, and requires the solution of a system of ordinary differential equations. Thus, the model consists of a PDE (in 1-, 2- or 3-dimensions) coupled to a system of ODEs, and it is very expensive to solve, especially in 2 and 3 dimensions. In order to solve this equation numerically, we develop an algorithm, extended from the Parareal Algorithm, to efficiently incorporate space-parallelized solvers into the framework of the Parareal algorithm, to achieve time-and-space parallelization. Numerical results and comparison of the performance of several serial, space-parallelized and time-and-space-parallelized time-stepping numerical schemes in one-dimension and in two-dimensions are also presented

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

Fast Simulations of Models of Cardiac Electrophysiology

Cardiac electrophysiology studies the electrical activity of the heart. Researchers from different fields are working together to model and simulate the electrical activity in the heart. Such simulations may lead to help cardiologists to treat a patient's heart condition with better techniques and diagnoses. Mathematical models of cardiac electrophysiology are described by a system of partial differential equations~(PDEs) and a non-linear system of ordinary differential equations~(ODEs). One way to reach real-time simulation of the electrical activity of the heart is through more efficient time integration of these ODEs. Larger time steps lead to more efficient computations provided the solutions remain accurate enough. Producing the largest stable time step for solving these ODEs is a daunting task. Usually, trial and error is used to get the largest stable time step, but it is time consuming. In this thesis, we propose a new efficient method to find the largest stable time step to solve these ODEs efficiently. In this thesis, we present thirty-seven cardiac cell models. We compare, the forward Euler method, the explicit midpoint method, the four-stage, fourth-order Runge--Kutta method, the two-stage, first-order Runge--Kutta--Chebyshev method, the three-stage, first-order Runge--Kutta--Chebyshev method, the three-stage, third-order strong-stability-preserving Runge--Kutta method, and the four-stage, third-order strong-stability-preserving Runge--Kutta method based on their largest stepsize, mixed root mean square errors, and CPU time to solve the cardiac cell models. From the theoretical largest stepsize results, the forward Euler method outperforms all the other methods considered on all thirty-seven cardiac cell models. Next to the forward Euler method, the two-stage, first-order Runge--Kutta--Chebyshev method outperforms all the other methods considered on thirty-seven cardiac cell models. From the experimental largest stepsize results, the forward Euler method outperforms all the other methods considered on all thirty-seven cardiac cell models. Next to the forward Euler method, the two-stage, first-order Runge--Kutta--Chebyshev and the explicit midpoint methods outperform all the other methods considered on thirty-six cardiac cell models and one cardiac cell model, respectively. Also, we solve the monodomain problem coupled with FitzHugh--Nagumo model using the forward Euler, the two-stage, first-order Runge--Kutta--Chebyshev method, and the three-stage, first-order Runge--Kutta--Chebyshev method in one, two, and three dimensions. We compare these methods based on their theoretical largest stepsize, experimental largest stepsize, CPU time, and mixed root mean square errors. From these results, the three-stage, first-order Runge--Kutta--Chebyshev method outperforms all the other methods over a one-dimensional problem, a two-dimensional problem, and a three-dimensional problem

Efficient Numerical Methods for Heart Simulation

The heart is one the most important organs in the human body and many other live creatures. The electrical activity in the heart controls the heart function, and many heart diseases are linked to the abnormalities in the electrical activity in the heart. Mathematical equations and computer simulation can be used to model the electrical activity in the heart. The heart models are challenging to solve because of the complexity of the models and the huge size of the problems. Several cell models have been proposed to model the electrical activity in a single heart cell. These models must be coupled with a heart model to model the electrical activity in the entire heart. The bidomain model is a popular model to simulate the propagation of electricity in myocardial tissue. It is a continuum-based model consisting of non-linear ordinary differential equations (ODEs) describing the electrical activity at the cellular scale and a system of partial differential equations (PDEs) describing propagation of electricity at the tissue scale. Because of this multi-scale, ODE/PDE structure of the model, splitting methods that treat the ODEs and PDEs in separate steps are natural candidates as numerical methods. First, we need to solve the problem at the cellular scale using ODE solvers. One of the most popular methods to solve the ODEs is known as the Rush-Larsen (RL) method. Its popularity stems from its improved stability over integrators such as the forward Euler (FE) method along with its easy implementation. The RL method partitions the ODEs into two sets: one for the gating variables, which are treated by an exponential integrator, and another for the remaining equations, which are treated by the FE method. The success of the RL method can be understood in terms of its relatively good stability when treating the gating variables. However, this feature would not be expected to be of benefit on cell models for which the stiffness is not captured by the gating equations. We demonstrate that this is indeed the case on a number of stiff cell models. We further propose a new partitioned method based on the combination of a first-order generalization of the RL method with the FE method. This new method leads to simulations of stiff cell models that are often one or two orders of magnitude faster than the original RL method. After solving the ODEs, we need to use bidomain solvers to solve the bidomain model. Two well-known, first-order time-integration methods for solving the bidomain model are the semi-implicit method and the Godunov operator-splitting method. Both methods decouple the numerical procedure at the cellular scale from that at the tissue scale but in slightly different ways. The methods are analyzed in terms of their accuracy, and their relative performance is compared on one-, two-, and three-dimensional test cases. As suggested by the analysis, the test cases show that the Godunov method is significantly faster than the semi-implicit method for the same level of accuracy, specifically, between 5 and 15 times in the cases presented. Second-order bidomain solvers can generally be expected to be more effective than first-order bidomain solvers under normal accuracy requirements. However, the simplest and the most commonly applied second-order method for the PDE step, the Crank-Nicolson (CN) method, may generate unphysical oscillations. We investigate the performance of a two-stage, L-stable singly diagonally implicit Runge-Kutta method for solving the PDEs of the bidomain model and present a stability analysis. Numerical experiments show that the enhanced stability property of this method leads to more physically realistic numerical simulations compared to both the CN and Backward Euler (BE) methods

GEMS: A Fully Integrated PETSc-Based Solver for Coupled Cardiac Electromechanics and Bidomain Simulations

Cardiac contraction is coordinated by a wave of electrical excitation which propagates through the heart. Combined modeling of electrical and mechanical function of the heart provides the most comprehensive description of cardiac function and is one of the latest trends in cardiac research. The effective numerical modeling of cardiac electromechanics remains a challenge, due to the stiffness of the electrical equations and the global coupling in the mechanical problem. Here we present a short review of the inherent assumptions made when deriving the electromechanical equations, including a general representation for deformation-dependent conduction tensors obeying orthotropic symmetry, and then present an implicit-explicit time-stepping approach that is tailored to solving the cardiac mono- or bidomain equations coupled to electromechanics of the cardiac wall. Our approach allows to find numerical solutions of the electromechanics equations using stable and higher order time integration. Our methods are implemented in a monolithic finite element code GEMS (Ghent Electromechanics Solver) using the PETSc library that is inherently parallelized for use on high-performance computing infrastructure. We tested GEMS on standard benchmark computations and discuss further development of our software

Efficient Implicit Runge-Kutta Methods for Fast-Responding Ligand-Gated Neuroreceptor Kinetic Models

Neurophysiological models of the brain typically utilize systems of ordinary differential equations to simulate single-cell electrodynamics. To accurately emulate neurological treatments and their physiological effects on neurodegenerative disease, models that incorporate biologically-inspired mechanisms, such as neurotransmitter signalling, are necessary. Additionally, applications that examine populations of neurons, such as multiscale models, can demand solving hundreds of millions of these systems at each simulation time step. Therefore, robust numerical solvers for biologically-inspired neuron models are vital. To address this requirement, we evaluate the numerical accuracy and computational efficiency of three L-stable implicit Runge-Kutta methods when solving kinetic models of the ligand-gated glutamate and gamma-aminobutyric acid (GABA) neurotransmitter receptors. Efficient implementations of each numerical method are discussed, and numerous performance metrics including accuracy, simulation time steps, execution speeds, Jacobian calculations, and LU factorizations are evaluated to identify appropriate strategies for solving these models. Comparisons to popular explicit methods are presented and highlight the advantages of the implicit methods. In addition, we show a machine-code compiled implicit Runge-Kutta method implementation that possesses exceptional accuracy and superior computational efficiency