Evaluation of a Rapid Anisotropic Model for ECG Simulation

Combination of personalized computational modeling and machine learning for optimization of left ventricular pacing site in cardiac resynchronization therapy

Introduction: The 30-50% non-response rate to cardiac resynchronization therapy (CRT) calls for improved patient selection and optimized pacing lead placement. The study aimed to develop a novel technique using patient-specific cardiac models and machine learning (ML) to predict an optimal left ventricular (LV) pacing site (ML-PS) that maximizes the likelihood of LV ejection fraction (LVEF) improvement in a given CRT candidate. To validate the approach, we evaluated whether the distance D_PS between the clinical LV pacing site (ref-PS) and ML-PS is associated with improved response rate and magnitude. Materials and methods: We reviewed retrospective data for 57 CRT recipients. A positive response was defined as a more than 10% LVEF improvement. Personalized models of ventricular activation and ECG were created from MRI and CT images. The characteristics of ventricular activation during intrinsic rhythm and biventricular (BiV) pacing with ref-PS were derived from the models and used in combination with clinical data to train supervised ML classifiers. The best logistic regression model classified CRT responders with a high accuracy of 0.77 (ROC AUC = 0.84). The LR classifier, model simulations and Bayesian optimization with Gaussian process regression were combined to identify an optimal ML-PS that maximizes the ML-score of CRT response over the LV surface in each patient. Results: The optimal ML-PS improved the ML-score by 17 ± 14% over the ref-PS. Twenty percent of the non-responders were reclassified as positive at ML-PS. Selection of positive patients with a max ML-score >0.5 demonstrated an improved clinical response rate. The distance D_PS was shorter in the responders. The max ML-score and D_PS were found to be strong predictors of CRT response (ROC AUC = 0.85). In the group with max ML-score > 0.5 and D_PS< 30 mm, the response rate was 83% compared to 14% in the rest of the cohort. LVEF improvement in this group was higher than in the other patients (16 ± 8% vs. 7 ± 8%). Conclusion: A new technique combining clinical data, personalized heart modelling and supervised ML demonstrates the potential for use in clinical practice to assist in optimizing patient selection and predicting optimal LV pacing lead position in HF candidates for CRT.

Impact of electrode orientation, myocardial wall thickness, and myofiber direction on intracardiac electrograms: numerical modeling and analytical solutions

Intracardiac electrograms (iEGMs) are time traces of the electrical potential recorded close to the heart muscle. We calculate unipolar and bipolar iEGMs analytically for a myocardial slab with parallel myofibers and validate them against numerical bidomain simulations. The analytical solution obtained via the method of mirrors is an infinite series of arctangents. It goes beyond the solid angle theory and is in good agreement with the simulations, even though bath loading effects were not accounted for in the analytical calculation. At a large distance from the myocardium, iEGMs decay as 1/R (unipolar), 1/R ² (bipolar and parallel), and 1/R ³ (bipolar and perpendicular to the endocardium). At the endocardial surface, there is a mathematical branch cut. Here, we show how a thicker myocardium generates iEGMs with larger amplitudes and how anisotropy affects the iEGM width and amplitude. If only the leading-order term of our expansion is retained, it can be determined how the conductivities of the bath, torso, myocardium, and myofiber direction together determine the iEGM amplitude. Our results will be useful in the quantitative interpretation of iEGMs, the selection of thresholds to characterize viable tissues, and for future inferences of tissue parameters.

Synchronization of repolarization after cardiac resynchronization therapy: a combined clinical and modeling study

Introduction

The changes in ventricular repolarization after cardiac resynchronization therapy (CRT) are poorly understood. This knowledge gap is addressed using a multimodality approach including electrocardiographic and echocardiographic measurements in patients and using patient‐specific computational modeling.

Methods

