Bidomain Equations Research Articles

Effects of deformation on transmural dispersion of repolarization using in silico models of human left ventricular wedge

Mechanical deformation affects the electrical activity of the heart through multiple feedback loops. The purpose of this work is to study the effect of deformation on transmural dispersion of repolarization and on surface electrograms using an in silico human ventricular wedge. To achieve this purpose, we developed a strongly coupled electromechanical cell model by coupling a human left ventricle electrophysiology model and an active contraction model reparameterized for human cells. This model was then embedded in tissue simulations on the basis of bidomain equations and nonlinear solid mechanics. The coupled model was used to evaluate effects of mechanical deformation on important features of repolarization and electrograms. Our results indicate an increase in the T-wave amplitude of the surface electrograms in simulations that account for the effects of cardiac deformation. This increased T-wave amplitude can be explained by changes to the coupling between neighboring myocytes, also known as electrotonic effect. The thickening of the ventricular wall during repolarization contributes to the decoupling of cells in the transmural direction, enhancing action potential heterogeneity and increasing both transmural repolarization dispersion and T-wave amplitude of surface electrograms. The simulations suggest that a considerable percentage of the T-wave amplitude (15%) may be related to cardiac deformation.

Optimal control of the bidomain system (III): Existence of minimizers and first-order optimality conditions

We consider optimal control problems for the bidomain equations of cardiac electrophysiology together with two-variable ionic models, e.g. the Rogers–McCulloch model. After ensuring the existence of global minimizers, we provide a rigorous proof for the system of first-order necessary optimality conditions. The proof is based on a stability estimate for the primal equations and an existence theorem for weak solutions of the adjoint system.

Modeling for cardiac excitation propagation based on the Nernst-Planck equation and homogenization

The bidomain model is a commonly used mathematical model of the electrical properties of the cardiac muscle that takes into account the anisotropy of both the intracellular and extracellular spaces. However, the equations contain self-contradiction such that the update of ion concentrations does not consider intracellular or extracellular ion movements due to the gradient of electric potential and the membrane charge as capacitive currents in spite of the fact that those currents are taken into account in forming Kirchhoff's first law. To overcome this problem, we start with the Nernst-Planck equation, the ionic conservation law, and the electroneutrality condition at the cellular level, and by introducing a homogenization method and assuming uniformity of variables at the microscopic scale, we derive rational bidomain equations at the macroscopic level.

Bidomain simulations of defibrillation: 20 years of progress

Twenty years ago, Natalia Trayanova, her student Lisa Malden, and I published our “spherical heart model.” We showed that an applied electric field produces a large transmembrane potential (nowadays called a “virtual electrode polarization”) at the heart surface, but also polarizes tissue throughout the heart wall caused by fiber curvature. This was the first application of the bidomain model—a mathematical model that describes the anisotropic electrical properties of cardiac tissue—to defibrillation. Since then, Trayanova and her colleagues have developed and refined these bidomain simulations, dramatically improving their realism and sophistication. In this issue of HeartRhythm, Rantner et al take these calculations to a new level of detail. They present bidomain simulations of the rabbit right ventricle having unprecedented resolution, featuring microstructures such as trabeculae and coronary vessels. They solve the bidomain equations using the finite element method with 29 million nodes (in our original calculation, we used 257 nodes), and an average element edge length in the tissue of only 59 μm. The heart is submerged in a perfusing bath, and an electric field is applied through distant electrodes. Rantner et al conclude that “virtual electrode polarizations ... were strongest at the trabecular grooves.” This result implies that much of the interesting behavior occurs on the endocardial surface. Optical mapping studies that image the epicardium might miss the crucial events triggered by the shock. In addition, mathematical models (such as our simple spherical heart) that lack the trabeculae will not predict the same behavior as will a model rich in structural detail. The conclusion that the trabecular length scale is important harkens back to a hypothesis proposed by Krassowska and Kumar, who suggested that excitation may not happen at the cellular scale (the “sawtooth” effect) or the whole heart scale, but at an intermediate scale of fiber bundles a few millimeters in size. The article by Rantner et al does not examine traditional defibrillation with a single strong shock. Instead, the authors apply a series of weak shocks timed just right to first convert the fibrillation to tachycardia and then to terminate the tachycardia. The idea of applying multiple weak stimuli

