Numerical quadrature and operator splitting in finite element methods for cardiac electrophysiology

Digital twinning of cardiac electrophysiology for congenital heart disease

In recent years, blending mechanistic knowledge with machine learning has had a major impact in digital healthcare. In this work, we introduce a computational pipeline to build certified digital replicas of cardiac electrophysiology in paediatric patients with congenital heart disease. We construct the patient-specific geometry by means of semi-automatic segmentation and meshing tools. We generate a dataset of electrophysiology simulations covering cell-to-organ level model parameters and using rigorous mathematical models based on differential equations. We previously proposed Branched Latent Neural Maps (BLNMs) as an accurate and efficient means to recapitulate complex physical processes in a neural network. Here, we employ BLNMs to encode the parametrized temporal dynamics of in silico 12-lead electrocardiograms (ECGs). BLNMs act as a geometry-specific surrogate model of cardiac function for fast and robust parameter estimation to match clinical ECGs in paediatric patients. Identifiability and trustworthiness of calibrated model parameters are assessed by sensitivity analysis and uncertainty quantification.

Effects of cardiac growth on electrical dyssynchrony in the single ventricle patient

Single ventricle patients, including those with hypoplastic left heart syndrome (HLHS), typically undergo three palliative heart surgeries culminating in the Fontan procedure. HLHS is associated with high rates of morbidity and mortality, and many patients develop arrhythmias, electrical dyssynchrony, and eventually ventricular failure. However, the correlation between ventricular enlargement and electrical dysfunction in HLHS physiology remains poorly understood. Here we characterize the relationship between growth and electrophysiology in HLHS using computational modeling. We integrate a personalized finite element model, a volumetric growth model, and a personalized electrophysiology model to perform controlled in silico experiments. We show that right ventricle enlargement negatively affects QRS duration and interventricular dyssynchrony. Conversely, left ventricle enlargement can partially compensate for this dyssynchrony. These findings have potential implications on our understanding of the origins of electrical dyssynchrony and, ultimately, the treatment of HLHS patients.

Branched Latent Neural Maps

We introduce Branched Latent Neural Maps (BLNMs) to learn finite dimensional input-output maps encoding complex physical processes. A BLNM is defined by a simple and compact feedforward partially-connected neural network that structurally disentangles inputs with different intrinsic roles, such as the time variable from model parameters of a differential equation, while transferring them into a generic field of interest. BLNMs leverage latent outputs to enhance the learned dynamics and break the curse of dimensionality by showing excellent in-distribution generalization properties with small training datasets and short training times on a single processor. Indeed, their in-distribution generalization error remains comparable regardless of the adopted discretization during the testing phase. Moreover, the partial connections, in place of a fully-connected structure, significantly reduce the number of tunable parameters. We show the capabilities of BLNMs in a challenging test case involving biophysically detailed electrophysiology simulations in a biventricular cardiac model of a pediatric patient with hypoplastic left heart syndrome. The model includes a 1D Purkinje network for fast conduction and a 3D heart-torso geometry. Specifically, we trained BLNMs on 150 in silico generated 12-lead electrocardiograms (ECGs) while spanning 7 model parameters, covering cell-scale, organ-level and electrical dyssynchrony. Although the 12-lead ECGs manifest very fast dynamics with sharp gradients, after automatic hyperparameter tuning the optimal BLNM, trained in less than 3 hours on a single CPU, retains just 7 hidden layers and 19 neurons per layer. The resulting mean square error is on the order of 10 - 4 on an independent test dataset comprised of 50 additional electrophysiology simulations. In the online phase, the BLNM allows for 5000x faster real-time simulations of cardiac electrophysiology on a single core standard computer and can be employed to solve inverse problems via global optimization in a few seconds of computational time. This paper provides a novel computational tool to build reliable and efficient reduced-order models for digital twinning in engineering applications. The Julia implementation is publicly available under MIT License at https://github.com/StanfordCBCL/BLNM.jl.

Digital twinning of cardiac electrophysiology for congenital heart disease