In 33 patients electrocardiographic and echocardiographic measurements were performed before and at various intervals after CRT, both during CRT‐ON and temporary CRT‐OFF. T‐wave area was calculated from vectorcardiograms, and reconstructed from the 12‐lead electrocardiography (ECG). Computer simulations were performed using a patient‐specific eikonal model of cardiac activation with spatially varying action potential duration (APD) and repolarization rate, fit to a patient's ECG.

Results

During CRT‐ON T‐wave area diminished within a day and remained stable thereafter, whereas QT‐interval did not change significantly. During CRT‐OFF T‐wave area doubled within 5 days of CRT, while QT‐interval and peak‐to‐end T‐wave interval hardly changed. Left ventricular (LV) ejection fraction only increased significantly increased after 1 month of CRT. Computer simulations indicated that the increase in T‐wave area during CRT‐OFF can be explained by changes in APD following chronic CRT that are opposite to the change in CRT‐induced activation time. These APD changes were associated with a reduction in LV dispersion in repolarization during chronic CRT.

Conclusion

T‐wave area during CRT‐OFF is a sensitive marker for adaptations in ventricular repolarization during chronic CRT that may include a reduction in LV dispersion of repolarization.

Machine Learning Prediction of Cardiac Resynchronisation Therapy Response From Combination of Clinical and Model-Driven Data

Background: Up to 30–50% of chronic heart failure patients who underwent cardiac resynchronization therapy (CRT) do not respond to the treatment. Therefore, patient stratification for CRT and optimization of CRT device settings remain a challenge.

Objective: The main goal of our study is to develop a predictive model of CRT outcome using a combination of clinical data recorded in patients before CRT and simulations of the response to biventricular (BiV) pacing in personalized computational models of the cardiac electrophysiology.

Materials and Methods: Retrospective data from 57 patients who underwent CRT device implantation was utilized. Positive response to CRT was defined by a 10% increase in the left ventricular ejection fraction in a year after implantation. For each patient, an anatomical model of the heart and torso was reconstructed from MRI and CT images and tailored to ECG recorded in the participant. The models were used to compute ventricular activation time, ECG duration and electrical dyssynchrony indices during intrinsic rhythm and BiV pacing from the sites of implanted leads. For building a predictive model of CRT response, we used clinical data recorded before CRT device implantation together with model-derived biomarkers of ventricular excitation in the left bundle branch block mode of activation and under BiV stimulation. Several Machine Learning (ML) classifiers and feature selection algorithms were tested on the hybrid dataset, and the quality of predictors was assessed using the area under receiver operating curve (ROC AUC). The classifiers on the hybrid data were compared with ML models built on clinical data only.

Results: The best ML classifier utilizing a hybrid set of clinical and model-driven data demonstrated ROC AUC of 0.82, an accuracy of 0.82, sensitivity of 0.85, and specificity of 0.78, improving quality over that of ML predictors built on clinical data from much larger datasets by more than 0.1. Distance from the LV pacing site to the post-infarction zone and ventricular activation characteristics under BiV pacing were shown as the most relevant model-driven features for CRT response classification.

Conclusion: Our results suggest that combination of clinical and model-driven data increases the accuracy of classification models for CRT outcomes.

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.

Reconstruction of three-dimensional biventricular activation based on the 12-lead electrocardiogram via patient-specific modelling

Aims

Non-invasive imaging of electrical activation requires high-density body surface potential mapping. The nine electrodes of the 12-lead electrocardiogram (ECG) are insufficient for a reliable reconstruction with standard inverse methods. Patient-specific modelling may offer an alternative route to physiologically constraint the reconstruction. The aim of the study was to assess the feasibility of reconstructing the fully 3D electrical activation map of the ventricles from the 12-lead ECG and cardiovascular magnetic resonance (CMR).

Methods and results