Chaste

The simulation of cardiac electrophysiology is a mature field in computational physiology. Recent advances in medical imaging, high-performance computing and numerical methods mean that computational models of electrical propagation in human heart tissue are ripe for use in patient-specific simulation for diagnosis, for prognosis and for selection of treatment methods. However, in order to move in this direction, it is necessary to make efficient use of modern petascale computing resources. This paper focuses on an existing open source simulation framework (Chaste) and documents work done to improve the parallel scaling on a small range of electrophysiology benchmark problems. These benchmarks involve the numerical solution of the monodomain or bidomain equations via the finite-element method. At the beginning of this study the electrophysiology libraries within Chaste were already enabled to run in parallel and were able to solve for electrical propagation using the monodomain or bidomain equations, but parallel efficiency dropped rapidly when run on more than about 64 processors. Throughout the course of the study, improvements were made to problem definition input; geometric mesh partitioning; finite-element assembly of large, sparse linear systems; problem-specific matrix preconditioning; numerical solution of the linear system; and output of the approximate solution. The consequence of these improvements is that, at the end of the study, Chaste is able to solve a monodomain benchmark problem in close to real time. While some of the improvements made to the parallel Chaste code are specific to cardiac electrophysiology, many of the techniques documented in this paper are generic to the parallel finite-element method in other scientific application areas.

Block Preconditioners for Coupled Physics Problems

Finite element discretizations of multiphysics problems frequently give rise to block-structured linear algebra problems that require effective preconditioners. We build two classes of preconditioners in the spirit of well-known block factorizations [M. F. Murphy, G. H. Golub, and A. J. Wathen, SIAM J. Sci. Comput., 21 (2000), pp. 1969--1972; I. C. F. Ipsen, SIAM J. Sci. Comput., 23 (2001), pp. 1050--1051] and apply these to the diffusive portion of the bidomain equations and the Bénard convection problem. An abstract generalized eigenvalue problem allows us to give application-specific bounds for the real parts of eigenvalues for these two problems. This analysis is accompanied by numerical calculations with several interesting features. One of our preconditioners for the bidomain equations converges in five iterations for a range of problem sizes. For Bénard convection, we observe mesh-independent convergence with reasonable robustness with respect to physical parameters, and offer some preliminary parallel scaling results on a multicore processor via message passing interface (MPI).

Dynamical modelling of cardiac electrical activity using bidomain approach: The effects of variation of ionic model parameters

This work presents the dynamical modelling of cardiac electrical activity using bidomain approach. It focuses on the effects of variation of the ionic model parameters on cardiac wave propagation. Cardiac electrical activity is governed by partial differential equations coupled to a system of ordinary differential equations. Numerical simulation of these equations is computationally expensive due to their non-linearity and stiffness. Nevertheless, we adopted the bidomain model due to its ability to reflect the actual cardiac wave propagation. The derived bidomain equations coupled with FitzHugh-Nagumo’s ionic equations were time-discretized using explicit forward Euler method and space-discretized using 2-D network modelling to obtain linearized equations for transmembrane potential Vm, extracellular potential φe and gating variable w. We implemented the discretized model and performed simulation experiments to study the effects of variation of ionic model parameters on the propagation of electrical wave across the cardiac tissue. Time characteristic of transmembrane potential, Vm, in the normal cardiac tissue was obtained by setting the values of ionic model parameters to 0.2, 0.2, 0.7 and 0.8 for excitation rate constant ε1, recovery rate constant ε2, recovery decay constant γ and excitation decay constant β respectively. Changing the values of ε1, ε2 to 0.04 and 0.28 respectively, the obtained Vm showed a time dilation at 0.04 indicating cardiac arrhythmia but no significant change to Vm was observed at 0.28. Also, changing β to 0.3 and 1.1 and γ to 0.4 and 1.2 sequentially, there was no significant change to the time characteristic of Vm. The obtained results revealed that only decrease in ε1, ε2 impacted significantly on the cardiac wave propagation.

