Accelerating cardiac bidomain simulations using graphics processing units

Simulation-free prediction of atrial fibrillation inducibility with the fibrotic kernel signature

Computational models of atrial fibrillation (AF) can help improve success rates of interventions, such as ablation. However, evaluating the efficacy of different treatments requires performing multiple costly simulations by pacing at different points and checking whether AF has been induced or not, hindering the clinical application of these models. In this work, we propose a classification method that can predict AF inducibility in patient-specific cardiac models without running additional simulations. Our methodology does not require re-training when changing atrial anatomy or fibrotic patterns. To achieve this, we develop a set of features given by a variant of the heat kernel signature that incorporates fibrotic pattern information and fiber orientations: the fibrotic kernel signature (FKS). The FKS is faster to compute than a single AF simulation, and when paired with machine learning classifiers, it can predict AF inducibility in the entire domain. To learn the relationship between the FKS and AF inducibility, we performed 2371 AF simulations comprising 6 different anatomies and various fibrotic patterns, which we split into training and a testing set. We obtain a median F1 score of 85.2% in test set and we can predict the overall inducibility with a mean absolute error of 2.76 percent points, which is lower than alternative methods. We think our method can significantly speed-up the calculations of AF inducibility, which is crucial to optimize therapies for AF within clinical timelines. An example of the FKS for an open source model is provided in https://github.com/tbanduc/FKS_AtrialModel_Ferrer.git.

An automated near-real time computational method for induction and treatment of scar-related ventricular tachycardias

Catheter ablation is currently the only curative treatment for scar-related ventricular tachycardias (VTs). However, not only are ablation procedures long, with relatively high risk, but success rates are punitively low, with frequent VT recurrence. Personalized in-silico approaches have the opportunity to address these limitations. However, state-of-the-art reaction diffusion (R-D) simulations of VT induction and subsequent circuits used for in-silico ablation target identification require long execution times, along with vast computational resources, which are incompatible with the clinical workflow. Here, we present the Virtual Induction and Treatment of Arrhythmias (VITA), a novel, rapid and fully automated computational approach that uses reaction-Eikonal methodology to induce VT and identify subsequent ablation targets. The rationale for VITA is based on finding isosurfaces associated with an activation wavefront that splits in the ventricles due to the presence of an isolated isthmus of conduction within the scar; once identified, each isthmus may be assessed for their vulnerability to sustain a reentrant circuit, and the corresponding exit site automatically identified for potential ablation targeting. VITA was tested on a virtual cohort of 7 post-infarcted porcine hearts and the results compared to R-D simulations. Using only a standard desktop machine, VITA could detect all scar-related VTs, simulating activation time maps and ECGs (for clinical comparison) as well as computing ablation targets in 48 minutes. The comparable VTs probed by the R-D simulations took 68.5 hours on 256 cores of high-performance computing infrastructure. The set of lesions computed by VITA was shown to render the ventricular model VT-free. VITA could be used in near real-time as a complementary modality aiding in clinical decision-making in the treatment of post-infarction VTs.

Resource-Efficient Use of Modern Processor Architectures For Numerically Solving Cardiac Ionic Cell Models

A central component in simulating cardiac electrophysiology is the numerical solution of nonlinear ordinary differential equations, also called cardiac ionic cell models, that describe cross-cell-membrane ion transport. Biophysically detailed cell models often require a considerable amount of computation, including calls to special mathematical functions. This paper systematically studies how to efficiently use modern multicore CPUs for this costly computational task. We start by investigating the code restructurings needed to effectively enable compiler-supported SIMD vectorisation, which is the most important performance booster in this context. It is found that suitable OpenMP directives are sufficient for achieving both vectorisation and parallelisation. We then continue with an evaluation of the performance optimisation technique of using lookup tables. Due to increased challenges for automated vectorisation, the obtainable benefits of lookup tables are dependent on the hardware platforms chosen. Throughout the study, we report detailed time measurements obtained on Intel Xeon, Xeon Phi, AMD Epyc and two ARM processors including Fujitsu A64FX, while attention is also paid to the impact of SIMD vectorisation and lookup tables on the computational accuracy. As a realistic example, the benefits of performance enhancement are demonstrated by a 109-run ensemble on the Oakforest-PACS system, where code restructurings and SIMD vectorisation yield an 84% reduction in computing time, corresponding to 63,270 node hours.

An Integrated Workflow for Building Digital Twins of Cardiac Electromechanics-A Multi-Fidelity Approach for Personalising Active Mechanics