Ventricular activation was estimated by iteratively optimizing the parameters (conduction velocity and sites of earliest activation) of a patient-specific model to fit the simulated to the recorded ECG. Chest and cardiac anatomy of 11 patients (QRS duration 126–180 ms, documented scar in two) were segmented from CMR images. Scar presence was assessed by magnetic resonance (MR) contrast enhancement. Activation sequences were modelled with a physiologically based propagation model and ECGs with lead field theory. Validation was performed by comparing reconstructed activation maps with those acquired by invasive electroanatomical mapping of coronary sinus/veins (CS) and right ventricular (RV) and left ventricular (LV) endocardium. The QRS complex was correctly reproduced by the model (Pearson’s correlation r = 0.923). Reconstructions accurately located the earliest and latest activated LV regions (median barycentre distance 8.2 mm, IQR 8.8 mm). Correlation of simulated with recorded activation time was very good at LV endocardium (r = 0.83) and good at CS (r = 0.68) and RV endocardium (r = 0.58).

Conclusion

Non-invasive assessment of biventricular 3D activation using the 12-lead ECG and MR imaging is feasible. Potential applications include patient-specific modelling and pre-/per-procedural evaluation of ventricular activation.

Atlas-based methods for efficient characterization of patient-specific ventricular activation patterns

AIMS

Ventricular activation patterns can aid clinical decision-making directly by providing spatial information on cardiac electrical activation or indirectly through derived clinical indices. The aim of this work was to derive an atlas of the major modes of variation of ventricular activation from model-predicted 3D bi-ventricular activation time distributions and to relate these modes to corresponding vectorcardiograms (VCGs). We investigated how the resulting dimensionality reduction can improve and accelerate the estimation of activation patterns from surface electrogram measurements.

METHODS AND RESULTS

Atlases of activation time (AT) and VCGs were derived using principal component analysis on a dataset of simulated electrophysiology simulations computed on eight patient-specific bi-ventricular geometries. The atlases provided significant dimensionality reduction, and the modes of variation in the two atlases described similar features. Utility of the atlases was assessed by resolving clinical waveforms against them and the VCG atlas was able to accurately reconstruct the patient VCGs with fewer than 10 modes. A sensitivity analysis between the two atlases was performed by calculating a compact Jacobian. Finally, VCGs generated by varying AT atlas modes were compared with clinical VCGs to estimate patient-specific activation maps, and the resulting errors between the clinical and atlas-based VCGs were less than those from more computationally expensive method.

CONCLUSION

Atlases of activation and VCGs represent a new method of identifying and relating the features of these high-dimensional signals that capture the major sources of variation between patients and may aid in identifying novel clinical indices of arrhythmia risk or therapeutic outcome.

GEASI: Geodesic-based earliest activation sites identification in cardiac models

The identification of the initial ventricular activation sequence is a critical step for the correct personalization of patient-specific cardiac models. In healthy conditions, the Purkinje network is the main source of the electrical activation, but under pathological conditions the so-called earliest activation sites (EASs) are possibly sparser and more localized. Yet, their number, location and timing may not be easily inferred from remote recordings, such as the epicardial activation or the 12-lead electrocardiogram (ECG), due to the underlying complexity of the model. In this work, we introduce GEASI (Geodesic-based Earliest Activation Sites Identification) as a novel approach to simultaneously identify all EASs. To this end, we start from the anisotropic eikonal equation modeling cardiac electrical activation and exploit its Hamilton-Jacobi formulation to minimize a given objective function, for example, the quadratic mismatch to given activation measurements. This versatile approach can be extended to estimate the number of activation sites by means of the topological gradient, or fitting a given ECG. We conducted various experiments in 2D and 3D for in-silico models and an in-vivo intracardiac recording collected from a patient undergoing cardiac resynchronization therapy. The results demonstrate the clinical applicability of GEASI for potential future personalized models and clinical intervention.

Computer Modeling of the Heart for ECG Interpretation—A Review