In recent years, blending mechanistic knowledge with machine learning has had a major impact in digital healthcare. In this work, we introduce a computational pipeline to build certified digital replicas of cardiac electrophysiology in pediatric patients with congenital heart disease. We construct the patient-specific geometry by means of semi-automatic segmentation and meshing tools. We generate a dataset of electrophysiology simulations covering cell-to-organ level model parameters and utilizing rigorous mathematical models based on differential equations. We previously proposed Branched Latent Neural Maps (BLNMs) as an accurate and efficient means to recapitulate complex physical processes in a neural network. Here, we employ BLNMs to encode the parametrized temporal dynamics of in silico 12-lead electrocardiograms (ECGs). BLNMs act as a geometry-specific surrogate model of cardiac function for fast and robust parameter estimation to match clinical ECGs in pediatric patients. Identifiability and trustworthiness of calibrated model parameters are assessed by sensitivity analysis and uncertainty quantification.

A mathematical model that integrates cardiac electrophysiology, mechanics, and fluid dynamics: Application to the human left heart

We propose a mathematical and numerical model for the simulation of the heart function that couples cardiac electrophysiology, active and passive mechanics and hemodynamics, and includes reduced models for cardiac valves and the circulatory system. Our model accounts for the major feedback effects among the different processes that characterize the heart function, including electro-mechanical and mechano-electrical feedback as well as force-strain and force-velocity relationships. Moreover, it provides a three-dimensional representation of both the cardiac muscle and the hemodynamics, coupled in a fluid-structure interaction (FSI) model. By leveraging the multiphysics nature of the problem, we discretize it in time with a segregated electrophysiology-force generation-FSI approach, allowing for efficiency and flexibility in the numerical solution. We employ a monolithic approach for the numerical discretization of the FSI problem. We use finite elements for the spatial discretization of partial differential equations. We carry out a numerical simulation on a realistic human left heart model, obtaining results that are qualitatively and quantitatively in agreement with physiological ranges and medical images.

A simple and efficient adaptive time stepping technique for low-order operator splitting schemes applied to cardiac electrophysiology

We present a simple, yet efficient adaptive time stepping scheme for cardiac electrophysiology (EP) simulations based on standard operator splitting techniques. The general idea is to exploit the relation between the splitting error and the reaction's magnitude-found in a previous one-dimensional analytical study by Spiteri and Ziaratgahi-to construct the new time step controller for three-dimensional problems. Accordingly, we propose to control the time step length of the operator splitting scheme as a function of the reaction magnitude, in addition to the common approach of adapting the reaction time step. This conforms with observations in numerical experiments supporting the need for a significantly smaller time step length during depolarization than during repolarization. The proposed scheme is compared with classical proportional-integral-differential controllers using state-of-the-art error estimators, which are also presented in details as they have not been previously applied in the context of cardiac EP with operator splitting techniques. Benchmarks show that choosing the time step as a sigmoidal function of the reaction magnitude is highly efficient and full cardiac cycles can be computed with precision even in a realistic biventricular setup. The proposed scheme outperforms common adaptive time stepping techniques, while depending on fewer tuning parameters.

Efficient time splitting schemes for the monodomain equation in cardiac electrophysiology

Approximating the fast dynamics of depolarization waves in the human heart described by the monodomain model is numerically challenging. Splitting methods for the PDE-ODE coupling enable the computation with very fine space and time discretizations. Here, we compare different splitting approaches regarding convergence, accuracy, and efficiency. Simulations were performed for a benchmark problem with the Beeler-Reuter cell model on a truncated ellipsoid approximating the left ventricle including a localized stimulation. For this configuration, we provide a reference solution for the transmembrane potential. We found a semi-implicit approach with state variable interpolation to be the most efficient scheme. The results are transferred to a more physiological setup using a bi-ventricular domain with a complex external stimulation pattern to evaluate the accuracy of the activation time for different resolutions in space and time.

Balancing conduction velocity error in cardiac electrophysiology using a modified quadrature approach

Conduction velocity error is often the main culprit behind the need for very fine spatial discretizations and high computational effort in cardiac electrophysiology problems. In light of this, a novel approach for simulating an accurate conduction velocity in coarse meshes with linear elements is suggested based on a modified quadrature approach. In this approach, the quadrature points are placed at arbitrary offsets of the isoparametric coordinates. A numerical study illustrates the dependence of the conduction velocity on the spatial discretization and the conductivity when using different quadrature rules and calculation approaches. Additionally, examples using the modified quadrature in coarse meshes for wave propagation demonstrate the improved accuracy of the conduction velocity with this method. This novel approach possesses great potential in reducing the computational effort required but remains limited to specific linear elements and experiences a reduction in accuracy for irregular meshes and heterogeneous conductivities. Further research can focus on developing an adaptive quadrature and extending the approach to other element formulations in order to make the approach more generally applicable.