Personalised computer models of cardiac function, referred to as cardiac digital twins, are envisioned to play an important role in clinical precision therapies of cardiovascular diseases. A major obstacle hampering clinical translation involves the significant computational costs involved in the personalisation of biophysically detailed mechanistic models that require the identification of high-dimensional parameter vectors. An important aspect to identify in electromechanics (EM) models are active mechanics parameters that govern cardiac contraction and relaxation. In this study, we present a novel, fully automated, and efficient approach for personalising biophysically detailed active mechanics models using a two-step multi-fidelity solution. In the first step, active mechanical behaviour in a given 3D EM model is represented by a purely phenomenological, low-fidelity model, which is personalised at the organ scale by calibration to clinical cavity pressure data. Then, in the second step, median traces of nodal cellular active stress, intracellular calcium concentration, and fibre stretch are generated and utilised to personalise the desired high-fidelity model at the cellular scale using a 0D model of cardiac EM. Our novel approach was tested on a cohort of seven human left ventricular (LV) EM models, created from patients treated for aortic coarctation (CoA). Goodness of fit, computational cost, and robustness of the algorithm against uncertainty in the clinical data and variations of initial guesses were evaluated. We demonstrate that our multi-fidelity approach facilitates the personalisation of a biophysically detailed active stress model within only a few (2 to 4) expensive 3D organ-scale simulations-a computational effort compatible with clinical model applications.

The openCARP simulation environment for cardiac electrophysiology

BACKGROUND AND OBJECTIVE

Cardiac electrophysiology is a medical specialty with a long and rich tradition of computational modeling. Nevertheless, no community standard for cardiac electrophysiology simulation software has evolved yet. Here, we present the openCARP simulation environment as one solution that could foster the needs of large parts of this community.

METHODS AND RESULTS

openCARP and the Python-based carputils framework allow developing and sharing simulation pipelines which automate in silico experiments including all modeling and simulation steps to increase reproducibility and productivity. The continuously expanding openCARP user community is supported by tailored infrastructure. Documentation and training material facilitate access to this complementary research tool for new users. After a brief historic review, this paper summarizes requirements for a high-usability electrophysiology simulator and describes how openCARP fulfills them. We introduce the openCARP modeling workflow in a multi-scale example of atrial fibrillation simulations on single cell, tissue, organ and body level and finally outline future development potential.

CONCLUSION

As an open simulator, openCARP can advance the computational cardiac electrophysiology field by making state-of-the-art simulations accessible. In combination with the carputils framework, it offers a tailored software solution for the scientific community and contributes towards increasing use, transparency, standardization and reproducibility of in silico experiments.

A memory optimization method combined with adaptive time-step method for cardiac cell simulation based on multi-GPU

Cardiac electrophysiological simulation is a very complex computational process, which can be run on graphics processing unit (GPU) to save computational cost greatly. The use of adaptive time-step can further effectively speed up the simulation of heart cells. However, if the adaptive time-step method applies to GPU, it suffers synchronization problem on GPU, weakening the acceleration of adaptive time-step method. The previous work ran on a single GPU with the adaptive time-step to get only 1.5 times (× 1.5) faster than the fixed time-step. This study proposes a memory allocation method, which can effectively implement the adaptive time-step method on GPU. The proposed method mainly focuses on the stimulus point and potential memory arrangement in order to achieve optimal memory storage efficiency. All calculation is implemented on GPU. Large matrices such as potential are arranged in column order, and the cells on the left are stimulated. The Luo-Rudy passive (LR1) and dynamic (LRd) ventricular action potential models are used with adaptive time-step methods, such as the traditional hybrid method (THM) and Chen-Chen-Luo's (CCL) "quadratic adaptive algorithm" method. As LR1 is solved by the THM or CCL on a single GPU, the acceleration is × 34 and × 75 respectively compared with the fixed time-step. With 2 or 4 GPUs, the acceleration of the THM and CCL is × 34 or × 35 and × 73 or × 75, but it would decrease to × 5 or × 3 and × 20 or × 15 without optimization. In an LRd model, the acceleration reaches × 27 or × 85 as solved by the THM or CCL compared with the fixed time-step on multi-GPU with linear speed up increase versus the number of GPU. However, with the increase of GPUs number, the acceleration of the THM and CCL is continuously weakened before optimization. The mixed root mean square error (MRMSE) lower than 5% is applied to ensure the accuracy of simulation. The result shows that the proposed memory arrangement method can save computational cost a lot to speed up the heart simulation greatly. Graphical abstract Acceleration ratio compared with CPU with fixed time-step (dt = 0.001 ms).

A Lead Field Two-Domain Model for Longitudinal Neural Tracts-Analytical Framework and Implications for Signal Bandwidth

Somatosensory evoked potentials are a well-established tool for assessing volley conduction in afferent neural pathways. However, from a clinical perspective, recording of spinal signals is still a demanding task due to the low amplitudes compared to relevant noise sources. Computer modeling is a powerful tool for gaining insight into signal genesis and, thus, for promoting future innovations in signal extraction. However, due to the complex structure of neural pathways, modeling is computationally demanding. We present a theoretical framework which allows computing the electric potential generated by a single axon in a body surface lead by the convolution of the neural lead field function with a propagating action potential term. The signal generated by a large cohort of axons was obtained by convoluting a single axonal signal with the statistical distribution of temporal dispersion of individual axonal signals. For establishing the framework, analysis was based on an analytical model. Our approach was further adopted for a numerical computation of body surface neuropotentials employing the lead field theory. Double convolution allowed straightforward analysis in the frequency domain. The highest frequency components occurred at the cellular membrane. A bandpass type spectral shape and a peak frequency of 1800 Hz was observed. The volume conductor transmitting the signal to the recording lead acted as an additional bandpass reducing the axonal peak frequency from 200 Hz to 500 Hz. The superposition of temporally dispersed axonal signals acted as an additional low-pass filter further reducing the compound action potential peak frequency from 90 Hz to 170 Hz. Our results suggest that the bandwidth of spinal evoked potentials might be narrower than the bandwidth requested by current clinical guidelines. The present findings will allow the optimization of noise suppression. Furthermore, our theoretical framework allows the adaptation in numerical methods and application in anatomically realistic geometries in future studies.

