Abstract
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.
1. Introduction
The electrocardiogram (ECG) is one of the most common tools in present-day medicine, yet its relation with the molecular biology of the heart is still poorly understood. The ECG witnesses the collective activity of about a million current-generating transmembrane proteins in each of the heart's muscle cells (Hille, 2001). Many of these proteins have been identified and their actions have been captured in mathematical models that predict their collective behavior on the scale of a cell (Noble and Rudy, 2001). By coupling millions of these membrane models one can create a model of whole-heart electrophysiology. Such models generate crucial insights in the functional effects of molecular-level changes, allowing for example to predict dangerous side effects of new drug designs (Passini et al., 2017) or to understand how cardiac ion-channel mutations influence cardiac rhythm disorders (Gima and Rudy, 2002). Moreover, from their results one can compute the corresponding ECG and predict how lab results on subcellular components would translate to everyday practice (Hoogendijk et al., 2010; Keller et al., 2012; Zemzemi et al., 2013).
Such realistic models are large and, when run on a single processor, would take days to simulate just one heartbeat. Fortunately the problem can be expressed in such a way that the work may be spread over many processors with little communication between them. Therefore, these computations are said to scale very well, meaning that they run almost twice as fast every time the number of processors is doubled (Vázquez et al., 2011). This makes them suitable for use on large-scale parallel computers, allowing models to run in nearly real time (Niederer et al., 2011b; Richards et al., 2013).
Simulation of a realistic ECG from the results of such a numerical heart model is much harder, because the electrical current generated by the heart meets a different conductivity at each point in the torso. As a result, each point influences the potential everywhere else, so to find the potential anywhere one must solve it everywhere at the same time.
Numerically this means that a large system of linear equations must be solved, one for each point in the torso model. These problems are harder when they are larger and require frequent communication between the processors in a parallel computer. This means that they cannot be solved much faster by using more processors. Therefore, ECG computation is becoming a bottleneck, limiting both the speed and the spatial resolution of our models.
To avoid this problem many researchers have used simplified torso models, resulting in a less accurate ECG. A solution that can avoid such a sacrifice is to simulate the ECG using an electrocardiographic concept named a lead field. This allows the problem to be split into a hard (poorly scaling) part and an easy (well scaling) part. The hard part is solved only once for each ECG lead, while the easy part is run repeatedly for each time step in a simulation and for multiple simulations on the same geometry. This approach has been used by several authors, but generally with simplified heart models (Pezzuto et al., 2017) or, again, with simplified torso models (Horacek, 1973; Miller and Geselowitz, 1978; Mailloux and Gulrajani, 1982; Aoki et al., 1987).
The purpose of this paper is to show that a lead-field approach can greatly improve scalability in a high-performance computing (HPC) context without sacrificing accuracy. This is not obvious, because the method requires a large set of transfer coefficients (the lead field) to be stored between the two phases of the computation. The efficiency of the method depends on the accuracy with which the lead field must be computed and the degree to which it can be downsampled without affecting the accuracy of the ECG too much. Finally, to provide answers to these questions an accurate reference solution is needed.
Using a reference solution computed on a full torso model at 0.2 mm resolution this study shows that the lead field can indeed be downsampled enough to achieve an efficient and scalable computation, providing roughly two orders of magnitude speedup with negligible loss in accuracy.
The results of this study make it possible to build more realistic heart models with higher spatial resolution, without spending much more time to compute the ECG.
2. Methods
2.1. Model equations
The methods in this study are based on the bidomain model of cardiac electrophysiology (Miller and Geselowitz, 1978; Tung, 1978), on which most of the current modeling work in this area is based (Niederer et al., 2011a; Henriquez, 2014). The bidomain model is a continuum approximation of the heart muscle, which in reality consists of a network of interconnected muscle cells embedded in an extracellular matrix and other structures such as fibroblasts and capillaries. The bidomain model approximates this as two co-located spaces: the intracellular domain, consisting of the interior of the cells and the gap junctions that connect them, and the extracellular domain, consisting of everything else.
The two domains are characterized by conductivity tensors Gi and Ge, respectively. Their values at each point in the model depend on the fiber direction and account for the partial volume occupation of the two domains. In addition the parameters Cm and β determine the capacitance of the cell membrane and the amount of membrane per unit volume, respectively. The state variables of the model are the potential fields ϕi in the intracellular and ϕe in the extracellular domain, and a set of variables describing the state of the membrane model at each location. Using the auxiliary variable Vm = ϕi − ϕe and agreeing that all variables are functions of time and position we can express the bidomain model compactly as
where the term Cm∂tVm represents the capacitive transmembrane current, the function Iion the density of ionic current flowing between the two domains, and F is a nonlinear vector-valued function describing how the membrane state evolves. The pair of functions Iion and F constitutes the membrane model. Suitable boundary conditions are
on the boundary ΩA of the cardiac muscle and
on the torso boundary ΩT (Tung, 1978; Krassowska and Neu, 1994).
The electrical activity of the heart can then be simulated by integrating Equations (1), (2), and (3) under the boundary conditions (4) and (5) (Vigmond et al., 2002). This is known as a bidomain reaction-diffusion model. In this study a simplified version, a “monodomain” reaction-diffusion model, was used. This model can be derived by assuming that Gi and Ge are proportional (Leon and Horácek, 1991). Although this is a gross simplification the effect of this assumption is negligible for most purposes if the model parameters are well chosen (Potse et al., 2006; Nielsen et al., 2007; Bishop and Plank, 2011; Coudière et al., 2014). The monodomain model reads
The “monodomain conductivity tensor” Gm was computed as the series conductivity of the two domains, Gm = GiGe/(Gi + Ge). With this choice the resistance encountered by a current loop through the cell membrane is the same as in a bidomain model, so that also the conduction velocity of a propagating activation wavefront is almost the same.
An ECG potential V(t) at time t is the difference in ϕe between two locations on the body surface or, more generally, a linear combination
where ci are the relative contributions of the two or more electrodes and are the potentials at the corresponding positions. The coefficients ci must fulfill charge conservation, ∑ci = 0.
To compute ϕe we must return to the bidomain model. Equations (1) and (2) can be combined and reorganized to yield
This equation can be solved for ϕe in the whole torso at once from a given distribution of Vm. However, for the ECG we need to know ϕe at a few locations only. Therefore, it can be more efficient to use a Green's function of the operator ∇ · ((Gi + Ge)∇.) for each of these locations. Since an ECG lead is a linear combination of ϕe at two or more points it can also be represented directly by a linear combination of Green's functions. In electrocardiology such linear combinations of Green's functions are named lead fields (McFee and Johnston, 1953; Geselowitz, 1989; Colli-Franzone et al., 2000). A lead field is computed once for each ECG lead. It is then used to evaluate the ECG at each time step of the reaction-diffusion model and, as long as the conductivity parameters are not changed, can be re-used for multiple simulations. In terms of a lead field the ECG potential V(t) at time t is
where the integration is over the myocardium. In contrast to the solution of the full system (8) this calculation is simple and a priori highly scalable. The lead field can be computed as the potential field resulting from a unit current applied at the electrode locations (Geselowitz, 1989):
where the coefficients ci are as in Equation (7) and δ is Dirac's delta function. To avoid a scaling factor in (9) the total injected current must be unitary, ∑|ci| = 2.
2.2. Model geometry
In order to run tests on a relevant geometry a model of the heart and torso was used that had been created for a previous study (Kania et al., 2017). The methods to build this geometry, only tersely described before, were as follows. High-resolution cardiac and thoracic computed tomography (CT) images were obtained from a female patient in her thirties. Images were segmented automatically using the MUSIC software (IHU Liryc, Université de Bordeaux and Inria Sophia Antipolis, France), under supervision of an expert operator. The boundaries of the segmented volumes were expressed as triangulated surfaces and meshing errors were manually corrected using Blender (The Blender Foundation, Amsterdam, The Netherlands). The resulting surface mesh defined the volumes of the ventricular myocardium, left and right cavities with parts of the great vessels, lungs, and the whole body. To define hexahedral meshes for the computations the surfaces were overlaid with a 3D cartesian mesh whose elements were assigned types according to the surfaces in which they were contained. The bones were also segmented and meshed but not included in the simulations. The atrial myocardium was not segmented.
The heart mesh was processed to define subendocardial and subepicardial layers and fiber directions using the rule proposed by Beyar and Sideman (1984), as previously described (Potse et al., 2006). The torso mesh was similarly processed to define a layer of 1 cm thickness directly under the skin as skeletal muscle and to define a sheet direction in this layer. Since the true fiber directions of the skeletal muscle layer are too complex to account for the model muscle simply had a low conductivity in the radial direction and a high conductivity in all circumferential directions (Table 1).
Table 1
| Material | Volume (mL) | σiL | σiT | σiC | σeL | σeT | σeC | β |
|---|---|---|---|---|---|---|---|---|
| Myocardium | 110 | 3.0 | 0.3 | 0.3 | 3.00 | 1.20 | 1.20 | 800 |
| Body | 16,482 | 0 | 0 | 0 | 2.00 | 2.00 | 2.00 | 0 |
| Blood | 236 | 0 | 0 | 0 | 6.00 | 6.00 | 6.00 | 0 |
| Lung | 4,352 | 0 | 0 | 0 | 0.50 | 0.50 | 0.50 | 0 |
| Muscle | 5,605 | 0 | 0 | 0 | 3.55 | 3.55 | 0.44 | 0 |
Tissues used in the simulations together with the volumes they occupy in the torso model, the conductivity parameters σ (in mS/cm), and β (cm−1); the subscript “i” stands for intracellular, “e” for extracellular, “L” for longitudinal, “T” for transverse (within a tissue sheet), and “C” for across-sheet.
During the thoracic scan the patient was wearing a vest with 252 embedded electrodes (Tilt et al., 2013; Cochet et al., 2014). The locations of these electrodes were extracted from the CT data using software provided by the manufacturer of the vest. In addition the locations of the 9 standard ECG electrodes were determined by referring to the bone mesh, and two electrode locations on the hips were chosen. The surface mesh with electrode positions is illustrated in Figure 1.
Figure 1
2.3. Spatial discretization
Spatial discretization was done using a finite-difference method. Differential operators of the form ∇·(G∇.), where G is any of the conductivity tensor fields employed, were computed using an expression proposed by Saleheen and Ng (1997). This expression assumes that G is constant on elements and that potentials are defined on the nodes of the mesh. It produces a 19-point stencil that takes anisotropy and inhomogeneities into account. The simulation code read its geometry in terms of elements, and created a node mesh, assigning node types such that all corners of a myocardial element would have myocardial nodes. In order to treat myocardial boundaries correctly, the β value of each node was the average of those associated with the 8 elements around it, which was zero for non-myocardium (Potse et al., 2006).
2.4. Simulation of cardiac activity
To prepare input data for ECG simulation propagating activation was simulated using the monodomain reaction-diffusion model (6) using the membrane model of Ten Tusscher and Panfilov (2006) for the functions F and Iion. A uniform time step of 10 μs was used. At each time step the code
evaluated the diffusion current β−1∇·(Gm∇VmΔVm),
communicated the diffusion current across domain boundaries,
integrated the membrane status variables ,
evaluated , and
integrated Vm.
After each 100 time steps results were written to file. Simulations were run on a heart mesh at 0.2 mm resolution. Tissue parameters determining Gm and β are listed in Table 1. Gating variables were integrated with the method of Rush and Larsen (1978) and all other variables with a forward Euler method.
Activation was started with a single stimulus at one location, at the beginning of the simulation. Seven simulations were run, each time with the stimulus at a different location. Simulations covered 500 ms to include the full depolarization and repolarization of the ventricles.
2.5. ECG simulation
The ECG was computed with several methods:
FSF, the fine-mesh full solution solved the full system (8) for given Vm on a heart-torso mesh with 0.2 mm resolution. This was an exceptionally large computation requiring 3.3·109 mesh nodes and 12 TB memory. It was combined in a single run with the integration of the monodomain reaction-diffusion model (6). Solutions for ϕe were computed after each 100 time steps.
FSC, the coarse-mesh full solution solved an alternate form of Equation (8) on a heart-torso mesh with 1 mm resolution (Potse and Kuijpers, 2010). In this case the equation read
where Iw is a projection of the term ∇·(Gi∇Vm) from a 0.2 mm resolution heart mesh onto a 1 mm resolution torso mesh. Each coarse-mesh node received contributions from a cube-shaped area including all fine-mesh nodes within the up to 8 coarse-mesh elements around it, with higher weights attributed to nearby nodes, as in a trilinear interpolation: Let Δx, Δy, Δz be the number of fine-mesh edges between a coarse-mesh node and a fine-mesh node along the x, y, and z axis, respectively. Then the contribution of the fine-mesh node to the coarse-mesh node was
The coarse mesh was constructed such that a myocardial fine-mesh node was always surrounded by 8 coarse-mesh nodes. Therefore, w added up to unity for each fine-mesh node and charge conservation was ensured.
For the FSC method the monodomain reaction-diffusion model (6) was integrated in a separate run which saved Iw to file. This method has been used routinely in several studies (Nguyên et al., 2015; Meijborg et al., 2016; Duchateau et al., 2017; Kania et al., 2017). The torso mesh in this case consisted of 2.7·107 nodes.
LF, the lead-field method evaluated the integral expression (9) in its discretized form. This took place during the reaction-diffusion simulation and on the same mesh, i.e., at 0.2 mm resolution, after each 100 time steps. Each component of ∇Vm was evaluated on model elements as an average of the differentials along 4 edges of the element. The conductivity tensor Gi was also evaluated on each element. For testing purposes the lead vector field ∇Z was evaluated at different resolutions. For this purpose the field was first downsampled by an external program, using a simple averaging of n × n × n elements, where n could be 2, 5, 10, or 25.
LFS, the lead-field method with selective downsampling was identical to the lead-field method except that the downsampling program took the tissue types of the elements into account. If any of the fine-mesh elements inside a coarse-mesh element E had a myocardial type, only fine-mesh elements with myocardial type were used in the average for E. The idea behind this was that ∇Z undergoes abrupt changes at the myocardial boundaries, and that it is more accurate to mix in a contribution from another myocardial area than, for example, one from the lung.
The notations LF(C, S) and LFS(C, S) will be used for the LF and LFS methods, respectively, with lead fields computed at a resolution of C millimeters and downsampled to a resolution of S millimeters.
2.6. Computation of lead fields
To prepare the lead fields Z for the ECG computation the system (10) was solved for each lead. This was done once with a torso model at 1 mm resolution and once with a torso model at 0.2 mm resolution. Like the FSF, the latter calculation was exceptionally large and was only intended to provide reference values, to test the hypothesis that 1 mm resolution suffices for such calculations.
In either case 266 lead fields were computed: the 12 standard ECG leads, and one lead for each of the 252 vest electrodes and 2 hip electrodes referenced against Wilson's central terminal (the average of the two arm electrodes and the left leg electrode).
The computed lead fields Z were stored in files. A dedicated program computed ∇Z and downsampled it using the two methods described in section 2.5, i.e., with and without consideration of the tissue types of the elements. The field computed at 0.2 mm resolution was downsampled by the factors 2, 5, 10, and 25 to obtain resolutions of 0.4, 1, 2, and 5 mm. The field computed at 1 mm resolution was downsampled by the factors 2 and 5 to obtain resolutions of 2 and 5 mm.
2.7. Testing protocol
ECGs were simulated using each of the 4 methods described in section 2.5 and, for the methods based on lead fields, at each of the resolutions mentioned in section 2.6.
The ECG potentials V were compared to a reference ECG Vref in terms of three measures: maximum, root-mean-square (RMS), and relative difference (RelDif) (van Oosterom, 2001; Tysler et al., 2007), defined as
where the index t ranges over all 500 samples and the index n ranges over all 266 leads. For the 252 vest leads the dependence of the error values on the position of the positive electrode was investigated.
The effect of the ECG computation on the run time of a reaction-diffusion model was investigated and the scalability of the 4 methods was investigated by running tests on 16, 32, …, 512 nodes of a Bull cluster. Each of these nodes was equipped with two 14-core Intel Xeon E5-2690 processors with 2.6 GHz clock frequency and 64 GB memory. Accuracy results are reported as averages over the 7 activation sequences. Performance tests were carried out 5 times to report average values and standard deviations of run time.
2.8. Numerical methods
Simulations were performed using the Propag-5 software (Krause et al., 2012), to which new code was added to compute a lead field-based ECG on the fly during a simulation of the heart, and to facilitate the computation of the lead fields themselves. Like its predecessor Propag-4 (Potse et al., 2006), the software uses a structured mesh, but stores information only for elements and nodes that are relevant for the computation: only myocardium for a monodomain model, and only conducting material for a bidomain model. As discussed by Krause et al. (2012) Propag-5 uses a hybrid MPI/OpenMP parallellization scheme. Using a naive temporary partitioning of the domain the code reads the geometry in terms of elements and creates a node mesh using rules that ensure consistency with the scheme discussed in section 2.3. It then uses the ParMetis library to partition this mesh in parallel and creates a definitive domain partitioning for the computations. This fully parallel workflow allowed it to load and partition a mesh with over 3 billion nodes.
Because in some of the computations the model size exceeded the maximum value of a signed 32-bit integer, Propag was compiled with a 64-bit integer type for global indices. The PetSC (Balay et al., 2017) and Parmetis libraries which Propag uses were compiled entirely with 64-bit integers because they do not have a distinct type for global indices.
The linear systems (8), (10), and (11) were solved with a biCGStab solver (van der Vorst, 1992) with a BoomerAMG preconditioner from the Hypre package (Henson and Meier Yang, 2002; Falgout et al., 2017). The solver terminated when the norm of the error term was 10−8 times smaller than the norm of the right-hand side. Multigrid preconditioners such as BoomerAMG are very powerful and well-suited for large bidomain problems (Sundnes et al., 2002; Weber dos Santos et al., 2004; Austin et al., 2006) so that the solver typically needs only a handful of iterations, in contrast to the problematic convergence observed on large models with an incomplete-LU preconditioner (Potse et al., 2006).
3. Results
An example of a computed lead field is shown in Figure 2. This field was computed and stored at 1-mm resolution. The figure shows how the field suddenly changes direction and magnitude at lung boundaries. There is a slight left-right asymmetry because the highly conductive cardiac cavities concentrate the field on the left side of the thorax.
Figure 2
The computed depolarization sequences of the 7 simulated heart beats that were used for ECG computation are shown in Figure 3.
Figure 3
Potentials computed with a full-torso solution from beat 5 are shown in Figure 4. They are about 10 times larger in the myocardium than near the body surface.
Figure 4
3.1. Lead-field ECG compared to full solution
To establish that the lead-field and full solution methods produce the same results, simulated ECGs were compared between the LF(1, 1) and FSC methods. Averaged over the 7 simulations, RelDif was 0.0016, RMS error 0.3 μV, and maximum error 4.6 μV, while ECG amplitudes were in the order of 1 mV.
Analogously, a single ECG was compared between the LF(0.2, 0.2) and FSF methods. In this case the differences were slightly smaller: RelDif was 0.0014, RMS error 0.2 μV, and maximum error 2.6 μV.
3.2. Effect of resolution
To determine the effect of lead-field resolution on ECG accuracy, 7 different activation sequences were simulated with a monodomain reaction-diffusion model and ECGs were simulated on the fly using a lead field. This was done for the lead fields computed at 0.2 and at 1.0 mm and all downsamplings thereof, both with the LF and with the LFS method. The resulting ECGs were compared to a reference ECG.
The results are shown in Figure 5. In Figure 5A errors are shown using the ECG computed with LF(0.2, 0.2) as the reference. For the fields subsampled from those computed at 0.2 mm resolution, differences are seen to increase roughly linearly with the stepsize of the lead field. The LFS method resulted in smaller differences. Results obtained with the field computed at 1.0 mm resolution and downsamplings differed from the reference solution with little dependence on the sampling level. Figure 5B shows that this dependency is recovered when ECGs computed with LF(1, 1) are used as the reference.
Figure 5
The relatively large influence of the spatial stepsize in the lead-field computation suggests that differences in model geometry dominate the error. Indeed, the difference between full solutions at 0.2 and 1.0 mm, computed only for one simulation, had a RelDif of 0.10, RMS error 12 μV, and maximum error 0.15 mV, which are very similar to the differences between LF(1, 1) and LF(0.2, 1) in Figure 5A.
To find out at which locations in the model the lead fields computed with LFS(0.2, 1) and LF(1, 1) differed, the L2 norm of the difference between the two vector fields was computed for all elements. Large differences were found to occur at locations where the fiber direction was highly variable. One such location, at the inferior septal junction, is illustrated in Figure 6. It is compared with a measure of variability in fiber direction in the underlying anatomy files, computed as
where is the fiber direction in the coarse-mesh element and are the fiber directions in the corresponding fine-mesh elements. The absolute value, denoted as |.|, was taken because the orientation of the direction vector is irrelevant.
Figure 6
In Figure 7 a few ECG leads are compared between different computation methods. In Figure 7A full solutions at 0.2 and 1.0 mm are compared. At the coarser resolution the ECG appears more fractionated; this is particularly visible in lead III. As discussed above, the RelDif between these ECGs was 0.10. In Figure 7B the same full solution at 0.2 mm is compared with an ECG computed with LFS(0.2, 2). Despite the 10-fold downsampling of the lead field the traces are visually identical; the RelDif was 0.02. Thus, an ECG computed with a lead field downsampled to 2 mm resolution is more faithful than a full solution at 1 mm resolution, when compared to a solution at 0.2 mm.
Figure 7
3.3. Performance
Table 2 shows how ECG computation with lead fields at different resolutions affects the run time of a typical simulation. The data in each row were obtained from 5 simulations of 500 ms activity with a reaction-diffusion model at 0.2 mm resolution, run on 32 compute nodes (896 cores). The table separates initialization time, ECG computation time, and simulation time (including ECG computation but excluding initialization). For lead fields at 0.2 and 0.4 mm resolution the initialization time is of the same order of magnitude as the simulation time, due to the time it takes to read the lead fields from file (141 and 53 GB in these cases). The time for ECG computation itself ranges between 4 and 5 % of the simulation time, slightly reducing with the lead-field resolution. At 1 mm resolution the memory accesses related to ∇Z (for 266 leads) are similar to those for Gi∇Vm so a further reduction would not be expected. At 0.2 mm resolution the ECG computation is faster than at 0.4 mm, likely because in this case the lead field has the same resolution as the reaction-diffusion model and the code then avoids an index conversion.
Table 2
| res | sim | ECG | init |
|---|---|---|---|
| 0.2 | 172.6 ± 0.7 | 8.3 ± 0.0 | 216.2 ± 3.6 |
| 0.4 | 179.2 ± 0.4 | 9.8 ± 0.1 | 103.3 ± 1.9 |
| 1.0 | 173.5 ± 1.8 | 8.1 ± 0.1 | 28.2 ± 1.5 |
| 2.0 | 171.8 ± 1.4 | 6.9 ± 0.3 | 25.8 ± 1.0 |
| 5.0 | 171.3 ± 2.2 | 6.2 ± 0.2 | 14.5 ± 0.7 |
Time required for LF-based ECG computation during a reaction-diffusion simulation of 500 ms.
res, lead-field resolution in mm; sim, total simulation time; init, initialization time. Time is given as average ± standard deviation over 5 simulations, in seconds.
Figure 8A, shows how the computation times scale with the number of cores used for a single lead-field resolution of 1.0 mm. The reaction-diffusion simulation and the ECG computation scale well. Initialization time increases with the number of cores, due to increasing communication for mesh distribution and data input. Tests with higher and lower lead-field resolutions, not presented in the figure, showed that the initialization time was highly variable and had no clear relation with the resolution (and thus the storage size) of the field. Rather, the number of collective read operations seemed to be determining.
Figure 8
The black trace in Figure 8A shows the scaling of a full solution (FSC method). It is over 2 orders of magnitude slower than the lead-field ECG and stops scaling at 7,168 cores.
Figure 8B shows how the ECG computation time scales with the number of nodes for all tested values of lead-field resolution. Lead-field resolution is seen not to affect the scaling with the number of cores. Generally the time decreases slightly with decreasing resolution but, as in Table 2, the computation at 0.2 mm was faster than the one at 0.4 mm.
4. Discussion
This study shows that a lead-field approach is an attractive solution for ECG simulation on (large) parallel computers whenever the number of ECG leads is smaller than the number of samples. It is about 100 times faster than a full solution, scalable to more than 104 cores, and does not cause a significant loss in accuracy. Lead fields can be stored at a resolution as low as 2 mm, meaning that they do not use excessive disk space even for a few hundred leads.
4.1. Previous work on lead fields
The concept of lead fields was initially proposed by McFee and Johnston (1953) as a method to understand how ECG leads “view” the heart. Their purpose was in the first place to design leads that would be better in the sense that their fields would be more uniform inside the heart muscle (McFee and Johnston, 1954). Later the idea has been adopted for the purpose of accurate numerical simulation of the ECG (Geselowitz, 1989) and even local electrograms inside the heart (Colli-Franzone et al., 2000; Western et al., 2015).
The idea to use lead-field methods for ECG simulation has been widely adopted. While the very earliest studies did not use them, for example because they computed only a small number of potential distributions (Gelernter and Swihart, 1964) or because a full solution required less memory (Barr et al., 1966; Barnard et al., 1967), numerous studies are based on some form of lead fields or transfer coefficients between Vm in the heart and ϕe on the body surface (Horacek, 1973; Miller and Geselowitz, 1978; Mailloux and Gulrajani, 1982; Aoki et al., 1987; Lorange and Gulrajani, 1993; Trudel et al., 2004).
Mailloux and Gulrajani (1982) and further work from the same group (Lorange and Gulrajani, 1993; Trudel et al., 2004) used transfer coefficients that are mathematically identical to lead fields. Their transfer coefficients were computed with a boundary element model (BEM) which accounted for heterogeneity of the torso, but not for anisotropy. They found that they needed <100 regions to define these coefficients, likely because their model was isotropic. In the anisotropic model used here the lead field changed considerably through the wall, requiring a much higher though not prohibitive resolution. Jacquemet (2015, 2017) evaluated the performance of the same (BEM-based) method on a reaction-diffusion model of the human atria and found that 1,000 regions sufficed for a 1% accuracy.
Boulakia et al. (2010) reported that an ECG simulation based on a transfer matrix was 60 times faster than solving a coupled heart-torso problem. They were using a finite-element model with about 1 million tetrahedra whose sizes gradually increased from the heart to the torso surface, and a serial code. Despite the obvious differences in methods the speedup was very similar to what was found in the current study.
Electrocardiographic inverse modeling studies that used volumetric transmembrane potentials or current dipoles as their source models have also used transfer coefficients that are similar to lead fields (Liu et al., 2006; Wang L. et al., 2013).
4.2. Other methods to compute the ECG
Many other studies have used full torso solutions to obtain the ECG from a reaction-diffusion model using finite-difference (Potse et al., 2009; Hoogendijk et al., 2010; Meijborg et al., 2016; Chamorro-Servent et al., 2017) or finite-element models (Lines et al., 2003; MacLachlan et al., 2005; Boulakia et al., 2010; Keller et al., 2010; Zemzemi et al., 2015; Janssen et al., 2017). In some cases this was done because intracardiac electrograms in a torso-coupled heart were also simulated (Hoogendijk et al., 2010; Meijborg et al., 2016). The ECG is then a free by-product.
An interesting alternative is a mixed approach in which anisotropic regions such as the heart and skeletal muscle are handled with finite elements and isotropic regions with boundary elements (Pullan and Bradley, 1996), resulting in fewer degrees of freedom than a complete volume discretization.
There is a considerable body of literature dedicated to the problem of solving body-surface potentials from epicardial (extracellular) potentials (Barr et al., 1977; Pilkington et al., 1987; Stenroos and Haueisen, 2008), which has found an application in cardiac inverse modeling (Greensite and Huiskamp, 1998; Ramanathan et al., 2004; Shou et al., 2008). A formulation in terms of transmembrane potentials on the (endocardial and epicardial) surface of the cardiac muscle is possible if equal anisotropy of the intracellular and extracellular domain is assumed (Geselowitz, 1989; van Oosterom and Jacquemet, 2005) and is also used to solve cardiac inverse problems (Oosterhoff et al., 2016).
4.3. Strengths and limitations
ECG simulation based on lead fields is very fast and as scalable as a monodomain reaction-diffusion model. This makes it suitable for inclusion in the same model run on a large-scale parallel computer or a GPGPU, in contrast to full solutions, which would limit the scalability of the entire computation. This advantage is present whenever the number of ECG samples to be simulated exceeds the number of leads.
Lead-field methods can also be used to compute local electrograms in the heart but this may require a higher spatial resolution at least near the electrode (Colli-Franzone et al., 2000).
For detailed spatial mapping of potentials, either in the heart or on the torso surface, lead-field methods are less advantageous, as the number of locations might exceed the number of samples and may even be so large that the storage of the lead fields becomes a performance bottleneck. In such cases full solutions remain the method of choice and a relatively long solution time will have to be accepted. Although new developments in scalable preconditioners may improve the situation somewhat (Munteanu et al., 2009; Ottino and Scacchi, 2015), it is unlikely that full solvers will ever scale as well as an ECG computation based on lead fields.
It would also be challenging to use a lead-field approach in an electromechanical, deforming heart model. A lead field that would be deformed with the mesh might be a reasonable approximation but this has not been tested here.
The results of this study also suggest further improvements, in the first place the use of non-uniform mesh density for lead-field computation. Comparison of ECGs computed at 0.2 and 1.0 mm resolution showed that the latter had artefactual notches of about 0.05 mV amplitude in the QRS complex, due to misrepresentation of fiber orientation at locations where this orientation changed rapidly. This applied to both full solutions and lead-field ECGs. To avoid such artifacts one could try to ensure a smooth fiber orientation throughout the model (Bayer et al., 2012), but this can be challenging at the interventricular junctions, or whenever measured fiber orientations rather than rule-based orientations are used. The only alternative seems to be computation of the lead field with a mesh at the same resolution as the reaction-diffusion model inside the heart, and for improved efficiency a lower resolution elsewhere in the torso (Pullan and Bradley, 1996; Boulakia et al., 2010). While the computations could still be hard on a mesh with a wide variation in element size, the memory requirements would be much lower than the 12 TB reported here for the reference torso model.
Another possible improvement that would be relevant for very accurate computations with high-resolution lead fields is to develop suitable compression methods for lead-field data. Very likely the regularity of the field could be exploited by using fixed-point numbers in combination with spatial differentiation and a variable-length encoding.
In Figure 8A, a particularly unfavorable scaling of the initialization phase was shown for the propagation model with lead-field ECG. This was probably due to an issue with the collective reading operation in the MPI library that was used, but also to the fact that for this feasibility study little care had been taken to organize this efficiently—after all the specifications for this code depended on the outcome of the study. With these results in hand it should be possible to avoid this problem by using a more efficient storage format and organizing the read operation in a different way. The figure also shows that the FSC method takes an order of magnitude more time than the reaction-diffusion model. This difference is partly due to the small solver tolerance that was chosen for this study.
4.4. Applications
The use of lead-field methods simplifies the workflow for large-scale cardiac simulations, as it allows the ECG to be computed on the fly with very little overhead during a reaction-diffusion simulation on a mesh of the heart alone. Moreover, its high scalability allows the resolution of the models to be increased without causing a disproportional increase in the time needed for ECG computation.
The results of this study are not only relevant for work on large-scale computers but also for simulations on general-purpose graphics processing units (GPGPU). Reaction-diffusion simulations on GPGPUs have been reported by several groups (e.g., Bartocci et al., 2011; Neic et al., 2012; Mena et al., 2015; Kudryashova et al., 2017), recently even for a whole human heart model run on a desktop computer (Vandersickel et al., 2016). The strength of a GPGPU is that it provides thousands of parallel processors for the price of a single CPU. However, communication between these processors is a distinct weakness. With a method based on lead fields it is nevertheless possible to add rapid ECG computation to a model running on a GPGPU. Pezzuto et al. (2017) have recently reported such a method, though in combination with an eikonal model rather than a reaction-diffusion model.
In the context of ECG inverse models and model personalization a variety of methods has been reported ranging from infinite-medium potentials (Giffard-Roisin et al., 2017; Neic et al., 2017) to full-torso bidomain solutions (Wang D. et al., 2013). A lead-field approach could offer a solution that combines the speed of the former (if the computation of the lead field itself is excluded) with the accuracy of the latter. Only methods based on equivalent double layers (Geselowitz, 1992; van Oosterom and Jacquemet, 2005) offer more efficiency as they need to evaluate only the surface of the heart, but the price for this efficiency is that these methods neglect anisotropy. A lead-field approach combined with an eikonal-diffusion model for cardiac propagation (Konukoglu et al., 2011; Jacquemet, 2012; Neic et al., 2017) could soon be a practical solution for ECG inverse problems with an accuracy very close to the state of the art in forward modeling of the ECG.
5. Conclusion
Lead fields are a practical alternative for full-torso solutions when the number of ECG leads that need to be simulated is smaller than the total number of samples that will be calculated. The method is fast and highly scalable. Lead fields can be stored at a resolution as low as 2 mm without unacceptable loss of accuracy.
Statements
Author contributions
The author confirms being the sole contributor of this work and approved it for publication.
Acknowledgments
This work was granted access to HPC resources of CINES under GENCI allocation 2018-A0030307379. This work was supported by the National Research Agency, grant reference ANR-10-IAHU04-LIRYC.
The author thanks Dr. Michael Leguèbe and Dr. Emmanuelle Saillard for proofreading the manuscript, and Dr. Hubert Cochet for providing geometry data and proofreading.
Conflict of interest
The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.00370/full#supplementary-material
References
1
AokiM.OkamotoY.MushaT.HarumiK.-I. (1987). Three-dimensional simulation of the ventricular depolarization and repolarization processes and body surface potentials: normal heart and bundle branch block. IEEE Trans. Biomed. Eng.34, 454–462.
2
AustinT. M.TrewM. L.PullanA. J. (2006). Solving the cardiac bidomain equations for discontinuous conductivities. IEEE Trans. Biomed. Eng.53, 1265–1272. 10.1109/TBME.2006.873750
3
BalayS.AbhyankarS.AdamsM. F.BrownJ.BruneP.BuschelmanK.et al. (2017). PETSc Web Page. Available online at: http://www.mcs.anl.gov/petsc
4
BarnardA. C.DuckI. M.LynnM. S.TimlakeW. P. (1967). The application of electromagnetic theory to electrocardiology; II. Numerical solution of the integral equations. Biophys. J.7, 463–490.
5
BarrR. C.PilkingtonT. C.BoineauJ. P.SpachM. S. (1966). Determining surface potentials from current dipoles, with application to electrocardiography. IEEE Trans. Biomed. Eng.13, 88–92.
6
BarrR. C.RamseyM.III.SpachM. S. (1977). Relating epicardial to body surface potential distributions by means of transfer coefficients based on geometry measurements. IEEE Trans. Biomed. Eng.24, 1–11.
7
BartocciE.CherryE. M.GlimmJ.GrosuR.SmolkaS. A.FentonF. H. (2011). Toward real-time simulation of cardiac dynamics, in CMSB 2011: Proceedings of the 9th ACM International Conference on Computational Methods in Systems Biology (Paris: ACM).
8
BayerJ. D.BlakeR. C.PlankG.TrayanovaN. A. (2012). A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Ann. Biomed. Eng.40, 2243–2254. 10.1007/s10439-012-0593-5
9
BeyarR.SidemanS. (1984). A computer study of the left ventricular performance based on fiber structure, sarcomere dynamics, and transmural electrical propagation velocity. Circ. Res.55, 358–375.
10
BishopM. J.PlankG. (2011). Bidomain ECG simulations using an augmented monodomain model for the cardiac source. IEEE Trans. Biomed. Eng.58, 2297–2307. 10.1109/TBME.2011.2148718
11
BoulakiaM.CazeauS.FernándezM. A.GerbeauJ.-F.ZemzemiN. (2010). Mathematical modeling of electrocardiograms: a numerical study. Ann. Biomed. Eng.38, 1071–1097. 10.1007/s10439-009-9873-0
12
Chamorro-ServentJ.DuboisR.PotseM.CoudièreY. (2017). Improving the spatial solution of electrocardiographic imaging: a new regularization parameter choice technique for the Tikhonov method, in Functional Imaging and Modeling of the Heart, eds PopM.WrightG. (Toronto, ON: Springer).
13
CochetH.DuboisR.SacherF.DervalN.SermesantM.HociniM.et al. (2014). Cardiac arrhythmias: multimodal assessment integrating body surface ECG mapping into cardiac imaging. Radiology271, 239–247. 10.1148/radiol.13131331
14
Colli-FranzoneP.PennacchioM.GuerriL. (2000). Accurate computation of electrograms in the left ventricular wall. Math. Mod. Methods Appl. Sci.10, 507–538. 10.1142/S0218202500000288
15
CoudièreY.BourgaultY.RiouxM. (2014). Optimal monodomain approximations of the bidomain equations used in cardiac electrophysiology. Math. Models Methods Appl. Sci.24, 1115–1140. 10.1142/S0218202513500784
16
DuchateauJ.PotseM.DuboisR. (2017). Spatially coherent activation maps for electrocardiographic imaging. IEEE Trans. Biomed. Eng.64, 1149–1156. 10.1109/TBME.2016.2593003
17
FalgoutR.BakerA.HensonV. E.YangU. M.KolevT.LeeB.et al. (2017). Hypre Web Page. Available online at: https://computation.llnl.gov/casc/linear_solvers
18
GelernterH. L.SwihartJ. C. (1964). A mathematical-physical model of the genesis of the electrocardiogram. Biophys. J.4, 285–301.
19
GeselowitzD. B. (1989). On the theory of the electrocardiogram. Proc. IEEE77, 857–876.
20
GeselowitzD. B. (1992). Description of cardiac sources in anisotropic cardiac muscle; application of the bidomain model. J. Electrocardiol.25(Suppl.) 65–67.
21
Giffard-RoisinS.JacksonT.FovargueL.LeeJ.DelingetteH.RazaviR.et al. (2017). Non-invasive personalisation of a cardiac electrophysiology model from body surface potential mapping. IEEE Trans. Biomed. Eng.64, 2206–2218. 10.1109/TBME.2016.2629849
22
GimaK.RudyY. (2002). Ionic current basis of electrocardiographic waveforms; A model study. Circ. Res.90, 889–896.
23
GreensiteF.HuiskampG. (1998). An improved method for estimating epicardial potentials from the body surface. IEEE Trans. Biomed. Eng.45, 98–104.
24
HenriquezC. S. (2014). A brief history of tissue models for cardiac electrophysiology. IEEE Trans. Biomed. Eng.61, 1457–1465. 10.1109/TBME.2014.2310515
25
HensonV. E.Meier YangU. (2002). BoomerAMG : a parallel algebraic multigrid solver and preconditioner. Appl. Numer. Math.41, 155–177. 10.1016/S0168-9274(01)00115-5
26
HilleB. (2001). Ion Channels of Excitable Membranes. Sunderland, MA: Sinauer Associates, Inc.
27
HoogendijkM. G.PotseM.LinnenbankA. C.VerkerkA. O.den RuijterH. M.van AmersfoorthS. C. M.et al. (2010). Mechanism of right precordial ST-segment elevation in structural heart disease: excitation failure by current-to-load mismatch. Heart Rhythm7, 238–248. 10.1016/j.hrthm.2009.10.007
28
HoracekB. M. (1973). Digital model for studies in magnetocardiography. IEEE Trans. Magn.3, 440–444.
29
JacquemetV. (2012). An eikonal-diffusion solver and its application to the interpolation and the simulation of reentrant cardiac activations. Comput. Methods Programs Biomed.108, 548–558. 10.1016/j.cmpb.2011.05.003
30
JacquemetV. (2015). Modeling left and right atrial contributions to the ECG: a dipole-current source approach. Comput. Biol. Med.65, 192–199. 10.1016/j.compbiomed.2015.06.007
31
JacquemetV. (2017). Equivalent dipole sources to estimate the influence of extracellular myocardial anisotropy in thin-walled cardiac forward models. Math. Biosci.286, 31–38. 10.1016/j.mbs.2017.01.008
32
JanssenA. M.PotyagayloD.DösselO.OostendorpT. F. (2017). Assessment of the equivalent dipole layer source model in the reconstruction of cardiac activation times on the basis of BSPMs produced by an anisotropic model of the heart. Med. Biol. Eng. Comput.10.1007/s11517-017-1715-x [Epub ahead of print].
33
KaniaM.CoudièreY.CochetH.JaïsP.PotseM. (2017). Prediction of the exit site of ventricular tachycardia based on different ECG lead systems, in Computing in Cardiology, eds PickettC.CorsiC.LagunaP.MacLeodR. (Rennes: Computing in Cardiology).
34
KellerD. U. J.WeberF. M.SeemannG.DösselO. (2010). Ranking the influence of tissue conductivities on forward-calculated ECGs. IEEE Trans. Biomed. Eng.57, 1568–1576. 10.1109/TBME.2010.2046485
35
KellerD. U. J.WeissD. L.DösselO.SeemannG. (2012). Influence of IKs heterogeneities on the genesis of the T-wave: a computational evaluation. IEEE Trans. Biomed. Eng.59, 311–322. 10.1109/TBME.2011.2168397
36
KonukogluE.RelanJ.CilingirU.MenzeB. H.ChinchapatnamP.JadidiA.et al. (2011). Efficient probabilistic model personalization integrating uncertainty on data and parameters: application to eikonal-diffusion models in cardiac electrophysiology. Prog. Biophys. Mol. Biol.107, 134–146. 10.1016/j.pbiomolbio.2011.07.002
37
KrassowskaW.NeuJ. C. (1994). Effective boundary conditions for syncytial tissues. IEEE Trans. Biomed. Eng.41, 143–150.
38
KrauseD.PotseM.DickopfT.KrauseR.AuricchioA.PrinzenF. W. (2012). Hybrid parallelization of a large-scale heart model, in Facing the Multicore-Challenge II, Vol. 7174 Lecture Notes in Computer Science, eds KellerR.KramerD.WeissJ.-P. (Berlin: Springer), 120–132.
39
KudryashovaN.TsvelayaV.AngaldzeK.PanfilovA. (2017). Virtual cardiac monolayers for electrical wave propagation. Sci. Rep.7:7887. 10.1038/s41598-017-07653-3
40
LeonL. J.HorácekB. M. (1991). Computer model of excitation and recovery in the anisotropic myocardium. I. Rectangular and cubic arrays of excitable elements. J. Electrocardiol.24, 1–15.
41
LinesG. T.GrøttumP.TveitoA. (2003). Modeling the electrical activity of the heart; A bidomain model of the ventricles embedded in a torso. Comput. Vis. Sci.5, 195–213. 10.1007/s00791-003-0100-5
42
LiuZ.LiuC.HeB. (2006). Noninvasive reconstruction of three-dimensional ventricular activation sequence from the inverse solution of distributed equivalent current density. IEEE Trans. Biomed. Eng.25, 1307–1318. 10.1109/TMI.2006.882140
43
LorangeM.GulrajaniR. M. (1993). A computer heart model incorporating anisotropic propagation: I. Model construction and simulation of normal activation. J. Electrocardiol.26, 245–261.
44
MacLachlanM. C.SundnesJ.LinesG. T. (2005). Simulation of ST segment changes during subendocardial ischemia using a realistic 3-D cardiac geometry. IEEE Trans. Biomed. Eng.52, 799–807. 10.1109/TBME.2005.844270
45
MaillouxG. E.GulrajaniR. M. (1982). Theoretical evaluation of the McFee and Frank vectorcardiographic lead systems using a numerical inhomogeneous torso model. IEEE Trans. Biomed. Eng.29, 322–332.
46
McFeeR.JohnstonF. D. (1953). Electrocardiographic leads; I. introduction. Circulation8, 554–568.
47
McFeeR.JohnstonF. D. (1954). Electrocardiographic leads; III. synthesis. Circulation9, 868–880.
48
MeijborgV. M. F.PotseM.ConrathC. E.BeltermanC. N. W.de BakkerJ. M. T.CoronelR. (2016). Reduced sodium current in the lateral ventricular wall induces inferolateral J-waves. Front. Physiol.7:365. 10.3389/fphys.2016.00365
49
MenaA.FerreroJ. M.MatasJ. F. R. (2015). GPU accelerated solver for nonlinear reaction-diffusion systems. Application to the electrophysiology problem. Comput. Phys. Commun.196, 280–289. 10.1016/j.cpc.2015.06.018
50
MillerW. T.III.GeselowitzD. B. (1978). Simulation studies of the electrocardiogram; I. The normal heart. Circ. Res.43, 301–315.
51
MunteanuM.PavarinoL. F.ScacchiS. (2009). A scalable Newton–Krylov–Schwarz method for the bidomain reaction-diffusion system. SIAM J. Sci. Comput.31, 3861–3883. 10.1137/08074355X
52
NeicA.CamposF. O.PrasslA. J.NiedererS. A.BishopM. J.VigmondE. J.et al. (2017). Efficient computation of electrograms and ECGs in human whole heart simulations using a reaction-eikonal model. J. Comput. Phys.346, 191–211. 10.1016/j.jcp.2017.06.020
53
NeicA.LiebmannM.HoetzlE.MitchellL.VigmondE. J.HaaseG.et al. (2012). Accelerating cardiac bidomain simulations using graphics processing units. IEEE Trans. Biomed. Eng.59, 2281–2290. 10.1109/TBME.2012.2202661
54
NguyênU. C.PotseM.RegoliF.CaputoM. L.ConteG.MurzilliR.et al. (2015). An in-silico analysis of the effect of heart position and orientation on the ECG morphology and vectorcardiogram parameters in patients with heart failure and intraventricular conduction defects. J. Electrocardiol.48, 617–625. 10.1016/j.jelectrocard.2015.05.004
55
NiedererS. A.KerfootE.BensonA. P.BernabeuM. O.BernusO.BradleyC.et al. (2011a). Verification of cardiac tissue electrophysiology simulators using an N-version benchmark. Philos. Trans. A Math. Phys. Eng. Sci.369, 4331–4351. 10.1098/rsta.2011.0139
56
NiedererS. A.MitchellL.SmithN.PlankG. (2011b). Simulating human cardiac electrophysiology on clinical time-scales. Front. Physiol.2:14. 10.3389/fphys.2011.00014
57
NielsenB. F.RuudT. S.LinesG. T.TveitoA. (2007). Optimal monodomain approximations of the bidomain equations. Appl. Math. Comput.184, 276–290. 10.1016/j.amc.2006.05.158
58
NobleD.RudyY. (2001). Models of cardiac ventricular action potentials: iterative interaction between experiment and simulation. Philos. Trans. R. Soc. A359, 1127–1142. 10.1098/rsta.2001.0820
59
OosterhoffP.MeijborgV.van DamP. M.van DesselP. F. H. M.BeltermanC. N. W.StreekstraG. J.et al. (2016). Experimental validation of noninvasive epicardial and endocardial activation imaging. Circ. Arrhythm. Electrophysiol.9:e004104. 10.1161/CIRCEP.116.004104
60
OttinoD.ScacchiS. (2015). BPX preconditioners for the bidomain model of electrocardiology. J. Comput. Appl. Math.285, 151–168. 10.1016/j.cam.2015.02.011
61
PassiniE.BrittonO. J.LuH. R.RohrbacherJ.HermansA. N.GallacherD. J.et al. (2017). Human in silico drug trials demonstrate higher accuracy than animal models in predicting clinical pro-arrhythmic cardiotoxicity. Front. Physiol.8:668. 10.3389/fphys.2017.00668
62
PezzutoS.KaľavskyP.PotseM.PrinzenF. W.AuricchioA.KrauseR. (2017). Evaluation of a rapid fully anisotropic model for ECG simulation. Front. Physiol.8:265. 10.3389/fphys.2017.00265
63
PilkingtonT. C.MorrowM. N.StanleyP. C. (1987). A comparison of finite element and integral equation formulations for the calculation of electrocardiographic potentials–II. IEEE Trans. Biomed. Eng.34, 258–260.
64
PotseM.DubéB.RicherJ.VinetA.GulrajaniR. M. (2006). A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart. IEEE Trans. Biomed. Eng.53, 2425–2435. 10.1109/TBME.2006.880875
65
PotseM.KuijpersN. H. L. (2010). Simulation of fractionated electrograms at low spatial resolution in large-scale heart models, in Computing in Cardiology, Vol. 37, eds MurrayA.DassenW.LagunaP.MoodyG.SimsA. (Belfast: Computing in Cardiology), 849–852.
66
PotseM.VinetA.OpthofT.CoronelR. (2009). Validation of a simple model for the morphology of the T wave in unipolar electrograms. Am. J. Physiol. Heart Circ. Physiol.297, H792–H801. 10.1152/ajpheart.00064.2009
67
PullanA. J.BradleyC. P. (1996). A coupled cubic Hermite finite element/boundary element procedure for electrocardiographic problems. Comput. Mech.18, 356–368.
68
RamanathanC.GhanemR.JiaP.RyuK.RudyY. (2004). Noninvasive electrocardiographic imaging for cardiac electrophysiology and arrhythmia. Nat. Med.10, 422–428. 10.1038/nm1011
69
RichardsD. F.GlosliJ. N.DraegerE. W.MirinA. A.ChanB.FattebertJ.-L.et al. (2013). Towards real-time simulation of cardiac electrophysiology in a human heart at high resolution. Comput. Meth. Biomech. Biomed. Eng.16, 802–805. 10.1080/10255842.2013.795556
70
RushS.LarsenH. (1978). A practical algorithm for solving dynamic membrane equations. IEEE Trans. Biomed. Eng.25, 389–392.
71
SaleheenH. I.NgK. T. (1997). New finite difference formulations for general inhomogeneous anisotropic bioelectric problems. IEEE Trans. Biomed. Eng.44, 800–809.
72
ShouG.XiaL.JiangM.WeiQ.LiuF.CrozierS. (2008). Truncated total least squares: a new regularization method for the solution of ECG inverse problems. IEEE Trans. Biomed. Eng.55, 1327–1355. 10.1109/TBME.2007.912404
73
StenroosM.HaueisenJ. (2008). Boundary element computations in the forward and inverse problems of electrocardiography: comparison of collocation and Galerkin weightings. IEEE Trans. Biomed. Eng.55, 2124–2133. 10.1109/TBME.2008.923913
74
SundnesJ.LinesG. T.MardalK. A.TveitoA. (2002). Multigrid block preconditioning for a coupled system of partial differential equations modeling the electrical activity in the heart. Comput. Meth. Biomech. Biomed. Eng.5, 397–409. 10.1080/1025584021000025023
75
ten TusscherK. H. W. J.PanfilovA. V. (2006). Alternans and spiral breakup in a human ventricular tissue model. Am. J. Physiol. Heart Circ. Physiol.291, H1088–H1100. 10.1152/ajpheart.00109.2006
76
TiltJ.TaggartJ. S.TuftsL.BaconD. L.DurantP.WodlingerH.et al. (2013). Multi-Layered Sensor Apparatus.Washington, DC: U.S. Patent and Trademark Office.
77
TrudelM.-C.DubéB.PotseM.GulrajaniR. M.LeonL. J. (2004). Simulation of propagation in a membrane-based computer heart model with parallel processing. IEEE Trans. Biomed. Eng.51, 1319–1329. 10.1109/TBME.2004.827934
78
TungL. (1978). A Bi-Domain Model for Describing Ischemic Myocardial D-C Potentials. Ph.D. thesis, MITx, Cambridge MA.
79
TyslerM.KneppoP.TurzováM.ŠvehlíkováJ.KarasS.HeblákováE.et al. (2007). Noninvasive assessment of local myocardium repolarization changes using high resolution surface ECG mapping. Physiol. Res.56(Suppl. 1), S133–S141.
80
van der VorstH. A. (1992). Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput.13, 631–644.
81
van OosteromA. (2001). Genesis of the T wave as based on an equivalent surface source model. J. Electrocardiol.34(Suppl.), 217–227. 10.1054/jelc.2001.28896
82
van OosteromA.JacquemetV. (2005). Genesis of the P wave: atrial signals as generated by the equivalent double layer source model. Europace7(Suppl. 2), S21–S29. 10.1016/j.eupc.2005.05.001
83
VandersickelN.de BoerT. P.VosM.PanfilovA. V. (2016). Perpetuation of torsade de pointes in heterogeneous hearts: competing foci or re-entry?J. Physiol.594, 6865–6878. 10.1113/JP271728
84
VázquezM.ArísR.HouzeauxG.AubryR.VillarP.Garcia-BarnésJ.et al. (2011). A massively parallel computational electrophysiology model of the heart. Int. J. Numer. Methods Biomed. Eng.27, 1911–1929. 10.1002/cnm.1443
85
VigmondE. J.AguelF.TrayanovaN. A. (2002). Computational techniques for solving the bidomain equations in three dimensions. IEEE Trans. Biomed. Eng.49, 1260–1269. 10.1109/TBME.2002.804597
86
WangD.KirbyR. M.MacLeodR. S.JohnsonC. R. (2013). Inverse electrocardiographic source localization of ischemia: an optimization framework and finite element solution. J. Comput. Phys.250, 403–424. 10.1016/j.jcp.2013.05.027
87
WangL.DawoudF.YeungS.-K.ShiP.WongK. C. L.LiuH.et al. (2013). Transmural imaging of ventricular action potentials and post-infarction scars in swine hearts. IEEE Trans. Med. Imaging32, 731–747. 10.1109/TMI.2012.2236567
88
Weber dos SantosR.PlankG.BauerS.VigmondE. J. (2004). Parallel multigrid preconditioner for the cardiac bidomain model. IEEE Trans. Biomed. Eng.51, 1960–1968. 10.1109/TBME.2004.834275
89
WesternD.HansonB.TaggartP. (2015). Measurement bias in activation-recovery intervals from unipolar electrograms. Am. J. Physiol. Heart Circ. Physiol.308, H331–H338. 10.1152/ajpheart.00478.2014
90
ZemzemiN.BernabeuM. O.SaizJ.CooperJ.PathmanathanP.MiramsG. R.et al. (2013). Computational assessment of drug-induced effects on the electrocardiogram: from ion channel to body surface potentials. Br. J. Pharmacol.168, 718–733. 10.1111/j.1476-5381.2012.02200.x
91
ZemzemiN.DobrzynskiC.BearL.PotseM.DalletC.CoudièreY.et al. (2015). Effect of the torso conductivity heterogeneities on the ECGI inverse problem solution, in Computing in Cardiology, Vol. 42, (Nice), 233–236.
Summary
Keywords
numerical modeling, electrocardiogram, high-performance computing, reaction-diffusion model, bidomain model, lead fields
Citation
Potse M (2018) Scalable and Accurate ECG Simulation for Reaction-Diffusion Models of the Human Heart. Front. Physiol. 9:370. doi: 10.3389/fphys.2018.00370
Received
14 January 2018
Accepted
27 March 2018
Published
20 April 2018
Volume
9 - 2018
Edited by
Mariano Vázquez, Barcelona Supercomputing Center, Spain
Reviewed by
Arun V. Holden, University of Leeds, United Kingdom; Mohammad Hasan Imam, American International University-Bangladesh, Bangladesh
Updates
Copyright
© 2018 Potse.
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) and the copyright owner 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.
*Correspondence: Mark Potse mark@potse.nl
This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.