A Computational Study of the Electrophysiological Substrate in Patients Suffering From Atrial Fibrillation

In the context of cardiac electrophysiology, we propose a novel computational approach to highlight and explain the long-debated mechanisms behind atrial fibrillation (AF) and to reliably numerically predict its induction and sustainment. A key role is played, in this respect, by a new way of setting a parametrization of electrophysiological mathematical models based on conduction velocities; these latter are estimated from high-density mapping data, which provide a detailed characterization of patients' electrophysiological substrate during sinus rhythm. We integrate numerically approximated conduction velocities into a mathematical model consisting of a coupled system of partial and ordinary differential equations, formed by the monodomain equation and the Courtemanche-Ramirez-Nattel model. Our new model parametrization is then adopted to predict the formation and self-sustainment of localized reentries characterizing atrial fibrillation, by numerically simulating the onset of ectopic beats from the pulmonary veins. We investigate the paroxysmal and the persistent form of AF starting from electro-anatomical maps of two patients. The model's response to stimulation shows how substrate characteristics play a key role in inducing and sustaining these arrhythmias. Localized reentries are less frequent and less stable in case of paroxysmal AF, while they tend to anchor themselves in areas affected by severe slow conduction in case of persistent AF.

A numerical study on the effects of spatial and temporal discretization in cardiac electrophysiology

Millions of degrees of freedom are often required to accurately represent the electrophysiology of the myocardium due to the presence of discretization effects. This study seeks to explore the influence of temporal and spatial discretization on the simulation of cardiac electrophysiology in conjunction with changes in modeling choices. Several finite element analyses are performed to examine how discretization affects solution time, conduction velocity and electrical excitation. Discretization effects are considered along with changes in the electrophysiology model and solution approach. Two action potential models are considered: the Aliev-Panfilov model and the ten Tusscher-Noble-Noble-Panfilov model. The solution approaches consist of two time integration schemes and different treatments for solving the local system of ordinary differential equations. The efficiency and stability of the calculation approaches are demonstrated to be dependent on the action potential model. The dependency of the conduction velocity on the element size and time step is shown to be different for changes in material parameters. Finally, the discrepancies between the wave propagation in coarse and fine meshes are analyzed based on the temporal evolution of the transmembrane potential at a node and its neighboring Gauss points. Insight obtained from this study can be used to suggest new methods to improve the efficiency of simulations in cardiac electrophysiology.

Efficiency of semi-implicit alternating direction implicit methods for solving cardiac monodomain model

It is well known that numerical simulations of the cardiac monodomain model require fine mesh resolution, which increases the computational resources required. In this paper, we construct three operator-splitting alternating direction implicit (ADI) schemes to efficiently solve the nonlinear cardiac monodomain model. The main objective of the proposed methods is to reduce the computational time and memory consumed for solving electrocardiology models, compared to standard numerical methods. The proposed methods have second-order accuracy in both space and time while evaluating the ionic model only once per time-step. Several examples using regular wave, spiral wave reentry, and nonsymmetrical scroll wave are conducted, and the efficiency of the proposed ADI methods is compared to the standard semi-implicit Crank-Nicolson/Adams-Bashforth method. Large-scale two- and three-dimensional simulations are performed.

Sensitivity analysis of an electrophysiology model for the left ventricle

Modelling the cardiac electrophysiology entails dealing with the uncertainties related to the input parameters such as the heart geometry and the electrical conductivities of the tissues, thus calling for an uncertainty quantification (UQ) of the results. Since the chambers of the heart have different shapes and tissues, in order to make the problem affordable, here we focus on the left ventricle with the aim of identifying which of the uncertain inputs mostly affect its electrophysiology. In a first phase, the uncertainty of the input parameters is evaluated using data available from the literature and the output quantities of interest (QoIs) of the problem are defined. According to the polynomial chaos expansion, a training dataset is then created by sampling the parameter space using a quasi-Monte Carlo method whereas a smaller independent dataset is used for the validation of the resulting metamodel. The latter is exploited to run a global sensitivity analysis with nonlinear variance-based indices and thus reduce the input parameter space accordingly. Thereafter, the uncertainty probability distribution of the QoIs are evaluated using a direct UQ strategy on a larger dataset and the results discussed in the light of the medical knowledge.