Personalization of electro-mechanical models of the pressure-overloaded left ventricle: fitting of Windkessel-type afterload models

Computer models of left ventricular (LV) electro-mechanics (EM) show promise as a tool for assessing the impact of increased afterload upon LV performance. However, the identification of unique afterload model parameters and the personalization of EM LV models remains challenging due to significant clinical input uncertainties. Here, we personalized a virtual cohort of N = 17 EM LV models under pressure overload conditions. A global-local optimizer was developed to uniquely identify parameters of a three-element Windkessel (Wk3) afterload model. The sensitivity of Wk3 parameters to input uncertainty and of the EM LV model to Wk3 parameter uncertainty was analysed. The optimizer uniquely identified Wk3 parameters, and outputs of the personalized EM LV models showed close agreement with clinical data in all cases. Sensitivity analysis revealed a strong dependence of Wk3 parameters on input uncertainty. However, this had limited impact on outputs of EM LV models. A unique identification of Wk3 parameters from clinical data appears feasible, but it is sensitive to input uncertainty, thus depending on accurate invasive measurements. By contrast, the EM LV model outputs were less sensitive, with errors of less than 8.14% for input data errors of 10%, which is within the bounds of clinical data uncertainty. This article is part of the theme issue 'Uncertainty quantification in cardiac and cardiovascular modelling and simulation'.

The impact of wall thickness and curvature on wall stress in patient-specific electromechanical models of the left atrium

The left atrium (LA) has a complex anatomy with heterogeneous wall thickness and curvature. The anatomy plays an important role in determining local wall stress; however, the relative contribution of wall thickness and curvature in determining wall stress in the LA is unknown. We have developed electromechanical finite element (FE) models of the LA using patient-specific anatomical FE meshes with rule-based myofiber directions. The models of the LA were passively inflated to 10mmHg followed by simulation of the contraction phase of the atrial cardiac cycle. The FE models predicted maximum LA volumes of 156.5 mL, 99.3 mL and 83.4 mL and ejection fractions of 36.9%, 32.0% and 25.2%. The median wall thickness in the 3 cases was calculated as [Formula: see text] mm, [Formula: see text] mm, and [Formula: see text] mm. The median curvature was determined as [Formula: see text] [Formula: see text], [Formula: see text], and [Formula: see text]. Following passive inflation, the correlation of wall stress with the inverse of wall thickness and curvature was 0.55-0.62 and 0.20-0.25, respectively. At peak contraction, the correlation of wall stress with the inverse of wall thickness and curvature was 0.38-0.44 and 0.16-0.34, respectively. In the LA, the 1st principal Cauchy stress is more dependent on wall thickness than curvature during passive inflation and both correlations decrease during active contraction. This emphasizes the importance of including the heterogeneous wall thickness in electromechanical FE simulations of the LA. Overall, simulation results and sensitivity analyses show that in complex atrial anatomy it is unlikely that a simple anatomical-based law can be used to estimate local wall stress, demonstrating the importance of FE analyses.

Delay differential equation-based models of cardiac tissue: Efficient implementation and effects on spiral-wave dynamics

Delay differential equations (DDEs) recently have been used in models of cardiac electrophysiology, particularly in studies focusing on electrical alternans, instabilities, and chaos. A number of processes within cardiac cells involve delays, and DDEs can potentially represent mechanisms that result in complex dynamics both at the cellular level and at the tissue level, including cardiac arrhythmias. However, DDE-based formulations introduce new computational challenges due to the need for storing and retrieving past values of variables at each spatial location. Cardiac tissue simulations that use DDEs may require over 28 GB of memory if the history of variables is not managed carefully. This paper addresses both computational and dynamical issues. First, we present new methods for the numerical solution of DDEs in tissue to mitigate the memory requirements associated with the history of variables. The new methods exploit the different time scales of an action potential to dynamically optimize history size. We find that the proposed methods decrease memory usage by up to 95% in cardiac tissue simulations compared to straightforward history-management algorithms. Second, we use the optimized methods to analyze for the first time the dynamics of wave propagation in two-dimensional cardiac tissue for models that include DDEs. In particular, we study the effects of DDEs on spiral-wave dynamics, including wave breakup and chaos, using a canine myocyte model. We find that by introducing delays to the gating variables governing the calcium current, DDEs can induce spiral-wave breakup in 2D cardiac tissue domains.

