Chapter 12. Electrocardiographic Imaging of Epicardial Potentials

Electrocardiograms (ECGs) recorded simultaneously at many body surface sites can be used to construct a sequence of body surface potential distributions that correspond to cardiac electrical activity; this technique is referred to as body surface potential mapping (BSPM).1,2 The noninvasively acquired BSPM data can be used, in turn, along with a mathematical model that accounts for geometry and electrical properties of the thorax, to reconstruct electrical potentials on the epicardial surface. The problem of calculating epicardial potentials from recorded BSPM data constitutes one of the formulations of the inverse problem of ECG.3,4

The calculation of epicardial potentials from noninvasively acquired electrical and anatomic data is also referred to, quite aptly, as electrocardiographic imaging.5,6 The procedure requires BSPM data, accurate representation of the patient's thoracic geometry, routines for calculating transfer coefficients relating epicardial and body surface potentials,7 and reliable inverse-solution techniques for estimating epicardial potentials from BSPM data.3,4 Patient-specific anatomic data can be obtained by computed tomography (CT) or magnetic resonance imaging (MRI). ECG imaging has shown its potential in theoretical and experimental studies,5,6 and it promises to find its way into routine clinical applications.8,9

The aim of this chapter is to illustrate the use of ECG imaging in clinical cardiac electrophysiology, in particular as an aid to radiofrequency catheter ablation of ventricular arrhythmias. In this chapter, we introduce our methodology of BSPM, with particular attention to applications in clinical cardiac electrophysiology; we deal with the mathematical formulation of the relationship between potentials on the epicardial surface and body surface and show how this relationship is used to solve the inverse problem of ECG in terms of epicardial potentials; we describe the application of ECG imaging during catheter-ablation procedures in the clinical cardiac electrophysiology laboratory; and we discusses results reported in this chapter in the context of current clinical practice.

To obtain the BSPM data required for ECG imaging, 120 disposable radiolucent Ag/AgCl electrodes (Ref. 31.8778.26; Covidien, Dublin, Ireland) were placed on the patient's torso in 18 strips according to the standard Dalhousie configuration (Fig. 12–1). The electrode array was connected to the 128-channel Mark-6 acquisition system (BioSemi Inc., Amsterdam, the Netherlands), which was coupled via a fiber optics cable with a laptop computer. The ECG signals were displayed and stored with MAPPER software (Dalhousie University, Halifax, Nova Scotia, Canada).

The acquisition system amplified and filtered (bandpass, 0.025-300 Hz) analog ECG signals and sampled each channel at 2000 Hz with 16-bit resolution. The raw data were recorded typically for 15 seconds during sinus rhythm, pacing, or ventricular tachycardia (VT) and stored on a hard disk for subsequent analysis.

All ECG tracings were inspected to mark faulty leads, such as leads affected by poor skin-electrode contact, motion artifact, or inaccessibility of the chest area where the electrode was supposed to be placed (eg, the area occupied by external defibrillator pads). Such faulty leads were estimated using a three-dimensional Laplacian interpolation scheme.13 Processing of the ECG data took into account ECG signals obtained during sinus rhythm, pacing, or VT.14 In each case, a different semi-automated routine was written in MATLAB (Mathworks Inc, Natick, MA) to determine fiducial marks (ie, onsets/offsets of P wave, QRS complex, and T wave). For ECG data obtained during sinus rhythm, the baseline was determined by using the second derivative of the lead potential with respect to time, and the fiducial marks of ECG waves were determined by algorithm; these marks were subsequently revised, if necessary. For paced data, the pacing artifact was identified by thresholding the time derivative of the sum of all body surface ECGs, and then the beginning and end of each pacing cycle was determined by algorithm. For ECG data recorded during VT, the onset of each cycle was determined from the narrowest point between the upper and lower envelope of all body surface ECG tracings preceding a rapid change in their amplitudes (Fig. 12–2).

Next, each lead was averaged over a number of beats in order to improve the signal-to-noise ratio. Dynamic beat averaging was performed to include in the averaging only those beats that were similar enough to the selected template beat.

An example of four instantaneous BSPM distributions corresponding to activation/repolarization sequence initiated by epicardial pacing is shown in Fig. 12–3. Note that the color scale is adjusted to the instantaneous values of extrema; thus, spatial patterns of low-level potential distributions can be discerned as clearly as those with much higher potentials at the peak activation.