Similarities and differences between electrocardiogram signs of left bundle-branch block and left-ventricular uncoupling

A left bundle-branch block (LBBB) electrocardiogram (ECG) type may be caused by either a block in the left branch of the ventricular conduction system or by uncoupling in the working myocardium. We used a realistic large-scale computer model to evaluate the effects of uncoupling with and without left-sided block and in combination with biventricular pacing. Action potential propagation was simulated using a reaction-diffusion model of the human ventricles. Electrocardiograms and cardiac electrograms were computed from the simulated action potentials by solving the bidomain equations. In all situations, diffuse uncoupling reduced QRS amplitude, prolonged QRS duration, and rotated the QRS axis leftward. Uncoupling by 50% increased QRS duration from 90 to 120 ms with a normal conduction system and from 140 to 190 ms when the left bundle branch was blocked. Biventricular pacing did not change QRS duration but reduced total ventricular activation time. Uncoupling in the working myocardium can mimic left-sided block in the ventricular conduction system and can explain an LBBB ECG pattern with relatively low amplitude. Biventricular pacing improves ventricular activation in true LBBB with or without uncoupling but not in case of 50% uncoupling alone.

Convergence order vs. parallelism in the numerical simulation of the bidomain equations

The propagation of electrical activity in the human heart can be modelled mathematically by the bidomain equations. The bidomain equations represent a multi-scale reaction-diffusion model that consists of a set of ordinary differential equations governing the dynamics at the cellular level coupled with a set of partial differential equations governing the dynamics at the tissue level. Significant computation is generally required to generate clinically useful data from the bidomain equations. Contemporary developments in computer architecture, in particular multi- and many-core computers and graphics processing units, have made such computations feasible. However, the zeal to take advantage to parallel architectures has typically caused another important aspect of numerical methods for the solution of differential equations to be overlooked, namely the convergence order. It is well known that higher-order methods are generally more efficient than lower-order ones when solutions are smooth and relatively high accuracy is desired. In these situations, serial implementations of high-order methods may remain surprisingly competitive with parallel implementations of low-order methods. In this paper, we examine the effect of order on the numerical solution of the bidomain equations in parallel. We find that high-order methods, in particular high-order time-integration methods with relatively better stability properties, tend to outperform their low-order counterparts, even when the latter are run in parallel. In other words, increasing integration order often trumps increasing available computational resources, especially when relatively high accuracy is desired.

Modelling the effect of gap junctions on tissue-level cardiac electrophysiology

When modelling tissue-level cardiac electrophysiology, continuum approximations to the discrete cell-level equations are used to maintain computational tractability. One of the most commonly used models is represented by the bidomain equations, the derivation of which relies on a homogenisation technique to construct a suitable approximation to the discrete model. This derivation does not explicitly account for the presence of gap junctions connecting one cell to another. It has been seen experimentally [Rohr, Cardiovasc. Res. 2004] that these gap junctions have a marked effect on the propagation of the action potential, specifically as the upstroke of the wave passes through the gap junction. In this paper we explicitly include gap junctions in a both a 2D discrete model of cardiac electrophysiology, and the corresponding continuum model, on a simplified cell geometry. Using these models we compare the results of simulations using both continuum and discrete systems. We see that the form of the action potential as it passes through gap junctions cannot be replicated using a continuum model, and that the underlying propagation speed of the action potential ceases to match up between models when gap junctions are introduced. In addition, the results of the discrete simulations match the characteristics of those shown in Rohr 2004. From this, we suggest that a hybrid model -- a discrete system following the upstroke of the action potential, and a continuum system elsewhere -- may give a more accurate description of cardiac electrophysiology.

Solving the cardiac bidomain equations using graphics processing units

