A New spatio-temporal regularization approach
to the electrocardiographic inverse problem

B. Messnarz^1,2, B. Tilg¹, R. Modre¹, G. Fischer¹, F. Hanser^1,2, P. Wach^2 
1University for Health Informatics and Technology Tyrol
Innrain 98, 6020 Innsbruck, AUSTRIA
²Institute of Biomedical Engineering, Graz University of Technology
Krenngasse 37, 8010 Graz, AUSTRIA

Abstract: In this study we propose a new spatio-temporal regularization approach for the inverse problem of electrocardiography. In a single beat potential based method we consider all time instants during atrial or ventricular depolarization at once. Combined spatial and temporal regularization is achieved with linear operators. This leads to a large scale constrained quadratic optimization problem which is solved with an interior point optimizer. The new method would be capable of localizing ischemic regions and infarcts.             

INTRODUCTION

Noninvasive imaging of the heart’s electrical activity from electrocardiographic mapping data (ECG) may become an important clinical diagnosis tool which should provide cardiologists to decide for a proper therapy (e.g. ablative, drug) for the individual patient.

A model based imaging method includes a model which relates the electrical sources within the human heart to the body surface potential, the ECG. This model constitutes the solution to the so-called forward problem. The source formulation is based on the bidomain theory [1]. Assuming electrical isotropy for the heart muscle’s conductivity the spatio-temporal transmembrane potential (TMP) distribution on the heart’s surface is the primary source. Therefore, the relationship between the TMP and the ECG is given by a surface integral equation, a so-called Fredholm integral equation of second kind. For a numerical treatment the boundary element method (BEM) [2] is applied for this type of equation, which leads to matrix equations.

The challenge of the inverse problem of electrocardiography is to reconstruct the electrical sources from the measured ECG employing the forward model. Unfortunately this problem is exponentially ill-posed and, therefore, straightforward solution procedures would not lead to physiologically meaningful results. Regularization methods have to be employed where additional a-priori information about the desired solution, e.g. the TMP, is implemented in the solution procedure.  Common regularization approaches which exploit only spatial constraints on the solution would show only weak performance. For an improvement also the time domain has to be considered resulting in a spatio-temporal regularization procedure. In the activation time imaging method [3] a default template function for the temporal shape of the TMP at each source point is impressed and only the onset (activation time) is determined in the inverse procedure. However, in some pathologic cases (e.g. ischemic regions, infarcts) this regularization is too restrictive. We propose a new potential based method where the temporal course of the TMP is less predefined. 

METHODS

Forward problem

Applying the BEM the resulting forward solution is a discrete model expressing the potential at each electrode position on the torso (m electrodes) as a linear combination of the transmembrane potentials at all nodes of the heart’s surface (n nodes):

(1)

L denotes the (m × n) lead-field or transfer matrix, x(t) is the (n × 1) vector of the transmembrane potential at time t and d(t) is the (m × 1) vector of the body surface potential.

Inverse Problem - Regularization

In the spatial domain of the TMP we require a certain smoothness of the spatial distribution on the heart’s surface. In the temporal domain we assume that during depolarization the TMP at each source point is monotonically increasing, which is a realistic approximation to the physiological behavior of action potentials. In addition, we restrict the amplitude of the TMP within lower (lb) and upper (ub) bounds.  To handle the time instants during atrial or ventricular depolarization in the interval [t_b, t_e] all together we stack the TMP and ECG vectors and rewrite the forward equation (1) in block-diagonal form:

(2)

We abbreviate the (mT × nT) lead-field block-diagonal matrix by , the (nT × 1) TMP vector by and the (mT × 1) ECG vector by . The number of time instants within the interval [t_b, t_e] is given by T.

To fulfill the above mentioned requirements we follow the style of Tikhonov regularization [4] and minimize the functional  

(3)

subject to

(4)

The first term in (3) is the model error while the second term is the spatial regularization term with denoting the block-diagonal matrix of an approximation to the surface Laplacian.

