Efficient integration of a realistic two-dimensional cardiac tissue model by domain decomposition

An adaptive hybridizable discontinuous Galerkin approach for cardiac electrophysiology

Cardiac electrophysiology simulations are numerically challenging because of the propagation of a steep electrochemical wave front and thus require discretizations with small mesh sizes to obtain accurate results. In this work, we present an approach based on the hybridizable discontinuous Galerkin method (HDG), which allows an efficient implementation of high-order discretizations into a computational framework. In particular, using the advantage of the discontinuous function space, we present an efficient p-adaptive strategy for accurately tracking the wave front. The HDG allows to reduce the overall degrees of freedom in the final linear system to those only on the element interfaces. Additionally, we propose a rule for a suitable integration accuracy for the ionic current term depending on the polynomial order and the cell model to handle high-order polynomials. Our results show that for the same number of degrees of freedom, coarse high-order elements provide more accurate results than fine low-order elements. Introducing p-adaptivity further reduces computational costs while maintaining accuracy by restricting the use of high-order elements to resolve the wave front. For a patient-specific simulation of a cardiac cycle, p-adaptivity reduces the average number of degrees of freedom by 95% compared to the nonadaptive model. In addition to reducing computational costs, using coarse meshes with our p-adaptive high-order HDG method also simplifies practical aspects of mesh generation and postprocessing.

A patient specific electro-mechanical model of the heart

This paper presents a patient specific deformable heart model that involves the known electrical and mechanical properties of the cardiac cells and tissue. The whole heart model comprises ten Tusscher's ventricular and Nygren's atrial cell models, the anatomical and electrophysiological model descriptions of the atria (introduced by Harrild et al.) and ventricle (given by Winslow et al.), and the mechanical model of the periodical cardiac contraction and resting phenomena proposed by Moireau et al. During the propagation of the depolarization wave, the kinetic, compositional and rotational anisotropy is handled by the tissue, organ and torso model. The applied patient specific parameters were determined by an evolutionary computation method. An intensive parameter reduction was performed using the abstract formulation of the searching space. This patient specific parameter representation enables the adjustment of deformable model parameters in real-time. The validation process was performed using simultaneously measured ECG and ultrasound image records that were compared with simulated signals and shapes using an abstract, parameterized evaluation criterion.

Models of cardiac tissue electrophysiology: progress, challenges and open questions

Models of cardiac tissue electrophysiology are an important component of the Cardiac Physiome Project, which is an international effort to build biophysically based multi-scale mathematical models of the heart. Models of tissue electrophysiology can provide a bridge between electrophysiological cell models at smaller scales, and tissue mechanics, metabolism and blood flow at larger scales. This paper is a critical review of cardiac tissue electrophysiology models, focussing on the micro-structure of cardiac tissue, generic behaviours of action potential propagation, different models of cardiac tissue electrophysiology, the choice of parameter values and tissue geometry, emergent properties in tissue models, numerical techniques and computational issues. We propose a tentative list of information that could be included in published descriptions of tissue electrophysiology models, and used to support interpretation and evaluation of simulation results. We conclude with a discussion of challenges and open questions.

CHASTE: incorporating a novel multi-scale spatial and temporal algorithm into a large-scale open source library

Recent work has described the software engineering and computational infrastructure that has been set up as part of the Cancer, Heart and Soft Tissue Environment (CHASTE) project. CHASTE is an open source software package that currently has heart and cancer modelling functionality. This software has been written using a programming paradigm imported from the commercial sector and has resulted in a code that has been subject to a far more rigorous testing procedure than that is usual in this field. In this paper, we explain how new functionality may be incorporated into CHASTE. Whiteley has developed a numerical algorithm for solving the bidomain equations that uses the multi-scale (MS) nature of the physiology modelled to enhance computational efficiency. Using a simple geometry in two dimensions and a purpose-built code, this algorithm was reported to give an increase in computational efficiency of more than two orders of magnitude. In this paper, we begin by reviewing numerical methods currently in use for solving the bidomain equations, explaining how these methods may be developed to use the MS algorithm discussed above. We then demonstrate the use of this algorithm within the CHASTE framework for solving the monodomain and bidomain equations in a three-dimensional realistic heart geometry. Finally, we discuss how CHASTE may be developed to include new physiological functionality--such as modelling a beating heart and fluid flow in the heart--and how new algorithms aimed at increasing the efficiency of the code may be incorporated.

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.

