Europace (2014) 16, iv120–iv128 doi:10.1093/europace/euu268
SUPPLEMENT: REVIEW
A comparison of electrocardiographic imaging based on two source types Adriaan van Oosterom* Radboud University Nijmegen, Nijmegen, The Netherlands Received 30 July 2014; accepted after revision 2 August 2014
----------------------------------------------------------------------------------------------------------------------------------------------------------Keywords
Electrocardiographic imaging † Inverse problem † Electric source descriptions
Introduction The dynamics of the contraction of the cardiac muscle is initiated and accompanied by the flow of electric current. These set up potential fields throughout the myocardium as well as in the body tissues surrounding the heart and on the surface bounding the current flow: the thorax. The observed time course of the potential difference between any two locations inside the thorax or on its surface is referred to as an electrocardiogram (ECG). In cardiology, the set of 12 simultaneously recorded signals, derived from nine strategically chosen locations on the body surface, forms one of the most frequently used diagnostic tools. As such, the interpretation of the ECG can be classed as a non-invasive method. Moreover, it can be classed as performing an inverse procedure, being the interpretation of a state, or phenomenon on the basis of observations made at a distance. Over the period of more than a century following the first recording of an ECG, the diagnostic accuracy of the ECG has steadily increased. Yet, in some types of cardiac disease the accuracy still needs to be improved. One of the later developments has been the increase in the number of electrodes sampling the time course of the potential.
The display of a sequence of instantaneous potential fields on the thorax by means of scalar maps, body surface potential maps (BSPMs),3 has inspired the search for methods which also display the interpretation of the signals close to the heart, on its bounding surface (epicardium and endocardium) or intramurally. These emphasize the spatial character of the cardiac electric sources. This procedure has recently been known as imaging. The image formation involved in biomedical imaging rests on the nature of the physics involved. In many of these, the image formation is based on the interaction of the applied, or induced, radiation encountered in the tissue while passing along it in straight lines. In contrast, the electric potential field in the bounded region of the thorax is formed along electric current density lines that form loops, starting at exit points within the cardiac tissue (the sources) and returning to nearby locations (the sinks). The charge carried along the current density lines is merely pumped around. This property necessitates the formulation of an electric source model (SM) as well as of the properties (electric conductivities) of the electrically passive tissues external to the heart: volume conductor model (VCM).
* Corresponding author. Tel: +31 485 316543. E-mail address: [email protected] Published on behalf of the European Society of Cardiology. All rights reserved. & The Author 2014. For permissions please email: [email protected].
Downloaded from by guest on February 17, 2015
The aim of the study to compare the performance of two major source types involved in the imaging of the electric activity of the heart on the basis of potential differences observed on the thorax. The images depict either the timing of activation and repolarization of the myocardium or the potential field on a surface closely encompassing the myocardium. The depolarization and repolarization timing on a closed surface bounding the ventricular myocardium was derived from measured body surface potentials (BSPs), an MRI-based electric volume conductor model comprising the geometry of thorax, lungs, heart surface Sh, and cavities. The solution was constrained by using a template of the local transmembrane potentials (TMPs). The latter serve as the strength of the Equivalent Double Layer (EDL) source model (SM), which was used to compute the potential field on Sh (epicardium and endocardium). The second SM is the potential distribution on the epicardium Sp, referred to here as the Equivalent Potential Distribution (EPD). Its strength was estimated directly from the BSPs. The inflection points of the estimated electrograms (ELGs) were taken as markers of the timing of local depolarization and repolarization. The endocardial potential fields estimated using both sources exhibited qualitative similarity, as did the ELGs. With reference to the one generated by the EDL source, the magnitude of the estimated endocardial EPD field was smaller, the downslopes of the ELGs were lower. The timing of depolarization estimated from the EDL-based ELGs was highly correlated with those of the TMP templates, the EPD-based correlation was lower. For the repolarization timing the corresponding test values indicated an insufficient similarity. The EDL- and EPD-based source variants deserve to be studied alongside each other in the future development of electrocardiographic imaging.
iv121
A comparison of electrocardiographic imaging
What’s new? † An analysis of the similarity between potential fields, electrograms (ELGs), and timing as estimated based on either the EDL or the EPD source model (SM). † An explicit demonstration of the potential field on the closed surface closing the myocardium as generated by the EDL SM. † An explicit demonstration of the ELGs at the closed surface closing the myocardium (endocardium and epicardium) as generated by the EDL SM. † A reconfirmation1 of the value of including the inhomogeneity of the lungs in EPD-based imaging. † Fresh food for doubt on the accuracy of estimating the timing of the repolarization at Sh, taken to be the timing of the inflection point of the transmembrane potential downslope, from the timings of the inflection points of the upstrokes of the nearby ELGs.2
F = A S,
(1)
f1 (r ′ , t) =
1 4ps
∇†(si (r )∇Vm (r , t)) dv(r ). R V
(2)
The integral is taken over the volume encompassing all myocytes, ∇ denotes the spatial gradient operator, ∇† the divergence operator, si the intracellular conductivity, and R the distance between source location r and external observation point r ′ . The term ∇†(si (r )∇ Vm (r, t)) represents the electric current source density within the myocardium, the strength of which is proportional to the varying source strength of the local TMP: Vm (r , t). By using the approximation si (r ) = si and applying Green’s second theorem (a standard result from differential calculus)11 Eq. 2 is equivalent to F1 (r ′ , t) =
−s i Vm (r , t) dv(r ′ ,r ), 4ps
(3)
S
in which F denotes a matrix (size: L,T) comprising the documentation of the potential differences between all observation points (field points, ℓ ¼ 1 . . . L) and the potential at a common reference (CR) at all instances (t ¼ 1 . . . T ) studied. Moreover, a matrix S (size: N,T ) expresses the source strength at all locations (n ¼ 1 . . . N ) of the source elements as a function of time. With the source strengths depicted as images, the elements of S may be viewed as the N pixels of the sequence of all T images, stills of a movie. Finally, matrix A (size: L,N ) represents the links (transfer) between all source elements and the potential at the observation points, i.e. the electrode locations. The computation of the source matrix S on the basis of Eq. 1, a set of observations F and a computed transfer A, is generally referred to as solving the inverse problem of electrocardiography. 5 The specific solution to inverse problems in which the solution S can be represented by images is referred to as electrocardiographic imaging (ECGI). Several other acronyms related to variations in the involved source description have appeared in the literature. In this paper two of the main source descriptions are discussed and their application to the data pertaining to one and the same subject is illustrated. Next, their similarities, differences, advantages, and disadvantages are discussed.
A tale of two source types The equivalent double layer Soon after the introduction of the electric current dipole,6 Wilson et al. 7 introduced an equivalent SM that is closely linked to bioelectricity: the uniform double layer (UDL). In the UDL the source is
with dv(r ′ ,r ) denoting the solid angle subtended by small segments of Sh, the surface bounding V. This indicates that within the limitation of the assumed homogeneity of the intracellular domain, all source activity within the source volume V (Eq. 2) may be specified wholly by the distribution of Vm(t) on the surface Sh of the source region. It forms the generalization of the classic EDL theory. The result was formulated by Geselowitz.8 In a later publication he showed that it also holds true for anisotropic tissues, on condition the anisotropy ratio is uniform.9 Equation 3 forms the basis for solving the inverse problem aimed at estimating the timing of depolarization and repolarization, but may also be used to compute the epicardial potential field as is discussed in the section ‘An evaluation’ (Figure 1). The numerical expression for the forward computation is F = BF1 = BV Vm = DVm ,
(4)
in which Vm is a matrix of size (N,T ), with N being the number of nodes of the mesh representing S. Matrix V of size (M,N ), expresses the numerical variants of dv(r ′ ,r ), scaled by the remaining factors in Eq. 2, with M being the number of field points external to Sh for which f1 (r ′ , t) is to be known. Matrix B, of size (L,M), accounts for the effects on the potential of any inhomogeneity in the electrical conductivity of the external domain. These M locations include those at the body surface Sb. If the medium between the thorax surface Sb and Sh is homogenous then M ¼ L. The elements of B can be computed, e.g. by means of the boundary element method (Section. 2.6.4.12).
Downloaded from by guest on February 17, 2015
Upon having specified both models, the electric potential field at any set of electrode locations can be computed. This computation is known as solving the ‘forward problem of electrocardiography’.4 A generic expression for the outcome of this computation is the matrix multiplication
restricted to the boundary of a propagating activation front, with its strength taken to be uniform over the wave front. The equivalent double layer (EDL) is the successor to the UDL; its elementary sources are located on the closed surface Sh bounding the myocardium (either ventricles or atria). It extends the effectiveness of the UDL concept to the repolarization phase.8 – 10 The EDL concept is based on the general expression linking transmembrane potential (TMP) Vm (r , t) of the active myocytes to the field potential f1 (r, t) in an infinite homeogeneous external medium, a generalization of the classic cable theory. It reads:
iv122
A. van Oosterom
with V being the matrix of elementary solid angles and G the matrix of the weights dS(r )/R(r ′ ,r ) involved in the integration. In both types of matrices, the respective scaling factors seen in Eq. 5 have been absorbed. The use of Eq. 5 requires the potential field at the boundary and its normal gradients to be known. Equation 5 may also be used to compute the link between the fields on the boundaries S1 and S2. As an example, the most commonly considered situation is where partial, but nonetheless useful information is available. Since the normal gradient of the potential at the isolated torso surface is zero, matrix J1 is zero. The field F2 acts as a potential source distribution.14 The resulting potential field at S1 is computed as follows. Equation 6 may be applied to the potential fields just inside the two bounding surfaces. By combining the resulting two equations and, after rearranging terms, it follows that F1 can be found by solving the system of linear equations
Geometry of torso, heart and body tissue conductivity
Measured BSPs F
Source models S1; Vm on Sh
S2; Y on epicard
Vm=Vm(d,r) TMP template d,r Vm=Vm(d,r)
Constraints priors
Y: norm; smoothness; temporal
Solutions Forward computation
Y: images electrograms
Timing inflection points
Figure 1 Flowchart of the procedures involved in the imaging based on the two different source descriptions discussed. On the left: pertaining to the EDL source, on the right: to the EPD source. Note the clockwise loop shown at the bottom.
1 1 1 Jn (r ) dS(r ) F1 (r )dv + F2 (r )dv + 4p 4p 4ps R 1 Jn (r ) dS(r ), + 4ps R
S2
S1
(5)
S2
∂F (unit: A/m2) (a scalar function) expresses the ∂n magnitude of the local gradient of the potential field normal to the surfaces scaled by the local electric conductivity. The potential field Fi, Eq. 6, at any field point i in V can be computed numerically from the expression where Jn (r) = s
Fi = Vi,1 F1 + Vi,2 F2 + Gi,1 J1 + Gi,2 J2 ,
− G1,2 G2,2
F1 J2
=
V1,2 I − V2,2
F2 ,
(7)
(6)
A F2 = F1 ,
(8)
which is the essential expression for computing the pericardial potential distribution F2 on S2 as estimated on the basis of observations F1 on S1. The derivation can be generalized for cases in which V comprises inhomogeneous, closed, non-intersecting source-free regions such as poorly conducting lungs.17
Numerical aspects The integrals as appear in Eqs. 3 and 5 need to be computed numerically. The functions involved (potentials, source densities) are integrated over small surface elements. The most convenient shape of such elements is the planar triangle. Closed form analytical expressions are available for the integral over triangular elements in which the function over the triangle is approximated to be bilinear, specified by the function values at the vertices. For the integration involving solid angles, the analytical expression can be found in Munck,18 the one for the distributed monolayer in a more recent paper.19 For the results involving the solid angles, the potential field is discontinuous, i.e. when crossing the surface element the solid angle changes by 4p, scaled by a factor comprising the two conductivities on either side of the surface. Hence, the value assigned to the integral at or close to the surface needs to be tuned to the application in hand. In contrast, the potential field on the triangle generated by a monolayer on the triangle itself is continuous whereas the component of its gradient along the normal of the surface is discontinuous.
Summarizing The basic, numerical expressions linking SM strength to observed field potentials as discussed in this section are Eqs. 4 and 8. The matrix formulations used, incorporate the temporal nature of source strengths Vm (EDL) and F2 (equivalent potential distribution,
Downloaded from by guest on February 17, 2015
The second source variant is the potential distribution on a closed surface encompassing Sh generated by its internal bioelectric sources, such as the pericardium. Its strength is to be estimated on the basis of potentials observed on the thorax.13 – 15 In experimental studies the surface involved is frequently taken to be a closed variant of the epicardium. The thorax and pericardium form two nested non-intersecting surfaces. No active bioelectric sources are assumed to be present in the domain bounded by these two surfaces. Below, the essence is shown of the computation of the forward transfer. It was first derived by Barr et al. 14 It is illustrated by considering a source-free domain V with uniform isotropic conductivity s bounded by closed surfaces S1 and S2, with S2 contained in S1. The derivation is shown here in some detail, motivated by the high relevance of the transfer. As shown in page 12 of the reference,16 the potential field within the source-free domain depends uniquely on the potential and the normal gradient of the potential at its boundaries. The general expression reads
S1
I − V1,1 V2,1
with I being the identity matrix having the size of the accompanying matrix. From the solution to this expression a forward transfer matrix A is isolated, yielding the linear system
The equivalent potential distribution
F(r ′ ) =
iv123
A comparison of electrocardiographic imaging
EPD), respectively. Matrices D and A, which express the volume conduction effects only, are taken to be static. The isolation of matrices D and A in the forward procedure increases the transparency of the handling of the inverse procedure.20 In Eq. 4, matrix D is seen as transferring the current double layer source density, whose strength is proportional to Vm at the heart surface Sh, to the potential field inside a bounded medium, viz. the thorax, as well as on its boundary. In Eq. 8, the transfer A expresses the effect of the potential distribution on the pericardium.
parameters p are estimated from Eq. 9 by treating it as a least-squares problem. The source strength Vm(p) depends non-linearly on the timing parameters p (mainly stretching and/or shifting the TMP shape). As a consequence, finding p from Eq. 9 poses a so-called non-linear least-squares (NLLS) problem. Such problems need to be solved iteratively, starting from a posed initial solution p0. In the example shown in the next section, the solution method used is the Marquardt algorithm.25 Details of this procedure may be found in references therein.32,33
The imaging problem
Including constraints
The equivalent double layer A paper by Salu26 provided the first introduction to the EDL positioned on the surface Sh bounding the entire myocardium, with unit strength on depolarized sections of Sh and zero strength on the remaining parts. The link between UDL and EDL is based on the properties of solid angles.27 The EDL as the source implied in an inverse procedure was introduced some years later.28 In the years thereafter, various studies on possible applications and improvements followed.29,30 The current state of the development is based on the forward formulation in Eq. 4, in which a template of the TMP is included. The shape of the TMP source strength at each node n (n ¼ 1 . . . N ) of the N nodes (vertices) of a triangulated representation of Sh is specified by an analytical function. Its parameters include the timing of depolarization (dn) and repolarization (rn). These may be denoted by the column vectors d and r. The nature of the analytical function is such that these two parameters relate to the inflection points of the upslope of the local TMP and of their final downslope, respectively.31 A full description of the procedure is included in van Oosterom.32 The estimation is that of finding the parameter set that minimizes the least-squares difference between the observed signals and those resulting from the model ‘prediction’. We first reformulate Eq. 4 as DVm = DVm (p) = F,
(9)
with p denoting the timing parameters (either d or r) setting the shape of the TMP, and F the CR signals observed on the thorax. The latter may be a subset of the potentials at all nodes involved in the numerical representation of the thorax boundary. The timing
Results of numerical experiments based on the estimation of p from Eq. 9 indicated that additional constraints are needed if realistic, physiologically sound solutions are to be found. These constraints may be combined in Eq. 9. Different constraints were tested; the one adopted29 is the Frobenius norm of the surface Laplacian,11 i.e. ||lL p||F with matrix L denoting the surface Laplacian operator34 and l a scalar variable setting the weight of this constraint. The norm involved is the square root of the sum of the squares of all matrix elements, the so-called Frobenius norm. The value of l may be set to drive the estimated p toward a spatial smoothness having a physiologically realistic value.29 The iterative procedure for solving an NLLS problem may end up at a local minimum of the parameter space rather than at the desired global minimum. The quality, i.e. realism, of the solution found depends greatly on the realism of the initial estimate p0. The initial estimate used in recent work33,35,36 is based on the fastest route algorithm based on Huygens’ principle applied to propagating activation wave fronts.
The equivalent potential transfer Equation 8 specifies the forward transfer A between two potential fields on two nested surfaces: the thorax surface Sb and the pericardial surface Sp. The inverse of matrix A may be used to estimate the pericardial potential C from F. In contrast to the inverse procedure for the EDL, the basic relationship between SM C and observations F is linear. It may seem as if this type of inverse merely requires computing the inverse of the forward transfer matrix A. However, the existence of a unique inverse solution C depends on the properties of A. This has a dominant role in these types of inverse computations, hence the following, brief recapitulation is in order.
Regularization Just as when dealing with the EDL source, in studies estimating the EPD it was found to be necessary to include additional constraints. This was necessary even if the number of leads L used in sampling the potential field on the body surface was greater than the number N of the nodes on Sp. Constraining the solutions is carried out by demanding the solution to approximate CC R, with C and R specifying the nature of the constraint and its value aimed for, respectively. If the constraint can be expressed by the matrix formulation CC R, the estimate C is found as the linear least squares (LLS) solution to
A
l C
C
F l R
(10)
Downloaded from by guest on February 17, 2015
In this section, some basic methods involved in solving the inverse problem are discussed with regard to their application to the imaging of either of the two source types introduced in the previous section. The basic, numerical expressions linking SM strength to observed field potentials are Eqs. 4 and 8. In the sequel, the thorax potentials are denoted by F, and the potential distribution on the pericardium by C. The inverse problem is treated from the perspective of source parameter estimation. The parameters specifying the involved VCM, such as geometry and passive tissue properties, are taken to be known. Fundamental methods for parameter estimation in a wide variety of scientific domains have been developed over the last 50 years. Major reference books that widely cover the theory are provided under references.21 – 25
iv124
A. van Oosterom
The constant l is the Lagrange multiplier as used in the classic, constrained optimization theory.37 The higher the values of l, the more the estimate is forced to comply with the constraint. The procedure is also referred to as regularizing the solution.38 The standard expression for the LLS solution to Eq. 10 is C(l) = (AT A + l2 CT C)−g (AT F + l2 CT R).
(11)
For each value of l, this produces a unique, minimum norm solution, which poses the question of which value l should be taken. For a recent discussion on this problem see Milanicˇ et al.39
Comparing both imaging variants
Below, an evaluation is presented of the theory exposed in Figure 1. It is based on a reference data set similar to the one included in the interactive forward simulation package41,42 available from www. ecgsim.org. These EDL derived data are taken as a gold standard against which the performances of the EPD estimates are compared. The specifications of the geometry of the relevant VCM surfaces, as well as the reference potential field F on the thorax used in this comparison, were identical. Material The reference data comprises the numerical representation of an inhomogeneous, multi-compartmental VCM. It includes MRI-derived geometry of the body surface, the lungs, the left cavities (combined ventricular and atrial), the right cavities (combined ventricular and atrial), and the surface Sh closing the ventricular myocardium. A three-fold higher electric conductivity than the surrounding space was assigned to the blood inside the cavities, and a five-fold lower value to the interiors of the surfaces representing the lungs. The timings (dn and rn) of the TMP templates at Nh ¼ 1500 nodes on Sh was derived from a dedicated inverse procedure The EDL SM was used to simulate body surface potentials at all L ¼ 300 nodes of the triangulated thorax geometry St. These formed the reference thorax potentials used in the sequel, referred to as F.35 The rows of this matrix F comprise T ¼ 401 subsequent samples at 1 ms intervals, starting at the onset of the QRS complex. A subset of these comprised the L ¼ 65 electrode locations specified by the Nijmegen-lead system.43 In addition, by including the nodes of Sh as field points in matrix B (Eq. 4), the potential field Ch on Sh was computed and taken as a reference field Cref. Methods The EPD estimates were derived as follows. A closed triangulated surface Sp was created comprising 807 nodes of the epicardium of Sh as well as 1090 Nc nodes specifying the surfaces of both atrial cavities. These in all Np ¼ 1897 nodes formed the vertices of Sp, which represent the pericardium. Its interior comprises all active (EDL) sources as well as the major inhomogeneity of the VCM: the blood inside the cavities. The potential field at the Np pericardial nodes was estimated by treating them as field points in the volume conductor specified in the previous section. The resulting field serves as the strength for the EPD source on Sp, referred to as Cest. Its size is Np,T. The domain bounded by Sp and Sb is source-free. In order to test the EPD-based imaging, two variants of the required transfer matrix A (Eq. 8) were computed. In one variant the effect of the inhomogeneity of the lungs is accounted for. In the other the conductivity of the domain is uniform. The estimation of the pericardial potential field was carried out on the basis of Eq. 11. The results based on the two types of regularization are documented: Tikhonov’ zero- and the second order. For the zero-order variant the constraint matrix C is a diagonal matrix of size N, for the second order it is the surface Laplacian operator L as specified previously. Its size is N,N. In both cases the term matrix R in Eq. 11 is zero and drops out of the expression. In applications to clinical data, where no gold standard is available, the regularization weight l needs to be found. A frequently used method is to scan a suitably selected range of l values and plot the norm of the regularized solution, ||CCest(l )||F, as a function of the relative residual
Downloaded from by guest on February 17, 2015
An overview of the two imaging recipes is displayed in the flow diagram shown in Figure 1. It summarizes and recapitulates the essence of the various topics discussed above. The ingredients indicated at the top of the figure are common in both procedures. The myocardium may be either that of the ventricles or that of the atria. The track on the left pertains to the EDL source. The parameters estimated are the timing of activation (dn) and recovery (rn) at Nh nodes on the closed surface Sh bounding the myocardium. The track on the right pertains to the EPD source. Here the parameters to be estimated are the potentials C(n,t) at the Np nodes on the pericardium at T time instances, relative to those present at the position of a CR electrode. The results may be displayed on the heart surface, colour-coded and/or supported by isofunction lines, for the EDL source: isochrones, for the EPD source: isopotential lines. As discussed, both modalities require the inclusion of some kind of prior knowledge that may be used for constraining the solution. For the EDL-based imaging method an important intrinsic constraint is highlighted: the template of a TMP wave form. An initial estimate is used based on the fastest route algorithm.40 For the EPD source various constraints have been tested, such as the norm of the solution, its spatial smoothness, or the inclusion of a temporal prior.5 The estimated timing parameters, the elements of d (depolarization) and r (repolarization) specify completely the time course of the TMP template and thus of the current source EDL strength at Sh. This source description may in turn be used to compute Ch on Sh, of which the epicardium is a part. In the lower part of Figure 1, this procedure is indicated by the arrows entering and leaving the box ‘forward computation’. At the lower right box the two imaging modalities meet. The results of this forward computation may be expected to be similar to the field Ce on the epicardial surface estimated from the EPD source description. The timing of depolarization and repolarization can be estimated from the inflection points of the electrograms (ELGs): Ce(n,t), the rows of matrix C. These may be compared with the corresponding elements of d and r as used in the EDL-based inverse procedure. This comparison is symbolized by the left pointing arrow at the bottom of the flow diagram (Figure 1). Note that this is a part of a clockwise loop.
An evaluation
iv125
A comparison of electrocardiographic imaging
Table 1 Descriptive statistics of the epicardial potential field Cest computed from F VCM
Order TIkh
L
RD (F)
0 0
300 300
1.85 e– 3 0.54 e– 3
STD (Cest) (mV)
P2P(Cest) (mV)
2.023
28.39
1.35 1.49
26.2 26.3
RE (Cest)
CC (Cest, Cref )
0.7076 0.6478
0.7071 0.7618
............................................................................................................................................................................... Reference lungs hom lungs hom
2
300
2.41 e– 3
1.44
25.0
0.6983
0.7152
lungs hom
2 0
300 65
0.53 e– 3 0.67 e– 3
1.60 1.12
25.7 19.4
0.6323 0.8174
0.7745 0.5760
lungs
0
65
0.44 e– 3
1.17
23.1
0.7982
0.6021
hom lungs
2 2
65 65
0.74 e– 4 0.47 e– 3
1.24 1.32
16.4 17.2
0.8070 0.7659
0.5912 0.6424
Upper row: some statistics of Cref, the epicardial field based on the EDL source. Next rows: estimation results based on the respective optimal l values. Columns: 1: VCM type; 2: order of applied Tikhonov regularizarion; 3: number of ECG signals; 4: relative residual difference between thorax potentials Fref and Fest ¼ A Cest; 5: STD(Cest); 6: extreme peak-to-peak values of Cest (downslope magnitude of the electrograms); 7: relative error between all elements of Cest and Cref; 8: linear correlation coefficient (CC) between all elements of Cest and Cref. Note the discrepancy between the magnitude of values of the relative residual differences: Column 4, and the relative errors: Column 7.
Figure 2 Reference timing of depolarization at the 1500 nodes of the ventricular surface; epicardium and endocardium. Left panel: natural anterior view, white circle marks the site of the epicardial breakthrough; right panel posterior view; white circle marks the locus of septal initiation of the activation. Colour code of the timing (ms). Supporting isochrones drawn at 10 ms intervals.
difference, RD ¼ ||F 2 ACest(l)||F/||F||F.44 Since the situation treated here does involve a pure gold standard, Cref, the scan can be carried out by using the error norm ||Cref 2 Cest(l)||F in a search to find its minimum. In Table 1 the resulting minimum relative errors were listed:, i.e. RE ¼ ||Cref 2 Cest||F/||Cref||F.
Results Figure 2 is an illustration of the timing of depolarization as included in the specification of the EDL source. The left panel depicts, in the natural anterior view, a colour-coded representation of the timing of depolarization on Sh, supported by isochrones lines drawn at 10 ms intervals. The pattern of the timing of repolarization, not shown here, has a similar complexity. The right panel depicts a posterior view, rotated additionally to offer a view of the left endocardium. Snapshots of the corresponding potential distributions on Sh at 30 ms from onset depolarization are presented in Figure 3. Note that this EDL-based procedure yields images both of the epicardium and the endocardium, of the timing as well as of the potential fields. In Table 1. the major descriptive statistics are presented of the test results for eight configurations, computed for the 807 epicardial
15 10 5 0 –5 –10 –15
Figure 3 Reference potential field Cref at the 1500 nodes of the ventricular surface at 30 ms from onset QRS. Left panel: natural anterior view, right panel posterior view, providing a view on the endocardium and the septum. Color code of the potential in mV. Supporting isopotential lines at 2 mV intervals.
reference data. The entries listed pertain to the respective l values yielding the smallest RE value. In addition to the RE values (Column 7), the standard deviation of all elements are documented, as well as the maximum over the N signals of the peak-to-peak value of each signal over time, the magnitude of the downslopes of the ELGs during depolarization and their upslopes during the final part of the repolarization. It shows that the smallest RE value is found when using an inhomogeneous VCM (including lungs), using L ¼ 300 and applying the second order Tikhonov regularization. All subsequent details illustrated below refer to this optimal case. A single snapshot of the Cest images is shown in Figure 4. Its timing is at t ¼ 30 ms, the same as is used for the reference image shown in the left panel of Figure 3. A similar comparison with the endocardial potential field is not included since the EPD based estimate was restricted to the pericardium. These results demonstrate the feasibility of the procedure indicated by the part pointing to the right of the loop shown at the bottom of Figure 1. A comparison of the ELGs drawn from Cest, some of its rows, and the corresponding ones from Cref is presented in Figure 5. Note the lower peak to peak values in the estimated ELGs as is also apparent in column 6 of Table 1.
Downloaded from by guest on February 17, 2015
100 90 80 70 60 50 40 30 20 10
iv126
A. van Oosterom
100 90 80 70 60 50 40 30 20 10
10 5 0 –5 –10
Figure 4 Estimated potential field Cest at the 807 nodes of the epicardial surface at 30 ms from onset QRS. Statistics of this case listed on the fourth row of Table 1. Natural anterior view. Colour code of the potential in mV. Supporting isopotential lines at 2 mV intervals.
Figure 6 Timing of depolarization at the 807 nodes of the epicardium. Left panel: as identified on the basis of Cref. Right panel: those identified on the basis of Cest. Compare with the left panel of Figure 2. Natural anterior view, white circle marks the locus of the earliest epicardial breakthrough. Colour code of the timing (ms). Supporting isochrones drawn at 10 ms intervals.
ELGS_ref and_est
522
569
205
324
273
56
78
107
15.000 mV 400 (ms)
Figure 5 Electrograms extracted from nine locations on the heart surface. The middle one is located at the position marked by the white circle in Figure 4, the other eight at locations of 3 cm from the central one. In blue: the signals extracted from Cref. The red ones, correspondingly, from Cest. (Figure 4).
The final segment of the loop, pointing to the left, deals with deducing the timing of depolarization and repolarization from the ELGs. As discussed, these are identified as the timings of the inflection points of the downstrokes of the ELGs (depolarization) and of their final upstrokes (repolarization), respectively. As can be seen by comparing the left panel of Figure 2 with those in Figure 6, the timing of depolarization derived from Cref as well as the timing derived from Cest closely correspond to the EDL reference d.
Discussion Of late, the inverse problem of electrocardiography, now commonly referred as ECGI, has attained clinical interest.45 The main interest has been the creation of images of the epicardial potential field as
estimated from recorded body surface potentials, the track illustrated on the right of Figure 1. Initially, the production of epicardial potential maps formed the main goal. Subsequently, the temporal aspect of the inversely estimated potential field, the ELGs, and their timing was included in the EPD-based method. This is the top-down track illustrated on the right of Figure 1, arguably the best- known track. The SM involved is the EPD. This paper aims at emphasizing the existence of an alternative track that may equally be followed to reach the same objectives, i.e. the EDL-based approach. The question then remains which route, if any, is optimal. The performances of both approaches have mostly been documented for each of the methods separately, and using different test conditions. Most of the results published are heralded as ‘validations’, rather than using just ‘evaluations’ or ‘test reports’, or ‘non-falsifications’.46 Early ‘comparing
Downloaded from by guest on February 17, 2015
381
A comparison of electrocardiographic imaging
Conclusions This paper shows how a reference distribution Cref that pertains to the complexity of the signals of a sinus beat may be used in the testing of imaging variants. Such reference data are unfortunately not available from clinical studies. Studies in which the references are drawn from stimulated activity tend to provoke too optimistic expectations regarding the performance of ECGI methods in more complex situations. The quality of the EPD-based Cest results documented here may possibly be inferior to other imaging methods that have been described in the literature. What remains is the conclusion that the EDL- and EPD-based source variants deserve to be studied alongside each other in the future development of electrocardiographic imaging. Conflict of interest: none declared.
References 1. van Oosterom A. The use of the spatial covariance in computing pericardial potentials. IEEE Trans Biomed Eng 1999;46:778–87.
2. van Oosterom A. Repolarization features as detectable from electrograms and electrocardiograms. J Electrocardiol 2013;46:55– 560. 3. Taccardi B. Distribution of heart potentials on the thoracic surface of normal human subjects. Circ Res 1963;00:341 –52. 4. MacLeod R, Buist M. The forward problem of elecrocardiography. In: Macfarlane PW, van Oosterom A, Janse MC, Kligfield P, Camm J, Pahlm O, eds. Basic electrocardiology. London: Springer, 2012. 5. Pullan AJ, Cheng LK, Nash MP, Ghodrati A, MacLeod R, Brooks DH. The inverse problem of elecrocardiography. In: Macfarlane PW, van Oosterom A, Janse MC, Kligfield P, Camm J, Pahlm O, eds. Basic electrocardiology. London: Springer, 2012. ´´ ber die Richtung und die manifeste Gro´´sse der 6. Einthoven W, Fahr G, de Waart A. U Potential Schwankungen im menschlichen Herzen und u´´ber den Einfluss der Herzlage auf die Form des Elektrokardiogramms. Pflugers Arch 1913;150:275 –315. 7. Wilson FN, Macleod AG, Barker PS. The distribution of action currents produced by the heart muscle and other excitable tissues immersed in conducting media. J Gen Physiol 1933;16:423 –56. 8. Geselowitz DB. On the theory of the electrocardiogram. Proc IEEE 1989;77:857 – 76. 9. Geselowitz DB. Description of cardiac sources in anisotropic cardiac muscle. J Electrocardiogr 1992;25:65– 7. 10. van Oosterom A. The equivalent double layer: source models for repolarization. In: Macfarlane PW, van Oosterom A, Pahlm O, Kligfield P, Janse MC, Camm J, eds. Basic electrocardiology. London: Springer, 2012. 227 –46. 11. van Bladel J. Electromagnetic fields. New York: McGraw-Hill, 1964. 12. Plonsey R, van Oosterom A. Introductory physics and mathematics. In: Macfarlane PW, van Oosterom A, Pahlm O, Kligfield P, Janse MC, Camm J, eds. Basic electrocardiology. London: Springer, 2012. 13. Martin RO. Inverse electrocardiography. PhD, Duke, NC: Duke University, 1970. 14. Barr RC, Ramsey M, Spach MS. Relating epicardial to body surface potentials by means of transfer coefficients based on geometry measurements. IEEE Trans Biomed Eng 1977;24:1 –11. 15. Colli-Franzone PC, Guerri L, Taccardi B, Viganotti C. The direct and inverse potential problems in electrocardiology. Numerical aspects of some regularization methods and application to data collected in dog heart experiments. Pavia: IAN-CNR, 1979. 16. Panofski WKH, Phillips M. Classical electricity and magnetism. London: AddisonWesley, 1962. 17. van Oosterom A, Oostendorp TF. On computing pericardial potentials and current densities in inverse electrocardiography. J Electrocardiol 1993;25S:102–6. 18. Munck JCd. A linear discretization of the volume conductor boundary integral equation using analytically integrated elements. IEEE Trans Biomed Eng 1992;39:986 – 90. 19. van Oosterom A. Closed-form analytical expressions for the potential fields generated by triangular monolayers with linearly distributed source strength. Med Biol Eng Comput 2012;55:1 –9. 20. Oostendorp TF, van Oosterom A. Source parameter estimation in inhomogeneous volume conductors of arbitrary shape. IEEE Trans Biomed Eng 1989;36:382–91. 21. Lawson CL, Hanson RJ. Solving least squares problems. Englewood Cliffs, NJ: PrenticeHall, 1974. 22. Beck JV, Arnold KJ. Parameter estimation in engineering and science. 1977. 23. Gill PE, Murray W, Wright MH. Practical optimization. New York: Academic Press, 1981. 24. Forsythe GE, Malcolm MA, Moler CB. Computer methods for mathematical computations. Englewood Cliffs, NJ: Prentice-Hall, 1977. 25. Marquardt DW. An algorithm for least-squares estimation of non-linear parameters. J Soc Indust Appl Math 1963;11:431–41. 26. Salu Y. Relating the multipole moments of the heart to activated parts of the epicardium and endocardium. Ann Biomed Eng 1978;6:492 –505. 27. van Oosterom A. Solidifying the solid angle. J Electrocardiol 2002;35S:181 –92. 28. Cuppen JJM, van Oosterom A. Model studies with the inversely calculated isochrones of ventricular depolarization. IEEE Trans Biomed Eng 1984;31:652 –9. 29. Huiskamp GJM, van Oosterom A. The depolarization sequence of the human heart surface computed from measured body surface potentials. IEEE Trans Biomed Eng 1989;35:1047 – 58. 30. Modre R, Tilg B, Fisher G, Hanser F, Messarz B, Seger M et al. Atrial noninvasive activation mapping of paced rhythm data. J Cardiovasc Electrophysiol 2003;14:712 –9. 31. van Oosterom A, Jacquemet V. A parameterized description of transmembrane potentials used in forward and inverse procedures. Folia Electrocardiol 2005;12: 111 –3. 32. van Oosterom A. Genesis of the T wave as based on an equivalent surface source model. J Electrocardiogr 2001;34:217 –27. 33. van Oosterom A, Jacquemet V. Genesis of the P wave: atrial signals as generated by the equivalent double layer source model. Europace 2005;7:S21 –29. 34. Huiskamp GJM. Difference formulas for the surface Laplacian on a triangulated surface. J Comp Physiol 1991;95:477 –96. 35. van Dam PM, Oostendorp TF, Linnenbank AC, van Oosterom A. Non-invasive imaging of cardiac activation and recovery. Ann Biomed Eng 2009;37:1739 –56.
Downloaded from by guest on February 17, 2015
notes’ studies were performed at the Nora Eccles CVRTI, Utah, where EDL-based inverse timings47,48 as well as potential fields were compared with those observed in torso tank measurements. The loop as outlined at the bottom of Figure 1 has not previously been explored. By starting the loop clockwise at the EDL-based timing (left to right) the ‘forward computation’ of the epicardial potential field is feasible and the similarity between Cref and Cest convincing (left panel of Figure 3 vs. Figure 4). The same applies to the comparison of the EDL-based ELGs and EPD-based versions (Figure 5). A difference in the magnitude of the ELG peak-to-peak values at activation and an accompanying decrease in the steepness of their slopes can be observed. These effects were more profound in the less optimal cases listed in Table 1. Finally, the return segment of the loop, the activation maps based on the timing of the downslopes of Cref and Cest, dref and dest respectively, corresponded well with the EDL reference timing (d) as illustrated in Figure 6. The linear correlation coefficients (CC) derived from the 806 epicardial data were CC(dref, d) ¼ 0.999 and CC(dest, d) ¼ 0.886. The corresponding CC values pertaining to the timing of repolarization (maximal upslopes of the ELGs) were CC(rref, r) ¼ 0.751 and CC(rest, r) ¼ 0.638. Although the correlation values listed may pass some significance level, the repolarization maps did not resemble well those of the reference data. The fact that even CC(rref, r) was low raises once more the question on the correct interpretation of such markers for the timing of repolarization in the presence of the full complexity of the diffuse myocardial repolarization process.2 The use of model-to-model studies, such as the one described here, is valuable only for providing a test under ideal conditions: an available gold standard, no VCM modeling errors, no noise. As such the results presented here may be taken as an upper limit of the quality of the estimates. Several such studies have been described in the literature, most recently one dedicated to the analysis of regularization variants.39 The relevance of such studies depends on the realism of both the VCM and the source description.
iv127
iv128 36. van Dam PM, van Oosterom A. Atrial excitation assuming uniform propagation. J Cardiovasc Electrophysiol 2003;14:S166 –71. 37. Walsh GR. Methods of optimization. London: John Wiley & Sons, 1975. 38. Gulrajani RM, Roberge FA, Savard P. The inverse problem of electrocardiography. Compr Electrocardiol 1989;I:237 –88. 39. Milanicˇ M, Jazbinsˇek V, MacLeod RS, Brooks DH, Hren R. Assessment of regularization techniques for electricardiographic imaging. J Electrocardiol 2014;47:20– 8. 40. van Dam P, Oostendorp T, van Oosterom A. Application of the fastest route algorithm in the interactive simulation of the effect of local ischemia on the ECG. Med Biol Eng Comput 2009;47:11 –20. 41. van Oosterom A, Oostendorp TF. ECGSIM: an interactive tool for studying the genesis of QRST waveforms. Heart 2004;90:165 – 8. 42. van Oosterom A, Oostendorp TF. Cardiac simulation for education: the electrocardiogram according to ECGSIM. In: Pahlm O, Wagner GS, eds. Cardiovascular multimodal image guided diagnosis and therapy. New York: McGraw Hill, 2010. 263 –80.
A. van Oosterom
43. Heringa A, Uijen GJH, van Dam RT. A 64-channel system for body surface potential mapping. In: Antalo´zcy Z, Pre´da I (eds). Electrocardiology 1981. Budapest: Hungary, 1982. 297 –301. 44. Hansen PC, Jensen TK, Rodriguez G. An adaptive pruning algorithm for the discrete L-curve criterion. J Comput Appl Math 2007;198:483 – 92. 45. Rudy Y. Noninvasive electrocardiographic imaging of arrhythmogenic substrates in humans. Circ Res 2013;112:834 –48. 46. Popper K. The logic of scientific discovery. London: Routledge, 1992. 47. Oostendorp TF, MacLeod RS, van Oosterom A. Non-invasive determination of the activation sequence of the heart: validation with invasive data. in Proc 19th Int Conf IEEE/EMBS, Chicago, 1997; pp 335 –8. doi:10.1109/IMBS.1997.754543. 48. Oster HS, Taccardi B, Lux RL, Ershler PR, Rudy Y. Noninvasive electrocardiographic imaging: reconstruction of epicardial potentials, electrograms, and isochrones and localization of single and multiple electrocardiographic events. Circulation 1997;96: 1012 –24.
Downloaded from by guest on February 17, 2015
SUPPLEMENT: REVIEW
A comparison of electrocardiographic imaging based on two source types Adriaan van Oosterom* Radboud University Nijmegen, Nijmegen, The Netherlands Received 30 July 2014; accepted after revision 2 August 2014
----------------------------------------------------------------------------------------------------------------------------------------------------------Keywords
Electrocardiographic imaging † Inverse problem † Electric source descriptions
Introduction The dynamics of the contraction of the cardiac muscle is initiated and accompanied by the flow of electric current. These set up potential fields throughout the myocardium as well as in the body tissues surrounding the heart and on the surface bounding the current flow: the thorax. The observed time course of the potential difference between any two locations inside the thorax or on its surface is referred to as an electrocardiogram (ECG). In cardiology, the set of 12 simultaneously recorded signals, derived from nine strategically chosen locations on the body surface, forms one of the most frequently used diagnostic tools. As such, the interpretation of the ECG can be classed as a non-invasive method. Moreover, it can be classed as performing an inverse procedure, being the interpretation of a state, or phenomenon on the basis of observations made at a distance. Over the period of more than a century following the first recording of an ECG, the diagnostic accuracy of the ECG has steadily increased. Yet, in some types of cardiac disease the accuracy still needs to be improved. One of the later developments has been the increase in the number of electrodes sampling the time course of the potential.
The display of a sequence of instantaneous potential fields on the thorax by means of scalar maps, body surface potential maps (BSPMs),3 has inspired the search for methods which also display the interpretation of the signals close to the heart, on its bounding surface (epicardium and endocardium) or intramurally. These emphasize the spatial character of the cardiac electric sources. This procedure has recently been known as imaging. The image formation involved in biomedical imaging rests on the nature of the physics involved. In many of these, the image formation is based on the interaction of the applied, or induced, radiation encountered in the tissue while passing along it in straight lines. In contrast, the electric potential field in the bounded region of the thorax is formed along electric current density lines that form loops, starting at exit points within the cardiac tissue (the sources) and returning to nearby locations (the sinks). The charge carried along the current density lines is merely pumped around. This property necessitates the formulation of an electric source model (SM) as well as of the properties (electric conductivities) of the electrically passive tissues external to the heart: volume conductor model (VCM).
* Corresponding author. Tel: +31 485 316543. E-mail address: [email protected] Published on behalf of the European Society of Cardiology. All rights reserved. & The Author 2014. For permissions please email: [email protected].
Downloaded from by guest on February 17, 2015
The aim of the study to compare the performance of two major source types involved in the imaging of the electric activity of the heart on the basis of potential differences observed on the thorax. The images depict either the timing of activation and repolarization of the myocardium or the potential field on a surface closely encompassing the myocardium. The depolarization and repolarization timing on a closed surface bounding the ventricular myocardium was derived from measured body surface potentials (BSPs), an MRI-based electric volume conductor model comprising the geometry of thorax, lungs, heart surface Sh, and cavities. The solution was constrained by using a template of the local transmembrane potentials (TMPs). The latter serve as the strength of the Equivalent Double Layer (EDL) source model (SM), which was used to compute the potential field on Sh (epicardium and endocardium). The second SM is the potential distribution on the epicardium Sp, referred to here as the Equivalent Potential Distribution (EPD). Its strength was estimated directly from the BSPs. The inflection points of the estimated electrograms (ELGs) were taken as markers of the timing of local depolarization and repolarization. The endocardial potential fields estimated using both sources exhibited qualitative similarity, as did the ELGs. With reference to the one generated by the EDL source, the magnitude of the estimated endocardial EPD field was smaller, the downslopes of the ELGs were lower. The timing of depolarization estimated from the EDL-based ELGs was highly correlated with those of the TMP templates, the EPD-based correlation was lower. For the repolarization timing the corresponding test values indicated an insufficient similarity. The EDL- and EPD-based source variants deserve to be studied alongside each other in the future development of electrocardiographic imaging.
iv121
A comparison of electrocardiographic imaging
What’s new? † An analysis of the similarity between potential fields, electrograms (ELGs), and timing as estimated based on either the EDL or the EPD source model (SM). † An explicit demonstration of the potential field on the closed surface closing the myocardium as generated by the EDL SM. † An explicit demonstration of the ELGs at the closed surface closing the myocardium (endocardium and epicardium) as generated by the EDL SM. † A reconfirmation1 of the value of including the inhomogeneity of the lungs in EPD-based imaging. † Fresh food for doubt on the accuracy of estimating the timing of the repolarization at Sh, taken to be the timing of the inflection point of the transmembrane potential downslope, from the timings of the inflection points of the upstrokes of the nearby ELGs.2
F = A S,
(1)
f1 (r ′ , t) =
1 4ps
∇†(si (r )∇Vm (r , t)) dv(r ). R V
(2)
The integral is taken over the volume encompassing all myocytes, ∇ denotes the spatial gradient operator, ∇† the divergence operator, si the intracellular conductivity, and R the distance between source location r and external observation point r ′ . The term ∇†(si (r )∇ Vm (r, t)) represents the electric current source density within the myocardium, the strength of which is proportional to the varying source strength of the local TMP: Vm (r , t). By using the approximation si (r ) = si and applying Green’s second theorem (a standard result from differential calculus)11 Eq. 2 is equivalent to F1 (r ′ , t) =
−s i Vm (r , t) dv(r ′ ,r ), 4ps
(3)
S
in which F denotes a matrix (size: L,T) comprising the documentation of the potential differences between all observation points (field points, ℓ ¼ 1 . . . L) and the potential at a common reference (CR) at all instances (t ¼ 1 . . . T ) studied. Moreover, a matrix S (size: N,T ) expresses the source strength at all locations (n ¼ 1 . . . N ) of the source elements as a function of time. With the source strengths depicted as images, the elements of S may be viewed as the N pixels of the sequence of all T images, stills of a movie. Finally, matrix A (size: L,N ) represents the links (transfer) between all source elements and the potential at the observation points, i.e. the electrode locations. The computation of the source matrix S on the basis of Eq. 1, a set of observations F and a computed transfer A, is generally referred to as solving the inverse problem of electrocardiography. 5 The specific solution to inverse problems in which the solution S can be represented by images is referred to as electrocardiographic imaging (ECGI). Several other acronyms related to variations in the involved source description have appeared in the literature. In this paper two of the main source descriptions are discussed and their application to the data pertaining to one and the same subject is illustrated. Next, their similarities, differences, advantages, and disadvantages are discussed.
A tale of two source types The equivalent double layer Soon after the introduction of the electric current dipole,6 Wilson et al. 7 introduced an equivalent SM that is closely linked to bioelectricity: the uniform double layer (UDL). In the UDL the source is
with dv(r ′ ,r ) denoting the solid angle subtended by small segments of Sh, the surface bounding V. This indicates that within the limitation of the assumed homogeneity of the intracellular domain, all source activity within the source volume V (Eq. 2) may be specified wholly by the distribution of Vm(t) on the surface Sh of the source region. It forms the generalization of the classic EDL theory. The result was formulated by Geselowitz.8 In a later publication he showed that it also holds true for anisotropic tissues, on condition the anisotropy ratio is uniform.9 Equation 3 forms the basis for solving the inverse problem aimed at estimating the timing of depolarization and repolarization, but may also be used to compute the epicardial potential field as is discussed in the section ‘An evaluation’ (Figure 1). The numerical expression for the forward computation is F = BF1 = BV Vm = DVm ,
(4)
in which Vm is a matrix of size (N,T ), with N being the number of nodes of the mesh representing S. Matrix V of size (M,N ), expresses the numerical variants of dv(r ′ ,r ), scaled by the remaining factors in Eq. 2, with M being the number of field points external to Sh for which f1 (r ′ , t) is to be known. Matrix B, of size (L,M), accounts for the effects on the potential of any inhomogeneity in the electrical conductivity of the external domain. These M locations include those at the body surface Sb. If the medium between the thorax surface Sb and Sh is homogenous then M ¼ L. The elements of B can be computed, e.g. by means of the boundary element method (Section. 2.6.4.12).
Downloaded from by guest on February 17, 2015
Upon having specified both models, the electric potential field at any set of electrode locations can be computed. This computation is known as solving the ‘forward problem of electrocardiography’.4 A generic expression for the outcome of this computation is the matrix multiplication
restricted to the boundary of a propagating activation front, with its strength taken to be uniform over the wave front. The equivalent double layer (EDL) is the successor to the UDL; its elementary sources are located on the closed surface Sh bounding the myocardium (either ventricles or atria). It extends the effectiveness of the UDL concept to the repolarization phase.8 – 10 The EDL concept is based on the general expression linking transmembrane potential (TMP) Vm (r , t) of the active myocytes to the field potential f1 (r, t) in an infinite homeogeneous external medium, a generalization of the classic cable theory. It reads:
iv122
A. van Oosterom
with V being the matrix of elementary solid angles and G the matrix of the weights dS(r )/R(r ′ ,r ) involved in the integration. In both types of matrices, the respective scaling factors seen in Eq. 5 have been absorbed. The use of Eq. 5 requires the potential field at the boundary and its normal gradients to be known. Equation 5 may also be used to compute the link between the fields on the boundaries S1 and S2. As an example, the most commonly considered situation is where partial, but nonetheless useful information is available. Since the normal gradient of the potential at the isolated torso surface is zero, matrix J1 is zero. The field F2 acts as a potential source distribution.14 The resulting potential field at S1 is computed as follows. Equation 6 may be applied to the potential fields just inside the two bounding surfaces. By combining the resulting two equations and, after rearranging terms, it follows that F1 can be found by solving the system of linear equations
Geometry of torso, heart and body tissue conductivity
Measured BSPs F
Source models S1; Vm on Sh
S2; Y on epicard
Vm=Vm(d,r) TMP template d,r Vm=Vm(d,r)
Constraints priors
Y: norm; smoothness; temporal
Solutions Forward computation
Y: images electrograms
Timing inflection points
Figure 1 Flowchart of the procedures involved in the imaging based on the two different source descriptions discussed. On the left: pertaining to the EDL source, on the right: to the EPD source. Note the clockwise loop shown at the bottom.
1 1 1 Jn (r ) dS(r ) F1 (r )dv + F2 (r )dv + 4p 4p 4ps R 1 Jn (r ) dS(r ), + 4ps R
S2
S1
(5)
S2
∂F (unit: A/m2) (a scalar function) expresses the ∂n magnitude of the local gradient of the potential field normal to the surfaces scaled by the local electric conductivity. The potential field Fi, Eq. 6, at any field point i in V can be computed numerically from the expression where Jn (r) = s
Fi = Vi,1 F1 + Vi,2 F2 + Gi,1 J1 + Gi,2 J2 ,
− G1,2 G2,2
F1 J2
=
V1,2 I − V2,2
F2 ,
(7)
(6)
A F2 = F1 ,
(8)
which is the essential expression for computing the pericardial potential distribution F2 on S2 as estimated on the basis of observations F1 on S1. The derivation can be generalized for cases in which V comprises inhomogeneous, closed, non-intersecting source-free regions such as poorly conducting lungs.17
Numerical aspects The integrals as appear in Eqs. 3 and 5 need to be computed numerically. The functions involved (potentials, source densities) are integrated over small surface elements. The most convenient shape of such elements is the planar triangle. Closed form analytical expressions are available for the integral over triangular elements in which the function over the triangle is approximated to be bilinear, specified by the function values at the vertices. For the integration involving solid angles, the analytical expression can be found in Munck,18 the one for the distributed monolayer in a more recent paper.19 For the results involving the solid angles, the potential field is discontinuous, i.e. when crossing the surface element the solid angle changes by 4p, scaled by a factor comprising the two conductivities on either side of the surface. Hence, the value assigned to the integral at or close to the surface needs to be tuned to the application in hand. In contrast, the potential field on the triangle generated by a monolayer on the triangle itself is continuous whereas the component of its gradient along the normal of the surface is discontinuous.
Summarizing The basic, numerical expressions linking SM strength to observed field potentials as discussed in this section are Eqs. 4 and 8. The matrix formulations used, incorporate the temporal nature of source strengths Vm (EDL) and F2 (equivalent potential distribution,
Downloaded from by guest on February 17, 2015
The second source variant is the potential distribution on a closed surface encompassing Sh generated by its internal bioelectric sources, such as the pericardium. Its strength is to be estimated on the basis of potentials observed on the thorax.13 – 15 In experimental studies the surface involved is frequently taken to be a closed variant of the epicardium. The thorax and pericardium form two nested non-intersecting surfaces. No active bioelectric sources are assumed to be present in the domain bounded by these two surfaces. Below, the essence is shown of the computation of the forward transfer. It was first derived by Barr et al. 14 It is illustrated by considering a source-free domain V with uniform isotropic conductivity s bounded by closed surfaces S1 and S2, with S2 contained in S1. The derivation is shown here in some detail, motivated by the high relevance of the transfer. As shown in page 12 of the reference,16 the potential field within the source-free domain depends uniquely on the potential and the normal gradient of the potential at its boundaries. The general expression reads
S1
I − V1,1 V2,1
with I being the identity matrix having the size of the accompanying matrix. From the solution to this expression a forward transfer matrix A is isolated, yielding the linear system
The equivalent potential distribution
F(r ′ ) =
iv123
A comparison of electrocardiographic imaging
EPD), respectively. Matrices D and A, which express the volume conduction effects only, are taken to be static. The isolation of matrices D and A in the forward procedure increases the transparency of the handling of the inverse procedure.20 In Eq. 4, matrix D is seen as transferring the current double layer source density, whose strength is proportional to Vm at the heart surface Sh, to the potential field inside a bounded medium, viz. the thorax, as well as on its boundary. In Eq. 8, the transfer A expresses the effect of the potential distribution on the pericardium.
parameters p are estimated from Eq. 9 by treating it as a least-squares problem. The source strength Vm(p) depends non-linearly on the timing parameters p (mainly stretching and/or shifting the TMP shape). As a consequence, finding p from Eq. 9 poses a so-called non-linear least-squares (NLLS) problem. Such problems need to be solved iteratively, starting from a posed initial solution p0. In the example shown in the next section, the solution method used is the Marquardt algorithm.25 Details of this procedure may be found in references therein.32,33
The imaging problem
Including constraints
The equivalent double layer A paper by Salu26 provided the first introduction to the EDL positioned on the surface Sh bounding the entire myocardium, with unit strength on depolarized sections of Sh and zero strength on the remaining parts. The link between UDL and EDL is based on the properties of solid angles.27 The EDL as the source implied in an inverse procedure was introduced some years later.28 In the years thereafter, various studies on possible applications and improvements followed.29,30 The current state of the development is based on the forward formulation in Eq. 4, in which a template of the TMP is included. The shape of the TMP source strength at each node n (n ¼ 1 . . . N ) of the N nodes (vertices) of a triangulated representation of Sh is specified by an analytical function. Its parameters include the timing of depolarization (dn) and repolarization (rn). These may be denoted by the column vectors d and r. The nature of the analytical function is such that these two parameters relate to the inflection points of the upslope of the local TMP and of their final downslope, respectively.31 A full description of the procedure is included in van Oosterom.32 The estimation is that of finding the parameter set that minimizes the least-squares difference between the observed signals and those resulting from the model ‘prediction’. We first reformulate Eq. 4 as DVm = DVm (p) = F,
(9)
with p denoting the timing parameters (either d or r) setting the shape of the TMP, and F the CR signals observed on the thorax. The latter may be a subset of the potentials at all nodes involved in the numerical representation of the thorax boundary. The timing
Results of numerical experiments based on the estimation of p from Eq. 9 indicated that additional constraints are needed if realistic, physiologically sound solutions are to be found. These constraints may be combined in Eq. 9. Different constraints were tested; the one adopted29 is the Frobenius norm of the surface Laplacian,11 i.e. ||lL p||F with matrix L denoting the surface Laplacian operator34 and l a scalar variable setting the weight of this constraint. The norm involved is the square root of the sum of the squares of all matrix elements, the so-called Frobenius norm. The value of l may be set to drive the estimated p toward a spatial smoothness having a physiologically realistic value.29 The iterative procedure for solving an NLLS problem may end up at a local minimum of the parameter space rather than at the desired global minimum. The quality, i.e. realism, of the solution found depends greatly on the realism of the initial estimate p0. The initial estimate used in recent work33,35,36 is based on the fastest route algorithm based on Huygens’ principle applied to propagating activation wave fronts.
The equivalent potential transfer Equation 8 specifies the forward transfer A between two potential fields on two nested surfaces: the thorax surface Sb and the pericardial surface Sp. The inverse of matrix A may be used to estimate the pericardial potential C from F. In contrast to the inverse procedure for the EDL, the basic relationship between SM C and observations F is linear. It may seem as if this type of inverse merely requires computing the inverse of the forward transfer matrix A. However, the existence of a unique inverse solution C depends on the properties of A. This has a dominant role in these types of inverse computations, hence the following, brief recapitulation is in order.
Regularization Just as when dealing with the EDL source, in studies estimating the EPD it was found to be necessary to include additional constraints. This was necessary even if the number of leads L used in sampling the potential field on the body surface was greater than the number N of the nodes on Sp. Constraining the solutions is carried out by demanding the solution to approximate CC R, with C and R specifying the nature of the constraint and its value aimed for, respectively. If the constraint can be expressed by the matrix formulation CC R, the estimate C is found as the linear least squares (LLS) solution to
A
l C
C
F l R
(10)
Downloaded from by guest on February 17, 2015
In this section, some basic methods involved in solving the inverse problem are discussed with regard to their application to the imaging of either of the two source types introduced in the previous section. The basic, numerical expressions linking SM strength to observed field potentials are Eqs. 4 and 8. In the sequel, the thorax potentials are denoted by F, and the potential distribution on the pericardium by C. The inverse problem is treated from the perspective of source parameter estimation. The parameters specifying the involved VCM, such as geometry and passive tissue properties, are taken to be known. Fundamental methods for parameter estimation in a wide variety of scientific domains have been developed over the last 50 years. Major reference books that widely cover the theory are provided under references.21 – 25
iv124
A. van Oosterom
The constant l is the Lagrange multiplier as used in the classic, constrained optimization theory.37 The higher the values of l, the more the estimate is forced to comply with the constraint. The procedure is also referred to as regularizing the solution.38 The standard expression for the LLS solution to Eq. 10 is C(l) = (AT A + l2 CT C)−g (AT F + l2 CT R).
(11)
For each value of l, this produces a unique, minimum norm solution, which poses the question of which value l should be taken. For a recent discussion on this problem see Milanicˇ et al.39
Comparing both imaging variants
Below, an evaluation is presented of the theory exposed in Figure 1. It is based on a reference data set similar to the one included in the interactive forward simulation package41,42 available from www. ecgsim.org. These EDL derived data are taken as a gold standard against which the performances of the EPD estimates are compared. The specifications of the geometry of the relevant VCM surfaces, as well as the reference potential field F on the thorax used in this comparison, were identical. Material The reference data comprises the numerical representation of an inhomogeneous, multi-compartmental VCM. It includes MRI-derived geometry of the body surface, the lungs, the left cavities (combined ventricular and atrial), the right cavities (combined ventricular and atrial), and the surface Sh closing the ventricular myocardium. A three-fold higher electric conductivity than the surrounding space was assigned to the blood inside the cavities, and a five-fold lower value to the interiors of the surfaces representing the lungs. The timings (dn and rn) of the TMP templates at Nh ¼ 1500 nodes on Sh was derived from a dedicated inverse procedure The EDL SM was used to simulate body surface potentials at all L ¼ 300 nodes of the triangulated thorax geometry St. These formed the reference thorax potentials used in the sequel, referred to as F.35 The rows of this matrix F comprise T ¼ 401 subsequent samples at 1 ms intervals, starting at the onset of the QRS complex. A subset of these comprised the L ¼ 65 electrode locations specified by the Nijmegen-lead system.43 In addition, by including the nodes of Sh as field points in matrix B (Eq. 4), the potential field Ch on Sh was computed and taken as a reference field Cref. Methods The EPD estimates were derived as follows. A closed triangulated surface Sp was created comprising 807 nodes of the epicardium of Sh as well as 1090 Nc nodes specifying the surfaces of both atrial cavities. These in all Np ¼ 1897 nodes formed the vertices of Sp, which represent the pericardium. Its interior comprises all active (EDL) sources as well as the major inhomogeneity of the VCM: the blood inside the cavities. The potential field at the Np pericardial nodes was estimated by treating them as field points in the volume conductor specified in the previous section. The resulting field serves as the strength for the EPD source on Sp, referred to as Cest. Its size is Np,T. The domain bounded by Sp and Sb is source-free. In order to test the EPD-based imaging, two variants of the required transfer matrix A (Eq. 8) were computed. In one variant the effect of the inhomogeneity of the lungs is accounted for. In the other the conductivity of the domain is uniform. The estimation of the pericardial potential field was carried out on the basis of Eq. 11. The results based on the two types of regularization are documented: Tikhonov’ zero- and the second order. For the zero-order variant the constraint matrix C is a diagonal matrix of size N, for the second order it is the surface Laplacian operator L as specified previously. Its size is N,N. In both cases the term matrix R in Eq. 11 is zero and drops out of the expression. In applications to clinical data, where no gold standard is available, the regularization weight l needs to be found. A frequently used method is to scan a suitably selected range of l values and plot the norm of the regularized solution, ||CCest(l )||F, as a function of the relative residual
Downloaded from by guest on February 17, 2015
An overview of the two imaging recipes is displayed in the flow diagram shown in Figure 1. It summarizes and recapitulates the essence of the various topics discussed above. The ingredients indicated at the top of the figure are common in both procedures. The myocardium may be either that of the ventricles or that of the atria. The track on the left pertains to the EDL source. The parameters estimated are the timing of activation (dn) and recovery (rn) at Nh nodes on the closed surface Sh bounding the myocardium. The track on the right pertains to the EPD source. Here the parameters to be estimated are the potentials C(n,t) at the Np nodes on the pericardium at T time instances, relative to those present at the position of a CR electrode. The results may be displayed on the heart surface, colour-coded and/or supported by isofunction lines, for the EDL source: isochrones, for the EPD source: isopotential lines. As discussed, both modalities require the inclusion of some kind of prior knowledge that may be used for constraining the solution. For the EDL-based imaging method an important intrinsic constraint is highlighted: the template of a TMP wave form. An initial estimate is used based on the fastest route algorithm.40 For the EPD source various constraints have been tested, such as the norm of the solution, its spatial smoothness, or the inclusion of a temporal prior.5 The estimated timing parameters, the elements of d (depolarization) and r (repolarization) specify completely the time course of the TMP template and thus of the current source EDL strength at Sh. This source description may in turn be used to compute Ch on Sh, of which the epicardium is a part. In the lower part of Figure 1, this procedure is indicated by the arrows entering and leaving the box ‘forward computation’. At the lower right box the two imaging modalities meet. The results of this forward computation may be expected to be similar to the field Ce on the epicardial surface estimated from the EPD source description. The timing of depolarization and repolarization can be estimated from the inflection points of the electrograms (ELGs): Ce(n,t), the rows of matrix C. These may be compared with the corresponding elements of d and r as used in the EDL-based inverse procedure. This comparison is symbolized by the left pointing arrow at the bottom of the flow diagram (Figure 1). Note that this is a part of a clockwise loop.
An evaluation
iv125
A comparison of electrocardiographic imaging
Table 1 Descriptive statistics of the epicardial potential field Cest computed from F VCM
Order TIkh
L
RD (F)
0 0
300 300
1.85 e– 3 0.54 e– 3
STD (Cest) (mV)
P2P(Cest) (mV)
2.023
28.39
1.35 1.49
26.2 26.3
RE (Cest)
CC (Cest, Cref )
0.7076 0.6478
0.7071 0.7618
............................................................................................................................................................................... Reference lungs hom lungs hom
2
300
2.41 e– 3
1.44
25.0
0.6983
0.7152
lungs hom
2 0
300 65
0.53 e– 3 0.67 e– 3
1.60 1.12
25.7 19.4
0.6323 0.8174
0.7745 0.5760
lungs
0
65
0.44 e– 3
1.17
23.1
0.7982
0.6021
hom lungs
2 2
65 65
0.74 e– 4 0.47 e– 3
1.24 1.32
16.4 17.2
0.8070 0.7659
0.5912 0.6424
Upper row: some statistics of Cref, the epicardial field based on the EDL source. Next rows: estimation results based on the respective optimal l values. Columns: 1: VCM type; 2: order of applied Tikhonov regularizarion; 3: number of ECG signals; 4: relative residual difference between thorax potentials Fref and Fest ¼ A Cest; 5: STD(Cest); 6: extreme peak-to-peak values of Cest (downslope magnitude of the electrograms); 7: relative error between all elements of Cest and Cref; 8: linear correlation coefficient (CC) between all elements of Cest and Cref. Note the discrepancy between the magnitude of values of the relative residual differences: Column 4, and the relative errors: Column 7.
Figure 2 Reference timing of depolarization at the 1500 nodes of the ventricular surface; epicardium and endocardium. Left panel: natural anterior view, white circle marks the site of the epicardial breakthrough; right panel posterior view; white circle marks the locus of septal initiation of the activation. Colour code of the timing (ms). Supporting isochrones drawn at 10 ms intervals.
difference, RD ¼ ||F 2 ACest(l)||F/||F||F.44 Since the situation treated here does involve a pure gold standard, Cref, the scan can be carried out by using the error norm ||Cref 2 Cest(l)||F in a search to find its minimum. In Table 1 the resulting minimum relative errors were listed:, i.e. RE ¼ ||Cref 2 Cest||F/||Cref||F.
Results Figure 2 is an illustration of the timing of depolarization as included in the specification of the EDL source. The left panel depicts, in the natural anterior view, a colour-coded representation of the timing of depolarization on Sh, supported by isochrones lines drawn at 10 ms intervals. The pattern of the timing of repolarization, not shown here, has a similar complexity. The right panel depicts a posterior view, rotated additionally to offer a view of the left endocardium. Snapshots of the corresponding potential distributions on Sh at 30 ms from onset depolarization are presented in Figure 3. Note that this EDL-based procedure yields images both of the epicardium and the endocardium, of the timing as well as of the potential fields. In Table 1. the major descriptive statistics are presented of the test results for eight configurations, computed for the 807 epicardial
15 10 5 0 –5 –10 –15
Figure 3 Reference potential field Cref at the 1500 nodes of the ventricular surface at 30 ms from onset QRS. Left panel: natural anterior view, right panel posterior view, providing a view on the endocardium and the septum. Color code of the potential in mV. Supporting isopotential lines at 2 mV intervals.
reference data. The entries listed pertain to the respective l values yielding the smallest RE value. In addition to the RE values (Column 7), the standard deviation of all elements are documented, as well as the maximum over the N signals of the peak-to-peak value of each signal over time, the magnitude of the downslopes of the ELGs during depolarization and their upslopes during the final part of the repolarization. It shows that the smallest RE value is found when using an inhomogeneous VCM (including lungs), using L ¼ 300 and applying the second order Tikhonov regularization. All subsequent details illustrated below refer to this optimal case. A single snapshot of the Cest images is shown in Figure 4. Its timing is at t ¼ 30 ms, the same as is used for the reference image shown in the left panel of Figure 3. A similar comparison with the endocardial potential field is not included since the EPD based estimate was restricted to the pericardium. These results demonstrate the feasibility of the procedure indicated by the part pointing to the right of the loop shown at the bottom of Figure 1. A comparison of the ELGs drawn from Cest, some of its rows, and the corresponding ones from Cref is presented in Figure 5. Note the lower peak to peak values in the estimated ELGs as is also apparent in column 6 of Table 1.
Downloaded from by guest on February 17, 2015
100 90 80 70 60 50 40 30 20 10
iv126
A. van Oosterom
100 90 80 70 60 50 40 30 20 10
10 5 0 –5 –10
Figure 4 Estimated potential field Cest at the 807 nodes of the epicardial surface at 30 ms from onset QRS. Statistics of this case listed on the fourth row of Table 1. Natural anterior view. Colour code of the potential in mV. Supporting isopotential lines at 2 mV intervals.
Figure 6 Timing of depolarization at the 807 nodes of the epicardium. Left panel: as identified on the basis of Cref. Right panel: those identified on the basis of Cest. Compare with the left panel of Figure 2. Natural anterior view, white circle marks the locus of the earliest epicardial breakthrough. Colour code of the timing (ms). Supporting isochrones drawn at 10 ms intervals.
ELGS_ref and_est
522
569
205
324
273
56
78
107
15.000 mV 400 (ms)
Figure 5 Electrograms extracted from nine locations on the heart surface. The middle one is located at the position marked by the white circle in Figure 4, the other eight at locations of 3 cm from the central one. In blue: the signals extracted from Cref. The red ones, correspondingly, from Cest. (Figure 4).
The final segment of the loop, pointing to the left, deals with deducing the timing of depolarization and repolarization from the ELGs. As discussed, these are identified as the timings of the inflection points of the downstrokes of the ELGs (depolarization) and of their final upstrokes (repolarization), respectively. As can be seen by comparing the left panel of Figure 2 with those in Figure 6, the timing of depolarization derived from Cref as well as the timing derived from Cest closely correspond to the EDL reference d.
Discussion Of late, the inverse problem of electrocardiography, now commonly referred as ECGI, has attained clinical interest.45 The main interest has been the creation of images of the epicardial potential field as
estimated from recorded body surface potentials, the track illustrated on the right of Figure 1. Initially, the production of epicardial potential maps formed the main goal. Subsequently, the temporal aspect of the inversely estimated potential field, the ELGs, and their timing was included in the EPD-based method. This is the top-down track illustrated on the right of Figure 1, arguably the best- known track. The SM involved is the EPD. This paper aims at emphasizing the existence of an alternative track that may equally be followed to reach the same objectives, i.e. the EDL-based approach. The question then remains which route, if any, is optimal. The performances of both approaches have mostly been documented for each of the methods separately, and using different test conditions. Most of the results published are heralded as ‘validations’, rather than using just ‘evaluations’ or ‘test reports’, or ‘non-falsifications’.46 Early ‘comparing
Downloaded from by guest on February 17, 2015
381
A comparison of electrocardiographic imaging
Conclusions This paper shows how a reference distribution Cref that pertains to the complexity of the signals of a sinus beat may be used in the testing of imaging variants. Such reference data are unfortunately not available from clinical studies. Studies in which the references are drawn from stimulated activity tend to provoke too optimistic expectations regarding the performance of ECGI methods in more complex situations. The quality of the EPD-based Cest results documented here may possibly be inferior to other imaging methods that have been described in the literature. What remains is the conclusion that the EDL- and EPD-based source variants deserve to be studied alongside each other in the future development of electrocardiographic imaging. Conflict of interest: none declared.
References 1. van Oosterom A. The use of the spatial covariance in computing pericardial potentials. IEEE Trans Biomed Eng 1999;46:778–87.
2. van Oosterom A. Repolarization features as detectable from electrograms and electrocardiograms. J Electrocardiol 2013;46:55– 560. 3. Taccardi B. Distribution of heart potentials on the thoracic surface of normal human subjects. Circ Res 1963;00:341 –52. 4. MacLeod R, Buist M. The forward problem of elecrocardiography. In: Macfarlane PW, van Oosterom A, Janse MC, Kligfield P, Camm J, Pahlm O, eds. Basic electrocardiology. London: Springer, 2012. 5. Pullan AJ, Cheng LK, Nash MP, Ghodrati A, MacLeod R, Brooks DH. The inverse problem of elecrocardiography. In: Macfarlane PW, van Oosterom A, Janse MC, Kligfield P, Camm J, Pahlm O, eds. Basic electrocardiology. London: Springer, 2012. ´´ ber die Richtung und die manifeste Gro´´sse der 6. Einthoven W, Fahr G, de Waart A. U Potential Schwankungen im menschlichen Herzen und u´´ber den Einfluss der Herzlage auf die Form des Elektrokardiogramms. Pflugers Arch 1913;150:275 –315. 7. Wilson FN, Macleod AG, Barker PS. The distribution of action currents produced by the heart muscle and other excitable tissues immersed in conducting media. J Gen Physiol 1933;16:423 –56. 8. Geselowitz DB. On the theory of the electrocardiogram. Proc IEEE 1989;77:857 – 76. 9. Geselowitz DB. Description of cardiac sources in anisotropic cardiac muscle. J Electrocardiogr 1992;25:65– 7. 10. van Oosterom A. The equivalent double layer: source models for repolarization. In: Macfarlane PW, van Oosterom A, Pahlm O, Kligfield P, Janse MC, Camm J, eds. Basic electrocardiology. London: Springer, 2012. 227 –46. 11. van Bladel J. Electromagnetic fields. New York: McGraw-Hill, 1964. 12. Plonsey R, van Oosterom A. Introductory physics and mathematics. In: Macfarlane PW, van Oosterom A, Pahlm O, Kligfield P, Janse MC, Camm J, eds. Basic electrocardiology. London: Springer, 2012. 13. Martin RO. Inverse electrocardiography. PhD, Duke, NC: Duke University, 1970. 14. Barr RC, Ramsey M, Spach MS. Relating epicardial to body surface potentials by means of transfer coefficients based on geometry measurements. IEEE Trans Biomed Eng 1977;24:1 –11. 15. Colli-Franzone PC, Guerri L, Taccardi B, Viganotti C. The direct and inverse potential problems in electrocardiology. Numerical aspects of some regularization methods and application to data collected in dog heart experiments. Pavia: IAN-CNR, 1979. 16. Panofski WKH, Phillips M. Classical electricity and magnetism. London: AddisonWesley, 1962. 17. van Oosterom A, Oostendorp TF. On computing pericardial potentials and current densities in inverse electrocardiography. J Electrocardiol 1993;25S:102–6. 18. Munck JCd. A linear discretization of the volume conductor boundary integral equation using analytically integrated elements. IEEE Trans Biomed Eng 1992;39:986 – 90. 19. van Oosterom A. Closed-form analytical expressions for the potential fields generated by triangular monolayers with linearly distributed source strength. Med Biol Eng Comput 2012;55:1 –9. 20. Oostendorp TF, van Oosterom A. Source parameter estimation in inhomogeneous volume conductors of arbitrary shape. IEEE Trans Biomed Eng 1989;36:382–91. 21. Lawson CL, Hanson RJ. Solving least squares problems. Englewood Cliffs, NJ: PrenticeHall, 1974. 22. Beck JV, Arnold KJ. Parameter estimation in engineering and science. 1977. 23. Gill PE, Murray W, Wright MH. Practical optimization. New York: Academic Press, 1981. 24. Forsythe GE, Malcolm MA, Moler CB. Computer methods for mathematical computations. Englewood Cliffs, NJ: Prentice-Hall, 1977. 25. Marquardt DW. An algorithm for least-squares estimation of non-linear parameters. J Soc Indust Appl Math 1963;11:431–41. 26. Salu Y. Relating the multipole moments of the heart to activated parts of the epicardium and endocardium. Ann Biomed Eng 1978;6:492 –505. 27. van Oosterom A. Solidifying the solid angle. J Electrocardiol 2002;35S:181 –92. 28. Cuppen JJM, van Oosterom A. Model studies with the inversely calculated isochrones of ventricular depolarization. IEEE Trans Biomed Eng 1984;31:652 –9. 29. Huiskamp GJM, van Oosterom A. The depolarization sequence of the human heart surface computed from measured body surface potentials. IEEE Trans Biomed Eng 1989;35:1047 – 58. 30. Modre R, Tilg B, Fisher G, Hanser F, Messarz B, Seger M et al. Atrial noninvasive activation mapping of paced rhythm data. J Cardiovasc Electrophysiol 2003;14:712 –9. 31. van Oosterom A, Jacquemet V. A parameterized description of transmembrane potentials used in forward and inverse procedures. Folia Electrocardiol 2005;12: 111 –3. 32. van Oosterom A. Genesis of the T wave as based on an equivalent surface source model. J Electrocardiogr 2001;34:217 –27. 33. van Oosterom A, Jacquemet V. Genesis of the P wave: atrial signals as generated by the equivalent double layer source model. Europace 2005;7:S21 –29. 34. Huiskamp GJM. Difference formulas for the surface Laplacian on a triangulated surface. J Comp Physiol 1991;95:477 –96. 35. van Dam PM, Oostendorp TF, Linnenbank AC, van Oosterom A. Non-invasive imaging of cardiac activation and recovery. Ann Biomed Eng 2009;37:1739 –56.
Downloaded from by guest on February 17, 2015
notes’ studies were performed at the Nora Eccles CVRTI, Utah, where EDL-based inverse timings47,48 as well as potential fields were compared with those observed in torso tank measurements. The loop as outlined at the bottom of Figure 1 has not previously been explored. By starting the loop clockwise at the EDL-based timing (left to right) the ‘forward computation’ of the epicardial potential field is feasible and the similarity between Cref and Cest convincing (left panel of Figure 3 vs. Figure 4). The same applies to the comparison of the EDL-based ELGs and EPD-based versions (Figure 5). A difference in the magnitude of the ELG peak-to-peak values at activation and an accompanying decrease in the steepness of their slopes can be observed. These effects were more profound in the less optimal cases listed in Table 1. Finally, the return segment of the loop, the activation maps based on the timing of the downslopes of Cref and Cest, dref and dest respectively, corresponded well with the EDL reference timing (d) as illustrated in Figure 6. The linear correlation coefficients (CC) derived from the 806 epicardial data were CC(dref, d) ¼ 0.999 and CC(dest, d) ¼ 0.886. The corresponding CC values pertaining to the timing of repolarization (maximal upslopes of the ELGs) were CC(rref, r) ¼ 0.751 and CC(rest, r) ¼ 0.638. Although the correlation values listed may pass some significance level, the repolarization maps did not resemble well those of the reference data. The fact that even CC(rref, r) was low raises once more the question on the correct interpretation of such markers for the timing of repolarization in the presence of the full complexity of the diffuse myocardial repolarization process.2 The use of model-to-model studies, such as the one described here, is valuable only for providing a test under ideal conditions: an available gold standard, no VCM modeling errors, no noise. As such the results presented here may be taken as an upper limit of the quality of the estimates. Several such studies have been described in the literature, most recently one dedicated to the analysis of regularization variants.39 The relevance of such studies depends on the realism of both the VCM and the source description.
iv127
iv128 36. van Dam PM, van Oosterom A. Atrial excitation assuming uniform propagation. J Cardiovasc Electrophysiol 2003;14:S166 –71. 37. Walsh GR. Methods of optimization. London: John Wiley & Sons, 1975. 38. Gulrajani RM, Roberge FA, Savard P. The inverse problem of electrocardiography. Compr Electrocardiol 1989;I:237 –88. 39. Milanicˇ M, Jazbinsˇek V, MacLeod RS, Brooks DH, Hren R. Assessment of regularization techniques for electricardiographic imaging. J Electrocardiol 2014;47:20– 8. 40. van Dam P, Oostendorp T, van Oosterom A. Application of the fastest route algorithm in the interactive simulation of the effect of local ischemia on the ECG. Med Biol Eng Comput 2009;47:11 –20. 41. van Oosterom A, Oostendorp TF. ECGSIM: an interactive tool for studying the genesis of QRST waveforms. Heart 2004;90:165 – 8. 42. van Oosterom A, Oostendorp TF. Cardiac simulation for education: the electrocardiogram according to ECGSIM. In: Pahlm O, Wagner GS, eds. Cardiovascular multimodal image guided diagnosis and therapy. New York: McGraw Hill, 2010. 263 –80.
A. van Oosterom
43. Heringa A, Uijen GJH, van Dam RT. A 64-channel system for body surface potential mapping. In: Antalo´zcy Z, Pre´da I (eds). Electrocardiology 1981. Budapest: Hungary, 1982. 297 –301. 44. Hansen PC, Jensen TK, Rodriguez G. An adaptive pruning algorithm for the discrete L-curve criterion. J Comput Appl Math 2007;198:483 – 92. 45. Rudy Y. Noninvasive electrocardiographic imaging of arrhythmogenic substrates in humans. Circ Res 2013;112:834 –48. 46. Popper K. The logic of scientific discovery. London: Routledge, 1992. 47. Oostendorp TF, MacLeod RS, van Oosterom A. Non-invasive determination of the activation sequence of the heart: validation with invasive data. in Proc 19th Int Conf IEEE/EMBS, Chicago, 1997; pp 335 –8. doi:10.1109/IMBS.1997.754543. 48. Oster HS, Taccardi B, Lux RL, Ershler PR, Rudy Y. Noninvasive electrocardiographic imaging: reconstruction of epicardial potentials, electrograms, and isochrones and localization of single and multiple electrocardiographic events. Circulation 1997;96: 1012 –24.
Downloaded from by guest on February 17, 2015
