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.     

An operator splitting method for solving the bidomain equations coupled to a volume conductor model for the torso

In this paper we present a numerical method for the bidomain model, which describes the electrical activity in the heart. The model consists of two partial differential equations (PDEs), which are coupled to systems of ordinary differential equations (ODEs) describing electrochemical reactions in the cardiac cells. Many applications require coupling these equations to a third PDE, describing the electrical fields in the torso surrounding the heart. The resulting system is challenging to solve numerically, because of its complexity and very strict resolution requirements in time and space. We propose a method based on operator splitting and a fully coupled discretization of the three PDEs. Numerical experiments show that for simple simulation cases and fine discretizations, the algorithm is second-order accurate in space and time.     

Solvers for the cardiac bidomain equations

The bidomain equations are widely used for the simulation of electrical activity in cardiac tissue. They are especially important for accurately modeling extracellular stimulation, as evidenced by their prediction of virtual electrode polarization before experimental verification. However, solution of the equations is computationally expensive due to the fine spatial and temporal discretization needed. This limits the size and duration of the problem which can be modeled. Regardless of the specific form into which they are cast, the computational bottleneck becomes the repeated solution of a large, linear system. The purpose of this review is to give an overview of the equations and the methods by which they have been solved. Of particular note are recent developments in multigrid methods, which have proven to be the most efficient.     

A linear system of partial differential equations modeling the resting potential of a heart with regional ischemia

Ischemic ST-segment shift has been modeled using scalar stationary approximations of the bidomain model. Here, we propose an alternative simplification of the bidomain equations: a linear system modeling the resting potential, to be used in determining ischemic TP shift. Results of 2D and 3D simulations show that the linear system model is much more accurate than the scalar model. This improved accuracy is important if the model is to be used for solving the inverse problem of determining the size and location of an ischemic region. Furthermore, the model can provide insight into how the resting potential depends on the variations in the extracellular potassium concentration that characterize ischemic regions.     

Modelling the Effect of Gap Junctions on Tissue-Level Cardiac Electrophysiology

When modelling tissue-level cardiac electrophysiology, a continuum approximation to the discrete cell-level equations, known as the bidomain equations, is often used to maintain computational tractability. Analysing the derivation of the bidomain equations allows us to investigate how microstructure, in particular gap junctions that electrically connect cells, affect tissue-level conductivity properties. Using a one-dimensional cable model, we derive a modified form of the bidomain equations that take gap junctions into account, and compare results of simulations using both the discrete and continuum models, finding that the underlying conduction velocity of the action potential ceases to match up between models when gap junctions are introduced at physiologically realistic coupling levels. We show that this effect is magnified by: (i) modelling gap junctions with reduced conductivity; (ii) increasing the conductance of the fast sodium channel; and (iii) an increase in myocyte length. From this, we conclude that the conduction velocity arising from the bidomain equations may not be an accurate representation of the underlying discrete system. In particular, the bidomain equations are less likely to be valid when modelling certain diseased states whose symptoms include a reduction in gap junction coupling or an increase in myocyte length.     

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.     

Wavefront propagation in an activation model of the anisotropic cardiac tissue: asymptotic analysis and numerical simulations

In this paper we present a macroscopic model of the excitation process in the myocardium. The composite and anisotropic structure of the cardiac tissue is represented by a bidomain, i.e. a set of two coupled anisotropic media. The model is characterized by a non linear system of two partial differential equations of parabolic and elliptic type. A singular perturbation analysis is carried out to investigate the cardiac potential field and the structure of the moving excitation wavefront. As a consequence the cardiac current sources are approximated by an oblique dipole layer structure and the motion of the wavefront is described by eikonal equations. Finally numerical simulations are carried out in order to analyze some complex phenomena related to the spreading of the wavefront, like the front-front or front-wall collision. The results yielded by the excitation model and the eikonal equations are compared.     

Virtual electrode polarization in the far field: implications for external defibrillation