Computer modeling of the electrophysiology of the heart has undergone significant progress. A healthy heart can be modeled starting from the ion channels via the spread of a depolarization wave on a realistic geometry of the human heart up to the potentials on the body surface and the ECG. Research is advancing regarding modeling diseases of the heart. This article reviews progress in calculating and analyzing the corresponding electrocardiogram (ECG) from simulated depolarization and repolarization waves. First, we describe modeling of the P-wave, the QRS complex and the T-wave of a healthy heart. Then, both the modeling and the corresponding ECGs of several important diseases and arrhythmias are delineated: ischemia and infarction, ectopic beats and extrasystoles, ventricular tachycardia, bundle branch blocks, atrial tachycardia, flutter and fibrillation, genetic diseases and channelopathies, imbalance of electrolytes and drug-induced changes. Finally, we outline the potential impact of computer modeling on ECG interpretation. Computer modeling can contribute to a better comprehension of the relation between features in the ECG and the underlying cardiac condition and disease. It can pave the way for a quantitative analysis of the ECG and can support the cardiologist in identifying events or non-invasively localizing diseased areas. Finally, it can deliver very large databases of reliably labeled ECGs as training data for machine learning.

An Inverse Problem Involving a Viscous Eikonal Equation with Applications in Electrophysiology

In this work we discuss the reconstruction of cardiac activation instants based on a viscous Eikonal equation from boundary observations. The problem is formulated as a least squares problem and solved by a projected version of the Levenberg-Marquardt method. Moreover, we analyze the well-posedness of the state equation and derive the gradient of the least squares functional with respect to the activation instants. In the numerical examples we also conduct an experiment in which the location of the activation sites and the activation instants are reconstructed jointly based on an adapted version of the shape gradient method from (J. Math. Biol. 79, 2033-2068, 2019). We are able to reconstruct the activation instants as well as the locations of the activations with high accuracy relative to the noise level.

Data integration for the numerical simulation of cardiac electrophysiology

The increasing availability of extensive and accurate clinical data is rapidly shaping cardiovascular care by improving the understanding of physiological and pathological mechanisms of the cardiovascular system and opening new frontiers in designing therapies and interventions. In this direction, mathematical and numerical models provide a complementary relevant tool, able not only to reproduce patient-specific clinical indicators but also to predict and explore unseen scenarios. With this goal, clinical data are processed and provided as inputs to the mathematical model, which quantitatively describes the physical processes that occur in the cardiac tissue. In this paper, the process of integration of clinical data and mathematical models is discussed. Some challenges and contributions in the field of cardiac electrophysiology are reported.

Short P-Wave Duration is a Marker of Higher Rate of Atrial Fibrillation Recurrences after Pulmonary Vein Isolation: New Insights into the Pathophysiological Mechanisms Through Computer Simulations

Background Short ECG P-wave duration has recently been demonstrated to be associated with higher risk of atrial fibrillation (AF). The aim of this study was to assess the rate of AF recurrence after pulmonary vein isolation in patients with a short P wave, and to mechanistically elucidate the observation by computer modeling. Methods and Results A total of 282 consecutive patients undergoing a first single-pulmonary vein isolation procedure for paroxysmal or persistent AF were included. Computational models studied the effect of adenosine and sodium conductance on action potential duration and P-wave duration (PWD). About 16% of the patients had a PWD of 110 ms or shorter (median PWD 126 ms, interquartile range, 115 ms-138 ms; range, 71 ms-180 ms). At Cox regression, PWD was significantly associated with AF recurrence (P=0.012). Patients with a PWD <110 ms (hazard ratio [HR], 2.20; 95% CI, 1.24-3.88; P=0.007) and patients with a PWD ≥140 (HR, 1.87, 95% CI, 1.06-3.30; P=0.031) had a nearly 2-fold increase in risk with respect to the other group. In the computational model, adenosine yielded a significant reduction of action potential duration 90 (52%) and PWD (7%). An increased sodium conductance (up to 200%) was robustly accompanied by an increase in conduction velocity (26%), a reduction in action potential duration 90 (28%), and PWD (22%). Conclusions One out of 5 patients referred for pulmonary vein isolation has a short PWD which was associated with a higher rate of AF after the index procedure. Computer simulations suggest that shortening of atrial action potential duration leading to a faster atrial conduction may be the cause of this clinical observation.

An Inverse Eikonal Method for Identifying Ventricular Activation Sequences from Epicardial Activation Maps

A key mechanism controlling cardiac function is the electrical activation sequence of the heart's main pumping chambers termed the ventricles. As such, personalization of the ventricular activation sequences is of pivotal importance for the clinical utility of computational models of cardiac electrophysiology. However, a direct observation of the activation sequence throughout the ventricular volume is virtually impossible. In this study, we report on a novel method for identification of activation sequences from activation maps measured at the outer surface of the heart termed the epicardium. Conceptually, the method attempts to identify the key factors governing the ventricular activation sequence - the timing of earliest activation sites (EAS) and the velocity tensor field within the ventricular walls - from sparse and noisy activation maps sampled from the epicardial surface and fits an Eikonal model to the observations. Regularization methods are first investigated to overcome the severe ill-posedness of the inverse problem in a simplified 2D example. These methods are then employed in an anatomically accurate biventricular model with two realistic activation models of varying complexity - a simplified trifascicular model (3F) and a topologically realistic model of the His-Purkinje system (HPS). Using epicardial activation maps at full resolution, we first demonstrate that reconstructing the volumetric activation sequence is, in principle, feasible under the assumption of known location of EAS and later evaluate robustness of the method against noise and reduced spatial resolution of observations. Our results suggest that the FIMIN algorithm is able to robustly recover the full 3D activation sequence using epicardial activation maps at a spatial resolution achievable with current mapping systems and in the presence of noise. Comparing the accuracy achieved in the reconstructed activation maps with clinical data uncertainties suggests that the FIMIN method may be suitable for the patient- specific parameterization of activation models.

Patient-Specific Cardiac Parametrization from Eikonal Simulations

Simulations in cardiac electrophysiology use the bidomain equations to describe the electrical potential in the heart. If only the electrical activation sequence in the heart is needed, then the full bidomain equations can be substituted by the Eikonal equation which allows much faster responses w.r.t. the changed material parameters in the equation. We use our Eikonal solver optimized for memory usage and parallelization. Patient-specific simulations in cardiac electrophysiology require patient-specific conductivity parameters which are not accurately available in vivo. One chance to improve the given conductivity parameters consists in comparing the computed activation sequence on the heart surface with the measured ECG on the torso mapped onto this surface. By minimizing the squared distance between the measured solution and the Eikonal computed solution we are able to determine the material parameters more accurately. To reduce the number of optimization parameters in this process, we group the material parameters and introduce a specific scaling parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _k$$\end{document} for each group. The minimization takes place w.r.t. the scaling \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\gamma }$$\end{document}. We solve the minimization problem by the BFGS method and adaptive step size control. The required gradient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla _\gamma f({\gamma })$$\end{document} is computed either via finite differences or algorithmic differentiation using in tangent as well as in adjoint mode. We present convergence behavior as well as runtime and scaling results.

Patient-Specific Identification of Atrial Flutter Vulnerability-A Computational Approach to Reveal Latent Reentry Pathways

