Introduction

The purpose of this paper is to investigate the computational performance for myocardial activation time imaging applying different electrocardiographic (ECG) and magnetocardiographic (MCG) mapping systems. Based on the bidomain model [1], myocardial activation time imaging from ECG and MCG mapping data is an inverse nonlinear ill posed problem. The primary electrical sources in the cardiac muscle are the spatio-temporal transmembrane potential distribution φ_m. Identifying diagnostically useful parameters such as the activation time τ, which describes the onset of φ_m, non linearity enters the problem, because of the nonlinear behavior between τ and φ_m (see  Figure 1). Thus a linearized iterative algorithm [2] which considers the nonlinear characteristic was used in estimating the activation time map. Additionally, an optimization method with a regularization scheme usually used in inverse electro- and magnetocardiography was also applied [3]. The numerical experiments presented were based on synthetic ECG and MCG mapping data for normal cardiac depolarization of the human ventricle.

Figure 1. The nonlinear ill posed electro- and magnetocardiographic inverse problem

Methods

A 43-year old male without structural heart disease underwent an electrophysiological study for paroxysmal atrial fibrillation. From this patient the 62-channel ECG map and the torso anatomy ((cine) MRI) was recorded at the University of California San Francisco, San Francisco, California, USA. The ECG mapping system (Amsterdam, The Netherlands) was an on-line portable dual-computer acquisition and analysis system with data transmission by optical fiber. A radiotransparent carbon electrode array was used to record unipolar ECG data from 62 torso sites - 41 electrodes on the anterior and 21 on the posterior chest - relative to the Wilson's terminal as a reference. The torso was imaged with a 1.5 Tesla MR scanner (Signa?, General Electrics). From the torso 40 axial and 10 coronal cross-sections were recorded with a slice thickness of 10 mm. Additionally, for modeling the surface of the heart the cardiac muscle was imaged in an ECG-gated (cine) breath-hold oblique imaging mode with 7 mm slice thickness and seven phases during the cardiac cycle with a phase delay of 35 ms. Vitamin E markers were used to determine the electrode positions from the axial MR scans. Additional markers were positioned on the sternum, xiphoid and on the left and right rib cage. A subset of these markers was used to couple the MR and ECG data. The actual coordinate system for the inverse formulation was equivalent with the MR machine coordinate system. A BEM with linear triangular elements including the heart, the lungs and the chest surface was the basis for modeling the patient's volume conductor (see Figure 2) [4]. The following conductivities of the torso compartments were used: 0.2 S/m (torso), 0.05 S/m (lungs). The entire volume conductor was meshed by 2168 BEs. The entire surface of the ventricle (es) was modeled with 1074 BEs (539 nodal points). The mesh of the sub surfaces (endocardium of the left ventricle (le) and of the right ventricle (re) and epicardium (ep)) consisted of 177, 152 and 210 nodal points. The cardiac geometry was collected during the end diastolic phase.

Different ECG and MCG mapping systems were used in the inverse problem. The ECG mapping systems (A) consisted of 62 (anterior and posterior) electrodes. Two synthetic 62 channel MCG mapping systems with coil positions close (20 mm (B) and 70 mm (C)) to the 62 electrode positions of (A) were created. As a combination of system (B) and (C), also a synthetic MCG mapping system (D) with 62 gradiometers with a baseline of 50 mm was designed. Furthermore, different MCG mapping systems with 37 (E) and 67 (F) channels were investigated (see Figure 2).

Figure 2. Boundary element torso model from an cranial (left panel) and a left lateral (right panel) view. The circles indicate the positions of the 62 electrodes of the ECG mapping system. The MCG mapping systems with 37 and 67 channels are indicated by plus signs and crosses, respectively.

Two different methods were used for myocardial activation time imaging: a linearized iterative algorithm [2] termed with method A and a conventional approach [3] (method B) usually used in inverse ECG and MCG problems.

Method A

