Electrocardiographic Imaging Using a Spatio-Temporal Basis of Body Surface Potentials—Application to Atrial Ectopic Activity

Electrocardiographic imaging (ECGI) strongly relies on a priori assumptions and additional information to overcome ill-posedness. The major challenge of obtaining good reconstructions consists in finding ways to add information that effectively restricts the solution space without violating properties of the sought solution. In this work, we attempt to address this problem by constructing a spatio-temporal basis of body surface potentials (BSP) from simulations of many focal excitations. Measured BSPs are projected onto this basis and reconstructions are expressed as linear combinations of corresponding transmembrane voltage (TMV) basis vectors. The novel method was applied to simulations of 100 atrial ectopic foci with three different conduction velocities. Three signal-to-noise ratios (SNR) and bases of six different temporal lengths were considered. Reconstruction quality was evaluated using the spatial correlation coefficient of TMVs as well as estimated local activation times (LAT). The focus localization error was assessed by computing the geodesic distance between true and reconstructed foci. Compared with an optimally parameterized Tikhonov-Greensite method, the BSP basis reconstruction increased the mean TMV correlation by up to 22, 24, and 32% for an SNR of 40, 20, and 0 dB, respectively. Mean LAT correlation could be improved by up to 5, 7, and 19% for the three SNRs. For 0 dB, the average localization error could be halved from 15.8 to 7.9 mm. For the largest basis length, the localization error was always below 34 mm. In conclusion, the new method improved reconstructions of atrial ectopic activity especially for low SNRs. Localization of ectopic foci turned out to be more robust and more accurate. Preliminary experiments indicate that the basis generalizes to some extent from the training data and may even be applied for reconstruction of non-ectopic activity.

1. Introduction

Reconstructing the heart's electrical activity from non-invasively measured body surface potentials (BSP) is known as the inverse problem of electrocardiography (Pullan et al., 2010). The ill-posedness of this problem can be overcome by introducing additional information—a technique called regularization. Classical regularization methods such as Tikhonov regularization add the “information” that the solution must be of small signal energy or smooth in space or time. However, they do not complement the measurements with physiological information about the spread of cardiac excitation, which could greatly improve the uniqueness of the solution. For example, ambiguities can arise between sources on two different regions oriented in parallel, when body surface potentials are projected back onto the heart. Incorporating information about the spatio-temporal coherence of excitation spread, i.e., taking into account that excitation waves can only gradually propagate across connected regions, is expected to help resolve ambiguities. Several approaches have been proposed to incorporate electrophysiological knowledge (van Oosterom, 1999; Messnarz et al., 2004; Ghodrati et al., 2006; Wang et al., 2010; Potyagaylo et al., 2014, 2016a; Cluitmans et al., 2018). Cluitmans et al. (2017) reconstructed potentials on the ventricular epicardium as sparse combinations of spatial source basis vectors generated from simulations of many paced beats. In this work, we use a related approach. However, we suggest to create a basis of body surface potential patterns instead of source patterns to condense the information to what can possibly be measured on the body surface. Corresponding basis patterns in source space are then obtained and combined to express reconstructions. Furthermore, we use a spatio-temporal basis instead of a spatial-only one, as we believe this reduces ambiguities and increases the robustness to noise. We demonstrate the method in an application to atrial ectopic activity.

2. Methods

The outline of this study is illustrated in Figure 1. First, fast marching simulations of 200 ectopic foci are performed and forward calculated to create basis vectors. Reconstructions are then performed using these basis vectors for another set of 100 ectopic foci simulations. Here, the monodomain model is used with three different conduction velocities (CV) and forward calculated BSPs are corrupted with three different levels of noise. Finally, reference reconstructions are obtained using the Tikhonov-Greensite method and the same metrics are calculated for both reconstruction methods.

FIGURE 1

Figure 1. Flowchart illustrating the outline of this study.

2.1. Geometries

Figure 2 shows the geometries used for forward and inverse calculations. They are the same as in Schuler et al. (2017), which in turn are based on Figuera et al. (2016). The surface meshes of the atria and torso consist of 4,800 and 844 nodes and have an average edge length of 3.4 and 27.0 mm, respectively. A subset of 173 torso nodes have been selected as electrodes (blue spheres). Much finer tetrahedral meshes of the atria were used for excitation simulations (142 k nodes, average edge length: 0.9 mm). Blue and red spheres on the atria mark evenly distributed pacing locations used for fast marching and monodomain simulations, respectively. Note that pacing locations for both sets of simulations do generally not coincide.

FIGURE 2

Figure 2. Meshes of the atria and torso used for forward and inverse calculations. Blue spheres on the atria indicate 200 pacing locations used to create basis vectors, while red spheres indicate 100 pacing locations used for monodomain simulations to be reconstructed. Torso electrodes are shown in blue. Red spheres on the torso mark reference nodes.