The computational modeling of the heart has been shown to be a very useful tool. The models, which become more realistic each day, provide a better understanding of the complex biophysical processes related to the electrical activity in the heart, e.g., in the case of cardiac arrhythmias. However, the increasing complexity of the models challenges high performance computing in many aspects. This work presents a cardiac simulator based on the bidomain equations that exploits the new parallel architecture of graphics processing units (GPUs). The initial results are promising. The use of the GPU accelerates the cardiac simulator by about 6 times compared to the best performance obtained in a general-purpose processor (CPU). In addition, the GPU implementation was compared to an efficient parallel implementation developed for cluster computing. A single desktop computer equipped with a GPU is shown to be 1.4 times faster than the parallel implementation of the bidomain equations running on a cluster composed of 16 processing cores.

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.

MagnetoHemoDynamics in the aorta and electrocardiograms

This paper addresses a complex multi-physical phenomenon involving cardiac electrophysiology and hemodynamics. The purpose is to model and simulate a phenomenon that has been observed in magnetic resonance imaging machines: in the presence of a strong magnetic field, the T-wave of the electrocardiogram (ECG) gets bigger, which may perturb ECG-gated imaging. This is due to a magnetohydrodynamic (MHD) effect occurring in the aorta. We reproduce this experimental observation through computer simulations on a realistic anatomy, and with a three-compartment model: inductionless MHD equations in the aorta, bi-domain equations in the heart and electrical diffusion in the rest of the body. These compartments are strongly coupled and solved using finite elements. Several benchmark tests are proposed to assess the numerical solutions and the validity of some modeling assumptions. Then, ECGs are simulated for a wide range of magnetic field intensities (from 0 to 20 T).

A Comparative Study of Graph-Based, Eikonal, and Monodomain Simulations for the Estimation of Cardiac Activation Times

The bidomain and monodomain equations are well established as the standard set of equations for the simulation of cardiac electrophysiological behavior. However, the computational cost of detailed bidomain/monodomain simulations limits their applicability in scenarios where a large number of simulations needs to be performed (e.g., parameter estimation). In this study, we present a graph-based method, which relies on point-to-point path finding to estimate activation times for single points in cardiac tissue with minimal computational costs. To validate our approach, activation times are compared to monodomain simulation results for an anatomically based rabbit ventricular model, incorporating realistic fiber orientation and conduction heterogeneities. Differences in activation times between the graph-based method and monodomain results are less than 10% of the total activation time, and computational performance is orders of magnitude faster with the proposed method when calculating activation times at single points. These results suggest that the graph-based method is well suited for estimating activation times when the need for fast performance justifies a limited loss of accuracy.

Reduced-order modeling for cardiac electrophysiology. Application to parameter identification.

A reduced-order model based on proper orthogonal decomposition (POD) is proposed for the bidomain equations of cardiac electrophysiology. Its accuracy is assessed through electrocardiograms in various configurations, including myocardium infarctions and long-time simulations. We show in particular that a restitution curve can efficiently be approximated by this approach. The reduced-order model is then used in an inverse problem solved by an evolutionary algorithm. Some attempts are presented to identify ionic parameters and infarction locations from synthetic electrocardiograms.

Ionic Mechanism for the Formation of Excitable Scroll Wave Filaments

Scroll waves revolving at high frequency in the heart are responsible for fatal arrhythmias. Some mechanisms of arrhythmias require scroll waves to revolve at high frequency around an excitable filament. Recent ventricular cell models including calcium dynamics cannot reproduce this phenomenon. We address this problem revising sodium current kinetics and key phase II inward currents. We perform a nonlinear analysis of the sodium current gathered in canine cardiac myocytes. Despite an extensive data set, our nonlinear analysis shows that several Hodgkin-Huxley formalisms, i.e., a model family, reproduce the voltage clamp data. We incorporate formalisms taken from this family in the latest version of the Luo and Rudy cell model (LRd) and explore the parameter space for scroll wave dynamics. The simulations are performed on a monolayer of cells (3cm x 3cm) and portions (about 2/3) of the left ventricular free wall with realistic representation of the microanatomy. The Bidomain equations are solved at a resolution of 100um in space and 100us in time and are carried out on a supercomputer of the Texas Advanced Computer Center. Our bifurcation analysis shows that sodium current formalisms associated with higher threshold and slightly slower rate of rise during early depolarization may be more realistic. When combined with relatively minor revision of phase II plateau inward currents, the LRd model can produce scroll waves revolving around an excitable filament. The revised model exhibits a rotation period significantly briefer than the original one, and even briefer than the refractory period measured on a plane wave. This mode of excitation allows us to study mechanisms of cardiac death. In conclusion, revision of sodium current kinetics and phase II inward currents of the LRd model allows to realistically reproduce scroll wave revolving around unexcited filament.