The non linear inverse ill posed problem

Ft = D

(1)

in which F  is a nonlinear operator mapping the activation time t onto the ECG or MCG mapping data D can be solved by an iterative method [2]. This method is based on a regularization scheme, which extents the regularization method - Tikhonov's approach of second order - from the linear to the nonlinear case. For growing k-values the iteration method

tk+1 = tk + dtk

(2)

with the descent direction dt_k

dtk = ( FTk Fk + lk2DTD )-1 [ FTk (D - Ftk) - lk2DTDtk]

converges to a regularized approximation of the activation time pattern. F_k represents the Frechet derivative of F.

Method B

Method B consists of an optimization routine and a regularization technique which is usually used in estimating the activation time map. The functional which has to be minimized is given by

|| Ljm - D ||2 + l2  || Dt ||2 = min

(3)

L stands for the lead field, which describes the relationship between the transmembrane potential jm and the ECG or MCG mapping data D. The minimizing in (3) can be done by an optimization routine, e.g., the Newton-type algorithm E04UCF in the NAG library (NAG Ltd. UK). E04UCF is designed to solve nonlinear programming problems and minimizing of smooth nonlinear functions subject to a set of constraints on the variables.

In both methods, a second order Tikhonov regularization approach with the surface laplacian operator is applied to stabilize the inverse solution in the face of measurement noise and errors in modeling the patient's geometry. D represents an approximation of the surface Laplacian on the surface of the ventricle [5]. The regularization parameter l specifies the relative weight between the residual norm and the smoothing norm. The determination of l in (2) as well in (3)  is based on the L-curve method [6]. This method involves a parametric log-log scale plot of the  smoothing norm on the ordinate against the residual norm on the abscissa, with l as parameter. The corner of the characteristic L-shaped curve, the point of maximum curvature, reflects the proper l.

A reference activation time map shown in Figure 3 was estimated from single-beat ECG recordings. This activation time map, which is in accordance with the activation time pattern of the sinus rhythm in humans, was used to simulate synthetic ECG and MCG mapping data. Gaussian noise was added to these synthetic mapping data with a signal-to-noise ratio of 50dB in order to take into account the unavoidable noise, which occurs in real recordings. The activation time map on the endo- and epicardium is estimated applying method A and B. In both cases a high-quality starting activation time map is required to achieve satisfying results. This starting vector is calculated applying the critical point theorem [3].

Figure 3. Reference activation time map t [ms] on the epicardium (top panels) and on the endocardium (bottom panels). Left panels show an anterior and right panels a posterior view. Isochrones are shown at steps of 5 ms.

Results

Firstly, the two different methods for noninvasive myocardial activation time imaging were investigated estimating the activation time map from synthetic ECG mapping data applying a 62-channel ECG setup. The effective rank of the data matrix was found to be 21. The relative root mean square error relRMS and the correlation coefficient CC of the activation time maps on the entire surface (es), on the left endocardium (le), on the right endocardium (re) and on the epicardium (ep) are summarized in  Table 1. These maps were estimated using the critical point theorem, method A (iterative algorithm) with two different numbers of iteration, and method B (optimization method). The optimization routine E04UCF took 169 iteration steps for finding a stable inverse solution. The L-curve applying method B as well as the individual L-curves of method A at iteration step 1, 5, 10 and 15 are depicted in Figure 4. As can be seen in Table 1 both methods result in nearly the same inverse solution. Using the 62 channel MCG setups (B) and (D), similar results were obtained. The effective rank of the data matrix was found to be 20 and 19, respectively. The results are summarized in Table 2 and Table 4. The MCG mapping setups (C), (E) and (F) showed a lower effective rank of 16, 9 and 14. Thus, the relRMS of these activation time maps estimated applying method A are higher than the relRMS of the activation time maps estimated using (A), (B), or (D). Applying method B better results were obtained. The relRMS and the CC of the inverse solutions using MCG setup (C), (E), and (F) can be found in Table 3, Table 5, and Table 6.