Factors Promoting Conduction Slowing as Substrates for Block and Reentry in Infarcted Hearts

The development of effective and safe therapies for scar-related ventricular tachycardias requires a detailed understanding of the mechanisms underlying the conduction block that initiates electrical re-entries associated with these arrhythmias. Conduction block has been often associated with electrophysiological changes that prolong action potential duration (APD) within the border zone (BZ) of chronically infarcted hearts. However, experimental evidence suggests that remodeling processes promoting conduction slowing as opposed to APD prolongation mark the chronic phase. In this context, the substrate for the initial block at the mouth of an isthmus/diastolic channel leading to ventricular tachycardia is unclear. The goal of this study was to determine whether electrophysiological parameters associated with conduction slowing can cause block and re-entry in the BZ. In silico experiments were conducted on two-dimensional idealized infarct tissue as well as on a cohort of postinfarction porcine left ventricular models constructed from ex vivo magnetic resonance imaging scans. Functional conduction slowing in the BZ was modeled by reducing sodium current density, whereas structural conduction slowing was represented by decreasing tissue conductivity and including fibrosis. The arrhythmogenic potential of APD prolongation was also tested as a basis for comparison. Within all models, the combination of reduced sodium current with structural remodeling more often degenerated into re-entry and, if so, was more likely to be sustained for more cycles. Although re-entries were also detected in experiments with prolonged APD, they were often not sustained because of the subsequent block caused by long-lasting repolarization. Functional and structural conditions associated with slow conduction rather than APD prolongation form a potent substrate for arrhythmogenesis at the isthmus/BZ of chronically infarcted hearts. Reduced excitability led to block while slow conduction shortened the wavelength of propagation, facilitating the sustenance of re-entries. These findings provide important insights for models of patient-specific risk stratification and therapy planning.

Large-scale Interactive Numerical Experiments of Chaos, Solitons and Fractals in Real Time via GPU in a Web Browser

The study of complex systems has emerged as an important field with many discoveries still to be made. Computer simulation and visualization provide important tools for studying complex dynamics including chaos, solitons, and fractals, but available computing power has been a limiting factor. In this work, we describe a novel and highly efficient computing and visualization paradigm using a Web Graphics Library (WebGL 2.0) methodology along with our newly developed library (Abubu.js). Our approach harnesses the power of widely available and highly parallel graphics cards while maintaining ease of use by simplifying programming through hiding implementation details, running in a web browser without the need for compilation, and avoiding the use of plugins. At the same time, it allows for interactivity, such as changing parameter values on the fly, and its computing is so fast that zooming in on a region of a fractal like the Mandelbrot set can incur no delay despite having to recalculate values for the entire plane. We demonstrate our approach using a wide range of complex systems that display dynamics from fractals to standing and propagating waves in 1, 2 and 3 dimensions. We also include some models with instabilities that can lead to chaotic dynamics. For all the examples shown here we provide links to the codes for anyone to use, modify and further develop with other models. Overall, the enhanced visualization and computation capabilities provided by WebGL together with Abubu.js have great potential to facilitate new discoveries about complex systems.

Assessment of wall stresses and mechanical heart power in the left ventricle: Finite element modeling versus Laplace analysis

INTRODUCTION

Stenotic aortic valve disease (AS) causes pressure overload of the left ventricle (LV) that may trigger adverse remodeling and precipitate progression towards heart failure (HF). As myocardial energetics can be impaired during AS, LV wall stresses and biomechanical power provide a complementary view of LV performance that may aide in better assessing the state of disease.

OBJECTIVES

Using a high-resolution electro-mechanical (EM) in silico model of the LV as a reference, we evaluated clinically feasible Laplace-based methods for assessing global LV wall stresses and biomechanical power.

METHODS

We used N = 4 in silico finite element (FE) EM models of LV and aorta of patients suffering from AS. All models were personalized with clinical data under pretreatment conditions. Left ventricle wall stresses and biomechanical power were computed accurately from FE kinematic data and compared with Laplace-based estimation methods, which were applied to the same FE model data.

RESULTS AND CONCLUSION

Laplace estimates of LV wall stress are able to provide a rough approximation of global mean stress in the circumferential-longitudinal plane of the LV. However, according to FE results, spatial heterogeneity of stresses in the LV wall is significant, leading to major discrepancies between local stresses and global mean stress. Assessment of mechanical power with Laplace methods is feasible, but these are inferior in accuracy compared with FE models. The accurate assessment of stress and power density distribution in the LV wall is only feasible based on patient-specific FE modeling.

Towards a Computational Framework for Modeling the Impact of Aortic Coarctations Upon Left Ventricular Load