2.2. Fast Marching Simulations

For creating basis vectors, fast marching simulations of 200 paced beats were performed. First, local activation times (LAT) t_a were computed by solving the eikonal equation with the fast marching method (Pernod et al., 2011):

$\begin{array}{l}\left|\right|\nabla t_a|{|}_2=\frac{1}{c}\end{array}$

The CV was homogeneously set to c = 0.8m/s. A TMV template was then aligned with LATs. Experiments with a template based on the Courtemanche et al. cell model and a step-function-like template (Figure 3) showed that basis vectors created as described in section 2.5.1 are mainly determined by the depolarization upstroke and it is not necessary to include the repolarization. Therefore, we decided to use the step-function-like template. The temporal sampling period was chosen to be 2ms.

FIGURE 3

Figure 3. TMV templates based on the Courtemanche et al. cell model (red) and the step function (blue).

2.3. Monodomain Simulations

As activity to be reconstructed, 100 paced beats were simulated using the monodomain model:

$\begin{array}{l}\nabla ·({\sigma }_{mono}\nabla V_m)=\beta \left(C_m\frac{\partial V_m}{\partial t}+I_{ion}\right)\end{array}$

The monodomain conductivity σ_mono was chosen homogeneous and isotropic and, together with the surface-to-volume ratio β, was adjusted to obtain the desired CV. In order to study what happens, if a wrong CV is assumed for the creation of basis vectors, we varied the CV by ±50% of its baseline value. Therefore, three different CVs were used for monodomain simulations: 0.4, 0.8, and 1.2m/s. Ionic currents I_ion across the cell membrane were defined according to Courtemanche et al. (1998). As for fast marching simulations, the temporal resolution is 2ms.

2.4. Forward Calculation

According to bidomain theory, extracellular potentials ϕ are related to volumetric TMVs V_m by:

$\begin{array}{ll}\nabla ·(({\sigma }_i+{\sigma }_e)\nabla \varphi )=-\nabla ·({\sigma }_i\nabla V_m)& (1)\end{array}$

In this work, we assume isotropic, i.e., scalar intra- and extra-cellular conductivities σ_i and σ_e, respectively. As this is a special case of equal intra- and extra-cellular anisotropy ratios, volumetric TMVs may be replaced by TMVs on the myocardial surface (Yamashita and Geselowitz, 1985). For a homogeneous torso with the same bulk conductivity as the heart σ, the potential ϕ at an observation point can be calculated from surface TMVs V_m by (Simms and Geselowitz, 1995):

$\begin{array}{ll}\varphi =-\frac{1}{4\pi \sigma }{\displaystyle {\int }_{S_H}}{\sigma }_i{V}_{m}d{\Omega }_H-\frac{1}{4\pi }{\displaystyle {\int }_{S_B}}{\varphi }_{B}d{\Omega }_B& (2)\end{array}$

S_H and S_B are the surfaces bounding the myocardium and torso, respectively. dΩ_H and dΩ_B are the solid angles subtended by an area element on the corresponding surface as seen from an observation point. By applying the boundary element method to solve (2) for BSPs ϕ = ϕ_B, a lead field matrix A is obtained that transforms surface TMVs to BSPs. Intracellular and bulk conductivities are set to σ_i = 0.05S/m, and σ = 0.2S/m and potentials are referenced to Wilson's central terminal. We assume the following linear forward model:

$\begin{array}{l}{\text{b}}_k=\text{A}{\text{x}}_k+{\epsilon }_k\end{array}$

where b_k, x_k and ε_k are BSPs, TMVs, and white Gaussian noise for all nodes at a time step k, respectively. For fast marching simulations used to create basis vectors, we do not add any noise. For monodomain simulations, we consider three different signal-to-noise ratios (SNRs): 40, 20, and 0dB. It is assumed that each electrode is affected by the same absolute noise power, which is set to the average signal power of all electrodes divided by the SNR. For solving the inverse problem, we assume perfect knowledge of A and thus neglect errors due to imperfect geometries and conductivities.

2.5. Reconstruction Using Spatio-Temporal BSP Basis

The following identifiers and terminology will be used in this section:

2.5.1. Creation of Basis Vectors

As in Cluitmans et al. (2017), we use the singular value decomposition (SVD) to create basis vectors. In order to get a spatio-temporal basis, we define an observation as the column-wise concatenation of values at all nodes for all time steps within a time window of length L. As we do not know the time delay between the activity to be reconstructed and the activities used to create the basis, we include all possible delays by continuously time-shifting the window by a single time step. This way, a total of K−L+1 observations are generated for each simulation. Row-wise concatenation of all observations (all time shifts of all simulations) then yields a data matrix D. Using x_k to denote TMVs for all nodes at a time step k, the TMV data matrix D_x is thus given as:

${\mathbf{\text{D}}}_{\mathbf{\text{x}}}=\left[\begin{array}{c}\left(\begin{array}{cccc}{\mathbf{\text{x}}}_1^{\top }& {\mathbf{\text{x}}}_2^{\top }& \cdots & {\mathbf{\text{x}}}_L^{\top }\\ {\mathbf{\text{x}}}_2^{\top }& {\mathbf{\text{x}}}_3^{\top }& \cdots & {\mathbf{\text{x}}}_{L+1}^{\top }\\ ⋮& ⋮& & ⋮\\ {\mathbf{\text{x}}}_{K-L+1}^{\top }& {\mathbf{\text{x}}}_{K-L+2}^{\top }& \cdots & {\mathbf{\text{x}}}_K^{\top }\end{array}\right)\\ ⋮\\ \left(\text{repeat for all simulations}\right)\end{array}\right]$

By replacing x with b, a BSP data matrix D_b can be constructed in exactly the same way. From Ax_k = b_k, it follows that the whole TMV data matrix can be forward calculated using a block diagonal lead field matrix $\tilde{A}$:

$\begin{array}{ll}\tilde{A}{\text{D}}_{\text{x}}^{\top }={\text{D}}_{\text{b}}^{\top }\iff {\text{D}}_{\text{b}}={\text{D}}_{\text{x}}{\tilde{A}}^{\top }\mathrm{\text{with }}\tilde{A}={\text{I}}_L\otimes \text{A}& (3)\end{array}$

I_L is the L×L identity matrix and ⊗ denotes the Kronecker product. We now perform an SVD of the BSP data matrix:

$\begin{array}{ll}{\text{D}}_{\text{b}}=\text{U}\text{S}{\text{V}}_{\text{b}}^{\top }& (4)\end{array}$

The columns of V_b are spatio-temporal basis vectors of BSPs. We now want to find the corresponding TMV basis vectors V_x, for which holds:

$\begin{array}{ll}{\text{V}}_{\text{b}}=\tilde{A}{\text{V}}_{\text{x}}& (5)\end{array}$

Substituting (3) and (5) in (4) yields:

$\begin{array}{l}{\text{D}}_{\text{x}}{\tilde{A}}^{\top }=\text{U}\text{S}{\text{V}}_{\text{x}}^{\top }{\tilde{A}}^{\top }\iff {\text{V}}_{\text{x}}^{\top }={(\text{U}\text{S})}^{+}{\text{D}}_{\text{x}}={\text{S}}^{+}{\text{U}}^{\top }{\text{D}}_{\text{x}}\end{array}$

(·)⁺ denotes the Moore–Penrose pseudoinverse. This shows that the TMV basis can directly be calculated from the TMV data matrix using an inversion of US, the scores matrix obtained from the SVD of the BSP data matrix. An inversion of $\tilde{A}$ is not required.

2.5.2. Reconstruction in Terms of Basis Vectors

Now the BSPs of a patient show up. They will be called B = [b₁, b₂, …, b_K] in the following. Here, K is the total number of time steps to be reconstructed. To reconstruct TMVs in terms of a reduced number P of basis vectors, we first perform a least-squares regression using BSP basis vectors. This results in the optimal basis vector weights W:

$\begin{array}{ll}\text{W}=arg{\displaystyle \underset{\text{W}}{min}}\left||{\text{V}}_{\text{b}}\left(:,1:P\right)\text{W}-\tilde{B}\right|{|}_F^2& (6)\end{array}$

||·||_F denotes the Frobenius norm and $\tilde{B}$ are the “measured” BSPs reshaped into the format of basis vectors. Using MATLAB notation, they are given by:

$\begin{array}{l}\tilde{B}\left(:,k\right)=\text{reshape}\left(\text{B}\left(:,k:k+L-1\right),LM,1\right)\mathrm{\text{ for }}k=1,2,\dots ,\\ K-L+1\end{array}$

Since the columns of V_b form an orthonormal basis, the solution to (6) is given by:

$\begin{array}{ll}\text{W}={\text{V}}_{\text{b}}{\left(:,1:P\right)}^{\top }\tilde{B}& (7)\end{array}$

This can be seen as filtering the BSPs by projecting them onto the P most important BSP basis vectors. The weights W are now used to obtain the reconstructed TMVs ${\tilde{X}}_r$ as linear combination of corresponding TMV basis vectors:

$\begin{array}{l}{\tilde{X}}_r={\text{V}}_{\text{x}}\left(:,1:P\right)\text{W}\end{array}$

Each column of ${\tilde{X}}_r$ contains the row-wise concatenation of all time windows with length L. As final solution X_r, we therefore extract the central time step of each window:

$\begin{array}{l}{\text{X}}_r={\text{V}}_{\text{x}}\left(\frac{L-1}{2}N+(1:N),1\mathrm{\text{:}}P\right)\text{W}\end{array}$

Corresponding BSPs can be obtained in the same manner:

$\begin{array}{l}{\text{B}}_r={\text{V}}_{\text{b}}\left(\frac{L-1}{2}M+(1:M),1\mathrm{\text{:}}P\right)\text{W}=\text{A}{\text{X}}_r\end{array}$

We would like to point out that creating the basis vectors in BSP space instead of in source space, i.e., calculating the SVD of D_b instead of D_x, is the key step which makes it possible to perform an unregularized least-squares regression without additional constraints. If the basis vectors were created in source space, they may still contain redundant information with respect to BSPs and therefore a regularized regression would be necessary: $\text{W}=arg{\displaystyle \underset{\text{W}}{min}}\left\{{\Vert \tilde{A}{\text{V}}_{\text{x}}\left(:,1:P\right)\text{W}-\tilde{B}\Vert }_F^2+\lambda {\Vert \text{R}{\text{V}}_{\text{x}}\left(:,1:P\right)\text{W}\Vert }_F^2\right\}$, where R is a regularization matrix. This would require inverting a P × P matrix for every combination (P, λ), while only a fast matrix multiplication (7) has to be computed for each P when using the BSP basis.

2.5.3. Choice of Basis Dimension

Following the concept of the L-curve, the “optimal” basis dimension is determined from the log-log graph of the residual norm ||B_r−B||_F versus the corresponding basis dimension P (left diagram in Figure 4). Instead of the maximal curvature, which is not very pronounced in the resulting “L-curves,” we found that the minimal absolute slope (blue circles) of a smoothing spline fit is a good criterion for selecting the basis dimension. To ensure that solutions are not underregularized, we set the basis dimension 10% lower than at the point of minimal absolute slope (red circles). The right diagram in Figure 4 shows that the resulting basis dimension for one specific SNR depends linearly on the basis length L.

FIGURE 4

Figure 4. (Left) Exemplary “L-curves” for a basis length of L = 17. Black dots: original data points. Red curves: smoothing spline fits. Blue circles: Points of min. absolute slope. Red circles: Points 10% left of min. slope points. (Right) Dependency of basis dimension on basis length.

2.6. Reference Reconstruction With Tikhonov-Greensite

For comparison, reconstructions with the Tikhonov-Greensite method (TikhGS) are performed (Greensite and Huiskamp, 1998). While standard Tikhonov methods regularize each time step individually, TikhGS performs Tikhonov regularization for the p most important temporal singular vectors of BSPs, which helps to eliminate noise and generally leads to better reconstructions, given the number of singular vectors used is chosen appropriately. We apply the epi-endo projection described in Schuler et al. (2017), which further improves the solution. TikhGS is then used with simultaneous zero- and second-order constraints:

$\begin{array}{l}B=US{V}^{\top },\overline{V}=V(:,1:p)\hfill \\ X\overline{V}=\arg \underset{(X\overline{V})}{\min }\left\{\Vert A(X\overline{V})-B\overline{V}{\Vert }_F^2+\text{λ}\Vert L(X\overline{V}){\Vert }_F^2+\eta \Vert X\overline{V}{\Vert }_F^2\right\}\hfill \\ X=(X\overline{V}){\overline{V}}^{\top }\hfill \end{array}$

In order to provide the reference method with the best possible parameters, p was varied between 3 and 13 and regularization parameters λ and η were optimized for each case individually to maximize the mean of spatial TMV correlation with the ground truth over the period of depolarization. The downhill simplex method was used for optimization.

2.7. Post-Processing and Metrics

2.7.1. Transmembrane Voltages

To assess the quality of reconstructed TMVs, the Pearson correlation coefficient between reconstructions and the ground truth is computed separately for each time step across all nodes (spatial CC). As we want to quantify, how good the depolarization is being reconstructed, we only calculate the spatial CC for time steps between the first and the last activation of each simulation, as defined by the ground truth.

2.7.2. Local Activation Times

Local activation times represent one of the most important characteristic of cardiac excitation spread and are therefore estimated from reconstructed TMVs. We use the “global activation time” approach described in Dubois et al. (2012), which is based on cross-correlating signals of nearby nodes to find their time delay. The method has further been advanced in Duchateau et al. (2017) to combine delay-based and deflection-based activation times. In this work, however, we stick with the delay-only formulation. The procedure is illustrated in Figure 5. First, TMVs are oversampled to allow for a precise alignment in time. As signal for cross-correlation, we then use a Gaussian filtered version of the magnitude of the surface gradient of TMVs: ||∇V_m(t)||₂. Using a lowpass filtered time derivative as described in Duchateau et al. (2017) yielded unsatisfactory results for both TikhGS and the BSP basis reconstruction. This might be explained by the fact that the spatial gradient of TMVs is the source of body surface potentials according to (1). LATs are finally estimated from the delays using least-squares regression of a linear model. As LAT metric, we calculate the spatial CC between LATs estimated for reconstructions and the ground truth.