Integration of activation maps of epicardial veins in computational cardiac electrophysiology

In this work we address the issue of validating the monodomain equation used in combination with the Bueno-Orovio ionic model for the prediction of the activation times in cardiac electro-physiology of the left ventricle. To this aim, we consider four patients who suffered from Left Bundle Branch Block (LBBB). We use activation maps performed at the septum as input data for the model and maps at the epicardial veins for the validation. In particular, a first set (half) of the latter are used to estimate the conductivities of the patient and a second set (the remaining half) to compute the errors of the numerical simulations. We find an excellent agreement between measures and numerical results. Our validated computational tool could be used to accurately predict activation times at the epicardial veins with a short mapping, i.e. by using only a part (the most proximal) of the standard acquisition points, thus reducing the invasive procedure and exposure to radiation.

Stabilized hybrid discontinuous Galerkin finite element method for the cardiac monodomain equation

Numerical methods for solving the cardiac electrophysiology model, which describes the electrical activity in the heart, are proposed. The model problem consists of a nonlinear reaction-diffusion partial differential equation coupled to systems of ordinary differential equations that describes electrochemical reactions in cardiac cells. The proposed methods combine an operator splitting technique for the reaction-diffusion equation with primal hybrid methods for spatial discretization considering continuous or discontinuous approximations for the Lagrange multiplier. A static condensation is adopted to form a reduced global system in terms of the multiplier only. Convergence studies exhibit optimal rates of convergence and numerical experiments show that the proposed schemes can be more efficient than standard numerical techniques commonly used in this context when preconditioned iterative methods are used for the solution of linear systems.

Non-ohmic tissue conduction in cardiac electrophysiology: Upscaling the non-linear voltage-dependent conductance of gap junctions

Gap junctions are key mediators of intercellular communication in cardiac tissue, and their function is vital to sustaining normal cardiac electrical activity. Conduction through gap junctions strongly depends on the hemichannel arrangement and transjunctional voltage, rendering the intercellular conductance highly non-Ohmic, particularly under steady-state regimes of conduction. Despite this marked non-linear behavior, current tissue-level models of cardiac conduction are rooted in the assumption that gap-junctions conductance is constant (Ohmic), which results in inaccurate predictions of electrical propagation, particularly in the low junctional-coupling regime observed under pathological conditions. In this work, we present a novel non-Ohmic homogenization model (NOHM) of cardiac conduction that is suitable to tissue-scale simulations. Using non-linear homogenization theory, we develop a conductivity model that seamlessly upscales the voltage-dependent conductance of gap junctions, without the need of explicitly modeling gap junctions. The NOHM model allows for the simulation of electrical propagation in tissue-level cardiac domains that accurately resemble that of cell-based microscopic models for a wide range of junctional coupling scenarios, recovering key conduction features at a fraction of the computational complexity. A unique feature of the NOHM model is the possibility of upscaling the response of non-symmetric gap-junction conductance distributions, which result in conduction velocities that strongly depend on the direction of propagation, thus allowing to model the normal and retrograde conduction observed in certain regions of the heart. We envision that the NOHM model will enable organ-level simulations that are informed by sub- and inter-cellular mechanisms, delivering an accurate and predictive in-silico tool for understanding the heart function. Codes are available for download at https://github.com/dehurtado/NonOhmicConduction.

The heart relies on the propagation of electrical impulses that are mediated gap junctions, whose conduction properties vary depending on the transjunctional voltage. Despite this non-linear feature, current mathematical models assume that cardiac tissue behaves like an Ohmic (linear) material, thus delivering inaccurate results when simulated in a computer. Here we present a novel mathematical multiscale model that explicitly includes the non-Ohmic response of gap junctions in its predictions. Our results show that the proposed model recovers important conduction features modulated by gap junctions at a fraction of the computational complexity. This contribution represents an important step towards constructing computer models of a whole heart that can predict organ-level behavior in reasonable computing times.

Semi-implicit Non-conforming Finite-Element Schemes for Cardiac Electrophysiology: A Framework for Mesh-Coarsening Heart Simulations