Atypical atrial flutter (AFlut) is a reentrant arrhythmia which patients frequently develop after ablation for atrial fibrillation (AF). Indeed, substrate modifications during AF ablation can increase the likelihood to develop AFlut and it is clinically not feasible to reliably and sensitively test if a patient is vulnerable to AFlut. Here, we present a novel method based on personalized computational models to identify pathways along which AFlut can be sustained in an individual patient. We build a personalized model of atrial excitation propagation considering the anatomy as well as the spatial distribution of anisotropic conduction velocity and repolarization characteristics based on a combination of a priori knowledge on the population level and information derived from measurements performed in the individual patient. The fast marching scheme is employed to compute activation times for stimuli from all parts of the atria. Potential flutter pathways are then identified by tracing loops from wave front collision sites and constricting them using a geometric snake approach under consideration of the heterogeneous wavelength condition. In this way, all pathways along which AFlut can be sustained are identified. Flutter pathways can be instantiated by using an eikonal-diffusion phase extrapolation approach and a dynamic multifront fast marching simulation. In these dynamic simulations, the initial pattern eventually turns into the one driven by the dominant pathway, which is the only pathway that can be observed clinically. We assessed the sensitivity of the flutter pathway maps with respect to conduction velocity and its anisotropy. Moreover, we demonstrate the application of tailored models considering disease-specific repolarization properties (healthy, AF-remodeled, potassium channel mutations) as well as applicabiltiy on a clinical dataset. Finally, we tested how AFlut vulnerability of these substrates is modulated by exemplary antiarrhythmic drugs (amiodarone, dronedarone). Our novel method allows to assess the vulnerability of an individual patient to develop AFlut based on the personal anatomical, electrophysiological, and pharmacological characteristics. In contrast to clinical electrophysiological studies, our computational approach provides the means to identify all possible AFlut pathways and not just the currently dominant one. This allows to consider all relevant AFlut pathways when tailoring clinical ablation therapy in order to reduce the development and recurrence of AFlut.

A general framework dedicated to computational morphogenesis Part I - Constitutive equations

In order to understand living organisms, considerable experimental efforts and resources have been devoted to correlate genes and their expressions with cell, tissue, organ and whole organisms' phenotypes. This data driven approach to knowledge discovery has led to many breakthrough in our understanding of healthy and diseased states, and is paving the way to improve the diagnosis and treatment of diseases. Complementary to this data-driven approach, computational models of biological systems based on first principles have been developed in order to deepen our understanding of the multi-scale dynamics that drives normal and pathological biological functions. In this paper we describe the biological, physical and mathematical concepts that led to the design of a Computational Morphogenesis (CM) platform baptized Generic Modeling and Simulating Platform (GMSP). Its role is to generate realistic 3D multi-scale biological tissues from virtual stem cells and the intended target applications include in virtuo studies of normal and abnormal tissue (re)generation as well as the development of complex diseases such as carcinogenesis. At all space-scales of interest, biological agents interact with each other via biochemical, bioelectrical, and mechanical fields that operate in concert during embryogenesis, growth and adult life. The spatio-temporal dependencies of these fields can be modeled by physics-based constitutive equations that we propose to examine in relation to the landmark biological events that occur during embryogenesis.

Fast uncertainty quantification of activation sequences in patient-specific cardiac electrophysiology meeting clinical time constraints

We present a fast, patient-specific methodology for uncertainty quantification in electrophysiology, aimed at meeting the time constraints of clinical practitioners. We focus on computing the statistics of the activation map, given the uncertainties associated with the conductivity tensor modeling the fiber orientation in the heart. We use a fast parallel solution method implemented on a graphics processing unit for the eikonal approximation, in order to compute the activation map and to sample the random fiber field with correlation on the basis of geodesic distances. While this enables to perform uncertainty quantification studies with a manageable computational effort, the required time frame still exceeds clinically suitable time expectations. In order to reduce it further by 2 orders of magnitude, we rely on Bayesian multifidelity methods. In particular, we propose a low-fidelity model that is patient-specific and free from the additional training cost associated with reduced models. This is achieved by a sound physics-based simplification of the full eikonal model. The low-fidelity output is then corrected by the standard multifidelity framework. In practice, the complete procedure only requires approximately 100 new runs of our eikonal graphics processing unit solver for producing the sought estimates and their associated credible intervals, enabling a full online analysis in less than 5 minutes.

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.