Joakim Sundnes

On the Computational Complexity of the Bidomain and the Monodomain Models of Electrophysiology

The bidomain model, coupled with accurate models of cell membrane kinetics, is generally believed to provide a reasonable basis for numerical simulations of cardiac electrophysiology. Because of changes occurring in very short time intervals and over small spatial domains, discretized versions of these models must be solved on fine computational grids, and small time-steps must be applied. This leads to huge computational challenges that have been addressed by several authors. One popular way of reducing the CPU demands is to approximate the bidomain model by the monodomain model, and thus reducing a two by two set of partial differential equations to one scalar partial differential equation; both of which are coupled to a set of ordinary differential equations modeling the cell membrane kinetics. A reduction in CPU time of two orders of magnitude has been reported. It is the purpose of the present paper to provide arguments that such a reduction is not present when order-optimal numerical methods are applied. Theoretical considerations and numerical experiments indicate that the reduction factor of the CPU requirements from bidomain to monodomain computations, using order-optimal methods, typically is about 10 for simple cell models and less than two for more complex cell models.

Efficient solution of ordinary differential equations modeling electrical activity in cardiac cells

The contraction of the heart is preceded and caused by a cellular electro-chemical reaction, causing an electrical field to be generated. Performing realistic computer simulations of this process involves solving a set of partial differential equations, as well as a large number of ordinary differential equations (ODEs) characterizing the reactive behavior of the cardiac tissue. Experiments have shown that the solution of the ODEs contribute significantly to the total work of a simulation, and there is thus a strong need to utilize efficient solution methods for this part of the problem. This paper presents how an efficient implicit Runge-Kutta method may be adapted to solve a complicated cardiac cell model consisting of 31 ODEs, and how this solver may be coupled to a set of PDE solvers to provide complete simulations of the electrical activity.

Multigrid block preconditioning for a coupled system of partial differential equations modeling the electrical activity in the heart.

The electrical activity of the heart may be modeled with a system of partial differential equations (PDEs) known as the bidomain model. Computer simulations based on these equations may become a helpful tool to understand the relationship between changes in the electrical field and various heart diseases. Because of the rapid variations in the electrical field, sufficiently accurate simulations require a fine-scale discretization of the equations. For realistic geometries this leads to a large number of grid points and consequently large linear systems to be solved for each time step. In this paper, we present a fully coupled discretization of the bidomain model, leading to a block structured linear system. We take advantage of the block structure to construct an efficient preconditioner for the linear system, by combining multigrid with an operator splitting technique.

Simulation of ST segment changes during subendocardial ischemia using a realistic 3-D cardiac geometry

The mechanisms underlying the ST segment shifts associated with subendocardial ischemia remain unclear. The aim of this paper is to shed further light on the subject through numerical simulations of these shifts. A realistic three-dimensional model of the ventricles, including fiber rotation and anisotropy, is embedded in a nonhomogeneous torso model. A simplification of the bidomain model is used to calculate only the ST segment shift, assuming known values of the transmembrane potential during the plateau and rest phases. A similar simulation is performed in two dimensions. The simulation results suggest that subendocardial ischemia can be located by ST segment shift on the epicardial and torso surfaces. It is shown that ST elevation is associated with the transmural ischemic boundary, while ST depression is associated with the lateral ischemic boundaries.

Numerical solution of the bidomain equations

Knowledge of cardiac electrophysiology is efficiently formulated in terms of mathematical models. However, most of these models are very complex and thus defeat direct mathematical reasoning founded on classical and analytical considerations. This is particularly so for the celebrated bidomain model that was developed almost 40 years ago for the concurrent analysis of extra- and intracellular electrical activity. Numerical simulations based on this model represent an indispensable tool for studying electrophysiology. However, complex mathematical models, steep gradients in the solutions and complicated geometries lead to extremely challenging computational problems. The greatest achievement in scientific computing over the past 50 years has been to enable the solving of linear systems of algebraic equations that arise from discretizations of partial differential equations in an optimal manner, i.e. such that the central processing unit (CPU) effort increases linearly with the number of computational nodes. Over the past decade, such optimal methods have been introduced in the simulation of electrophysiology. This development, together with the development of affordable parallel computers, has enabled the solution of the bidomain model combined with accurate cellular models, on geometries resembling a human heart. However, in spite of recent progress, the full potential of modern computational methods has yet to be exploited for the solution of the bidomain model. This paper reviews the development of numerical methods for solving the bidomain model. However, the field is huge and we thus restrict our focus to developments that have been made since the year 2000.

A Second-Order Algorithm for Solving Dynamic Cell Membrane Equations