We have shown previously15,16 that BSPM data acquired during catheter ablation of scar-related VT can assist with the identification of VT exit sites. Patients had BSPM data recorded during episodes of VT and endocardial pacing, and the data were analyzed to assess the similarity between the BSPM sequences occurring during VT and those induced by pacing. We found that, in general, the nearer the pacing site was to the exit site, the greater was the similarity between the VT and paced BSPM sequences. This suggested that BSPM can provide a useful adjunct to standard pace mapping. We realized, however, that additional processing—in particular, the inverse calculation of epicardial potentials/isochrones—would enhance the usefulness of the data gathered by BSPM. The rest of this chapter briefly reviews our methodology of inverse calculations and gives examples of their possible clinical use.

The forward problem of ECG involves computation of body surface potentials from specified cardiac generators that represent the electrical activity of the heart. The problem can be formulated in different ways, depending on the assumptions regarding the equivalent sources for the cardiac electrical activity and a representation of geometric and electrical properties of the volume conductor in which these sources are embedded.7,17,18

Numerical methods for solving the forward problem of ECG for arbitrarily shaped domain-wise homogeneous volume conductors have been in use since the 1960s; initially, the forward solution was performed to calculate potentials on the torso surface from dipolar sources.19-21 An alternative approach is to consider the distributed double-layer source on the epicardial surface as the equivalent generator.22 Barr et al7 introduced the mathematical formulation for computing linear transfer coefficients relating electric potentials on the epicardial surface and those on the body surface. This model—consisting of a homogeneous volume conductor bounded by realistically shaped surfaces of the heart and outer torso—was adopted in calculations presented in this chapter.

Forward Problem in Terms of Epicardial Potentials