Computational fluid dynamics (CFD) models of blood flow in the left ventricle (LV) and aorta are important tools for analyzing the mechanistic links between myocardial deformation and flow patterns. Typically, the use of image-based kinematic CFD models prevails in applications such as predicting the acute response to interventions which alter LV afterload conditions. However, such models are limited in their ability to analyze any impacts upon LV load or key biomarkers known to be implicated in driving remodeling processes as LV function is not accounted for in a mechanistic sense. This study addresses these limitations by reporting on progress made toward a novel electro-mechano-fluidic (EMF) model that represents the entire physics of LV electromechanics (EM) based on first principles. A biophysically detailed finite element (FE) model of LV EM was coupled with a FE-based CFD solver for moving domains using an arbitrary Eulerian-Lagrangian (ALE) formulation. Two clinical cases of patients suffering from aortic coarctations (CoA) were built and parameterized based on clinical data under pre-treatment conditions. For one patient case simulations under post-treatment conditions after geometric repair of CoA by a virtual stenting procedure were compared against pre-treatment results. Numerical stability of the approach was demonstrated by analyzing mesh quality and solver performance under the significantly large deformations of the LV blood pool. Further, computational tractability and compatibility with clinical time scales were investigated by performing strong scaling benchmarks up to 1536 compute cores. The overall cost of the entire workflow for building, fitting and executing EMF simulations was comparable to those reported for image-based kinematic models, suggesting that EMF models show potential of evolving into a viable clinical research tool.

Scalable and Accurate ECG Simulation for Reaction-Diffusion Models of the Human Heart

Realistic electrocardiogram (ECG) simulation with numerical models is important for research linking cellular and molecular physiology to clinically observable signals, and crucial for patient tailoring of numerical heart models. However, ECG simulation with a realistic torso model is computationally much harder than simulation of cardiac activity itself, so that many studies with sophisticated heart models have resorted to crude approximations of the ECG. This paper shows how the classical concept of electrocardiographic lead fields can be used for an ECG simulation method that matches the realism of modern heart models. The accuracy and resource requirements were compared to those of a full-torso solution for the potential and scaling was tested up to 14,336 cores with a heart model consisting of 11 million nodes. Reference ECGs were computed on a 3.3 billion-node heart-torso mesh at 0.2 mm resolution. The results show that the lead-field method is more efficient than a full-torso solution when the number of simulated samples is larger than the number of computed ECG leads. While the initial computation of the lead fields remains a hard and poorly scalable problem, the ECG computation itself scales almost perfectly and, even for several hundreds of ECG leads, takes much less time than the underlying simulation of cardiac activity.

Performance evaluation of GPU parallelization, space-time adaptive algorithms, and their combination for simulating cardiac electrophysiology

The use of computer models as a tool for the study and understanding of the complex phenomena of cardiac electrophysiology has attained increased importance nowadays. At the same time, the increased complexity of the biophysical processes translates into complex computational and mathematical models. To speed up cardiac simulations and to allow more precise and realistic uses, 2 different techniques have been traditionally exploited: parallel computing and sophisticated numerical methods. In this work, we combine a modern parallel computing technique based on multicore and graphics processing units (GPUs) and a sophisticated numerical method based on a new space-time adaptive algorithm. We evaluate each technique alone and in different combinations: multicore and GPU, multicore and GPU and space adaptivity, multicore and GPU and space adaptivity and time adaptivity. All the techniques and combinations were evaluated under different scenarios: 3D simulations on slabs, 3D simulations on a ventricular mouse mesh, ie, complex geometry, sinus-rhythm, and arrhythmic conditions. Our results suggest that multicore and GPU accelerate the simulations by an approximate factor of 33×, whereas the speedups attained by the space-time adaptive algorithms were approximately 48. Nevertheless, by combining all the techniques, we obtained speedups that ranged between 165 and 498. The tested methods were able to reduce the execution time of a simulation by more than 498× for a complex cellular model in a slab geometry and by 165× in a realistic heart geometry simulating spiral waves. The proposed methods will allow faster and more realistic simulations in a feasible time with no significant loss of accuracy.

Efficient computation of electrograms and ECGs in human whole heart simulations using a reaction-eikonal model

Anatomically accurate and biophysically detailed bidomain models of the human heart have proven a powerful tool for gaining quantitative insight into the links between electrical sources in the myocardium and the concomitant current flow in the surrounding medium as they represent their relationship mechanistically based on first principles. Such models are increasingly considered as a clinical research tool with the perspective of being used, ultimately, as a complementary diagnostic modality. An important prerequisite in many clinical modeling applications is the ability of models to faithfully replicate potential maps and electrograms recorded from a given patient. However, while the personalization of electrophysiology models based on the gold standard bidomain formulation is in principle feasible, the associated computational expenses are significant, rendering their use incompatible with clinical time frames. In this study we report on the development of a novel computationally efficient reaction-eikonal (R-E) model for modeling extracellular potential maps and electrograms. Using a biventricular human electrophysiology model, which incorporates a topologically realistic His-Purkinje system (HPS), we demonstrate by comparing against a high-resolution reaction-diffusion (R-D) bidomain model that the R-E model predicts extracellular potential fields, electrograms as well as ECGs at the body surface with high fidelity and offers vast computational savings greater than three orders of magnitude. Due to their efficiency R-E models are ideally suitable for forward simulations in clinical modeling studies which attempt to personalize electrophysiological model features.

Improved hybrid/GPU algorithm for solving cardiac electrophysiology problems on Purkinje networks