This paper describes an extension of the so-called Rush-Larsen scheme, which is a widely used numerical method for solving dynamic models of cardiac cell electrophysiology. The proposed method applies a local linearization of nonlinear terms in combination with the analytical solution of linear ordinary differential equations to obtain a second-order accurate numerical scheme. We compare the error and computational load of the second-order scheme to the original Rush-Larsen method and a second-order Runge-Kutta (RK) method. The numerical results indicate that the new method outperforms the original Rush-Larsen scheme for all the test cases. The comparison with the RK solver reveals that the new method is more efficient for stiff problems.

Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour

Models of cardiac mechanics are increasingly used to investigate cardiac physiology. These models are characterized by a high level of complexity, including the particular anisotropic material properties of biological tissue and the actively contracting material. A large number of independent simulation codes have been developed, but a consistent way of verifying the accuracy and replicability of simulations is lacking. To aid in the verification of current and future cardiac mechanics solvers, this study provides three benchmark problems for cardiac mechanics. These benchmark problems test the ability to accurately simulate pressure-type forces that depend on the deformed objects geometry, anisotropic and spatially varying material properties similar to those seen in the left ventricle and active contractile forces. The benchmark was solved by 11 different groups to generate consensus solutions, with typical differences in higher-resolution solutions at approximately 0.5%, and consistent results between linear, quadratic and cubic finite elements as well as different approaches to simulating incompressible materials. Online tools and solutions are made available to allow these tests to be effectively used in verification of future cardiac mechanics software.

A comparison of non-standard solvers for ODEs describing cellular reactions in the heart

Mathematical models for the electrical activity in cardiac cells are normally formulated as systems of ordinary differential equations (ODEs). The equations are nonlinear and describe processes occurring on a wide range of time scales. Under normal accuracy requirements, this makes the systems stiff and therefore challenging to solve numerically. As standard implicit solvers are difficult to implement, explicit solvers such as the forward Euler method are commonly used, despite their poor efficiency. Non-standard formulations of the forward Euler method, derived from the analytical solution of linear ODEs, can give significantly improved performance while maintaining simplicity of implementation. In this paper we study the performance of three non-standard methods on two different cell models with comparable complexity but very different stiffness characteristics.

Improved discretisation and linearisation of active tension in strongly coupled cardiac electro-mechanics simulations

Mathematical models of cardiac electro-mechanics typically consist of three tightly coupled parts: systems of ordinary differential equations describing electro-chemical reactions and cross-bridge dynamics in the muscle cells, a system of partial differential equations modelling the propagation of the electrical activation through the tissue and a nonlinear elasticity problem describing the mechanical deformations of the heart muscle. The complexity of the mathematical model motivates numerical methods based on operator splitting, but simple explicit splitting schemes have been shown to give severe stability problems for realistic models of cardiac electro-mechanical coupling. The stability may be improved by adopting semi-implicit schemes, but these give rise to challenges in updating and linearising the active tension. In this paper we present an operator splitting framework for strongly coupled electro-mechanical simulations and discuss alternative strategies for updating and linearising the active stress component. Numerical experiments demonstrate considerable performance increases from an update method based on a generalised Rush–Larsen scheme and a consistent linearisation of active stress based on the first elasticity tensor.

Electromechanical feedback with reduced cellular connectivity alters electrical activity in an infarct injured left ventricle: a finite element model study

Myocardial infarction (MI) significantly alters the structure and function of the heart. As abnormal strain may drive heart failure and the generation of arrhythmias, we used computational methods to simulate a left ventricle with an MI over the course of a heartbeat to investigate strains and their potential implications to electrophysiology. We created a fully coupled finite element model of myocardial electromechanics consisting of a cellular physiological model, a bidomain electrical diffusion solver, and a nonlinear mechanics solver. A geometric mesh built from magnetic resonance imaging (MRI) measurements of an ovine left ventricle suffering from a surgically induced anteroapical infarct was used in the model, cycled through the cardiac loop of inflation, isovolumic contraction, ejection, and isovolumic relaxation. Stretch-activated currents were added as a mechanism of mechanoelectric feedback. Elevated fiber and cross fiber strains were observed in the area immediately adjacent to the aneurysm throughout the cardiac cycle, with a more dramatic increase in cross fiber strain than fiber strain. Stretch-activated channels decreased action potential (AP) dispersion in the remote myocardium while increasing it in the border zone. Decreases in electrical connectivity dramatically increased the changes in AP dispersion. The role of cross fiber strain in MI-injured hearts should be investigated more closely, since results indicate that these are more highly elevated than fiber strain in the border of the infarct. Decreases in connectivity may play an important role in the development of altered electrophysiology in the high-stretch regions of the heart.