We recently suggested that failure of implantable defibrillation therapy may be explained by the virtual electrode-induced phase singularity mechanism. The goal of this study was to identify possible mechanisms of vulnerability and defibrillation by externally applied shocks in vitro. We used bidomain simulations of realistic rabbit heart fibrous geometry to predict the passive polarization throughout the heart induced by external shocks. We also used optical mapping to assess anterior epicardium electrical activity during shocks in Langendorff-perfused rabbit hearts (n = 7). Monophasic shocks of either polarity (10-260 V, 8 ms, 150 microF) were applied during the T wave from a pair of mesh electrodes. Postshock epicardial virtual electrode polarization was observed after all 162 applied shocks, with positive polarization facing the cathode and negative polarization facing the anode, as predicted by the bidomain simulations. During arrhythmogenesis, a new wave front was induced at the boundary between the two regions near the apex but not at the base. It spread across the negatively polarized area toward the base of the heart and reentered on the other side while simultaneously spreading into the depth of the wall. Thus a scroll wave with a ribbon-shaped filament was formed during external shock-induced arrhythmia. Fluorescent imaging and passive bidomain simulations demonstrated that virtual electrode polarization-induced scroll waves underlie mechanisms of shock-induced vulnerability and failure of external defibrillation.     

Mathematical modeling of the excitation process in myocardial tissue: influence of fiber rotation on wavefront propagation and potential field

In our macroscopic model the heart tissue is represented as a bidomain coupling the intra- and extracellular media. Owing to the fiber structure of the myocardium, these media are anisotropic, and their conductivity tensors have a principal axis parallel to the local fiber direction. A reaction-diffusion system is derived that governs the distribution and evolution of the extracellular and transmembrane potentials during the depolarization phase of the heart beat. To investigate frontlike solutions, the system is rescaled and transformed into a system dependent on a small parameter. Subsequently a perturbation analysis is carried out that yields zero- and first-order approximations called eikonal equations. The effects of the transmural fiber rotation on wavefront propagation and the corresponding potential field, elicited by point stimulations, are investigated by means of numerical simulations.     

Probing field-induced tissue polarization using transillumination fluorescent imaging

Despite major successes of biophysical theories in predicting the effects of electrical shocks within the heart, recent optical mapping studies have revealed two major discrepancies between theory and experiment: 1), the presence of negative bulk polarization recorded during strong shocks; and 2), the unexpectedly small surface polarization under shock electrodes. There is little consensus as to whether these differences result from deficiencies of experimental techniques, artifacts of tissue damage, or deficiencies of existing theories. Here, we take advantage of recently developed near-infrared voltage-sensitive dyes and transillumination optical imaging to perform, for the first time that we know of, noninvasive probing of field effects deep inside the intact ventricular wall. This technique removes some of the limitations encountered in previous experimental studies. We explicitly demonstrate that deep inside intact myocardial tissue preparations, strong electrical shocks do produce considerable negative bulk polarization previously inferred from surface recordings. We also demonstrate that near-threshold diastolic field stimulation produces activation of deep myocardial layers 2-6 mm away from the cathodal surface, contrary to theory. Using bidomain simulations we explore factors that may improve the agreement between theory and experiment. We show that the inclusion of negative asymmetric current can qualitatively explain negative bulk polarization in a discontinuous bidomain model.     

Nonlinear simulations of solid tumor growth using a mixture model: invasion and branching

We develop a thermodynamically consistent mixture model for avascular solid tumor growth which takes into account the effects of cell-to-cell adhesion, and taxis inducing chemical and molecular species. The mixture model is well-posed and the governing equations are of Cahn-Hilliard type. When there are only two phases, our asymptotic analysis shows that earlier single-phase models may be recovered as limiting cases of a two-phase model. To solve the governing equations, we develop a numerical algorithm based on an adaptive Cartesian block-structured mesh refinement scheme. A centered-difference approximation is used for the space discretization so that the scheme is second order accurate in space. An implicit discretization in time is used which results in nonlinear equations at implicit time levels. We further employ a gradient stable discretization scheme so that the nonlinear equations are solvable for very large time steps. To solve those equations we use a nonlinear multilevel/multigrid method which is of an optimal order O(N) where N is the number of grid points. Spherically symmetric and fully two dimensional nonlinear numerical simulations are performed. We investigate tumor evolution in nutrient-rich and nutrient-poor tissues. A number of important results have been uncovered. For example, we demonstrate that the tumor may suffer from taxis-driven fingering instabilities which are most dramatic when cell proliferation is low, as predicted by linear stability theory. This is also observed in experiments. This work shows that taxis may play a role in tumor invasion and that when nutrient plays the role of a chemoattractant, the diffusional instability is exacerbated by nutrient gradients. Accordingly, we believe this model is capable of describing complex invasive patterns observed in experiments.     