The weight of both the model error and the regularization term is determined by the regularization parameter λ. The L-curve criterion [5] yields the optimal value of λ and thus the final solution.  Simultaneous temporal regularization is introduced by the side constraint where G is a discrete approximation of the first time derivative operator. This constraint ensures increasing of the TMP within the lower and upper bound during depolarization as is normally the case. We mention that the TMP need not to increase from the lower to the upper bound but can start increasing from a higher level as is the case in ischemic regions. In addition, the TMP can also remain at a constant level, especially at the zero level which is important for infarcted regions.     

Minimizing the functional (3) with respect to the constraints (4) is a constrained convex (more precisely quadratic) optimization problem which is, due to the enormous number of unknowns (nT is in the range of 5∙10⁴), computationally very expensive. The optimization task is achieved by an interior point optimizer [6] which is capable of exploiting the sparse nature of the underlying matrices.

RESULTS

The performance of our new inverse method is demonstrated by comparing reconstructed TMP-patterns with the underlying reference patterns. Within the scope of this work we give just an example for a reference depolarization sequence simulating a WPW-syndrome in applying a realistic anisotropic FEM heart model [1].

Fig. 1: Reference (upper panel) and reconstructed (lower panel) TMP patterns for a WPW-syndrome at time 48ms (left) and 84ms (right) from a cranial view. The dark region indicates activated tissue. The excitation starts from the right ventricular outflow tract.

For 159 discrete time samples with a sampling interval of 1ms (ventricular depolarization) ECG-data were calculated at 62 electrode locations at the torso’s surface utilizing the lead-field matrix. Before the inverse algorithm was applied to the simulated data Gaussian noise was added in order to obtain a signal-to-noise ratio of 35 dB. In addition, for reducing computational costs we increased the sampling interval to 3ms leading to 53 time instants (T=53). According to 628 nodes of the ventricles surface and 53 time instants we have 33284 unknowns in the optimization procedure. The computing time on a PC AMD K7 1.2 GHz with 1 GB RAM was about an hour for each value of the regularization parameter λ.

The reference and reconstructed patterns at two different time instants are shown in Fig. 1. We can observe that the reconstructed pattern is qualitatively in good accordance with its reference. Numerically we obtained a correlation coefficient of 0.94 and a relative error of 0.26. The relatively high value of the relative error is explained by the smearing-out-effect due to linearity of the spatial regularization operator . Similar values of the correlation coefficient were obtained with other reference patterns, too.   

DISCUSSION

The proposed linear method showed good results for single-beat reconstruction of TMP-patterns. The less strict temporal regularization in comparison with the activation time approach would enable localization of ischemic regions and infarcts. The latter will be subject of further research. 

Acknowledgments: This work was supported by the START Y144-INF program funded by the Federal Ministry of Education, Science and Culture, Vienna, Austria.

REFERENCES

[1]  G. Fischer, B. Tilg, R. Modre, G.J.M. Huiskamp, J. Fetzer, W. Rucker, P. Wach, ”A bidomain model based BEM-FEM formulation for anisotropic cardiac tissue,“ Ann. Biomed. Eng., vol. 28, pp.1229-1243, 2000

[2]  G. Fischer, B. Tilg, P. Wach, R. Modre, U. Leder, H. Nowak, “Application of high-order boundary elements to the electrocardiographic inverse problem,” Comput. Meth. Prog. Biomed., vol. 58(2), pp. 119-131, 1999

[3]  R. Modre, B. Tilg, G. Fischer, P. Wach, “An iterative algorithm for myocardial activation time imaging,” Comput. Meth. Prog. Biomed., vol. 64, pp. 1-7, 2001

[4]  F. Greensite, “The mathematical basis for imaging cardiac electrical function,” Crit. Rev. Biomed. Eng., vol. 22, pp. 347-399, 1994

[5]  P.C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev., vol. 34, pp. 561-580, 1992

[6]  E. D. Andersen, Y. Ye, “A computational study of the homogeneous algorithm for large-scale convex optimization,” Computational Optimization and Applications, vol. 10, pp. 243-269, 1998