Cardiac Purkinje fibers provide an important pathway to the coordinated contraction of the heart. We present a numerical algorithm for the solution of electrophysiology problems across the Purkinje network that is efficient enough to be used in in silico studies on realistic Purkinje networks with physiologically detailed models of ion exchange at the cell membrane. The algorithm is on the basis of operator splitting and is provided with 3 different implementations: pure CPU, hybrid CPU/GPU, and pure GPU. Compared to our previous work, we modify the explicit gap junction term at network bifurcations to improve its mathematical consistency. Due to this improved consistency of the model, we are able to perform an empirical convergence study against analytical solutions. The study verified that all 3 implementations produce equivalent convergence rates, and shows that the algorithm produces equivalent result across different hardware platforms. Finally, we compare the efficiency of all 3 implementations on Purkinje networks of increasing spatial resolution using membrane models of increasing complexity. Both hybrid and pure GPU implementations outperform the pure CPU implementation, but their relative performance difference depends on the size of the Purkinje network and the complexity of the membrane model used.

PDE constrained optimization of electrical defibrillation in a 3D ventricular slice geometry

A computational study of an optimal control approach for cardiac defibrillation in a 3D geometry is presented. The cardiac bioelectric activity at the tissue and bath volumes is modeled by the bidomain model equations. The model includes intramural fiber rotation, axially symmetric around the fiber direction, and anisotropic conductivity coefficients, which are extracted from a histological image. The dynamics of the ionic currents are based on the regularized Mitchell-Schaeffer model. The controls enter in the form of electrodes, which are placed at the boundary of the bath volume with the goal of dampening undesired arrhythmias. The numerical optimization is based on Newton techniques. We demonstrated the parallel architecture environment for the computation of potentials on multidomains and for the higher order optimization techniques.

Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation

Electromechanical (EM) models of the heart have been used successfully to study fundamental mechanisms underlying a heart beat in health and disease. However, in all modeling studies reported so far numerous simplifications were made in terms of representing biophysical details of cellular function and its heterogeneity, gross anatomy and tissue microstructure, as well as the bidirectional coupling between electrophysiology (EP) and tissue distension. One limiting factor is the employed spatial discretization methods which are not sufficiently flexible to accommodate complex geometries or resolve heterogeneities, but, even more importantly, the limited efficiency of the prevailing solver techniques which are not sufficiently scalable to deal with the incurring increase in degrees of freedom (DOF) when modeling cardiac electromechanics at high spatio-temporal resolution. This study reports on the development of a novel methodology for solving the nonlinear equation of finite elasticity using human whole organ models of cardiac electromechanics, discretized at a high para-cellular resolution. Three patient-specific, anatomically accurate, whole heart EM models were reconstructed from magnetic resonance (MR) scans at resolutions of 220 μm, 440 μm and 880 μm, yielding meshes of approximately 184.6, 24.4 and 3.7 million tetrahedral elements and 95.9, 13.2 and 2.1 million displacement DOF, respectively. The same mesh was used for discretizing the governing equations of both electrophysiology (EP) and nonlinear elasticity. A novel algebraic multigrid (AMG) preconditioner for an iterative Krylov solver was developed to deal with the resulting computational load. The AMG preconditioner was designed under the primary objective of achieving favorable strong scaling characteristics for both setup and solution runtimes, as this is key for exploiting current high performance computing hardware. Benchmark results using the 220 μm, 440 μm and 880 μm meshes demonstrate efficient scaling up to 1024, 4096 and 8192 compute cores which allowed the simulation of a single heart beat in 44.3, 87.8 and 235.3 minutes, respectively. The efficiency of the method allows fast simulation cycles without compromising anatomical or biophysical detail.

Image-Based Personalization of Cardiac Anatomy for Coupled Electromechanical Modeling

Computational models of cardiac electromechanics (EM) are increasingly being applied to clinical problems, with patient-specific models being generated from high fidelity imaging and used to simulate patient physiology, pathophysiology and response to treatment. Current structured meshes are limited in their ability to fully represent the detailed anatomical data available from clinical images and capture complex and varied anatomy with limited geometric accuracy. In this paper, we review the state of the art in image-based personalization of cardiac anatomy for biophysically detailed, strongly coupled EM modeling, and present our own tools for the automatic building of anatomically and structurally accurate patient-specific models. Our method relies on using high resolution unstructured meshes for discretizing both physics, electrophysiology and mechanics, in combination with efficient, strongly scalable solvers necessary to deal with the computational load imposed by the large number of degrees of freedom of these meshes. These tools permit automated anatomical model generation and strongly coupled EM simulations at an unprecedented level of anatomical and biophysical detail.

Deranged sodium to sudden death