Soft Tissue Modelling of Cardiac Fibres for Use in Coupled Mechano-Electric Simulations

The numerical solution of the coupled system of partial differential and ordinary differential equations that model the whole heart in three dimensions is a considerable computational challenge. As a consequence, it is not computationally practical—either in terms of memory or time—to repeat simulations on a finer computational mesh to ensure that convergence of the solution has been attained. In an attempt to avoid this problem while retaining mathematical rigour, we derive a one dimensional model of a cardiac fibre that takes account of elasticity properties in three structurally defined axes within the myocardial tissue. This model of a cardiac fibre is then coupled with an electrophysiological cell model and a model of cellular electromechanics to allow us to simulate the coupling of the electrical and mechanical activity of the heart. We demonstrate that currently used numerical methods for coupling electrical and mechanical activity do not work in this case, and identify appropriate numerical techniques that may be used when solving the governing equations. This allows us to perform a series of simulations that: (i) investigate the effect of some of the assumptions inherent in other models; and (ii) reproduce qualitatively some experimental observations.     

Current injection into a two-dimensional anisotropic bidomain.

A two-dimensional sheet of anisotropic cardiac tissue is represented with the bidomain model, and the finite element method is used to solve the bidomain equations. When the anisotropy ratios of the intracellular and extracellular spaces are not equal, the injection of current into the tissue induces a transmembrane potential that has a complicated spatial dependence, including adjacent regions of depolarized and hyperpolarized tissue. This behavior may have important implications for the electrical stimulation of cardiac tissue and for defibrillation.     

Linear algebraic transformations of the bidomain equations: Implications for numerical methods

A mathematical framework is presented for the treatment of the bidomain equations used to model propagation in cardiac tissue. This framework is independent of the model used to represent membrane ionic currents and incorporates boundary conditions and other constraints. By representing the bidomain equations in the operator notation , various algebraic transformations can be expressed as , where P and Q are linear operators. The authors show how previous work fits into this framework and discuss the implications of various transformation for numerical methods of solution. Although such transformations allow many choices of independent variable, these results emphasize the fundamental importance of the transmembrane potential.     

Optimal control approach to termination of re-entry waves in cardiac electrophysiology

This work proposes an optimal control approach for the termination of re-entry waves in cardiac electrophysiology. The control enters as an extracellular current density into the bidomain equations which are well established model equations in the literature to describe the electrical behavior of the cardiac tissue. The optimal control formulation is inspired, in part, by the dynamical systems behavior of the underlying system of differential equations. Existence of optimal controls is established and the optimality system is derived formally. The numerical realization is described in detail and numerical experiments, which demonstrate the capability of influencing and terminating reentry phenomena, are presented.     

The effect of simplifying assumptions in the bidomain model of cardiac tissue: application to ST segment shifts during partial ischaemia

In this study various electrical conductivity approximations used in bidomain models of cardiac tissue are considered. Comparisons are based on epicardial surface potential distributions arising from regions of subendocardial ischaemia situated within the cardiac tissue. Approximations studied are a single conductivity bidomain model, an isotropic bidomain model and equal and reciprocal anisotropy ratios both with and without fibre rotation. It is demonstrated both analytically and numerically that the approximations involving a single conductivity bidomain, an isotropic bidomain or equal anisotropy ratios (ignoring fibre rotation) results in identical epicardial potential distributions for all degrees of subendocardial ischaemia. This result is contrary to experimental observations. It is further shown that by assuming reciprocal anisotropy ratios, epicardial potential distributions vary with the degree of subendocardial ischaemia. However, it is concluded that unequal anisotropy ratios must be used to obtain the true character of experimental observations.     

On the possibility for computing the transmembrane potential in the heart with a one shot method: an inverse problem