From mitochondrial ion channels to arrhythmias in the heart: computational techniques to bridge the spatio-temporal scales

Computer simulations of electrical behaviour in the whole ventricles have become commonplace during the last few years. The goals of this article are (i) to review the techniques that are currently employed to model cardiac electrical activity in the heart, discussing the strengths and weaknesses of the various approaches, and (ii) to implement a novel modelling approach, based on physiological reasoning, that lifts some of the restrictions imposed by current state-of-the-art ionic models. To illustrate the latter approach, the present study uses a recently developed ionic model of the ventricular myocyte that incorporates an excitation-contraction coupling and mitochondrial energetics model. A paradigm to bridge the vastly disparate spatial and temporal scales, from subcellular processes to the entire organ, and from sub-microseconds to minutes, is presented. Achieving sufficient computational efficiency is the key to success in the quest to develop multiscale realistic models that are expected to lead to better understanding of the mechanisms of arrhythmia induction following failure at the organelle level, and ultimately to the development of novel therapeutic applications.

Chaste: using agile programming techniques to develop computational biology software

Cardiac modelling is the area of physiome modelling where the available simulation software is perhaps most mature, and it therefore provides an excellent starting point for considering the software requirements for the wider physiome community. In this paper, we will begin by introducing some of the most advanced existing software packages for simulating cardiac electrical activity. We consider the software development methods used in producing codes of this type, and discuss their use of numerical algorithms, relative computational efficiency, usability, robustness and extensibility. We then go on to describe a class of software development methodologies known as test-driven agile methods and argue that such methods are more suitable for scientific software development than the traditional academic approaches. A case study is a project of our own, Cancer, Heart and Soft Tissue Environment, which is a library of computational biology software that began as an experiment in the use of agile programming methods. We present our experiences with a review of our progress thus far, focusing on the advantages and disadvantages of this new approach compared with the development methods used in some existing packages. We conclude by considering whether the likely wider needs of the cardiac modelling community are currently being met and suggest that, in order to respond effectively to changing requirements, it is essential that these codes should be more malleable. Such codes will allow for reliable extensions to include both detailed mathematical models--of the heart and other organs--and more efficient numerical techniques that are currently being developed by many research groups worldwide.

A guide to modelling cardiac electrical activity in anatomically detailed ventricles

One of the most recent trends in cardiac electrophysiology is the development of integrative anatomically accurate models of the heart, which include description of cardiac activity from sub-cellular and cellular level to the level of the whole organ. In order to construct this type of model, a researcher needs to collect a wide range of information from books and journal articles on various aspects of biology, physiology, electrophysiology, numerical mathematics and computer programming. The aim of this methodological article is to survey recent developments in integrative modelling of electrical activity in the ventricles of the heart, and to provide a practical guide to the resources and tools that are available for work in this exciting and challenging area.

Physiology driven adaptivity for the numerical solution of the bidomain equations

Previous work [Whiteley, J. P. IEEE Trans. Biomed. Eng. 53:2139-2147, 2006] derived a stable, semi-implicit numerical scheme for solving the bidomain equations. This scheme allows the timestep used when solving the bidomain equations numerically to be chosen by accuracy considerations rather than stability considerations. In this study we modify this scheme to allow an adaptive numerical solution in both time and space. The spatial mesh size is determined by the gradient of the transmembrane and extracellular potentials while the timestep is determined by the values of: (i) the fast sodium current; and (ii) the calcium release from junctional sarcoplasmic reticulum to myoplasm current. For two-dimensional simulations presented here, combining the numerical algorithm in the paper cited above with the adaptive algorithm presented here leads to an increase in computational efficiency by a factor of around 250 over previous work, together with significantly less computational memory being required. The speedup for three-dimensional simulations is likely to be more impressive.

3D Heart Simulation And Recognition Of Various Events

