Non-invasive reconstruction of dynamic myocardial transmembrane potential with graph-based total variation constraints

Non-invasive reconstruction of electrophysiological activity in the heart is of great significance for clinical disease prevention and surgical treatment. The distribution of transmembrane potential (TMP) in three-dimensional myocardium can help us diagnose heart diseases such as myocardial ischemia and ectopic pacing. However, the problem of solving TMP is ill-posed, and appropriate constraints need to be added. The existing state-of-art method total variation minimisation only takes advantage of the local similarity in space, which has the problem of over-smoothing, and fails to take into account the relationship among frames in the dynamic TMP sequence. In this work, the authors introduce a novel regularisation method called graph-based total variation to make up for the above shortcomings. The graph structure takes the TMP value of a time sequence on each heart node as the criterion to establish the similarity relationship among the heart. Two sets of phantom experiments were set to verify the superiority of the proposed method over the traditional constraints: infarct scar reconstruction and activation wavefront reconstruction. In addition, experiments with ten real premature ventricular contractions patient data were used to demonstrate the accuracy of the authors’ method in clinical applications.

1. Introduction: Using electrophysiological imaging (ECGI) to depict the electrophysiological information within the heart has become a research hotspot for the diagnosis and therapy of heart disease. By selecting appropriate source and volume conductor model, one can reconstruct the potential distribution or the excitation propagation on the surface of the heart or in the three-dimensional (3D) ventricles, and then determine the location and extent of the lesion based on the anomalous characteristics that appeared [1,2].
Common sources include transmembrane potential (TMP) [3], endo-and epicardial potential (EEP) [4], and activation time [5]. Among them, TMP was selected as the research target of this paper since it provides the most abundant physiological information and allows stronger prior information constraints [6]. However, it is ill-posed to reconstruct myocardial TMP from body surface potential (BSP) due to the following two points: (i) the dimensions of known quantities do not match those of unknown quantities (for accuracy, the number of unknown points on the heart is much larger than the number of surface leads). (ii) According to the Helmholtz's equivalent double-layer principle of the electromagnetic field, an infinite number of intramural solutions fit the same electrical field on the surface. Hence, researchers have made great efforts to add appropriate constraints on the solution space, called regularisation [7].
The most classic method is Tikhonov regularisation [8], which imposes a neighbourhood smoothing constraint and provides a solution with compromise accuracy. In recent years, combined with the sparsity and piece-wise smoothness of the potential on the heart, reconstruction methods based on sparse expression have been proposed, such as total variation (TV) minimisation [9]. These methods impose a constrain term with L1 norm form, which performs better than the L2 norm form in maintaining the steeply changing area of potential. However, they can only process one frame of data at the same time. If the whole sequence is to be solved, the computational time will increase linearly with the length of the sequence. The low rank and sparse constraint assume that the solution is composed of low-rank background and sparse foreground [10]. This method treats the solution of the whole sequence as a matrix, which takes into account the time dependence of the solution. However, singular value decomposition (SVD) in high-dimensional solution space in each iteration also makes the method computationally expensive.
With the rapid development of deep learning in the medical field, researchers try to solve the inverse problem of ECG by data-driven methods. Ghimire et al. trained a large number of simulation data to learn the general transformations of BSPs to TMP [11], and discussed the influence of different autoencoder architectures on the results [12]. However, there is no clinical method for directly measuring cardiac TMP values. Therefore, no real patient data participates in the training.
In order to make use of time dependence while avoiding frequent singular value decomposition, and taking into account the piecewise smoothness and sparsity in space, we propose a novel method, the graph-based TV reconstruction. This method improves the classical TV minimisation. We regard the records on each point in the sequence as a graph signal, making full use of the temporal and spatial correlations of the myocardial TMP distribution. Specific contributions are as follows: (i) We construct a graph structure on the preliminary estimation of the solution which is based on the spatial surface smoothing hypothesis. This structure constructs a similarity relationship among the whole 3D myocardium, making full use of the underlying nonlocal features. (ii) Different penalty weights are given to nodes with different similarities, which enhance intra-class consistency and inter-class differences. (iii) The criterion for judging similarity is the Euclidean distance between TMP value vectors over a period of time at each node. Therefore, the temporal correlation among dynamic TMP sequence frames is taken into account. This makes the method more robust than that based on a single frame, and less affected by noise.
(iv) The mentioned method solves the whole sequence at the same time, avoiding the computational time increases linearly with the length of the sequence.
(v) The mentioned method is based on the spatio-temporal distribution characteristics of the solution, without the need for training based on large amounts of data.
2. Methodology 2.1. TMP imaging model: The forward relationship between TMPs and BSPs can be modelled as the following linear equation by applying the quasi-static approximation of electromagnetic theory [13]: in which F [ R m×t represents the m-lead BSPs of a length t, H [ R m×n is the transform matrix which contains the geometry and conductivity information in the heart-torso structures, obtained by finite element method (FEM) or boundary FEM (BEM). U [ R n×t is the n-dimensional TMPs of a length t. t represents the information of time dimension. The equation denotes the forward mapping relationship between dynamic TMP sequence and dynamic BSP sequence. To overcome the ill-posedness of this problem, we start from the characteristics of the spatio-temporal distribution of dynamic TMP sequence to find reasonable constraints for the solution. First, we considered the distribution of TMP on the myocardium during infarction or myocardial ischaemia cases. At the ST segment of the electrocardiogram (ECG), the ventricular cells are in a plateau phase in which all depolarisations are completed and sharp repolarisation has not yet begun. The potential distribution is flat, which corresponds to the condition of the normal heart. For infarcted or ischaemic cells, their action potentials are reduced due to the decreased or absent cell activity. This is prominent in a flat TMP distribution. There is a clear boundary between the normal area and the lesion area, as shown in Fig. 1.
For the case of the activation sequence, there are many nodes activated at the same time, but they are not necessarily concentrated near the same area. The TMP value changes for nodes activated at the same time are similar for a period of time. Fig. 1c shows this feature. This example is the spatial distribution of TMP at about 120 ms after the onset of pacing rhythm. The red part indicates high potential. The white line is the active wavefront and is annular in the myocardium. We can see that the nodes j and k have approximate values and may be considered adjacent in the graph, although they are separated in space.
The piecewise smoothing and gradient sparse feature of the TMP distribution prompts us to find a constraint to enhance the similarity between nodes in similar classes, and to distinguish between different classes of nodes. This makes the border clearer and can accurately indicate the lesion and the activation wavefronts.

Constraint regularisation:
The TV method is the most classic sparse constraint method. Comparing Figs. 1a and b, we can see that in matrix U, the same values are not simply clustered together. This is caused by the way the FEM models the heart. The adjacent nodes in the 3D space are not spatially adjacent in the 2D matrix U expanded by the serial number. The classic TV method uses adjacency matrix in 3D space to find spatial adjacent nodes and minimise the TV among them. However, this brings several problems that limit the improvement of precision: (i) The adjacency matrix of 3D space is the only prior to this problem, so it puts high requirement on the accuracy of cardiac modelling. (ii) It can only use the information between spatial adjacent nodes, similarities between those nodes which are not adjacent in space, but still belong to the same class are not available. (iii) It only considers the situation at a certain moment, without using the time correlation of dynamic sequences. This inspired us to propose a graph-based TV constraint.
2.2.1 Neighbourhood-smoothness estimation: First, we apply a spatial neighbourhood-smoothing constraint to overcome the mathematical ill-posedness of the problem. This low-resolution solutionÛ will be used as the initial input for the subsequent optimisation algorithm.
L is the surface Laplacian operator, which utilises the spatial connection information of the nodes in the 3D heart model.

2.2.2
Graph-based TV constraint: In recent years, graph-based regularisation methods have attracted much attention in image and manifold processing, for the graph structure extracts and enhance the self-similarity of the signal [14]. Similar ideas have also been applied to medical image processing [15,16]. A graph structure contains three core constituent factors: node, edge and weight matrix W. Adjacent nodes are connected by edges. Here judging whether adjacent or not is based on the similarity of the values on the nodes, not just the spatial connection which classical total variation takes as a priori. The degree of similarity of two nodes i, j is measured by the Euclidean distance: The vector u represents the value on one node, which is, in our case, a record of the TMP values overtime on a single node of the heart. A smaller ℓ indicates a high consistency between two nodes, meaning that they are more likely to be in the same state, e.g. both infarct cells or cells activated at close time. The k-nearest neighbour (KNN) search algorithm finds the distance for each pair of nodes, then returns the nearest k ones for each node, which are considered to be adjacent. The graph TV constraint can be written as where n denotes the total number of nodes representing the heart. n i is a collection of all the neighbours of node i, the result from KNN. For a pair of adjacent nodes with a short edge (small ℓ), we apply a large weight to minimise the gradient between them. For those nodes that have a large dissimilarity, we reduce the weight between them. Based on this, we choose the following form of weight matrix definition with parameter σ being the average distance of the connected nodes: 2.3. Optimisation algorithm: Our reconstruction problem can be written as the following constraint minimisation form: where F represents the BSPs of a length t. H is the transform matrix. U is the transmembrane potentials of a length t. m is a positive weighting parameter used to balance data fidelity term with regularisation term. The problem (6) can be solved using forward-backward primal-dual method [17,18]. This algorithm minimises a nondifferentiable function by combining a gradient descent step (forward) with a proximal point step (backward). We split the target problem (6) into two sub-functions: h U ( ) is the data fitting term with a Lipschitz continuous gradient As for g LU ( ), L is a linear operator, which in our case the graph gradient operator ∇ G . Since the L1 norm is non-differentiable, we calculate its proximal operator by the soft-thresholding method [19] prox sg x By alternately solving the forward gradient and backward proximal problems, the algorithm finally converges to the final optimal solution. To ensure the convergence of the algorithm, the step size parameters should be chosen rigorously. The Lipschitz constant b of h U ( ) equals to 2 H 2 . We set t = 1/b, g = 0.99. The algorithm converges under the condition that (see Fig. 2) We have tested the reconstruction quality of different m values within a certain range, for there are no mature ways to decide the best choice for this regularisation parameter. Specific details will be discussed in Section 4.1.

Experiments:
We designed two sets of simulation experiments and one set of real patient experiments to verify the superiority and accuracy of our method. The ground truth and heart model of the simulation experiments are from the 'Experimental Data and Geometric Analysis Repository' (EDGAR) database [20]. This is an Internet-based archive of curated data that is freely distributed to the international research community for the application and validation of electrocardiographic imaging (ECGI) techniques.
3.1. Infarct scar reconstruction: In this section, we set the infarct nucleus locating at eight different segments to explore the performance of our method in reconstructing different infract locations. To test the robustness of the method to different noises, we add 5, 10, 15, 20, 25, 30 dB noise to the simulated BSP. We focus on comparing the performance between the graph-TV method and the classic TV method. The TV method we choose is Iteratively Re-weighted Minimisation of TV (IRTV) [10], which has been used to solve the myocardial TMP inverse problem.
We set the TMP value of the infarct site to −84 mV and the healthy site to 27 mV as the reference ground truth, with sequence length to 100 ms. It can be found from Fig. 3 that the infarct area reconstructed by Graph-TV method has a higher consistency with the ground truth, with the scar shape being closer to the reference and the edge representation being more accurate. In contrast, the Tikhonov method provides an overly smooth solution, and its location has some artefacts. The IRTV method eliminates artefacts and provides a higher quality solution. Fig. 4 shows the numerical analysis results for 8 locations, total 48 cases when the body surface potentials were disturbed by noise in the range of 5-30 dB. Since our method considers the TMP distribution over the entire sequence, it provides better performance than the frame-by-frame calculation methods.
The accuracy of the method was evaluated by comparing the correlation coefficient (CC) and the relative error (RE) between   the ground truth (subscript t) and the reconstructed TMP value (subscript r).
3.2. Activation wavefronts reconstruction: This part of the experiment was designed to verify the accuracy of our method in reconstructing the activation propagation sequence. Eight ventricular pacing transmembrane potentials were interpolated on the coarse tetrahedral mesh. The depolarisation phase of pacing can reveal the abnormality in the conduction pathway, which is helpful in locating the ectopic site or the underlying pathology.
In the initial phase of pacing, the Tikhonov regularisation method is almost impossible to indicate the location of the pacing site. The L1-norm based sparse constraint makes the IRTV method perform better at the initial and final phases of the pacing. However, the position of the propagation wavefront cannot be accurately captured. The result of Graph-TV method has the highest concordance to reference ground truth during the whole process of propagation, with the closest wavefront shape. Fig. 5 shows the statistical results of the reconstruction at different pacing sites. Each sequence length is 150-250 ms. Outliers often appear at the beginning and end of the sequence. In this time period, the number of activated nodes (or inactivated nodes) is small, and accurate reconstruction is more difficult than others. The results of reconstruction at the septum are worse than those at other locations, for this location hides the deepest from the body surface.
As another property describing the process of excitement propagation, the activation time map is used to visually indicate the order in which cardiomyocytes are activated. The earliest point of the activation time map is considered as the beginning point of pacing. In the case of Fig. 6, the location errors of listed three methods are 15.63 mm (Tikhonov-2), 9.95 mm (IRTV) and 7.44 mm (Graph-TV).
3.3. Real patient experiment: In this section, we collected the data of ten premature ventricular contractions (PVCs) patients, and to non-invasively locate the lesion by the mentioned three methods to explore the potential of the methods in clinical application. We locate the pacing site by reconstructing the high-potential area in the early stages of pacing. The 64-lead ECG recorded the distribution of body surface potentials overtime during arrhythmia with a sampling frequency of 2 kHz, i.e. F, in our problem. The patient wears lead electrodes for a CT scan, whereby the position of the lead on the torso can be obtained to establish a torso model. We sliced the CT image containing the relative position of the heart-torso along the short axis of the heart, then labelled the contours of the epicardium, ventricular and right ventricular outflow tracts to obtain a 3D cardiac model. Finally, the two models are registered in the same coordinate system, and the anisotropic conduction information is combined to obtain a personalised TMP-BSP transfer matrix H. The blue part of Fig. 7 shows this process. The Ensite3000 diagnostic results of invasive surgery records are considered as the gold standard to verify the accuracy of non-invasive reconstruction results. Nine cases of the ectopic pacing site were located in the right ventricular outflow tract and one in the left ventricular apex.
The occurrence of PVC is usually accompanied by a large and deformed QRS wave, hence we select this segment for reconstruction. Table 1 lists the clinical diagnosis and reconstruction results of these ten patients. Fig. 8 shows the 37th lead ECG, the Ensite3000 diagnostic screenshot and the distribution of TMP on the 3D heart model of cases 1, 2 and 3. We reconstructed the QRS part (about 140-180 ms) and selected the early time nodes for 3D display. The high level (red) indicates the location of the pacing. Combined with Table 1, it can be seen that for patients with complete data collection (with CT scan data with 64 lead positions, so that a personalised cardiac torso model can be established), the proposed method can accurately determine the location of the pacemaker. For cases where CT images were not acquired, we used a general model, which made the reconstruction results of the two cases slightly offset (see cases 7 and 10). The personalised model contains patient-specific geometric location information and has an important impact on the quality of the reconstruction results. Experiments have shown that our method can clearly indicate the boundary between the active site and the resting site, so as to accurately locate the lesion.  final result. The k value affects the accuracy of the graph structure, while μ is an important factor that weighs the data fidelity and sparse constraints.
We first determine the value of k and then test the effect of different μ values based on a fixed k. After test we found that the value of k has little effect on the result compared to the value of μ. Fig. 9a shows that the gap between CC and RE is less than 0.005 when k changes from 0 to 20. Considering that a large k value will produce a significant oscillation near the optimal solution, while a too small k increases the calculation time, we finally choose k = 5 as our experimental parameter. Fig. 9b shows that the optimal value of μ increases as the noise level increases. This is in line with the general denoising problem we are familiar with.

Computing time:
Here we discuss the computational time cost of the three methods mentioned. Table 2 shows the results of one set of activation wavefronts reconstruction experiments. The pacing site was located at LV-APEX, with a sequence length of 221 ms. The number of mesh nodes of the heart model is 2223. The experimental environment is MATLAB R2016a, with a 3.2 GHz processor and 16 GB RAM.
The Tikhonov method takes the shortest time because it does not require iterative computation. The calculation of both single frame and sequence involve one singular-value decomposition. IRTV is an iterative algorithm, and the regularisation parameter is adjusted in each iteration to converge to the optimal solution. Each frame is computed independently, so the time-consuming increases almost linearly. Graph TV is an iterative algorithm, which converges to the optimal solution by solving the primal and dual problems alternately. Because the algorithm can solve the whole sequence at the same time, it has a great advantage in speed compared with the IRTV algorithm.

Conclusion:
This paper proposes a graph-based TV constraint to solve the inverse problem of ECG imaging. Simulation experiments show that the proposed method has higher reconstruction quality than the Tikhonov method based on L2-norm and the classical TV method based on spatial adjacency constraints. Real patient experiments also illustrate the accuracy of the method of clinical application. Compared to the frame-by-frame reconstruction method, the sequence-based solution overcomes the disadvantage that the computation time increases linearly with the sequence length. Due to the lack of sufficient training data, this work has not been compared with the deep learning method for the time being. In the future, we will consider combining deep learning and traditional optimisation methods to improve the performance of TMP reconstruction in a data-driven manner.