In February 2014, a group of scientists convened as part of the University of California Davis Cardiovascular Symposium to bring together experimental and mathematical modelling perspectives and discuss points of consensus and controversy on the topic of sodium in the heart. This paper summarizes the topics of presentation and discussion from the symposium, with a focus on the role of aberrant sodium channels and abnormal sodium homeostasis in cardiac arrhythmias and pharmacotherapy from the subcellular scale to the whole heart. Two following papers focus on Na(+) channel structure, function and regulation, and Na(+)/Ca(2+) exchange and Na(+)/K(+) ATPase. The UC Davis Cardiovascular Symposium is a biannual event that aims to bring together leading experts in subfields of cardiovascular biomedicine to focus on topics of importance to the field. The focus on Na(+) in the 2014 symposium stemmed from the multitude of recent studies that point to the importance of maintaining Na(+) homeostasis in the heart, as disruption of homeostatic processes are increasingly identified in cardiac disease states. Understanding how disruption in cardiac Na(+)-based processes leads to derangement in multiple cardiac components at the level of the cell and to then connect these perturbations to emergent behaviour in the heart to cause disease is a critical area of research. The ubiquity of disruption of Na(+) channels and Na(+) homeostasis in cardiac disorders of excitability and mechanics emphasizes the importance of a fundamental understanding of the associated mechanisms and disease processes to ultimately reveal new targets for human therapy.

Parallel Optimization of 3D Cardiac Electrophysiological Model Using GPU

Large-scale 3D virtual heart model simulations are highly demanding in computational resources. This imposes a big challenge to the traditional computation resources based on CPU environment, which already cannot meet the requirement of the whole computation demands or are not easily available due to expensive costs. GPU as a parallel computing environment therefore provides an alternative to solve the large-scale computational problems of whole heart modeling. In this study, using a 3D sheep atrial model as a test bed, we developed a GPU-based simulation algorithm to simulate the conduction of electrical excitation waves in the 3D atria. In the GPU algorithm, a multicellular tissue model was split into two components: one is the single cell model (ordinary differential equation) and the other is the diffusion term of the monodomain model (partial differential equation). Such a decoupling enabled realization of the GPU parallel algorithm. Furthermore, several optimization strategies were proposed based on the features of the virtual heart model, which enabled a 200-fold speedup as compared to a CPU implementation. In conclusion, an optimized GPU algorithm has been developed that provides an economic and powerful platform for 3D whole heart simulations.

Image-Based Personalization of Cardiac Anatomy for Coupled Electromechanical Modeling

Computational models of cardiac electromechanics (EM) are increasingly being applied to clinical problems, with patient-specific models being generated from high fidelity imaging and used to simulate patient physiology, pathophysiology and response to treatment. Current structured meshes are limited in their ability to fully represent the detailed anatomical data available from clinical images and capture complex and varied anatomy with limited geometric accuracy. In this paper, we review the state of the art in image-based personalization of cardiac anatomy for biophysically detailed, strongly coupled EM modeling, and present our own tools for the automatic building of anatomically and structurally accurate patient-specific models. Our method relies on using high resolution unstructured meshes for discretizing both physics, electrophysiology and mechanics, in combination with efficient, strongly scalable solvers necessary to deal with the computational load imposed by the large number of degrees of freedom of these meshes. These tools permit automated anatomical model generation and strongly coupled EM simulations at an unprecedented level of anatomical and biophysical detail.

OpenCOR: a modular and interoperable approach to computational biology

Computational biologists have been developing standards and formats for nearly two decades, with the aim of easing the description and exchange of experimental data, mathematical models, simulation experiments, etc. One of those efforts is CellML (cellml.org), an XML-based markup language for the encoding of mathematical models. Early CellML-based environments include COR and OpenCell. However, both of those tools have limitations and were eventually replaced with OpenCOR (opencor.ws). OpenCOR is an open source modeling environment that is supported on Windows, Linux and OS X. It relies on a modular approach, which means that all of its features come in the form of plugins. Those plugins can be used to organize, edit, simulate and analyze models encoded in the CellML format. We start with an introduction to CellML and two of its early adopters, which limitations eventually led to the development of OpenCOR. We then go onto describing the general philosophy behind OpenCOR, as well as describing its openness and its development process. Next, we illustrate various aspects of OpenCOR, such as its user interface and some of the plugins that come bundled with it (e.g., its editing and simulation plugins). Finally, we discuss some of the advantages and limitations of OpenCOR before drawing some concluding remarks.

A brief history of tissue models for cardiac electrophysiology

The last four decades have produced a number of significant advances in the developments of computer models to simulate and investigate the electrical activity of cardiac tissue. The tissue descriptions that underlie these simulations have been built from a combination of clever insight and careful comparison with measured data at multiple scales. Tissue models have not only led to greater insights into the mechanisms of life-threatening arrhythmias but have been used to engineer new therapies to treat the consequences of cardiac disease. This paper is a look back at the early years in the cardiac modeling and the challenges facing the field as models move toward the clinic.

An efficient finite element approach for modeling fibrotic clefts in the heart