We analyze the possibility for using body surface potential maps (BSPMs), a priori information about the voltage distribution in the heart and the bidomain equations to compute the transmembrane potential throughout the myocardium. Our approach is defined in terms of an inverse problem for elliptic partial differential equations (PDEs). More precisely, we formulate it in terms of an output least squares framework in which a goal functional is minimized subject to suitable PDE constraints. The problem is highly unstable and, even under optimal recording conditions, it does not have a unique solution. We propose a methodology for stabilizing and enforcing uniqueness for this inverse problem. Moreover, a fully implicit method for solving the involved minimization problem is presented. In other words, we show how one may solve it in terms of a system consisting of three linear elliptic PDEs, i.e. we derive a so-called one shot method (also commonly referred to as an all-at-once method). Finally, our theoretical findings are illuminated by a series of numerical experiments. These examples indicate that, in the presence of regional ischemia, it might be possible to approximately recover the transmembrane potential during the resting and plateau phases of the heart cycle. This is probably due to the fact that rather accurate a priori information is available during these time intervals. The problem of computing the transmembrane potential at an arbitrary time instance during a heart beat is still an open problem.     

A cylindrical model for studying subendocardial ischaemia in the left ventricle

In this paper a mathematical model of a left ventricle with a cylindrical geometry is presented with the aim of gaining a better understanding of the relationship between subendocardial ischaemia and ST depression. The model is formulated as an infinite cylinder and takes into account the full bidomain nature of cardiac tissue, as well as fibre rotation. A detailed solution method (based on Fourier series, Fourier transforms and a one dimensional finite difference scheme) for the governing equations for electric potential in the tissue and the blood is also presented. The model presented is used to study the effect increasing subendocardial ischaemia has on the epicardial potential distribution as well as the effects of changing the bidomain conductivity values. The epicardial potential distributions obtained with this cylindrical geometry are compared with results obtained using a previously published slab model. Results of the simulations presented show that the morphologies of the epicardial potential distributions are similar between the two geometries, with the main difference being that the cylindrical model predicts slightly higher potentials.     

Electrically silent magnetic fields.

There has been a significant controversy over the past decade regarding the relative information content of bioelectric and biomagnetic signals. In this paper we present a new, theoretical example of an electrically-silent magnetic field, based on a bidomain model of a cylindrical strand of tissue generalized to include off-diagonal components in the conductivity tensors. The physical interpretation of the off-diagonal components is explained, and analytic expressions for the electrical potential and the magnetic field are found. These expressions show that information not obtainable from electrical potential measurements can be obtained from measurements of the magnetic field in systems with conductivity tensors more complicated than those previously examined.     

Fast Simulation of Mechanical Heterogeneity in the Electrically Asynchronous Heart Using the MultiPatch Module

Cardiac electrical asynchrony occurs as a result of cardiac pacing or conduction disorders such as left bundle-branch block (LBBB). Electrically asynchronous activation causes myocardial contraction heterogeneity that can be detrimental for cardiac function. Computational models provide a tool for understanding pathological consequences of dyssynchronous contraction. Simulations of mechanical dyssynchrony within the heart are typically performed using the finite element method, whose computational intensity may present an obstacle to clinical deployment of patient-specific models. We present an alternative based on the CircAdapt lumped-parameter model of the heart and circulatory system, called the MultiPatch module. Cardiac walls are subdivided into an arbitrary number of patches of homogeneous tissue. Tissue properties and activation time can differ between patches. All patches within a wall share a common wall tension and curvature. Consequently, spatial location within the wall is not required to calculate deformation in a patch. We test the hypothesis that activation time is more important than tissue location for determining mechanical deformation in asynchronous hearts. We perform simulations representing an experimental study of myocardial deformation induced by ventricular pacing, and a patient with LBBB and heart failure using endocardial recordings of electrical activation, wall volumes, and end-diastolic volumes. Direct comparison between simulated and experimental strain patterns shows both qualitative and quantitative agreement between model fibre strain and experimental circumferential strain in terms of shortening and rebound stretch during ejection. Local myofibre strain in the patient simulation shows qualitative agreement with circumferential strain patterns observed in the patient using tagged MRI. We conclude that the MultiPatch module produces realistic regional deformation patterns in the asynchronous heart and that activation time is more important than tissue location within a wall for determining myocardial deformation. The CircAdapt model is therefore capable of fast and realistic simulations of dyssynchronous myocardial deformation embedded within the closed-loop cardiovascular system.