A Bidomain Model of the Ventricular Specialized Conduction System of the Heart

An efficient bidomain model of electrical activity in cardiac specialized conduction system fibers is developed and applied to a geometric model of the specialized conduction system and the ventricles. The bidomain model allows the impact of externally applied electric fields on the specialized conduction system to be studied. To model this system, the three-dimensional bidomain equations for a fiber are reduced to one-dimensional equations governing electrical propagation by averaging over the fiber cross section. The one-dimensional equations are coupled to the surrounding three-dimensional extracellular electrical field to allow their use in defibrillation studies. A finite element method, with semi-implicit time stepping, is developed to numerically solve the equations. Current flow through fiber branch points is governed by Kirchhoff's law and imposed directly in a weak formulation of the equations for use with the finite element method. Coupling between intracellular potential in distal system Purkinje fibers and the myocardium at discrete Purkinje-ventricular junctions is modelled using a resistor model. Simulations of electrical activation through a geometric model of the specialized conduction system are presented. The simulations demonstrate physiologically realistic activation sequences under sinus rhythm and the impact of a large externally applied extracellular field on the specialized conduction system.

Optimal control of the bidomain system (I): The monodomain approximation with the Rogers–McCulloch model

For the monodomain approximation of the bidomain equations, a weak solution concept is proposed. We analyze it for the FitzHugh–Nagumo and the Rogers–McCulloch ionic models, obtaining existence and uniqueness theorems. Subsequently, we investigate optimal control problems subject to the monodomain equations. After proving the existence of global minimizers, the system of the first-order necessary optimality conditions is rigorously characterized. For the adjoint system, we prove an existence and regularity theorem as well.

An active strain electromechanical model for cardiac tissue.

We propose a finite element approximation of a system of partial differential equations describing the coupling between the propagation of electrical potential and large deformations of the cardiac tissue. The underlying mathematical model is based on the active strain assumption, in which it is assumed that there is a multiplicative decomposition of the deformation tensor into a passive and active part holds, the latter carrying the information of the electrical potential propagation and anisotropy of the cardiac tissue into the equations of either incompressible or compressible nonlinear elasticity, governing the mechanical response of the biological material. In addition, by changing from a Eulerian to a Lagrangian configuration, the bidomain or monodomain equations modeling the evolution of the electrical propagation exhibit a nonlinear diffusion term. Piecewise quadratic finite elements are employed to approximate the displacements field, whereas for pressure, electrical potentials and ionic variables are approximated by piecewise linear elements. Various numerical tests performed with a parallel finite element code illustrate that the proposed model can capture some important features of the electromechanical coupling and show that our numerical scheme is efficient and accurate.

A sensitivity study of conductivity values in the passive bidomain equation

There is a complex interplay between the four conductivity values used in the bidomain equation and the resulting electric potential distribution in cardiac tissue arising from subendocardial ischaemia. Based on the three commonly used experimentally derived conductivity data sets, a non-dimensional formulation of the passive bidomain equation is derived, which gives rise naturally to several dimensionless conductivity ratios. The data sets are then used to define a parameter space of these ratios, which is studied by considering the correlation coefficients between different epicardial potential distributions. From this study, it is shown that the ratio of the intracellular longitudinal conductivity to the intracellular transverse conductivity is the key parameter in explaining the differences between the epicardial potential distributions observed with these three data sets.