This paper presents a new way to solve the inverse problem of electrocardiography in terms of heart model parameters. The developed event estimation and recognition method is based on an optimization system of heart model parameters. An ANN-based preliminary ECG analyzer system has been created to reduce the searching space of the optimization algorithm. The optimal model parameters were determined by minimizing the objective functions, as relations of the observed and model-generated body surface ECGs. The final evaluation results, validated by physicians were about 86% correct. Starting from the fact that input ECGs contained various malfunction cases, such as Wolff-Parkinson-White (WPW) syndrome, atrial and ventricular fibrillation, these results suggest that this approach provides a robust inverse solution, circumventing most of the difficulties of the ECG inverse problem.

An efficient numerical technique for the solution of the monodomain and bidomain equations

Most numerical schemes for solving the monodomain or bidomain equations use a forward approximation to some or all of the time derivatives. This approach, however, constrains the maximum timestep that may be used by stability considerations as well as accuracy considerations. Stability may be ensured by using a backward approximation to all time derivatives, although this approach requires the solution of a very large system of nonlinear equations at each timestep which is computationally prohibitive. In this paper we propose a semi-implicit algorithm that ensures stability. A linear system is solved on each timestep to update the transmembrane potential and, if the bidomain equations are being used, the extracellular potential. The remainder of the equations to be solved uncouple into small systems of ordinary differential equations. The backward Euler method may be used to solve these systems and guarantee numerical stability: as these systems are small, only the solution of small nonlinear systems are required. Simulations are carried out to show that the use of this algorithm allows much larger timesteps to be used with only a minimal loss of accuracy. As a result of using these longer timesteps the computation time may be reduced substantially.

Asynchronous adaptive time step in quantitative cellular automata modeling

BACKGROUND

The behaviors of cells in metazoans are context dependent, thus large-scale multi-cellular modeling is often necessary, for which cellular automata are natural candidates. Two related issues are involved in cellular automata based multi-cellular modeling: how to introduce differential equation based quantitative computing to precisely describe cellular activity, and upon it, how to solve the heavy time consumption issue in simulation.

RESULTS

Based on a modified, language based cellular automata system we extended that allows ordinary differential equations in models, we introduce a method implementing asynchronous adaptive time step in simulation that can considerably improve efficiency yet without a significant sacrifice of accuracy. An average speedup rate of 4-5 is achieved in the given example.

CONCLUSIONS

Strategies for reducing time consumption in simulation are indispensable for large-scale, quantitative multi-cellular models, because even a small 100 x 100 x 100 tissue slab contains one million cells. Distributed and adaptive time step is a practical solution in cellular automata environment.

Three-dimensional pseudospectral modelling of cardiac propagation in an inhomogeneous anisotropic tissue

Various investigators have used the monodomain model to study cardiac propagation behaviour. In many cases, the governing non-linear parabolic equation is solved using the finite-difference method. An adequate discretisation of cardiac tissue with realistic dimensions, however, often leads to a large model size that is computationally demanding. Recently, it has been demonstrated, for a two-dimensional homogeneous monodomain, that the Chebyshev pseudospectral method can offer higher computational efficiency than the finite-difference technique. Here, an extension of the pseudospectral approach to a three-dimensional inhomogeneous case with fibre rotation is presented. The unknown transmembrane potential is expanded in terms of Chebyshev polynomial trial functions, and the monodomain equation is enforced at the Gauss-Lobatto node points. The forward Euler technique is used to advance the solution in time. Numerical results are presented that demonstrate that the Chebyshev pseudospectral method offered an even larger improvement in computational performance over the finite-difference method in the three-dimensional case. Specifically, the pseudospectral method allowed the number of nodes to be reduced by approximately 85 times, while the same solution accuracy was maintained. Depending on the model size, simulations were performed with approximately 18-41 times less memory and approximately 99-169 times less CPU time.

Efficient simulation of three-dimensional anisotropic cardiac tissue using an adaptive mesh refinement method