FIGURE 5

Figure 5. Delay-based LAT estimation. D is a difference matrix to obtain pair-wise delays from LATs.

2.7.3. Ectopic Focus Localization

BSPs are directly proportional to the solid angle at a measurement point subtended by the depolarization wavefront, which separates regions of low and high TMVs (see Equation 2). For a given noise level, the instantaneous SNR therefore rises with the size of the depolarization wavefront and the excitation origin cannot reliably be determined from the small signal at the very onset of excitation. If the depolarization wavefront spreads too far from the origin, however, the uncertainty of localizing the origin within the depolarized region increases as well. As the CV is not known beforehand, we therefore do not base the focus localization on one specific time step, but use the temporal mean of many time steps after excitation onset. If the TMV waveform was a Heaviside step function, the TMV time integral (and thus also its temporal mean) would be proportional to reversed activation times. For other TMV waveforms, the temporal mean still yields a valid activation “sequence,” as long as the TMV time integral is increasing (Schulze, 2015). Based on the TMV waveform of the Courtemanche et. al. model (Figure 3), we chose to calculate the mean over 200 ms after excitation onset. In general, this time should not be chosen larger than the effective refractory period. In this work, the time of excitation onset is assumed to be known. In practice, it would have to be determined as the P wave onset. Having obtained activation sequences, a template matching approach is employed to detect the focus location: For every mesh node, the zero-mean normalized correlation (ZNC) of the temporal mean of TMVs $\overline{V_m}$ and the reversed geodesic distance field, originating from the respective node and truncated at 3 cm, is calculated. As the ZNC only measures similarity in shape, not magnitude, the result is further weighted with $(\overline{V_m}-min\left\{\overline{V_m}\right\})$. The maximum of this “focus measure” is finally detected as focus. The method is illustrated in Figure 6 for both BSP basis and TikhGS reconstructions. Correlating with the geodesic distance fields can also be seen as a transformation to find the center of mass of the TMV distribution on a curved surface.

FIGURE 6

Figure 6. Ectopic focus localization illustrated for exemplary reconstructions with the BSP basis (L = 33) and TikhGS. CV = 0.8m/s. SNR = 0dB. Black spheres: true focus, white spheres: reconstructed focus. Geodesic distance fields shown on the left are reversed and truncated at 3 cm (dark blue). Gray areas are not taken into account for correlation.

As localization error, the geodesic distance between the true and reconstructed focus is evaluated. In contrast to the Euclidean distance, this metric correctly yields large errors for nearby points that are not directly connected via the geometry, such as two points on the opposite side of the interatrial region.

3. Results

3.1. Basis Vectors and Singular Values

Figure 7 (left) depicts 3 time steps of exemplary TMV basis vectors for a basis length of L = 17. It can be seen that the spatial and/or temporal frequency increases with the basis vector number and that spatial patterns evolve over time. The diagram on the right shows that the larger the basis length, the more basis vectors are needed to represent the same proportion of information contained within all basis vectors.

FIGURE 7

Figure 7. (Left) Examples of TMV basis vectors for L = 17. p: basis vector number. k: time step. Colors only visualize spatial morphology. (Right) Cumulative sum of singular values for different basis lengths L.

3.2. Robustness to Noise

Metrics for reconstructions of monodomain simulations with the same CV as used for basis creation are shown in Figure 8. For all noise levels considered, BSP basis reconstructions perform consistently better than TikhGS. Even a spatial-only basis (L = 1) leads to an improvement. Increasing the basis length further increases the correlation coefficients for both TMVs and LATs, especially for low SNRs. This can also be seen from Figure 9, where LATs for an ectopic focus near the left inferior pulmonary vein are shown for different SNRs and basis lengths. While the LAT map for L = 1 and 40dB already looks much like the ground truth, basis lengths of at least L = 9 and L = 17 are needed for 20 and 0dB, respectively.

FIGURE 8

Figure 8. Metrics for a CV of 0.8m/s. Colors represent different SNRs. Blue: 40dB, red: 20dB, yellow: 0dB. Boxes: 25–75th percentile. Whiskers: 1.5 inter-quartile range. Filled circles and lines represent the mean. (Left) Spatial CC of TMVs. (Right) Spatial CC of LATs.

FIGURE 9

Figure 9. LATs for different SNRs and basis lengths L. CV = 0.8m/s.