The field of computational cardiology has steadily progressed toward reliable and accurate simulations of the heart, showing great potential in clinical applications such as the optimization of cardiac interventions and the study of pro-arrhythmic effects of drugs in humans, among others. However, the computational effort demanded by in-silico studies of the heart remains challenging, highlighting the need of novel numerical methods that can improve the efficiency of simulations while targeting an acceptable accuracy. In this work, we propose a semi-implicit non-conforming finite-element scheme (SINCFES) suitable for cardiac electrophysiology simulations. The accuracy and efficiency of the proposed scheme are assessed by means of numerical simulations of the electrical excitation and propagation in regular and biventricular geometries. We show that the SINCFES allows for coarse-mesh simulations that reduce the computation time when compared to fine-mesh models while delivering wavefront shapes and conduction velocities that are more accurate than those predicted by traditional finite-element formulations based on the same coarse mesh, thus improving the accuracy-efficiency trade-off of cardiac simulations.

Predicting drug-induced arrhythmias by multiscale modeling

Drugs often have undesired side effects. In the heart, they can induce lethal arrhythmias such as torsades de pointes. The risk evaluation of a new compound is costly and can take a long time, which often hinders the development of new drugs. Here, we establish a high-resolution, multiscale computational model to quickly assess the cardiac toxicity of new and existing drugs. The input of the model is the drug-specific current block from single cell electrophysiology; the output is the spatio-temporal activation profile and the associated electrocardiogram. We demonstrate the potential of our model for a low-risk drug, ranolazine, and a high-risk drug, quinidine: For ranolazine, our model predicts a prolonged QT interval of 19.4% compared with baseline and a regular sinus rhythm at 60.15 beats per minute. For quinidine, our model predicts a prolonged QT interval of 78.4% and a spontaneous development of torsades de pointes both in the activation profile and in the electrocardiogram. Our model reveals the mechanisms by which electrophysiological abnormalities propagate across the spatio-temporal scales, from specific channel blockage, via altered single cell action potentials and prolonged QT intervals, to the spontaneous emergence of ventricular tachycardia in the form of torsades de pointes. Our model could have important implications for researchers, regulatory agencies, and pharmaceutical companies on rationalizing safe drug development and reducing the time-to-market of new drugs.

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.

Incorporating inductances in tissue-scale models of cardiac electrophysiology

In standard models of cardiac electrophysiology, including the bidomain and monodomain models, local perturbations can propagate at infinite speed. We address this unrealistic property by developing a hyperbolic bidomain model that is based on a generalization of Ohm's law with a Cattaneo-type model for the fluxes. Further, we obtain a hyperbolic monodomain model in the case that the intracellular and extracellular conductivity tensors have the same anisotropy ratio. In one spatial dimension, the hyperbolic monodomain model is equivalent to a cable model that includes axial inductances, and the relaxation times of the Cattaneo fluxes are strictly related to these inductances. A purely linear analysis shows that the inductances are negligible, but models of cardiac electrophysiology are highly nonlinear, and linear predictions may not capture the fully nonlinear dynamics. In fact, contrary to the linear analysis, we show that for simple nonlinear ionic models, an increase in conduction velocity is obtained for small and moderate values of the relaxation time. A similar behavior is also demonstrated with biophysically detailed ionic models. Using the Fenton-Karma model along with a low-order finite element spatial discretization, we numerically analyze differences between the standard monodomain model and the hyperbolic monodomain model. In a simple benchmark test, we show that the propagation of the action potential is strongly influenced by the alignment of the fibers with respect to the mesh in both the parabolic and hyperbolic models when using relatively coarse spatial discretizations. Accurate predictions of the conduction velocity require computational mesh spacings on the order of a single cardiac cell. We also compare the two formulations in the case of spiral break up and atrial fibrillation in an anatomically detailed model of the left atrium, and we examine the effect of intracellular and extracellular inductances on the virtual electrode phenomenon.

The importance of mechano-electrical feedback and inertia in cardiac electromechanics

In the past years, a number cardiac electromechanics models have been developed to better understand the excitation-contraction behavior of the heart. However, there is no agreement on whether inertial forces play a role in this system. In this study, we assess the influence of mass in electromechanical simulations, using a fully coupled finite element model. We include the effect of mechano-electrical feedback via stretch activated currents. We compare five different models: electrophysiology, electromechanics, electromechanics with mechano-electrical feedback, electromechanics with mass, and electromechanics with mass and mechano-electrical feedback. We simulate normal conduction to study conduction velocity and spiral waves to study fibrillation. During normal conduction, mass in conjunction with mechano-electrical feedback increased the conduction velocity by 8.12% in comparison to the plain electrophysiology case. During the generation of a spiral wave, mass and mechano-electrical feedback generated secondary wavefronts, which were not present in any other model. These secondary wavefronts were initiated in tensile stretch regions that induced electrical currents. We expect that this study will help the research community to better understand the importance of mechanoelectrical feedback and inertia in cardiac electromechanics.