A recently developed space-time adaptive mesh refinement algorithm (AMRA) for simulating isotropic one- and two-dimensional excitable media is generalized to simulate three-dimensional anisotropic media. The accuracy and efficiency of the algorithm is investigated for anisotropic and inhomogeneous 2D and 3D domains using the Luo-Rudy 1 (LR1) and FitzHugh-Nagumo models. For a propagating wave in a 3D slab of tissue with LR1 membrane kinetics and rotational anisotropy comparable to that found in the human heart, factors of 50 and 30 are found, respectively, for the speedup and for the savings in memory compared to an algorithm using a uniform space-time mesh at the finest resolution of the AMRA method. For anisotropic 2D and 3D media, we find no reduction in accuracy compared to a uniform space-time mesh. These results suggest that the AMRA will be able to simulate the 3D electrical dynamics of canine ventricles quantitatively for 1 s using 32 1-GHz Alpha processors in approximately 9 h.

Computational techniques for solving the bidomain equations in three dimensions

The bidomain equations are the most complete description of cardiac electrical activity. Their numerical solution is, however, computationally demanding, especially in three dimensions, because of the fine temporal and spatial sampling required. This paper methodically examines computational performance when solving the bidomain equations. Several techniques to speed up this computation are examined in this paper. The first step was to recast the equations into a parabolic part and an elliptic part. The parabolic part was solved by either the finite-element method (FEM) or the interconnected cable model model (ICCM). The elliptic equation was solved by FEM on a coarser grid than the parabolic problem and at a reduced frequency. The performance of iterative and direct linear equation system solvers was analyzed as well as the scalability and parallelizability of each method. Results indicate that the ICCM was twice as fast as the FEM for solving the parabolic problem, but when the total problem was considered, this resulted in only a 20% decrease in computation time. The elliptic problem could be solved on a coarser grid at one-quarter of the frequency at which the parabolic problem was solved and still maintain reasonable accuracy. Direct methods were faster than iterative methods by at least 50% when a good estimate of the extracellular potential was required. Parallelization over four processors was efficient only when the model comprised at least 500,000 nodes. Thus, it was possible to speed up solution of the bidomain equations by an order of magnitude with a slight decrease in accuracy.

A computer model of human atria with reasonable computation load and realistic anatomical properties

Atrial fibrillation is the most frequent arrhythmia, provoking discomfort, heart failure and arterial embolisms. The aim of this work is to develop a simplified anatomical computer model of human atria for the study of atrial arrhythmias and the understanding of electrical propagation mechanisms. With the model we propose, up to 40 s of real-time propagation have been simulated on a single-processor computer. The size and the electrophysiological properties of the simulated atria are within realistic values and information about anatomy has been taken into account in a three-dimensional structure. Besides normal sinus beat, pathological phenomena such as flutter and fibrillation have been induced using a programmed stimulation protocol. One important observation in our model is that atrial arrhythmias are a combination of functional and anatomical reentries and that the geometry plays an important role. This virtual atrium can reproduce electrophysiological observations made in humans but with the advantage of showing in great detail how arrhythmias are initiated and sustained. Such details are difficult or impossible to study in humans. This model will serve us as a tool to evaluate the impact of new therapeutic strategies and to improve them.

A theoretical analysis of acute ischemia and infarction using ECG reconstruction on a 2-D model of myocardium

We developed a two-dimensional ventricular tissue model in order to probe the determinants of electrocardiographic (ECG) morphology during acute and chronic ischemia. Hyperkalemia was simulated by step changes in [K+]out, while acidosis was induced by reducing Na+ and Ca2+ conductances. Hypoxia was introduced by its effect on potassium activity. During the initial moments of ischemia, ECG changes were characterized by increases in QRS amplitude and ST segment shortening, followed in the advanced phase by ST baseline elevation, T conformation changes, widening of the QRS and significant decreases in QRS amplitude in spite of an enlarged Q. During each phase, potential proarrhythmic mechanisms were investigated. The presence of unexcitable regions of simulated myocardial infarction led to polymorphic ECG. We also observed a nonuniform deflection of the ST segment from beat to beat. We used similar protocols to explore the responses of infarcted myocardium after impairment resolving. We found that despite irreversible uncoupling of the necrotic region, the restored normal ionic concentrations produced an isopotential ST segment and monomorphic ECG complexes, while an enlarged Q wave was still visible. In summary, these numerical experiments indicate the possibility to track in the ECG pathologic changes following the altered electrophysiology of the ischemic heart.