3.3. Influence of Conduction Velocity

In theory, errors in the spatial dimension, i.e., the position of the wavefront, increase linearly with time for a mismatch of CV. Therefore it is expected that the reconstruction quality deteriorates for large basis lengths if the CV used to create basis vectors deviates from the actual CV. Figures 10A,B show the metrics for a CV of 0.4 and 1.2m/s, respectively, while the CV assumed for basis vectors remains at 0.8m/s. It can be seen that there is now indeed an upper limit for the improvement with increasing basis lengths. Results show a clear tradeoff between the error due to a wrong CV for large basis lengths and the error due to lower robustness to noise for small basis lengths. Although the BSP basis reconstruction still outperforms TikhGS in every case, it can be seen from comparison of both figures that overestimating the true CV during basis creation seems to be less problematic than underestimating it. Figure 11 shows another interesting effect of wrong CVs. If the CV is overestimated (top row), LATs estimated from reconstructions with large basis lengths are smaller than true LATs, suggesting a larger than actual CV. To a lesser extent, the opposite effect can be seen for an underestimation (bottom row).

FIGURE 10

Figure 10. (A) Metrics for a CV of 0.4m/s. (B) Metrics for a CV of 1.2m/s. See caption of Figure 8 for details.

FIGURE 11

Figure 11. LATs for different CVs and basis lengths L. SNR = 40dB.

3.4. Ectopic Focus Localization

Localization errors are shown in Figure 12. The mean error decreases with increasing basis length, even for CVs different from 0.8m/s. For L = 41, it is approximately halved compared to TikhGS (see Table 1). Another important improvement is the reduction of the maximum localization error. While errors for TikhGS range up to 90mm, the maximum error for L = 41 is 34mm. For a CV of 0.4m/s, however, there is an outlier for L = 25 and L = 33. This case is illustrated on the right of Figure 12C. For this ectopic focus at the orifice of the inferior vena cava, the reconstruction first shows a false activity on the nearby left atrium. Only after the excitation has further increased in size, the reconstruction continues to show the activity at the correct location. A better localization thus would have been obtained for a later (or longer) time window used for calculating the temporal mean.

FIGURE 12

Figure 12. Localization errors for all 100 ectopic foci. See caption of Figure 8 for colors and markers. Outliers are shown as empty circles. (A) CV = 0.8m/s. (B) 1.2m/s. (C) 0.4m/s. The outlier for CV = 0.4m/s and L = 33 is depicted on the bottom right. Black sphere: true focus, white sphere: reconstructed focus.

TABLE 1

Table 1. Localization errors in mm.

In general, the BSP basis reconstruction largely resolves ambiguities. In Figure 13A, this is demonstrated for a focus on the anterior-septal wall of the left atrium. TikhGS and a purely spatial basis fail to recover the activity at the correct spot and show multiple activities on the right and left atrium instead. Using a spatio-temporal basis, however, leads to a reconstruction at the correct spot. Increasing the basis length progressively increases the uniqueness of the solution. An even better impression on where ambiguities arise between sources can be obtained by taking a look at the distributions of localization errors across the atria, as shown in Figure 13B. For TikhGS, the largest errors occur at the interatrial region, where left and right atrial surfaces are very close to each other and oriented in parallel. Using a BSP basis of sufficient length greatly reduces these errors.

FIGURE 13

Figure 13. (A) Focus measure for a focus on the anterior-septal wall of the left atrium (CV = 0.8m/s, SNR = 0dB). Black sphere: true focus, white sphere: reconstructed focus. (B) Distribution of localization errors across the atria for TikhGS and the BSP basis reconstruction (CV = 0.8m/s, SNR = 20dB, results for 0dB are qualitatively similar). Values were interpolated by minimizing the Laplacian at all nodes (Oostendorp et al., 1989).

3.5. Effect of Non-conducting Region

In order to test, how much the reconstruction relies on the excitation patterns in the training data, we added a non-conducting scar region to one simulation. For that purpose, the monodomain conductivity was set to zero in a circular region with a diameter of 4cm on the right atrium. The results are shown in Figure 14. It can be seen that the reconstructed wavefront propagates around the non-conducting region on the right atrium, even though the training data did not include such a pattern. The basis created by including all time shifts of excitation patterns in the data matrix therefore generalizes from individual patterns.

FIGURE 14

Figure 14. TMVs for a simulation containing a non-conducting scar region on the right atrium. SNR = 20dB. CV = 0.8m/s. Each time step was normalized by subtracting the spatial mean and dividing by the spatial standard deviation (SD).

4. Discussion