Space-discretization error analysis and stabilization schemes for conduction velocity in cardiac electrophysiology

In cardiac electrophysiology, the propagation of the action potential may be described by a set of reaction-diffusion equations known as the bidomain model. The shape of the solution is determined by a balance of a strong reaction and a relatively weak diffusion, which leads to steep variations in space and time. From a numerical point of view, the sharp spatial gradients may be seen as particularly problematic, because computational grid resolution on the order of 0.1 mm or less is required, yielding considerable computational efforts on human geometries. In this paper, we discuss a number of well-known numerical schemes for the bidomain equation and show how the quality of the solution is affected by the spatial discretization. In particular, we study in detail the effect of discretization on the conduction velocity (CV), which is an important quantity from a physiological point of view. We show that commonly applied finite element techniques tend to overestimate the CV on coarse grids, while it tends to be underestimated by finite difference schemes. Furthermore, the choice of interpolation and discretization scheme for the nonlinear reaction term has a strong impact on the CV. Finally, we exploit the results of the error analysis to propose improved numerical methods, including a stabilized scheme that tends to correct the CV on coarse grids but converges to the correct solution as the grid is refined. Copyright © 2016 John Wiley & Sons, Ltd.

Electrophysiology of Heart Failure Using a Rabbit Model: From the Failing Myocyte to Ventricular Fibrillation

Heart failure is a leading cause of death, yet its underlying electrophysiological (EP) mechanisms are not well understood. In this study, we use a multiscale approach to analyze a model of heart failure and connect its results to features of the electrocardiogram (ECG). The heart failure model is derived by modifying a previously validated electrophysiology model for a healthy rabbit heart. Specifically, in accordance with the heart failure literature, we modified the cell EP by changing both membrane currents and calcium handling. At the tissue level, we modeled the increased gap junction lateralization and lower conduction velocity due to downregulation of Connexin 43. At the biventricular level, we reduced the apex-to-base and transmural gradients of action potential duration (APD). The failing cell model was first validated by reproducing the longer action potential, slower and lower calcium transient, and earlier alternans characteristic of heart failure EP. Subsequently, we compared the electrical wave propagation in one dimensional cables of healthy and failing cells. The validated cell model was then used to simulate the EP of heart failure in an anatomically accurate biventricular rabbit model. As pacing cycle length decreases, both the normal and failing heart develop T-wave alternans, but only the failing heart shows QRS alternans (although moderate) at rapid pacing. Moreover, T-wave alternans is significantly more pronounced in the failing heart. At rapid pacing, APD maps show areas of conduction block in the failing heart. Finally, accelerated pacing initiated wave reentry and breakup in the failing heart. Further, the onset of VF was not observed with an upregulation of SERCA, a potential drug therapy, using the same protocol. The changes introduced at the cell and tissue level have increased the failing heart’s susceptibility to dynamic instabilities and arrhythmias under rapid pacing. However, the observed increase in arrhythmogenic potential is not due to a steepening of the restitution curve (not present in our model), but rather to a novel blocking mechanism.

Ventricular fibrillation (VF) is one of the leading causes of sudden death. During VF, the electrical wave of activation in the heart breaks up chaotically. Consequently, the heart is unable to contract synchronously and pump blood to the rest of the body. In our work we formulate and validate a model of heart failure (HF) that allows us to evaluate the arrhythmogenic potential of individual and combined electrophysiological changes. In diagnostic cardiology, the electrocardiogram (ECG) is one of the most commonly used tools for detecting abnormalities in the heart electrophysiology. One of our goals is to use our numerical model to link changes at the cellular and tissue level in a failing heart to a numerically computed ECG. This allows us to characterize the precursor to and the risk of VF. In order to understand the mechanisms underlying VF in HF, we design a test that simulates a HF patient performing physical exercise. We show that under fast heart rates with changes in pacing, HF patients are more prone to VF due to a new conduction blocking mechanism. In the long term, our mathematical model is suitable for investigating the effect of drug therapies in HF.

High-order finite element methods for cardiac monodomain simulations