Asynchronous algorithm for integration of reaction-diffusion equations for inhomogeneous excitable media

An asynchronous algorithm for the integration of reaction-diffusion equations for inhomogeneous excitable media is described. Since many physical systems are inhomogeneous where either the local kinetics or the diffusion or conduction properties vary significantly in space, integration schemes must be able to account for wide variations in the temporal and spatial scales of the solutions. The asynchronous algorithm utilizes a fixed spatial grid and automatically adjusts the time step locally to achieve an efficient simulation where the errors in the solution are controlled. The scheme does not depend on the specific form of the local kinetics and is easily applied to systems with complex geometries. (c) 2000 American Institute of Physics.

Distributed computing for membrane-based modeling of action potential propagation

Action potential propagation simulations with physiologic membrane currents and macroscopic tissue dimensions are computationally expensive. We, therefore, analyzed distributed computing schemes to reduce execution time in workstation clusters by parallelizing solutions with message passing. Four schemes were considered in two-dimensional monodomain simulations with the Beeler-Reuter membrane equations. Parallel speedups measured with each scheme were compared to theoretical speedups, recognizing the relationship between speedup and code portions that executed serially. A data decomposition scheme based on total ionic current provided the best performance. Analysis of communication latencies in that scheme led to a load-balancing algorithm in which measured speedups at 89 +/- 2% and 75 +/- 8% of theoretical speedups were achieved in homogeneous and heterogeneous clusters of workstations. Speedups in this scheme with the Luo-Rudy dynamic membrane equations exceeded 3.0 with eight distributed workstations. Cluster speedups were comparable to those measured during parallel execution on a shared memory machine.

Scroll wave dynamics in a three-dimensional cardiac tissue model: roles of restitution, thickness, and fiber rotation

Scroll wave (vortex) breakup is hypothesized to underlie ventricular fibrillation, the leading cause of sudden cardiac death. We simulated scroll wave behaviors in a three-dimensional cardiac tissue model, using phase I of the Luo-Rudy (LR1) action potential model. The effects of action potential duration (APD) restitution, tissue thickness, filament twist, and fiber rotation were studied. We found that APD restitution is the major determinant of scroll wave behavior and that instabilities arising from APD restitution are the main determinants of scroll wave breakup in this cardiac model. We did not see a "thickness-induced instability" in the LR1 model, but a minimum thickness is required for scroll breakup in the presence of fiber rotation. The major effect of fiber rotation is to maintain twist in a scroll wave, promoting filament bending and thus scroll breakup. In addition, fiber rotation induces curvature in the scroll wave, which weakens conduction and further facilitates wave break.

A space-time adaptive method for simulating complex cardiac dynamics

For plane-wave and many-spiral states of the experimentally based Luo-Rudy 1 model of heart tissue in large (8 cm square) domains, we show that a space-time-adaptive time-integration algorithm can achieve a factor of 5 reduction in computational effort and memory-but without a reduction in accuracy-when compared to an algorithm using a uniform space-time mesh at the finest resolution. Our results indicate that such an algorithm can be extended straightforwardly to simulate quantitatively three-dimensional electrical dynamics over the whole human heart.

A computer model of cardiac electrical activity for the simulation of arrhythmias

Modern computer power allows development of models of the heart that may be helpful for the understanding of arrhythmia mechanisms if, based on realistic physiological parameters, such models can display phenomena difficult to study in nature. Therefore, a two-dimensional model of the cardiac tissue has been implemented, where the modeling of each cell is based on membrane ionic channels (Beeler-Reuter and Luo-Rudy models). In addition, an ECG was computed based on the ionic currents simulated. This model allows us to observe the propagation of the action potentials Vm across the cardiac tissue, the evolution of Vm for any of the cardiac cells, and the underlying ionic currents. The computation of the ECG makes it possible to relate this information with an often-used diagnostic tool. Simulations of normal and pathological phenomena such as functional and anatomic reentry have been performed. Our simulation results show that the applied computer model based on ionic currents seems accurate and realistic when compared with biological models and offers a new approach to study the origin, prevention, and termination of arrhythmias.