We demonstrated that using a spatio-temporal basis of BSPs to reconstruct TMVs improves the robustness to noise, resolves ambiguities between sources and leads to better localization of atrial ectopic foci than TikhGS. Of all possible solutions fitting to measured BSPs, the proposed method selects the one which is most probable with regard to the training data used to create basis vectors. This approach allows us to reconstruct hidden sources as well, given they occur in conjunction with other, visible sources. Compared to simply correlating measured BSPs with BSPs of many simulated beats (Potyagaylo et al., 2016b), the new method has two main advantages: First, the basis generalizes to some extent from the individual activities used as input and thus also allows us to represent excitation patterns that were not among the training data. Second, measured and simulated BSPs do not have to be aligned in time, as the basis contains all time shifts. Localization results obtained by correlating with all BSP patterns in the training data and finding the one with the maximum correlation coefficient are included in the Supplementary Material. Although this approach works comparably well if the CV matches the CV in the training data, it performs considerably worse than the BSP basis reconstruction for non-matching CVs.

One remaining question regarding the creation of basis vectors is how to best set up the training data, so that each basic activation pattern occurs equally often in the data matrix. In the current geometry, the left and right atrium are only connected by one bridge representing the Bachmann's bundle. This results in almost the same activation pattern on the opposite atrium, once the excitation has passed the bridge. This activation pattern (corresponding to pacing at the bridge) occurs disproportionately often in the data matrix and will therefore best be represented by the resulting basis vectors. Reconstructions might therefore be biased toward this specific pattern. In geometries with multiple bridges, this effect is not as pronounced. Limiting the information in the data matrix to a time window around the beginning of excitations is expected to optimize the basis for better reconstructions near the excitation origin.

4.1. Limitations

There are several limitations of this study that need to be acknowledged. The most important one is that atrial conductivities were assumed to be homogeneous and isotropic. Although we think this is a logical first step in evaluating the new reconstruction method systematically, its performance has to be studied with atrial anisotropy and several atrial geometries. This has to be done in a future work. However, exemplary reconstructions for four ectopic foci in a highly anisotropic model of the atria are included in the Supplementary Material. These results indicate that the method improves reconstructions over TikhGS for anisotropic spread of excitation as well, although the reconstruction quality does decrease compared to the isotropic case. Another aspect not considered is fibrosis, which may hamper the localization of atrial ectopic foci (Godoy et al., 2018). Finally, further studies are needed to evaluate the sensitivity to imperfect geometries and conductivities.

4.2. Outlook

We are planning to study the new reconstruction method with clinical measurements during (focal) ventricular tachycardia. We hope to make use of scar-related information from late gadolinium enhancement MRI during the creation of basis vectors. For non-focal arrhythmias, it would be interesting to see whether including non-focal activities in the basis is beneficial.

Author Contributions

SS developed the BSP basis reconstruction method, designed and conducted the simulation study and wrote the manuscript. AW created atrial fiber orientations for the simulations in the Supplementary Material. OD provided technical expertise and contributed during study design and manuscript preparation.

Funding

This study was supported by the German Research Foundation (DFG) under the grant DO 637/21-1.

Conflict of Interest Statement

The 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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01126/full#supplementary-material

Movies of atrial ectopic foci reconstructions and exemplary reconstructions of rotational activity (simple atrial fibrillation from Figuera et al., 2016) can be found in the Supplementary Material.

References

Cluitmans, M., Karel, J., Bonizzi, P., Volders, P., Westra, R., and Peeters, R. (2018). Wavelet-promoted sparsity for non-invasive reconstruction of electrical activity of the heart. Med. Biol. Eng. Comput. doi: 10.1007/s11517-018-1831-2. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Cluitmans, M. J. M., Clerx, M., Vandersickel, N., Peeters, R. L. M., Volders, P. G. A., and Westra, R. L. (2017). Physiology-based regularization of the electrocardiographic inverse problem. Med. Biol. Eng. Comput. 55, 1353–1365. doi: 10.1007/s11517-016-1595-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Courtemanche, M., Ramirez, R. J., and Nattel, S. (1998). Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am. J. Physiol. 275(1 Pt 2), H301–H321.

PubMed Abstract | Google Scholar

Dubois, R., Labarthe, S., Coudière, Y., Hocini, M., and Haïssaguerre, M. (2012). “Global and directional activation maps for cardiac mapping in electrophysiology,” in Computing in Cardiology (Kraków), 349–352.

Google Scholar

Duchateau, J., Potse, M., and Dubois, R. (2017). Spatially coherent activation maps for electrocardiographic imaging. IEEE Trans. Biomed. Eng. 64, 1149–1156. doi: 10.1109/TBME.2016.2593003

PubMed Abstract | CrossRef Full Text | Google Scholar

Figuera, C., Suárez-Gutiérrez, V., Hernández-Romero, I., Rodrigo, M., Liberos, A., Atienza, F., et al. (2016). Regularization techniques for ECG imaging during atrial fibrillation: a computational study. Front. Physiol. 7:466. doi: 10.3389/fphys.2016.00466