Figure 4. L-curve plots at iteration step k = 1, 5, 10 and 15 applying method A and L-curve plot applying method B (747 iterations). The short solid line indicates the regularization parameter corresponding with the point of maximal curvature of the individual L-curve.

TABLE 2. The relative root mean square error relRMS and the correlation coefficient CC of the activation time maps on the entire surface (es), the left endocardium (le), the right endocardium (re) and the epicardium (ep) estimated applying mapping setup (B). The number in brackets indicates the number of iterations.
 

TABLE 3. The relative root mean square error relRMS and the correlation coefficient CC of the activation time maps on the entire surface (es), the left endocardium (le), the right endocardium (re) and the epicardium (ep) estimated applying mapping setup (C). The number in brackets indicates the number of iterations.
 

TABLE 4. The relative root mean square error relRMS and the correlation coefficient CC of the activation time maps on the entire surface (es), the left endocardium (le), the right endocardium (re) and the epicardium (ep) estimated applying mapping setup (D). The number in brackets indicates the number of iterations.
 

TABLE 5. The relative root mean square error relRMS and the correlation coefficient CC of the activation time maps on the entire surface (es), the left endocardium (le), the right endocardium (re) and the epicardium (ep) estimated applying mapping setup (E). The number in brackets indicates the number of iterations.
 

TABLE 6. The relative root mean square error relRMS and the correlation coefficient CC of the activation time maps on the entire surface (es), the left endocardium (le), the right endocardium (re) and the epicardium (ep) estimated applying mapping setup (F). The number in brackets indicates the number of iterations.
 

Conclusions

An activation time map reconstructed from measured 62-channel ECG mapping data, which is in accordance with the activation time pattern of the sinus rhythm in humans, was used as a reference shown in Figure 3. The effective rank of the data matrix was found to be 20. In order to get a similar effective rank from synthetic ECG and MCG mapping data, Gaussian noise with a signal-to-noise ratio of 50dB was added. As can be seen in  Figure 5, the singular values of the lead field matrix of the ECG mapping setup (A) and the MCG mapping setups (B) and (D) similarly decrease, which results in a similar effective rank of the mapping data matrix. A high-quality inverse solution was found even after 15 iterations applying method A. Using the mapping setups (A), (B) and (D), about 200 iterations were needed to achieve the final solution applying method B. The reconstruction from ECG and MCG mapping data results in nearly the same inverse solution, because of the similar effective rank of the mapping data matrix.

A larger distance between the MCG sensors and the primary sources as is the case in (C), results in a faster decrease of the singular values (see Figure 5) and thus in an effective rank of 16. As can be seen in Figure 2 the MCG setup (E) and (F) have a smaller observation space, which results in some sense in a higher spatial dependency on the individual sensors. Thus, the effective rank of the data matrix was further reduced to 14 (setup (F)) and 9 (setup (E)), respectively. The solution of the critical point theorem showed less quality. Because of the less-quality of the starting vector, the quality of the activation time pattern estimated applying method A could not be adequately improved after 15 iterations. Using MCG setup (F), more than 700 iterations in method B were necessary to achieve a useful inverse solution. Applying (E) the effective rank was too low resulting in a starting vector, which could not be further improved even using method B with more than 800 iterations.

We can conclude that in contrast to mapping systems with sensors rather uniformly distributed over the chest surface, the MCG mapping systems (E, F) show lower effective ranks and thus resulting in a less quality of the activation time map.
 

Figure 5. Normalized singular values of the lead field matrix for different ECG and MCG setups.

Acknowledgements

The authors are indebted to the Section of Cardiac Electrophysiology, Department of Medicine, University of California San Francisco (UCSF), California for making available ECG mapping and MRI patient data to this analysis. This project was supported by the Austrian Science Fund grant J1506-MED, grant P11448-MED and grant P14207-MED.

References