Computational modeling of tissue-scale cardiac electrophysiology requires numerically converged solutions to avoid spurious artifacts. The steep gradients inherent to cardiac action potential propagation necessitate fine spatial scales and therefore a substantial computational burden. The use of high-order interpolation methods has previously been proposed for these simulations due to their theoretical convergence advantage. In this study, we compare the convergence behavior of linear Lagrange, cubic Hermite, and the newly proposed cubic Hermite-style serendipity interpolation methods for finite element simulations of the cardiac monodomain equation. The high-order methods reach converged solutions with fewer degrees of freedom and longer element edge lengths than traditional linear elements. Additionally, we propose a dimensionless number, the cell Thiele modulus, as a more useful metric for determining solution convergence than element size alone. Finally, we use the cell Thiele modulus to examine convergence criteria for obtaining clinically useful activation patterns for applications such as patient-specific modeling where the total activation time is known a priori.

Regional segmentation of ventricular models to achieve repolarization dispersion in cardiac electrophysiology modeling

The electrocardiogram (ECG) is one of the most significant outputs of a computational model of cardiac electrophysiology because it relates the numerical results to clinical data and is a universal tool for diagnosing heart diseases. One key features of the ECG is the T-wave, which is caused by longitudinal and transmural heterogeneity of the action potential duration (APD). Thus, in order to model a correct wave of repolarization, different cell properties resulting in different APDs must be assigned across the ventricular wall and longitudinally from apex to base. To achieve this requirement, a regional parametrization of the heart is necessary. We propose a robust approach to obtain the transmural and longitudinal segmentation in a general heart geometry without relying on ad hoc procedures. Our approach is based on auxiliary harmonic lifting analyses, already used in the literature to generate myocardial fiber orientations. Specifically, the solution of a sequence of Laplace boundary value problems allows parametrically controlled segmentation of both heart ventricles. The flexibility and simplicity of the proposed method is demonstrated through several representative examples, varying the locations and extents of the epicardial, midwall, and endocardial layers. Effects of the control parameters on the T-wave morphology are illustrated via computed ECGs.

Simulation Methods and Validation Criteria for Modeling Cardiac Ventricular Electrophysiology

We describe a sequence of methods to produce a partial differential equation model of the electrical activation of the ventricles. In our framework, we incorporate the anatomy and cardiac microstructure obtained from magnetic resonance imaging and diffusion tensor imaging of a New Zealand White rabbit, the Purkinje structure and the Purkinje-muscle junctions, and an electrophysiologically accurate model of the ventricular myocytes and tissue, which includes transmural and apex-to-base gradients of action potential characteristics. We solve the electrophysiology governing equations using the finite element method and compute both a 6-lead precordial electrocardiogram (ECG) and the activation wavefronts over time. We are particularly concerned with the validation of the various methods used in our model and, in this regard, propose a series of validation criteria that we consider essential. These include producing a physiologically accurate ECG, a correct ventricular activation sequence, and the inducibility of ventricular fibrillation. Among other components, we conclude that a Purkinje geometry with a high density of Purkinje muscle junctions covering the right and left ventricular endocardial surfaces as well as transmural and apex-to-base gradients in action potential characteristics are necessary to produce ECGs and time activation plots that agree with physiological observations.

The Living Heart Project: A robust and integrative simulator for human heart function

The heart is not only our most vital, but also our most complex organ: Precisely controlled by the interplay of electrical and mechanical fields, it consists of four chambers and four valves, which act in concert to regulate its filling, ejection, and overall pump function. While numerous computational models exist to study either the electrical or the mechanical response of its individual chambers, the integrative electro-mechanical response of the whole heart remains poorly understood. Here we present a proof-of-concept simulator for a four-chamber human heart model created from computer topography and magnetic resonance images. We illustrate the governing equations of excitation-contraction coupling and discretize them using a single, unified finite element environment. To illustrate the basic features of our model, we visualize the electrical potential and the mechanical deformation across the human heart throughout its cardiac cycle. To compare our simulation against common metrics of cardiac function, we extract the pressure-volume relationship and show that it agrees well with clinical observations. Our prototype model allows us to explore and understand the key features, physics, and technologies to create an integrative, predictive model of the living human heart. Ultimately, our simulator will open opportunities to probe landscapes of clinical parameters, and guide device design and treatment planning in cardiac diseases such as stenosis, regurgitation, or prolapse of the aortic, pulmonary, tricuspid, or mitral valve.