Advanced medical imaging technologies provide a wealth of information on cardiac anatomy and structure at a paracellular resolution, allowing to identify microstructural discontinuities which disrupt the intracellular matrix. Current state-of-the-art computer models built upon such datasets account for increasingly finer anatomical details, however, structural discontinuities at the paracellular level are typically discarded in the model generation process, owing to the significant costs which incur when using high resolutions for explicit representation. In this study, a novel discontinuous finite element (dFE) approach for discretizing the bidomain equations is presented, which accounts for fine-scale structures in a computer model without the need to increase spatial resolution. In the dFE method, this is achieved by imposing infinitely thin lines of electrical insulation along edges of finite elements which approximate the geometry of discontinuities in the intracellular matrix. Simulation results demonstrate that the dFE approach accounts for effects induced by microscopic size scale discontinuities, such as the formation of microscopic virtual electrodes, with vast computational savings as compared to high resolution continuous finite element models. Moreover, the method can be implemented in any standard continuous finite element code with minor effort.

Fast acceleration of 2D wave propagation simulations using modern computational accelerators

Recent developments in modern computational accelerators like Graphics Processing Units (GPUs) and coprocessors provide great opportunities for making scientific applications run faster than ever before. However, efficient parallelization of scientific code using new programming tools like CUDA requires a high level of expertise that is not available to many scientists. This, plus the fact that parallelized code is usually not portable to different architectures, creates major challenges for exploiting the full capabilities of modern computational accelerators. In this work, we sought to overcome these challenges by studying how to achieve both automated parallelization using OpenACC and enhanced portability using OpenCL. We applied our parallelization schemes using GPUs as well as Intel Many Integrated Core (MIC) coprocessor to reduce the run time of wave propagation simulations. We used a well-established 2D cardiac action potential model as a specific case-study. To the best of our knowledge, we are the first to study auto-parallelization of 2D cardiac wave propagation simulations using OpenACC. Our results identify several approaches that provide substantial speedups. The OpenACC-generated GPU code achieved more than speedup above the sequential implementation and required the addition of only a few OpenACC pragmas to the code. An OpenCL implementation provided speedups on GPUs of at least faster than the sequential implementation and faster than a parallelized OpenMP implementation. An implementation of OpenMP on Intel MIC coprocessor provided speedups of with only a few code changes to the sequential implementation. We highlight that OpenACC provides an automatic, efficient, and portable approach to achieve parallelization of 2D cardiac wave simulations on GPUs. Our approach of using OpenACC, OpenCL, and OpenMP to parallelize this particular model on modern computational accelerators should be applicable to other computational models of wave propagation in multi-dimensional media.

Toward GPGPU accelerated human electromechanical cardiac simulations

In this paper, we look at the acceleration of weakly coupled electromechanics using the graphics processing unit (GPU). Specifically, we port to the GPU a number of components of CHeart--a CPU-based finite element code developed for simulating multi-physics problems. On the basis of a criterion of computational cost, we implemented on the GPU the ODE and PDE solution steps for the electrophysiology problem and the Jacobian and residual evaluation for the mechanics problem. Performance of the GPU implementation is then compared with single core CPU (SC) execution as well as multi-core CPU (MC) computations with equivalent theoretical performance. Results show that for a human scale left ventricle mesh, GPU acceleration of the electrophysiology problem provided speedups of 164 × compared with SC and 5.5 times compared with MC for the solution of the ODE model. Speedup of up to 72 × compared with SC and 2.6 × compared with MC was also observed for the PDE solve. Using the same human geometry, the GPU implementation of mechanics residual/Jacobian computation provided speedups of up to 44 × compared with SC and 2.0 × compared with MC.

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.

Simulating the Mechanics of Myocardial Tissue Using Strongly Scalable Parallel Algorithms

Due to preferential orientations of fibers, such as collagen or myocytes, the modeling of the mechanics of myocardial tissue leads to anisotropic and highly nonlinear material models. For micro-anatomically realistic geometries the computational effort to handle these sophisticated models is very challenging and demands the usage of strongly scalable parallel algorithms. In this context we mention two possible approaches and show numerical examples.

Simulations of complex and microscopic models of cardiac electrophysiology powered by multi-GPU platforms

Key aspects of cardiac electrophysiology, such as slow conduction, conduction block, and saltatory effects have been the research topic of many studies since they are strongly related to cardiac arrhythmia, reentry, fibrillation, or defibrillation. However, to reproduce these phenomena the numerical models need to use subcellular discretization for the solution of the PDEs and nonuniform, heterogeneous tissue electric conductivity. Due to the high computational costs of simulations that reproduce the fine microstructure of cardiac tissue, previous studies have considered tissue experiments of small or moderate sizes and used simple cardiac cell models. In this paper, we develop a cardiac electrophysiology model that captures the microstructure of cardiac tissue by using a very fine spatial discretization (8 μm) and uses a very modern and complex cell model based on Markov chains for the characterization of ion channel's structure and dynamics. To cope with the computational challenges, the model was parallelized using a hybrid approach: cluster computing and GPGPUs (general-purpose computing on graphics processing units). Our parallel implementation of this model using a multi-GPU platform was able to reduce the execution times of the simulations from more than 6 days (on a single processor) to 21 minutes (on a small 8-node cluster equipped with 16 GPUs, i.e., 2 GPUs per node).