CrossRef Full Text | Google Scholar

Ghodrati, A., Brooks, D., Tadmor, G., and MacLeod, R. (2006). Wavefront-based models for inverse electrocardiography. IEEE Trans. Biomed. Eng. 53, 1821–1831. doi: 10.1109/TBME.2006.878117

PubMed Abstract | CrossRef Full Text | Google Scholar

Godoy, E. J., Lozano, M., García-Fernández, I., Ferrer-Albero, A., MacLeod, R., Saiz, J., et al. (2018). Atrial fibrosis hampers non-invasive localization of atrial ectopic foci from multi-electrode signals: a 3D simulation study. Front. Physiol. 9:404. doi: 10.3389/fphys.2018.00404

PubMed Abstract | CrossRef Full Text | Google Scholar

Greensite, F., and Huiskamp, G. (1998). An improved method for estimating epicardial potentials from the body surface. IEEE Trans. Biomed. Eng. 45, 98–104.

PubMed Abstract | Google Scholar

Messnarz, B., Tilg, B., Modre, R., Fischer, G., and Hanser, F. (2004). A new spatiotemporal regularization approach for reconstruction of cardiac transmembrane potential patterns. IEEE Trans. Biomed. Eng. 51, 273–281. doi: 10.1109/TBME.2003.820394

PubMed Abstract | CrossRef Full Text | Google Scholar

Oostendorp, T. F., van Oosterom, A., and Huiskamp, G. (1989). Interpolation on a triangulated 3d surface. J. Comput. Phys 80, 331–343.

Google Scholar

Pernod, E., Sermesant, M., Konukoglu, E., Relan, J., Delingette, H., and Ayache, N. (2011). A multi-front eikonal model of cardiac electrophysiology for interactive simulation of radio-frequency ablation. Comput. Graph. 35, 431–440. doi: 10.1016/j.cag.2011.01.008

CrossRef Full Text | Google Scholar

Potyagaylo, D., Cortés, E. G., Schulze, W. H. W., and Dössel, O. (2014). Binary optimization for source localization in the inverse problem of ECG. Med. Biol. Eng. Comput. 52, 717–728. doi: 10.1007/s11517-014-1176-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Potyagaylo, D., Dossel, O., and van Dam, P. (2016a). Influence of modeling errors on the initial estimate for nonlinear myocardial activation times imaging calculated with fastest route algorithm. IEEE Trans. Biomed. Eng. 63, 2576–2584. doi: 10.1109/TBME.2016.2561973

PubMed Abstract | CrossRef Full Text | Google Scholar

Potyagaylo, D., Loewe, A., van Dam, P., and Dössel, O. (2016b). “ECG imaging of focal atrial excitation: evaluation in a realistic simulation setup,” in Computing in Cardiology, Vol. 43 (Vancouver, BC), 113–116.

Google Scholar

Pullan, A. J., Cheng, L. K., Nash, M. P., Ghodrati, A., MacLeod, R., and Brooks, D. H. (2010). “The inverse problem of electrocardiography,” in Comprehensive Electrocardiology, eds P. MacFarlane, A. van Oosterom, O. Pahlm, P. Kligfield, M. Janse, and J. Camm (London: Springer), 299–344.

Google Scholar

Schuler, S., Potyagaylo, D., and Dssel, O. (2017). “ECG imaging of simulated atrial fibrillation: imposing epi-endocardial similarity facilitates the reconstruction of transmembrane voltages,” in Computing in Cardiology, Vol. 44 (Rennes).

Schulze, W. H. W. (2015). ECG Imaging of Ventricular Activity in Clinical Applications. Ph.D. thesis, Institute of Biomedical Engineering, Karlsruhe Institute of Technology1.

Simms, H. D. J., and Geselowitz, D. B. (1995). Computation of heart surface potentials using the surface source model. J. Cardiovasc. Electrophysiol. 6, 522–531.

PubMed Abstract | Google Scholar

van Oosterom, A. (1999). The use of the spatial covariance in computing pericardial potentials. IEEE Trans. Biomed. Eng. 46, 778–787.

PubMed Abstract | Google Scholar

Wang, L., Zhang, H., Wong, K. C., Liu, H., and Shi, P. (2010). Physiological-model-constrained noninvasive reconstruction of volumetric myocardial transmembrane potentials. IEEE Trans. Biomed. Eng. 57, 296–315. doi: 10.1109/TBME.2009.2024531

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamashita, Y., and Geselowitz, D. B. (1985). Source-field relationships for cardiac generators on the heart surface based on their transfer coefficients. IEEE Trans. Biomed. Eng. 32, 964–970.

PubMed Abstract | Google Scholar