Let us consider the forward problem of ECG formulated as a calculation of the electric potentials ΦB = φ1B, … φmB at m area elements (or discrete points) on the body surface, SB, from the observed electric potentials ΦH = (φ1H, … φnH at n area elements (or discrete points) on the epicardial surface, SH. This leads to the boundary-value problem for Laplace equation in a homogeneous volume conductor that is bounded by smooth surfaces SB and SH.4,7 The final result, in matrix notation, is:

which expresses body surface potentials, ΦB, as a linear combination of epicardial potentials, ΦH, through the transfer matrix A. This forward relation is the basis for inverse-problem formulation in terms of epicardial potentials (see Equation 12–2).

Generation of Patient-Specific Torso Model

Calculation of the transfer matrix A for Equation 12–1 requires a discretized representation of surfaces SB and SH of the torso model; this is accomplished by defining a set of points connected to form a three-dimensional mesh that represents each surface as a tessellation of planar triangles. Over the last decade, advanced imaging modalities—such as CT and MRI—have made it possible to acquire anatomic scans with high resolution, and thus, accurate patient-specific torso models can be readily generated.

We used axial CT scans (Siemens Sonata, Erlangen, Germany) acquired in supine position to generate patient-specific torso models. The location of body surface electrodes was estimated from known constraints: in-strip (5 cm) and between-strip (>2.5 cm) interelectrode distance, and correspondence with anatomic landmarks (location of precordial leads with reference to intercostal spaces). Figure 12–4 shows an example of the estimated distribution of 120 electrodes for BSPM on the patient-specific torso surface.

The 120 patient-specific locations of electrodes were then used as vertices to create a triangulated closed surface of the outer torso. A method was developed to make this patient-specific triangulation compatible with the Dalhousie standard 352-node torso geometry. (The Dalhousie standard torso consists of 352 nodes, which are vertices of 700 triangular area elements of the body surface.) Figure 12–5 illustrates the generation of the customized 352-node torso.

Generation of the discretized epicardial surface involved two steps: segmentation and triangulation. CT images were imported into the Amira package (Mercury Computer Systems, Chelmsford, MA), and an oblique slicing plane parallel to the coronary sinus was defined to facilitate segmentation of the ventricles down to the left ventricular apex. The labeled ventricles were then triangulated. Figure 12–6 shows the discretized epicardial surface, the section of the torso, and the oblique cut of the CT data through the coronary sinus to facilitate segmentation and discretization. Multiple passes of spatial filtering were applied to render a smooth surface.

Inverse Problem in Terms of Epicardial Potentials

The objective of ECG imaging is to solve the inverse problem defined as the calculation of epicardial potential distributions from potential distributions acquired on the body surface (BSPM data). Whereas the mathematical solution of the forward problem of ECG is unique and depends continuously on the data (ie, the problem is "well posed"), the solution of the inverse problem may be neither unique nor depending continuously on the data (in which case, the problem is said to be "ill posed"), and small errors in the input data can generate unbound errors in the inverse solution. Thus, any attempted inverse solution usually requires that auxiliary constraints be imposed. The "ill-posedness" of the problem can be determined by examining the singular values of the transfer matrix A of Equation 12–1. Methods of mathematical regularization of ill-posed problems, most notably Tikhonov regularization,23 provide a key to their solution. Tikhonov regularization, applied to solving the inverse problem of ECG in terms of epicardial potentials, was our method of choice.

The inverse problem is to estimate in the Equation 12–1 epicardial potentials ΦH from the measured body-surface potentials ΦB and the transfer matrix A. Because both ΦB and ΦH are functions of time (body surface ECGs and heart surface electrograms, respectively), the inverse problem has to be solved at a series of time instants. At each instant, the solution for ΦH can be estimated by using Tikhonov regularization, which aims to solve a perturbed version of the least squares problem,

where ||.|| denotes the l2 norm, λ is the regularization parameter, and B is the regularizing operator. For a given A and appropriately selected λ > 0, a solution ΦH (λ) of the perturbed minimization problem in Equation 12–2 can be obtained from the generalized form of normal equations, which is derived by means of the generalized singular value decomposition.

Given a sequence of sampled BSPM data ΦB for a specific patient, together with the generalized singular value decomposition of the transfer matrix A and regularizing operator B for that patient (Appendix A), the solution for the epicardial potential distribution ΦH at each time step is calculated in two steps. First the regularization parameter λ is determined by one of several possible methods (Appendix B), and then, using Equation 12–A2 in Appendix A, the corresponding epicardial potential distribution ΦH is evaluated.

The ECG inverse solution can be applied to aid catheter ablation of scar-related VT, which is still one of the most challenging procedures in clinical electrophysiology. The majority of these tachycardias are poorly tolerated or difficult to induce, and they frequently transform to other tachycardia morphologies.24 Successful ablation of scar-related VT requires an understanding of its mechanism and of the underlying electroanatomic substrate. Such insight can be partly achieved by percutaneous endocardial or epicardial mapping with a single catheter, which is steered to multiple sites to define the substrate for arrhythmia during sinus rhythm.25-27 However, this alone is often not sufficient for predicting successful ablation sites; some tachycardias are "unmappable" with point-by-point mapping, and alternative methods for identifying exit sites for these tachycardias have to be sought.28 One of the frequently used methods is pace mapping in the peri-infarct zone.29 Our objective here was to investigate how BSPM and the ECG inverse solution, in conjunction with pace mapping, can be used to facilitate radiofrequency ablation of scar-related VT.

Three-Dimensional Electroanatomic Mapping

Catheter mapping and ablation were revolutionized by electroanatomic mapping. The CARTO electroanatomic mapping system (Biosense Webster, Inc., Diamond Bar, CA) employed in this study uses a magnetic sensor in the catheter tip to detect its location in the magnetic field created by magnets placed beneath the patient, and a triangulated ventricular surface is displayed by a system (see Fig. 12–8), allowing visualization of various electrical measurements in real time.24,30 This system can be used to perform sinus rhythm voltage mapping, which involves sampling of the endocardial or epicardial bipolar electrograms recorded from an electrode catheter at multiple locations.25 Low-amplitude bipolar electrograms have been shown to correspond to infarcted myocardium, with amplitudes typically less than 0.5 mV over the dense scar.27 The CARTO system can also perform activation mapping, which annotates the anatomic construct with local activation times.25

The limitations of point-by-point mapping of activation sequences have led to the development of pace mapping methods for assessing the location of the mapping catheter relative to the reentry circuit.31-33 The procedure involves successive stimulation at various sites and comparing the paced QRS patterns with the template pattern of the clinical VT. Distinct patterns of BSPM distributions were observed for ectopic ventricular activation sequences initiated at various endocardial pacing sites for patients with idiopathic VT,34 as well as for patients who have had a prior myocardial infarction.33,35

Propagation of electrical activation in the diseased myocardium is discontinuous, and small differences in the pacing site can induce different propagation patterns and resulting QRS complexes.31,33 Nevertheless, the pacing from a catheter located near the exit of a reentry circuit of VT usually produces a QRS morphology similar to that of VT, with a short stimulus-to-QRS delay.36 Pacing that produces a QRS matching VT after a delay (stimulus-to-QRS interval >40 ms) has been shown to mark a site in a reentry channel within an infarct scar that is considered a good candidate for ablation.36 Trying to approach the exit site using sequential pace mapping can be challenging and requires an intuitive interpretation of the 12-lead ECG. We hypothesized that this process could be aided by BSPM and ECG imaging.

Some intraoperative studies involving endocardial and epicardial mapping have shown that scar-related VT does not always arise in the subendocardium.37-40 This evidence was further supported by studies using hearts explanted from patients who have had a myocardial infarction.41 VTs in which epicardial breakthrough preceded the earliest endocardial activity were reported,39 and some were successfully terminated by an epicardially directed procedure.38,42,43 Sosa et al44 described subxiphoid access to the pericardial space for mapping and ablation; this approach was used in the study reported here.

Case Study: Methods

The case presented here was referred for ablation of scar-related VT to the Cardiac Electrophysiology Laboratory of the Queen Elizabeth II Health Sciences Centre in Halifax. The patient (age, 64 years) had a history of severe ischemic cardiomyopathy. Previous VT storms were resistant to drug therapy and were treated with VT ablation in 2003 and 2005. Several morphologies of VT from the patient's study in 2005 suggested an epicardial origin of arrhythmia. In 2007, the patient presented with incessant slow VT (with two distinctly different morphologies, denoted VT A and VT B). A cardiac resynchronization therapy defibrillator (CONTAK RENEWAL 3 RF; Boston Scientific, Natick, MA) had been implanted, and the patient was chronically paced. After written informed consent and detailed explanation of the procedure, the patient underwent a successful ablation of both VTs in June 2007.

The detections on the implantable cardioverter defibrillator (ICD) were deactivated prior to the ablation procedure, which was performed under general anesthesia. The groins were infiltrated with 2% lidocaine without epinephrine; 5- and 6-French sheaths were inserted percutaneously into the right femoral vein, and an 8-French sheath was inserted percutaneously into the right femoral artery. To gain access to the pericardial space, a Weiss needle was advanced below the xiphoid process to the pericardium under fluoroscopic guidance. A wire was then inserted into the pericardial space, a tract was progressively dilated with 6- and 8-French dilators, and an 8.5-French convoy sheath was advanced to the pericardial space. A 6-French hexapolar catheter was advanced to the right ventricular apex. A Navistar Thermocool F-curve catheter (Biosense Webster) was introduced into the pericardial space via the convoy sheath and used for epicardial mapping. The same catheter was later introduced into the left ventricle via the retrograde aortic approach and used for endocardial mapping and ablation.

BSPM data were collected during sinus rhythm, ICD pacing, roving catheter pacing, and VT. Multipoint catheter mapping and pacing were performed by the CARTO electroanatomic mapping system (Biosense Webster) using the Navistar Thermocool catheter. Bipolar and unipolar electrograms were recorded both by the CARTO system and by the GE Cardiolab system (GE Healthcare, Piscataway, NJ). The time delay was determined from the 12-lead ECG on the GE Cardiolab system; a stimulus-QRS interval >40 ms was referred to as delay.45

ECGs were processed as described earlier in the "Body Surface Potential Mapping" section. The patient-specific geometry of the torso and epicardial surface was extracted from the CT data, the transfer matrix, A, was calculated, and the inverse solution was obtained by means of a second-order Tikhonov regularization with the L-curve method used to determine the regularization parameter as described earlier in the "Forward and Inverse Problems of ECG" section and in Appendices A and B. MAP3D visualization software46 was used to display potential distributions on epicardial and torso surfaces.

The gold standard for assessing the accuracy of the inverse solution was the information collected by the CARTO electroanatomic mapping system during the ablation procedure. The patient-specific CARTO geometry was registered manually and fused with the CT data, as illustrated in Fig. 12–7. The locations of the CARTO sites where pacing was delivered and the estimated locations of these sites obtained by the inverse solution were compared to determine localization accuracy. CARTO voltage mapping was performed to delineate the scar and scar margin.

Case Study: Results

Reference Maps Delineating Scar Region

Figure 12–8 shows the epicardial voltage substrate maps obtained with the CARTO system. These maps delineate electroanatomic substrate by using color mapping to annotate bipolar signal amplitude (30-400 Hz) that distinguishes between normal myocardium, infarct scar, and border zone around scarred tissue. Regions of low-amplitude signal appear on the anterolateral left ventricle and over the lateral right ventricular wall. Moderate-size patches of confluent scar cover the midsection of the lateral left ventricle.

Localization of Epicardial Pacing Sites by Inverse Solution

Pacing on the epicardial surface initially produces epicardial potential distributions with an area of negative potentials surrounding the pacing site, along with two areas of low-level positive potentials due to conduction anisotropy at the stimulation site, as demonstrated in vivo using dense electrode arrays.47 By using BSPM recordings obtained during epicardial pace mapping (n = 22), we tested the accuracy of the inverse solution in localizing the epicardial stimulation site by comparing the calculated location of the potential minimum with the known location of epicardial pacing determined by the CARTO system. Figure 12–9 shows, for an epicardial pacing site at the basal inferior left ventricle, computed epicardial potential maps at four instants of time, with corresponding input data (recorded body surface potential maps, also shown for both the anterior and posterior torso in Fig. 12–3). The initial minimum on the epicardial surface (of –0.96 mV, near the pacing site) with an area of negative potentials that surround it is discernible after a long delay at 124 ms after pacing stimulus. At peak depolarization, the potential minimum drifts far from the pacing site. During early repolarization, at 275 ms, a low-amplitude maximum (0.347 mV) of epicardial potentials appears again in close proximity to the pacing site. With the insight provided by epicardial maps, measured BSPM distributions (see both Fig. 12–3 and top of Fig. 12–9) seem to be compatible with epicardial-to-endocardial direction of activation initiated on the inferior wall of the left ventricle (as evidenced by the equivalent dipole pointing straight up at 124 ms).

The results for pacing site localization by means of the inverse solution for all (n = 22) epicardial pacing sites of this patient are summarized in Fig. 12–10, which shows, in four views (corresponding to electroanatomic substrate maps in Fig. 12–8), the actual pacing sites, plotted from coordinates provided by the CARTO system, and the estimates of the location of these sites obtained by the inverse solution. The electroanatomic substrate around each pacing site was characterized by the bipolar voltage detected by the CARTO system and by the delay in response to pacing.

We found, not surprisingly, that the electroanatomic substrate around the pacing site affected localization accuracy. For the group of pacing sites in structurally normal myocardium with no delay (n = 7), the median Euclidean distance between the actual and estimated pacing site was just 11 mm, and the potential minimum appeared after a mean interval of 29 ms following the pacing spike. For pacing in scar and scar margin—which implied involvement of the scar substrate in modifying early activation sequence—the discrepancy was larger, but it did not exceed 35 mm.

Localization of ICD Pacing Site by Inverse Solution

Localizing the paced activation initiated from the leads of the ICD device provides another approach to testing the localization accuracy of the inverse solution because the tip of the ICD lead can be accurately determined from the CT images. The BSPM recording was done during pacing from an ICD device implanted in another patient, also recruited under scar-related VT study, who had an ICD pacing electrode in the endocardial right ventricular apex.

The epicardial potential maps obtained for early activation by the inverse solution (Fig. 12–11) provided a clue for estimating the site of pacing. Because the ICD pacing site was on the right ventricular endocardial surface, the activation wave front propagated first across the right ventricular wall before it broke through on the epicardial surface. Thus, an area of positive potentials appeared initially near the pacing site, reflecting propagation of the wave front toward the epicardial surface, and then, after the breakthrough, an area of negative potentials emerged and intensified throughout depolarization and was replaced by a distribution with opposite polarity during repolarization. The Euclidean distance between the actual location of an endocardial pacing electrode and the location of the early minimum on the epicardial surface was 12 mm. This small discrepancy can be attributed either to the involvement of right ventricular conduction system or to the anisotropic intramural propagation, or to both.

Activation Times Obtained by Inverse Solution

To assess the inverse solution's ability to preserve timing information, the activation time was determined by the CARTO system from bipolar electrograms at multiple points on the epicardial surface by the roving catheter, whereas the ICD device paced at the endocardial site near the right ventricular apex. The epicardial surface reconstructed from measurements made by the CARTO system was aligned with the epicardial surface reconstructed from the CT scan, and CARTO points were projected onto the nearest nodes of this patient-specific epicardial surface. Figure 12–12 shows the activation times at all nodes of this surface as estimated by three-dimensional interpolation from data provided by the CARTO system. The activation time obtained via the inverse solution was determined in each calculated epicardial electrogram by the steepest-descent criterion during depolarization; the isochronal maps obtained by this method are shown in Fig. 12–13.

The activation times yielded by the CARTO system show early depolarization of the right ventricular apex (red) and the activation wave spreading over the right ventricular anterior wall in approximately 100 ms (see Fig. 12–12, left); the latest area to activate (purple) is the basal region of the posterior left ventricle (see Fig. 12–12, right). The inverse-solution isochrones in Fig. 12–13 show a qualitatively similar spread of activation; the early activation of the right ventricular apical region, the progression of the activation wave over the right ventricular anterior wall and then the left ventricular posterior wall, and the latest activation of the posterobasal wall are all correctly estimated by the inverse solution. However, there is some discrepancy in the time scale, which reflects inability of the steepest-descent criterion to detect the activation wave front from low-amplitude electrograms in the region with depressed conduction.

Localization of the Reentrant Circuit's Exit Site by Inverse Solution

The earliest region to depolarize during scar-related VT is the site where the reentry circuit exits the scar to reexcite normal myocardium and sustain the arrhythmia. Localizing these exit sites is essential during the ablation procedure to guide delivery of radiofrequency energy that interrupts the reentrant pathway. To evaluate the potential usefulness of inverse ECG in aiding radiofrequency ablation, we applied the methodology tested during paced activation to BSPM data recorded during VT.

The onset of reentrant activation can usually be determined from body surface ECGs (see Fig. 12–2), and the BSPM distribution observed during this phase can be used to estimate the locus of exit site from calculated epicardial potential distributions or activation isochrones.

Two separate morphologies of VT (denoted VT A and VT B) were recorded for this patient, and both tachycardias were successfully ablated. Figure 12–14 shows (in 20-ms increments) a sequence of body surface potential maps for VT A. Initial distributions (during the first 10 ms) have very low amplitudes; the first discernible pattern of body surface potentials with a maximum near precordial site of lead V4 develops at 11 ms, and this distribution remains stationary until nearly 70 ms; subsequently, the maximum migrates toward the right chest. The repolarization minimum near the site of V2 appears at approximately 250 ms and remains stationary until approximately 410 ms.

Epicardial potential distributions obtained by the inverse solution for VT A (Fig. 12–15A) show early negative potentials on the high basal inferior region (at 69 ms) near scar border zone; this area of negative potentials remains stationary during the depolarization phase (106 ms). A map of the isochrones of activation (Fig. 12–15B) shows earliest activation starting at the high inferobasal left ventricle in the region where the earliest potential minimum appeared. The isochronal map shows activation spreading in a counterclockwise direction around the area of block on the inferior epicardial surface, starting at the inferior left ventricular base and ending on the basal inferior and inferolateral right ventricle while the superior view of the isochronal map shows activation spreading globally from apex to base. The ablation site was endocardial; its projection on the epicardial surface (marked by asterisk) was located near the early potential minimum and the region of early activation.

Figure 12–16 shows (in 20-ms increments) a sequence of body surface potential maps for pacing from the endocardial site that can help to locate the site of successful ablation of VT A.

The aim of this chapter was to illustrate potential applications of ECG imaging as a noninvasive investigative tool for clinical cardiac electrophysiology applications. ECG imaging can provide detailed spatial and temporal information on the electrical activity of the heart from body surface measurements, and thus, it can become an important aid in investigating mechanisms of cardiac arrhythmias and in making their treatment more effective.

We first demonstrated the localization accuracy of ECG imaging by comparing calculated epicardial electrograms, and measures derived from them, with a gold standard provided by CARTO electroanatomic mapping. Next, we used these tested methods in recording of VT itself to investigate their usefulness in facilitating catheter ablation of this arrhythmia.

Comparing the estimated location of the early epicardial potential minimum with the known site of epicardial pacing has been used previously as a measure of the inverse-solution accuracy in experimental studies using canine heart preparation in a torso tank, as reported by Oster et al,5 and more recently in studies involving patients with ICDs using a right ventricular apical pacing lead, as reported by Ramanathan et al.8 The median of localization accuracy in the present study, for epicardial pacing in the structurally normal myocardium with stimulus-to-QRS delay of less than 40 ms, was 11 mm, which compares well with the localization distance of within 11 mm reported by Ramanathan et al8 and the 6- to 29-mm distance reported by Ghanem et al48 in their clinical study performed during open-chest surgery.

In the pace mapping recordings with stimulus-to-QRS delay of greater than 40 ms and/or those for stimuli delivered inside scar/scar margin, the effect of pacing from multiple sites around the epicardium caused clustering of early potential minima around one or two regions of the ventricular surface, in particular, in the patient presented here, in the high basal inferolateral region (five sites). A possible explanation of this behavior is that in cases where pacing is delivered inside scar, the activation wave front travels with depressed velocity, taking longer to emerge out of the scar, at which point a substantial bulk of myocardial tissue activates, generating the multiple potential minima. Thus, the location of such clustering of potential minima may indicate a part of the tachycardia circuit close to the exit site that should be investigated further by the catheter mapping. In such cases, the early potential minima may be far from the pacing site. This argument is further supported by the fact that the location of the early potential minimum during recordings of VT of the same patient lies very close to the region where the early potential minima induced by pace mapping tend to cluster; these regions are also near the location where radiofrequency ablation was performed.

The display of isochrones of activation is an effective tool for visualizing propagation of the activation wave front in the ventricular myocardium. Producing the global isochronal map over the entire epicardial/endocardial surface by direct measurement requires maneuvering of the catheter tip to collect electrograms point by point and determining a time of activation for each point. This elaborate process is used in clinical electrophysiology to map sustained arrhythmias; however, the same approach cannot be used when arrhythmia is not sustained or when it is not hemodynamically tolerated. Therefore, one of our aims was to assess accuracy of activation times derived from the epicardial electrograms obtained by the inverse solution in comparison to those derived from directly measured epicardial electrograms collected invasively through the CARTO system. We used peak negative slope on calculated unipolar epicardial electrograms to determine the local activation time. However, in diseased hearts, there are many situations in which the peak negative derivative may not be a reliable predictor of activation time.49 Thus, the presence of diseased tissue complicated interpretation of electrograms obtained by the inverse solution in the case presented here. The method of steepest downslope applied to electrograms obtained by the inverse solution relies on derivatives calculated from successive samples, which makes it vulnerable to measurement noise in body surface potentials or to changes in the regularization parameter that controls the amount of smoothing in the inverse calculation. In cases of electrograms with low peak amplitude, there were multiple points with very similar steepest slope values, which caused spurious assignment of activation times. This problem was partially mitigated by spatial smoothing of inverse electrograms.

During endocardial pacing at the right ventricular apex, delivered by ICD, regions of early and late activation determined from calculated epicardial electrograms agreed well with those determined from electrograms measured directly by the CARTO system. However, the isochrones determined from calculated electrograms for paced activation (see Fig. 12–13) and for VT (see Fig. 12–15) show irregularities that include compression of the total activation time range, near-instant (bulk) activation of large regions, and crowding of isochronal lines at multiple locations. These flaws are due to limitations of the inverse-solution method used here, which can only produce unipolar electrograms that inherently contain superposition of both near- and distant-field potentials. On the other hand, the reference activation times (yielded by the CARTO system) were obtained from directly measured bipolar electrograms, which reflect only intrinsic near-field activity50; they exhibit a sharp peak when the activation wave front passes by the catheter tip.

In inverse potential maps calculated from body surface potentials recorded during VT, an early potential minimum was found on the epicardial surface near the endocardial site where radiofrequency ablation successfully terminated the arrhythmia. These observations support the argument that the early potential minimum corresponds to the location where the tachycardia reentrant pathway exits dense scar to reexcite normal myocardium to continue the arrhythmia cycle. However, because the VT reported here was not hemodynamically tolerated, it was not possible to exactly map the tachycardia's circuit using either activation mapping or prolonged entrainment mapping. Thus, the location of arrhythmia exit site was inferred from the site where radiofrequency ablation was successfully performed, aided during the procedure by goodness of match between the standard 12-lead ECG recorded during pace mapping and that recorded during VT. Comparison of Figs. 12–14 and 12–16 demonstrates how matching sequences captured by BSPM can be potentially used in real time to guide the ablation procedure.

A number of investigators have pointed out that geometric inaccuracies of the representation of the heart surface relative to the body surface contribute substantial errors in the estimates of epicardial potentials.51-53 Recently, Cheng et al54 showed, in an experimental study, that there are minor effects on the reconstructed epicardial potential, as measured by relative error, when body surface potentials were corrupted by measurement noise or correlated electrode displacement. However, they found that heart surface scaling, translation in three orthogonal directions, and rotation around coronal or sagittal planes resulted in significantly poorer epicardial solutions. They also noted that the L-curve method was generally more robust and stable in the presence of higher levels of geometric errors compared with other methods of selecting the regularization parameter.

The ultimate goal of studies conducted using inverse ECG is to provide adjunct electrical information to the electrophysiologist about local electrical events in the heart that aids diagnosis and decision making in treating the arrhythmia. Two important attributes for deployment of such an inverse ECG imaging system in the clinical setting are computational time and any human interaction required to obtain the inverse solution. The process of calculating the inverse epicardial solution involves (1) anatomic body imaging and processing of images to obtain transfer matrix; (2) BSPM and processing of recorded signals; (3) inverse-solution procedure; and (4) visualization. Currently, MRI and CT are routinely used to obtain anatomic images, and two advanced mapping systems (CARTO electroanatomic mapping and EnSite NavX mapping [Endocardial Solutions, St. Paul, MN]) allow importing of these images and segmentation of heart chamber of interest to be incorporated within the same display of data collected by these systems. The discretization of ventricular geometries and electrode locations from segmented images online is possible with current mesh generation algorithms.

Overall, it is well within the capability of current technology to perform the functions required for epicardial inverse ECG imaging during the ablation procedure. What remains to be achieved is rigorous validation of the inverse solution using physiologic data collected under clinical conditions against reference data provided by established clinical imaging tools.

The generalized singular value decomposition (GSVD) of Equation 12–2 is defined as follows.4 Let A be any m × n real matrix m > n, let B be a p × n(p ≤ n) real matrix of rank q ≤ p, and let k + q be the numerical effective rank of (

Here C = diag(α(k+1),..., αn and S = diag(β(k+1),..., βn are diagonal matrices, O denotes the additive identity, Um×m and Vp×p are orthogonal and Yn×n is a nonsingular matrix. The generalized singular values μi = αi/βi satisfy 0≤αi≤1, 1≥βi>0, and α2i+β2i = 1 for i = k + 1, …, n.

The unique solution to Equation 12–2 for each λ > 0 is given by

where ui are the columns of U, Yi are the columns of Y, and b stands for ΦB.

The generalized normal equations associated with the perturbed least-square problem in Equation 12–2 can be stated as4

For each λ > 0 the residual ‖Ax(λ)–b‖ can be expressed in the form

and for each λ > 0, the seminorm ‖BΦH(λ)‖ can be expressed in the form

Thus, the residual ‖AΦH(λ)–b‖ and solution seminorm ‖BΦH(λ)‖ can be expressed simply in terms of λ, the singular values μi,i=k+1,...,n, and the scalar products ui·b, i=k+1,...,m.

The Tikhonov regularization scheme represented by Equation 12–2 arrives at the solution by minimizing two quantities; the first is the least squares solution of Equation 12–1, and the second is a penalty function that imposes conditions on the smoothness of the solution by controlling the weighting regularization parameter λ and the regularization matrix B. Typically, the regularization matrix was chosen to be the identity matrix (I), which is called zero-order Tikhonov regularization, or a discrete approximation of the surface Laplacian matrix (L), which is referred to as second-order Tikhonov regularization. An important step in Tikhonov regularization is the choice of the regularization parameter (λ). As λ goes to zero, the solution to Equation 12–A2 tends to the least squares represented by the first term, which produces unstable epicardial potentials (underregularized solution). However, if λ is too large, the solution becomes overly smooth because the second term of Equation 12–A2 dominates the solution (overregularized solution). The optimal value of the regularization parameter provides a balance between instability and smoothness.

Several methods have been suggested in the literature for choosing the regularization parameter λ. These methods include the generalized cross-validation method,55 the composite residual and smoothing operator (CRESO) criterion,56 the L-curve method,57 and the zero-crossing method.58 We used the L-curve method of Hansen,57 which determines the λ parameter from the plot of the seminorm of the solution ‖BΦH‖ versus the norm of the residual ‖AΦH–ΦB‖. Choosing the regularization parameter in the vertical segment of the L-curve puts more weight on minimizing the residual, thus underregularizing the solution and making it unstable, whereas on the horizontal part of the L-curve, more weight is put on minimizing the seminorm of the solution, thus overregularizing it and smoothing out fine details. The corner of the L-curve provides a good balance between the two quantities in Equation 12–2. The corner can be defined mathematically as the point of maximum curvature and can be computed as the maximum of the function κ(λ):

where x(λ)=‖AΦH(λ)–ΦB‖, y(λ)=‖BΦH(λ)‖, with prime and double prime denoting first and second derivatives with respect to λ, respectively. Functions x(λ) and y(λ) can be expressed directly in terms of the regularization parameter (λ), the singular values (μ), and vector quantities (u and b) by means of Equations 12–A4 and 12–A5.

Various regularization methods have been implemented in this laboratory, making use of LAPACK and BLAS linear algebra packages.59 In this study, we used the implementation in MATLAB, making use of the Hansen's regularization toolbox.60