Edited by: Alexander Panfilov, Ghent University, Belgium
Reviewed by: Simone Scacchi, University of Milan, Italy; Chris Patrick Bradley, University of Auckland, New Zealand
*Correspondence: Simone Pezzuto
This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
State-of-the-art cardiac electrophysiology models that are able to deliver physiologically motivated activation maps and electrocardiograms (ECGs) can only be solved on high-performance computing architectures. This makes it nearly impossible to adopt such models in clinical practice. ECG imaging tools typically rely on simplified models, but these neglect the anisotropic electric conductivity of the tissue in the forward problem. Moreover, their results are often confined to the heart-torso interface. We propose a forward model that fully accounts for the anisotropic tissue conductivity and produces the standard 12-lead ECG in a few seconds. The activation sequence is approximated with an eikonal model in the 3d myocardium, while the ECG is computed with the lead-field approach. Both solvers were implemented on graphics processing units and massively parallelized. We studied the numerical convergence and scalability of the approach. We also compared the method to the bidomain model in terms of ECGs and activation maps, using a simplified but physiologically motivated geometry and 6 patient-specific anatomies. The proposed methods provided a good approximation of activation maps and ECGs computed with a bidomain model, in only a few seconds. Both solvers scaled very well to high-end hardware. These methods are suitable for use in ECG imaging methods, and may soon become fast enough for use in interactive simulation tools.
The cardiac muscle or myocardium consists of electrically active muscle cells (myocytes) connected to each other by gap junctions, embedded in extracellular fluid and other structures. The current state of the art in forward simulation of the electrocardiogram (ECG) is to use a bidomain reaction-diffusion model to simulate propagating electrical activation and repolarization as well as the potential field in the torso (Vigmond et al.,
where
The functions
Much simpler models are presently being used inside inverse models, which aim at recovering cardiac events from non-invasive signals (Ramanathan et al.,
The conductivities play three roles in the ECG:
They determine the conduction velocity, which will also be anisotropic;
They determine the current dipole density resulting from a voltage gradient; and
They determine how the resulting current flows through the heart and torso and thus what voltages will be measured by the ECG electrodes.
While virtually all models nowadays account for anisotropic propagation, rapid ECG models typically neglect the second and third role of the anisotropic conductivities. The major reason for this is that the conductivities are non-uniform: the fiber orientation rotates around the cavities of the heart and also changes transmurally. To account for it, therefore, a volumetric discretization of the heart and torso is necessary. By assuming isotropy, in contrast, existing methods simplify the problem to one on the surface of the muscle (Geselowitz,
One possible way to reduce the computational intensity of the bidomain model is to use an eikonal model (Colli Franzone and Guerri,
The method we propose here uses a volumetric anisotropic eikonal propagation model, pre-computed action potentials, and an ECG simulation based on dipole sources and a set of transfer functions known historically as lead fields. The volumetric eikonal model, in contrast to the surface models used by others (van Dam et al.,
Precomputed action potentials are assigned and shifted according to the computed activation times, resulting in a field of transmembrane potentials at any desired time instant. For each time instant, the current dipoles are computed with an anisotropic formula, and the ECG potentials are computed from them using a set of transfer coefficients. The transfer coefficients are computed once for each model by solving a bidomain problem in the full torso, again taking all anisotropic conductivities into account. This method of ECG simulation is mathematically equivalent with full bidomain solutions of the torso potential based on monopole or dipole sources (Potse et al.,
The eikonal model allows for the use of a much coarser spatial resolution of 1 mm (Colli Franzone and Guerri,
The method was entirely implemented in C, C++ and CUDA so it could run efficiently on general-purpose graphics processing units (GPGPU).
Using a reaction-diffusion and bidomain model as a reference, we evaluated the proposed model in terms of the accuracy of the simulated activation sequence and ECG.
Evaluations were performed using six patient-specific heart-torso models. Finally, we evaluated the parallel scalability of the methods on GPGPUs.
A “monodomain” reaction-diffusion model of the heart coupled to a bidomain heart-torso model was used as reference. The monodomain model can be derived from the bidomain model by assuming that the two domains have the same anisotropy ratio,
where
To compute the ECG, we solved the equation
in a 3d finite-difference full torso model at 1-mm resolution, including a downsampled heart model (Potse et al.,
All reference computations were performed with the propag-5 software (Potse et al.,
To compute the activation sequence of the heart we used a non-linear first-order Hamilton-Jacobi eikonal equation in the form
where Ψ(
The Equation (6), which assigns a locally constant front velocity, is one of the simplest eikonal models. More sophisticated versions of the eikonal equation are also available, where the velocity of the front depends also on the curvature of the front itself and higher-order curvature terms, or where the propagation is defined by a Finsler metric rather than a Riemannian metric (Colli Franzone et al.,
In general, the scaling parameter α is the velocity that the front would have if the conductivity and β were set to one. Indeed, given a propagation direction
which is exactly the conduction velocity given by the eikonal model for α = 1. α depends only on the membrane capacitance and the ionic model, in particular on the sodium-channel conductivity and gating dynamics.
The implementation of our eikonal solver is based on a local variational principle proposed by Bornemann and Rasch (
for
The quantity δ(
In order to compute it, we discretize the metric
The minimization problem in Equation (8) is solved by taking the minimum of the minima for each triangular face t of ∂ω(
The minimization over the single triangular face is a convex optimization problem, hence the minimum is unique.
It is possible to show that the map Λ has a unique fixed point (Bornemann and Rasch,
Each ECG lead potential
where Ω is the heart domain and
where
To obtain
However, we may also introduce position-dependent action potentials
that is, the convolution of the first derivative of
where Ψ−1(
We implemented two methods for the computation of the ECGs: one based on the general formula (11), termed “simple method” (SM), and one based on formulas (14) and (15), named “fast method” (FM). The former can account for action potential heterogeneity, which is important for the T-wave in the ECG. When focusing on the QRS complex, the main differences between the proposed methods are the way the cardiac sources are computed and resulting execution speed of the computation of the ECGs.
The SM computes the ECG by direct evaluation of Equation (11). The heart domain Ω is discretized into voxels v ∈ Ω
The convergence rate with respect to
The FM, based on Equation (14), consists of two steps. First, the piecewise linear interpolants in time of
The interpolated function
where
It is important to note at this stage that Ψ−1(τ
Finally, the second step is the convolution, approximated by
which results from a central finite difference approximation of
The proposed computational methods for activation times and ECG were designed to run on GPGPU architectures. The power of GPGPU computation lies in a massive parallelism that can hide long-latency operations such as memory accesses. The latency hiding is possible thanks to the warp schedulers, which continuously switch among warps whose instructions are ready for execution.
The solvers were implemented in C and C++ with CUDA extension (Nickolls et al.,
The eikonal Equation (6) is solved by means of the Fast Iterative Method (FIM), proposed by Jeong and Whitaker (
Briefly, the FIM algorithm proceeds as follows: first, the activation map is initialized to +∞ except at the boundary points; then, we iterate until a list of
The implementation of ECG computation with the SM and FM is similar. Each cube is assigned to an individual thread, and each thread updates the local contribution of one cube to the ECG. For the SM, the kernel iterates over time and computes
The FM has two kernels: one for the computation of
The computational complexity in time of the SM is
Since the computations of individual ECG leads are independent from each other, we also utilized parallel CUDA streams to support simultaneous execution of multiple kernels computing multiple ECG leads. All memory transfers between the CPU and GPU are asynchronous, in order to support also overlapping of the data transfers with the kernel executions.
Patient-specific heart and torso geometries were created as described previously (Nguyên et al.,
The model geometry was described in terms of a set of Catmull-Clark subdivision surfaces. The simulations were performed on semi-structured hexahedral meshes, i.e., while the basis of the mesh was regular, we did not store information for mesh nodes that did not play a role in the computation. These meshes were created by overlying the mesh base on the set of surfaces and assigning the properties of the mesh elements based on the surfaces in which they were included. Passive tissue characteristics were properties of elements, while active parameters and variables applied to nodes. The node properties followed from the element properties using a set of rules that ensured model consistency (Potse et al.,
In order to check the convergence of the method, we performed two numerical experiments. In both of them, we solved the eikonal equation on a tissue slab of size 1.5 × 2 × 2 cm, with a single pacing site at x0 = (0, 1, 1) cm.
In the first experiment, parameters were selected such that an analytical solution for both the activation times and the function
where the definitions and numerical values of the parameters are given in Table
α | Conduction velocity scaling | 2.0 cm ms−1 mS−1/2 |
β | Surface-to-volume ratio | 1000 cm−1 |
σel | Longitudinal extra-cellular conductivity | 3 mS cm−1 |
σet | Transverse extra-cellular conductivity | 1.2 mS cm−1 |
σec | Cross extra-cellular conductivity | 1.2 mS cm−1 |
σil | Longitudinal intra-cellular conductivity | 3 mS cm−1 |
σit | Transverse intra-cellular conductivity | 0.3 mS cm−1 |
σic | Cross intra-cellular conductivity | 0.3 mS cm−1 |
Vref | Resting potential | −80 mV |
Vdep | Depolarization potential | 30 mV |
εdep | Depolarization time-scale | 1 ms |
εrep | Repolarization time-scale | 100 ms |
APd | Action potential duration | 200 ms |
The exact solution for the activation times reads:
The ECG solution is the convolution between the first derivative of the action potential (20) and the function
The activated region on these two faces of the slab was always the intersection between an axis-aligned ellipse and a square. The corresponding area was derived from elementary geometry. See the
The convergence rate for the activation map was
The SM and the FM perform similarly, with FM being slightly more accurate than SM at coarser grids. The high accuracy of the SM is surprising because the template action potentials had very steep gradients in space (0.2 mm thick), clearly not captured by the coarse grid. But a closer inspection of the method showed that while the current per unit volume was not correctly approximated, the total current once integrated over the voxel was accurate.
The function
The second numerical experiment mimicked a transmural specimen of left ventricle. Fibers were parallel to the
with negative terminal x− = (3, 1, 1) cm and positive terminal x+ = (−1, 1, 1) cm.
In contrast to the first experiment, an analytical solution for the activation times, and thus for the ECG, was not available. The error was estimated using a solution obtained at a resolution of 2−8 ≈ 0.004 cm as reference. The convergence rate was roughly
From these experiments we concluded that a resolution of 1 or 0.5 mm for the eikonal solver provides reasonably accurate solutions, with a relative error lower than 5%.
The quantitative evaluation of the modeling error introduced by the proposed model with respect to the reference reaction-diffusion (R-D) model was conducted on a simple geometry but with physiologically motivated heterogeneities in the parameters. For this comparison we focused on the depolarization phase.
We considered a tissue slab shaped as in the previous section and embedded at the center of an electrically conductive bath of dimension 2.5 × 3 × 3 cm. The fiber organization in the slab was exactly as in the second experiment above, but we also introduced a 0.5 mm thick rapidly conducting layer at the bottom face of the slab. The propagation was initiated in a 0.5 × 1 × 1 mm volume at the center of the rapidly conducting layer. The bath was composed of two homogeneous and isotropic media, one 1 mm thick placed on top of the bath, which we named “skeletal muscle,” and the other, called “fluid,” filling the remaining part of the volume. The activation time was measured only in the tissue slab. A bipolar ECG was also produced by taking the potential difference between two electrodes placed at the centers of the top and bottom faces of the bath. An overview of the experimental setting is depicted in Figure
The parameters in the R-D model and the proposed model, summarized in Table
Tissue | 3.0 | 0.3 | 3.0 | 1.2 | 800 |
Fast | 3.0 | 0.3 | 3.0 | 1.2 | 356 |
Muscle | – | – | 0.44 | 0.44 | – |
Fluid | – | – | 6.0 | 6.0 | – |
The value for α in the eikonal model was obtained from the conduction velocity θ observed in a 1d preparation of the R-D model (Equation 4) through the relation
For instance, with β = 1,000 cm−1, σ = 1.5 mS/cm we observed θ = 0.075935 cm/ms, the latter being 4-digit accurate, hence α = 1.961 cm ms−1 mS−1/2. From the same 1d preparation we also extracted the template action potential
The propagation was triggered differently in the two models. In the eikonal model, we assumed that the pacing region was activated with an initial time of 1.12 ms, the latter obtained from the monodomain simulation. This delay is associated to the time required by the current injection to start the propagation.
In order to minimize the numerical error with respect to the modeling error, we performed the comparison on a uniform grid with voxels of side 0.05 mm. The eikonal model and the monodomain equation were solved only in the tissue slab. The forward bidomain problem was eventually solved including also the bath, and injecting the transmembrane currents in the tissue slab. The lead field, depicted in the right panel of Figure
The results are reported in Figure
We computed three ECGs in total: one with the R-D model, one with the proposed model, and one using the lead-field approach with the activation map obtained from the R-D model. The third was computed to evaluate the modeling error introduced by assuming a template action potential. The ECGs were very similar in shape and amplitude. In the terminal part of the QRS complex from the R-D model we observed a deeper inflection than in the other cases, and this seemed to be due to the assumption (Equation 13), since this was absent also in the case of the R-D model with lead field.
Finally, we performed the same comparison for a coarser grid resolution, observing a longer QRS duration in the R-D model. This is consistent with the fact that the finite-difference scheme underestimates the conduction velocity at coarser resolution (Pezzuto et al.,
Similarly to the previous section, we conducted an extensive comparison between the proposed model and the reference model on six patient-tailored anatomies, assuming the same parameters for both the models. Patients included in this study had a clinical indication for CRT; the ECG, cardiac resonance imaging, and an electrophysiological study, including 3D electro-anatomical mapping were clinically indicated. They were performed as part of the work-up for device selection, i.e., to evaluate scar location and extent, inducibility of ventricular tachycardia, and pacing site selection. The study was performed in compliance with the Declaration of Helsinki. The institutional review board approved the study protocol, and all patients gave oral and written informed consent for each investigation.
Characteristics of the heart meshes are listed in Table
1 | 165 × 91 × 119 | 191,734 | 151,959 | 3 |
2 | 205 × 107 × 151 | 292,585 | 235,502 | 1 |
3 | 186 × 123 × 126 | 289,252 | 230,968 | 1 |
4 | 178 × 117 × 147 | 350,584 | 290,166 | 2 |
5 | 156 × 91 × 108 | 222,951 | 187,534 | 1 |
6 | 238 × 145 × 162 | 493,006 | 413,252 | 4 |
The following electric conductivities were used for all the patients: torso 2 mS/cm; blood 6 mS/cm; lungs 0.5 mS/cm; skeletal muscle 3.55 mS/cm in the tangent plane and 0.44 mS/cm in the radial direction. The conductivities in the myocardium differed per patient as they were tailored to fit the measured ECG and activation sequence. A tuning was necessary for the activation time at the early activation sites (EASs) and the CV scaling parameter α, not present in the reaction-diffusion model. Left bundle branch block (LBBB) was assumed for all patients, and was modeled by placing the EASs only in the RV. The number of EASs for each patient is reported in Table
The initiation of the activation in the R-D model was performed by injecting a transmembrane source current of 200 mA/cm3 for 2 ms on 1 mm3 of tissue centered at the corresponding EAS location. The currents can be applied at different times. In the eikonal solver, the activation time of the EASs was enforced as boundary condition and shifted in time in order to mimic the delay observed in the R-D model between current injection and depolarization.
The activation times were compared in terms of maximum absolute error quantiles at the nodes and specific timing markers of cardiological interest at the LV endocardium. In particular, we measured the Trans-Septal Time (TST) and the Total Activation Time (TAT), which are the time of initial and the latest time of activation of the LV endocardium, respectively, and the QRS complex duration (QRSd), defined as the last activation time. Results are reported in Table
1 | 155.48 (157.39) | 45.93 (39.10) | 109.54 (118.29) | 7.11 | 10.77 | 12.91 |
2 | 184.20 (225.03) | 20.15 (19.76) | 136.58 (136.66) | 5.02 | 8.95 | 11.99 |
3 | 137.95 (137.92) | 42.04 (40.16) | 87.98 (88.67) | 2.58 | 5.06 | 8.06 |
4 | 143.15 (145.41) | 38.47 (36.85) | 104.67 (107.94) | 3.76 | 6.15 | 9.15 |
5 | 144.64 (146.22) | 32.01 (31.42) | 105.84 (109.17) | 2.64 | 4.52 | 6.83 |
6 | 165.60 (165.58) | 30.04 (28.30) | 128.40 (128.99) | 2.88 | 5.02 | 7.72 |
The overall propagation direction was captured for all the patients. The total activation time differed by 3 to 5 ms, being higher in the reaction-diffusion model. Thus, in the latter the excitation front was slower. The scaling parameter α was roughly 2 cm ms−1 mS−1/2 in all cases, with a variability within 1 %.
The isochrones were very similar in shape. The absolute error was well-distributed in the anatomy, although we observed a slightly larger error at the apex and in the LV. The EASs regions also differed, with the eikonal solution more ellipsoidal than the R-D model.
On average, we found that 90% of the nodes had an absolute error lower than 10 ms. QRSd, TAT and TST were very similar in all the cases, with the exception of patient 2, where QRSd with the R-D model was very high. The difference was due to the lower resolution geometry employed in the eikonal model (1 vs. 0.2 mm of the R-D model), which may introduce artifacts in the propagation. In patient 2, the anatomy included a thin, late-activated ventricular bundle stemming from the posterior LV free-wall, observed in MRI. In the downsampled geometry, this bundle was poorly reproduced resulting in a shorter activation time.
For the sake of completeness, in Figure
The ECGs were computed according to the standard 12-lead ECG definition (Malmivuo and Plonsey,
The ECGs for all the patients were computed from the activation time provided by the eikonal model and then compared against the bidomain solution. We adopted the action potential template in Equation (20) with uniform parameters in space with values from Table
The error analysis was limited to the QRS complex, since T-waves cannot be properly modeled by the eikonal approach. We compared the ECGs computed with the R-D model and the eikonal model in terms of: difference in maximum amplitude, area and positivity range, correlation, and
where
The ECGs computed with the SM and FM were very similar (not shown), with a small difference in amplitude in the QRS complex (FM was lower). The QRS complex was correctly captured in most of the leads for all the patients except for patient 6, where we observed large deviations in limb leads and moderate differences in precordial leads, whilst the activation map was fairly similar between the two approaches. This patient has a scar in the LV free wall that was modeled in the R-D model as purely passive. We suspect that in the downsampled geometry adopted by the eikonal model the scar was particularly jagged, hence affecting the ECG. In some cases our model provided a slightly larger amplitude in the signal than the bidomain model. Only in one case, for lead V5 of patient 5, we observed discordant QRS complexes. The discrepancies were particularly marked in leads with small amplitude. For patient 2, the discrepancy in the QRS duration was not present when calculated according to the ECG rather then activation map. As reported above, this was due to a ventricular bundle that did not affect the ECG. In general, the correlation between the bidomain model and the proposed model was very good in most of the leads, and the error lower than 0.5 mV.
The runtimes of the proposed eikonal solver and ECG solvers on the GPU were tested for six patient-specific geometries with 1-mm resolution.
We tested our code on two different Nvidia GPUs: a low-end (LE) GeForce GT 650M on a laptop (384 cores, 950 MHz clock, 0.73 Tflops of theoretical peak performance), and a high-end (HE) GeForce GTX 1080 (2,560 cores, 1,607 MHz clock, 8.2 Tflops of peak performance) on a local cluster node.
In Figure
HE vs. LE GPU | FIM | 5.1 | 5.9 | 5.3 | 6.1 | 5.2 | 6.2 | 5.6 |
SM | 17.3 | 17.6 | 17.8 | 17.9 | 18.1 | 19.4 | 18.0 | |
FM | 1.4 | 2.5 | 1.6 | 1.7 | 1.5 | 3.0 | 1.9 | |
FM vs. SM | HE GPU | 4.2 | 7.3 | 6.1 | 7.2 | 4.8 | 12.2 | 7.0 |
LE GPU | 53.8 | 51.3 | 67.7 | 76.3 | 57.1 | 79.5 | 64.3 |
A critical issue in computational cardiac electrophysiology is the complexity of the reaction-diffusion model. It is nearly impossible to consider this model for inverse problems, unless appropriate model reduction techniques are employed (Gerbeau et al.,
The proposed method compares adequately to the bidomain model, in terms of accuracy of the delivered solution. The average modeling error in the patient-specific context was roughly 15 % for the activation map (Table
The T-wave was almost never captured by our model. This was expected, since our model assumed that the repolarization follows exactly the depolarization, which was not the case in the reaction-diffusion model.
Finally, scalability was excellent for SM, and very good for the FIM and FM. Moreover, the SM and the FM provided very similar QRS complexes for all the patients and with different resolutions, but the FM is two orders of magnitude faster than the SM.
The eikonal model can be much faster than a reaction-diffusion equations for three reasons. First, the computational domain is reduced by one dimension, because the dependent variable is no longer a function of time. Second, the mesh can be coarser (1 vs. 0.1 mm), which results in much shorter computation times and much lower memory requirements. Last, it does not require computation of the ionic currents, which takes most of the computational time in the reaction-diffusion solution. Pullan et al. (
The numerical comparison on a simple geometry has shown that the eikonal equation is able to approximate the activation times of a bidomain model very accurately, and that there is no advantage in using the full bidomain formulation when modeling the excitation wavefront propagation if the conduction velocities are given.
The disadvantages of the eikonal model are obvious: the effects of ionic currents on conduction, for example partial refractoriness, cannot be simulated at all, simulation of reentrant arrhythmia is very difficult, and there is no calcium transient available for a realistic coupling to a mechanical model. Although several alternative eikonal-like models have been proposed in the literature aiming to improve the accuracy in such situations, such as fibrillation (Herlin and Jacquemet,
We believe that most of the discrepancy in the activation map between the R-D model and the proposed model is explained by the numerical error, and only minimally by the physiological phenomena that the eikonal model neglects.
The bidomain equation requires a mesh resolution that is comparable to the excitation front thickness, which depends on the conductivity in the propagation direction (Pezzuto et al.,
The eikonal solver was affected by a 5 % error at 1 mm resolution on a simple geometry. Moreover, the activation time was generally overestimated by the solver, similarly to the slower propagation observed in the reaction-diffusion model when solved with a finite-difference scheme. However, it is hard to tell whether the differences between the reaction-diffusion model and the proposed model are due to numerical errors or modeling errors. Very likely it is a mix of the two.
In our opinion, the eikonal model may have been more accurate than the reaction-diffusion model in the presented experiments. A comparison of the two models in Section 3.2 reported a modeling error of about 6 %. A small modeling error is also reported by others for an idealized LV geometry with physiologically motivated fiber architecture (Colli Franzone et al.,
The FIM ran on the HE-GPU on average only 5.6 times faster than on the LE-GPU. This is about half the ideal scaling. This may be explained by the fact that, at each iteration, only a fraction of the threads are active. In the FIM, only the activation times at the active nodes are evaluated, and hence only a few thousands of threads are running concurrently. A lower number of threads implies a lower number of active CUDA blocks and warps, and this may cause a paucity of warps prepared to hide the latencies of the long-lasting operations. The critical parts are the initial phase of the excitation process when the front starts to propagate from the EASs and the final phase when the last cells are activated. In these two phases only tens or hundreds of threads evaluate the solution and the GPU resources are not fully exploited.
Scalability of the FM was also not optimal. The FM ran on the HE-GPU on average only 4.5 times faster than on the LE-GPU. This means that the performance of the FM is more limited by the internal structure of the kernel than by the number of available GPU cores. The kernel generates one thread per cube in the computational domain. Warp divergence arises because individual threads compute their contribution only when the excitation front intersects the cube. With respect to the scaling, not the warp divergence itself is problematic, but rather the fact that different warps may require different numbers of clock cycles to finish their execution. The higher the number of time instants for which the warp thread evaluates the solution, the higher the warp execution time. Due to unbalanced execution times of different warps in the computational domain the linear scaling of the FM was broken.
We proposed a combination of an eikonal model for action potential propagation and an ECG simulation method based on lead fields. It can simulate an activation sequence and ECG in about 3 s on a GPGPU desktop platform (less than a second on HE-GPGPU). This method is suitable as a component of an inverse electrocardiographic model in an HPC context, but may in the near future also become practical within interactive ECG simulation tools. While we based our work largely on classical ideas, we believe that the following aspects are novel.
We proposed a new method to compute the ECG, based on the marching cubes method, which allows for several orders of magnitude speedup.
We compared the results of the eikonal model to those of a monodomain reaction-diffusion model, pointing out the difficulties of such a comparison.
We compared the simulated ECGs with those of a bidomain torso model and identified the most important causes of differences.
We discussed how the proposed methods can be implemented on a GPGPU.
We showed that the proposed methods can simulate a highly realistic ECG in a few seconds.
The ideas behind this work were initially proposed by MP, who also provided the numerical experiments with propag-5 and the patient's anatomies. SP designed the numerical method, drafted the implementation and performed the simulations. PK implemented the GPGPU code and optimized it. He wrote the initial version of the manuscript, then extended and revised by SP and MP to the present state. FP contributed on several occasions to the manuscript by reviewing it. AA and RK respectively reviewed the clinical outcome and numerical aspects of the manuscript.
This work was supported by grants from the Swiss National Supercomputing Centre (CSCS) under project IDs s397, s598, and s669. PK was supported by the Sciex-NMS fellowship “CATION,” project code 14.158. The authors gratefully acknowledge financial support by the Theo Rossi di Montelera Foundation, the Metis Foundation Sergio Mantegazza, the Fidinam Foundation, and the Horten Foundation to the Center for Computational Medicine in Cardiology. The authors also thank the Swiss PASC initiative (Platform for Advanced Scientific Computing) for their support. This work was also partially supported by a restricted grant of Biologic Delivery Systems, Division of Biosense Webster a Johnson & Johnson Company.
FP has received research grants from Medtronic, St. Jude Medical, Sorin, MSD, and Biotronik. AA has been a consultant to Medtronic, Boston Scientific, LivaNova and St. Jude, and has received speakers' fees from Medtronic, Boston Scientific and LivaNova. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The authors sincerely thank Hanspeter Fischer of Biosense Webster for his contribution in data acquisition and for making the data available for offline analysis.
The Supplementary Material for this article can be found online at: