Multipole Approach to the Inverse Problem in Electrocardiology: Convergence of the Multipole Equivalent Generator on the Inhomogeneous Body Conductor

Bulletin of Mathematical Biology (2000) 62, 543–583 doi:10.1006/bulm.2000.0169 Available online at http://www.idealibrary.com on

Multipole Approach to the Inverse Problem in Electrocardiology: Convergence of the Multipole Equivalent Generator on the Inhomogeneous Body Conductor P. DE MUYNCK AND J. CORNELIS Vrije Universiteit Brussel, Electronics and Information Processing, ETRO, Pleinlaan 2, 1050 Brussel, Belgium L. I. TITOMIR Russian Academy of Sciences, Institute for Problems of Information Tranmission, B. Karetny 19, 101447 Moscow, Russia The multipole approach to the inverse electrocardiological problem consists of estimating the multipole components of the cardiac electric generator, starting from the measured body surface potential. This paper presents a critical investigation of the basic premise for the applicability of the multipole approach, namely the convergence of the multipole equivalent generator for the heart on the surface of an inhomogeneous body conductor. As an extension to multipole theory, a criterion for the convergence is derived. Based on realistic models for the body conductor and the cardiac electric generator, we observe that the criterion is not strictly satisfied in realistic conditions. Numerical simulations with the same models point out that the multipole equivalent generator is indeed not convergent in the strict mathematical sense. On the other hand, we show that the multipole equivalent generator yields a rather close approximation of the electrocardiological potential for intermediate values of the order of the multipole generator. A discussion is given on how to explain the apparently ambiguous results for the estimation of cardiac multipole components. c 2000 Society for Mathematical Biology

1.

I NTRODUCTION : G ENERAL

The inverse problem of electrocardiology consists of determining information about the electric activity of the heart from measurements of the potential on the human torso. Whereas an ECG is obtained via the measurement of the potential in 12 distinct points, the number of measurement electrodes that are needed for 0092-8240/00/030543 + 41

$35.00/0

c 2000 Society for Mathematical Biology

544

P. De Muynck et al.

the inverse problem is much higher, typically 64 or more electrodes spread over the anterior and posterior parts of the human torso. Since the electric activity of the heart results from individual heart cells that behave as current generators when they are activated, the electric activity of the heart is quantitatively described at a macroscopic level by a current source density distribution I , which varies in time and is defined over the heart region. Due to the nature of the cardiac current sources and the body volume conductor, the electromagnetic field equations which determine the relation between I and the potential reduce to a quasi-stationary system (Plonsey and Heppner, 1967), so that the cardiac potential is uniquely defined by I at any given time in the heart cycle, i.e., the forward problem is uniquely defined. With respect to the inverse problem, however, some basic mathematical considerations show that it is fundamentally impossible to determine I from the measured body surface potential in a unique way (Plonsey, 1966). Since the real cardiac current generator cannot be estimated, equivalent generators were defined instead, i.e., models for the cardiac current sources that include the necessary constraints so that the model parameters can be estimated from the body surface potential in a unique way. The multipole equivalent generator, which was introduced in 1960 (Geselowitz, 1960), satisfies this criterion. It consists of an infinite superposition of multipole current sources (dipole, quadrupole, etc.), which are located in a single physical point near the centre of the heart region, the multipole origin. The parameters of the multipole equivalent generator are the superposition (intensity) coefficients of individual multipole sources. The validity of the multipole equivalent generator is based on the theoretical consideration that a convergent representation of the external electrocardiological potential is obtained, provided that we choose superposition coefficients which are equal to the multipole components of the cardiac current source distribution I (Geselowitz, 1960). The cardiac multipole components in turn have a direct relation to the cardiac current source distribution, since they are mathematically defined as weighted integrals of I (see Section 2). Evidently, the multipole approach to the inverse electrocardiological problem consists of estimating the multipole components of the cardiac electric generator, starting from the measured body surface potential. The cardiac equivalent multipole generator has been used in this context for over 30 years, and it is still being used [see recent publications: Takano et al. (1996), Nolte and Curio (1997), Arthur (1999)]. As a consequence, it is abundantly present in the scientific literature and is also included in many books on bioelectricity, e.g., the recent books Malmivuo and Plonsey (1995) and Geselowitz (1997). The main characteristic of the multipole approach is its universal validity, independent of the physiological model for the cardiac current sources. Furthermore, it can be shown that the set of all multipole components represents the maximum information about the heart current sources that can be obtained from external measurements (Plonsey, 1966). Moreover, this information is in a form that depends solely on the heart current sources themselves and not on the specificities of the surrounding torso conductor: this means that by estimating the multipole components, one also tries to eliminate the effect on mea-

Multipole Approach to the Inverse Problem in Electrocardiology

545

sured electrocardiographic data of interindividual variations of the torso geometry. Later, in Section 4.1, this principle statement will have to be somewhat relativated, because the effect of intracavitary bloodmasses cannot be removed from the measured electrocardiographic data. The results of the multipole approach may be presented in a diagnostically meaningful way, by means of an estimate of the cardiac potential on a surface around the heart (quasi-epicardial potential) (Pilkington and Morrow, 1982; Titomir and Kneppo, 1994). This potential is obtained by substitution of the estimated multipole components into the multipole series expansion of the infinite or bounded medium potential. A general discussion on the usefulness of multipoles in electrocardiology is given by Pilkington and Morrow (1982). The multipole approach was also efficiently applied to the problem of joint analysis of the electric and magnetic fields generated by the same bioelectric generator (Titomir and Kneppo, 1985, 1994). Although the multipole approach has been extensively used for various purposes in electrocardiology, a number of precautions should be taken when it is applied in the case of the inhomogeneous conductor. Since the validity constraints for the multipole approach can often not be verified experimentally, this paper describes a careful physical-mathematical analysis which can be used as an explanatory model for the origin of several apparent ambiguities in the results of the multipole studies. Moreover, numerical simulations reveal the extent to which the multipole approach can be trusted. There are basically two methods for estimating the multipole components, namely the integration method (Geselowitz, 1960) and the regression method (Arthur et al., 1972). Initially these methods were mainly applied using finite homogeneous torso conductor models. For practical application of these methods, it is then only required to measure the external geometry of the torso (Kneppo and Titomir, 1979). The results of Dub´e et al. (1985), however, clearly indicate that one cannot expect to obtain correct estimates of the multipole components unless the torso inhomogeneities are taken into account. This implies that the integration method has to be ruled out, because it is intrinsically based on the assumption of a homogeneous torso conductor. The regression method is not inherently incompatible with the inhomogeneous torso model. However, the validity of the method depends on the convergence of the multipole equivalent generator on the torso surface. Although research has been done on the regression method for the inhomogeneous torso conductor, not much attention has been paid to the problem of how the internal inhomogeneities affect the convergence. Some essential considerations can be found in Dub´e et al. (1985), de Guise et al. (1985), De Muynck et al. (1986), Nyssen et al. (1986) and Gulrajani et al. (1988). In this paper, we present an explicit convergence criterion which is proved mathematically. The constraints imposed by the criterion are defined in relation to the geometries and the relative positions of the internal inhomogeneities on the one hand, and of the cardiac electric current source region on the other hand. A major

546

P. De Muynck et al.

improvement compared to previous works is that the mathematical proof provides an insight into the factors which affect the rate of convergence of the generator. Moreover, numerical simulations give an idea of the rate of convergence of the multipole equivalent generator under realistic conditions. The set of quantitative results on the rate of convergence is a new contribution to the literature. The organization of the paper is as follows. In Section 2, preliminary considerations are made concerning the relation between multipole convergence and the validity of the regression method. In Section 3, two extensions to the multipole theory in electrocardiology are presented: the convergence criterion with its mathematical proof, and the multipole shift equations for arbitrarily high multipole order. In Section 4 we describe the models for the torso volume conductor and for the heart current sources that are used in our study, in particular, to verify the multipole convergence criterion. In Section 5, we use the models for numerical simulations to obtain an idea of the rate of multipole convergence, and we present additional theoretical considerations. A final discussion is given in Section 6.

2.

I NTRODUCTION : T HEORETICAL

The principle of the multipole equivalent generator is based on the following series expansion of the potential 8 (Geselowitz, 1960): 8=

+∞ X l X

qlm 8lm

(2.1)

l=0 m=−l

In (2.1), 8lm denotes the potentials of unit multipole sources of the order lm, which are located in the multipole origin O. The superposition coefficients qlm are the multipole components which are defined as weighted integrals of the cardiac current source density distribution I : qlm =

Z

ylm I d V.

(2.2)

The functions ylm are harmonic polynomial functions of degree l, which we will specify later on. Considering a maximum degree L, corresponding to the highest order of the multipoles, the potential of the multipole equivalent generator is given by the partial sum of (2.1) with l ≤ L. We denote this partial sum as 6 L 8 and we denote the residual error as δ L 8 = 8 − 6 L 8. In Section 3.2, we will use the similar notations 6 L and δ L for series expansions of other functions (e.g., 6 L 9 and 6 L 8∞ ). In the regression method, the multipole components are obtained by fitting the potential of the multipole generator to the measured body potential. The multir pole components are now considered as unknown parameters qlm which have to be

Multipole Approach to the Inverse Problem in Electrocardiology

547

estimated: 6 rL 8 =

L X l X

r qlm 8lm .

l=0 m=−l

The index r is used, here and also further in the text, to refer to the regression method. The fitting of the potential is done by minimizing the sum square residual error between 6 rL 8 and the measured potential 80 : Res =

M X

6 rL 8(Pi )

− 8 (Pi ) 0

!2

,

(2.3)

i=1

where Pi , i = 1, . . . , M, are the positions of the points at which the potential is measured. The relationship between 80 and 8 is given by 80 = 8 + ν, where ν is the measurement noise. The multipole components estimated by the regression method are obtained by solving the set of linear equations 6 rL 8(Pi ) = 80 (Pi ), i = 1, . . . , M, in the least squares sense. Ideally, the right-hand sides 80 (Pi ) should equal 6 L 8(Pi ) but due to the noise and to the finite speed of convergence of the multipole equivalent generator the values 80 (Pi ) are mere approximations of 6 L 8(Pi ). The effect of the two disturbing factors is clarified by the following triangular inequality: r r 0 0 6 8 − 6 8 ≤ 6 8 − 8 + |8 − 8| + 6 8 − 8 L L L L = Res + |ν| + |δ L 8|

(2.4)

with the quadratic norm | · | defined as | f |2 =

M X ( f (Pi ))2 .

(2.5)

i=1

The regression method minimizes Res whereas actually 6 rL 8 − 6 L 8 should be minimized. In view of (2.5), the minimization of Res makes sense only if at the same time |ν| and |δ L 8| can be made small enough. For |δ L 8| this can be achieved by increasing the value of L, provided that the multipole equivalent generator converges. However, if it does not converge, then it is clear that a small residual Res does not provide any guarantee at all that the estimated multipole components are sufficiently accurate. Pilkington and Morrow (1982) observed that the regression method can lead to very large errors for the multipole components, while at the same time the obtained residuals are very small. Interestingly, there are cases where the multipole equivalent generator does not converge, but where at the same time it can be proven mathematically that the residual Res converges to zero for

548

P. De Muynck et al.

L → +∞ because certain conditions are satisfied regarding the arrangement of the current sources and the conductor inhomogeneities. This is described in Section 5.2. Even more interestingly, it appears from the same theoretical considerations of Section 5.2 that in these cases the discrepancy between the true multipole components and those obtained from the regression method can be large. In this paper, we investigate the convergence of the multipole equivalent generator, as a basic condition for the validity of the regression method. Hereby we do not pretend that the convergence is the only determining factor: even in the case of convergence, the presence of noise [i.e., the term |ν| in (2.4)] and the inherent instability of estimating the multipole components of a high order (which is a manifestation of the ill-posedness of the inverse problem in electrocardiology) impose a limitation on the number and the accuracy of the multipole components which can be estimated. The mentioned instability also implies that in practice the value of L cannot be chosen arbitrarily high so that in the analysis presented above, convergence has to be considered as obtaining small residual errors for moderate values of L (see Section 5.1). Beyond these moderate values, the behaviour of the residual error in function of L is not really relevant anymore, at least not in the context of the inverse problem in electrocardiology (see Section 5.1 for explanations).

3.

M ULTIPOLE T HEORY

3.1. The multipole lead field theory. The starting point for the multipole series expansion of the electrocardiological potential 8 is the reciprocity principle (Helmholtz, 1853), which expresses the measurement µ of an electrocardiographic lead in terms of the cardiac current source density distribution I and the lead field potential 9. We apply it for a single-electrode lead at a point P on the body surface St : µ = 8(P) =

Z

8Is d S =

Z

9 P I d V,

(3.1)

St

where Is denotes the weighting function of the lead, i.e., Is = δ P (δ P denotes the Dirac delta function around P), and 9 P is the R resulting lead field potential for this lead. By applying the integration operator I d V to the Taylor series expansion of 9 P around a suitable origin O with coordinates x, y, z, one obtains the series expansion of 8(P) in terms of the heart and lead tensors of Brody et al. (1961). In this series expansion, the effect of I on 8(P) isR described via the heart tensors (which are the cartesian moments of I of the form x i y j z k I d V ), whereas the lead tensors (the partial derivatives of 9 P at O) represent the effect of the volume conductor. The (l+2)(l+1) partial derivatives of the order l can be expressed in terms of 2 only 2l + 1 independent derivative operators ∂lm , m = −l, . . . , 1 (this only holds

Multipole Approach to the Inverse Problem in Electrocardiology

549

for harmonic arguments, which is of course the case for 9 P at O): √   ∂ m ∂ l−m 4π 2 − δm0 ∂ +i , ∂lm + i(1 − δm0 )∂l,−m = (−1) √ 2l + 1 (l − m)!(l + m)! ∂ x ∂y ∂z l−m (3.2) where δi j is the Kronecker delta. The Taylor series expansion of 9 P becomes as follows: +∞ X l X 9P = (∂lm 9 P )(O)ylm . (3.3) m

r

l=0 m=−l

In this expansion, ylm are the harmonic polynomials of degree l in x, y, z given by: ylm = Ylm (θ, φ)r I ,

(3.4)

where Ylm are the well-known surface spherical harmonics specified as a function of the spherical coordinates r, θ, φ: Ylm + i(1 − δm0 )Yl,−m

r p 2l + 1 = (−1) 2 − δm0 4π s (l − m)! m × P (cos θ)eimφ (l + m)! l m

(m ≥ 0).

R Then, by applying the integral over I d V to (3.3) and using (3.1), we obtain the multipole series expansion of 8(P): 8(P) =

+∞ X l X

qlm · 8lm (P).

(3.5)

l=0 m=−l

Indeed it is obvious that the integration of the ylm yields the value of the multipole components (2.2), whereas the derivatives in (3.3) can be identified as the unit multipole potentials 8lm at P: 8lm (P) = (∂lm 9 P )(O).

(3.6)

It can be verified by the reciprocity principle that the 8lm defined by (3.6) are the finite medium potentials of the unit multipole current sources located at O [the sources are the same as in Morse and Feshbach (1953, p. 1281)]. It is noteworthy that other derivative operators ∂lm than (3.2) could be proposed (any linear combination will do), which results in an alternative multipole series expansion (3.5) with multipole components that are different from qlm . In Appendix A it is proven that (3.2) is the correct expression for the ∂lm for our particular choice of the multipole components qlm [or, in other words, for our particular choice (3.4)

550

P. De Muynck et al.

for ylm ]. Our choice of the multipole components corresponds with the conventional normalization of the qlm (Jackson, 1975). In electrocardiology, the most oftenly used multipole moments are alm , blm (Geselowitz, 1960) which are normalized differently: they are related to qlm as alm = (−1)m flm qlm for m ≥ 0 and blm = (−1)m flm ql,−m for m > 0, with the factor flm defined as flm =

p

2 − δm0

r

4π 2l + 1

s

(l − m)! . (l + m)!

It is evident that (3.5) is the transformation of the heart-lead tensor series expansion of 8(P) into a non-redundant form. The existence of such a transformation was recognized soon after the elaboration of the lead field theory, and expressions for the non-redundant derivative operators ∂lm matching the multipole moments alm , blm were derived for l up to l = 4 (Brody et al., 1961; Arzbaecher and Brody, 1976). However, no general expression for arbitrary values of l such as (3.2) has been derived, and we assume that this is the reason why in the literature on lead field theory one has also failed to find general expressions for the multipole shift equations. Indeed, (3.2) provides the basis for deriving such equations in Section 3.3. 3.2. Multipole convergence criterion. Besides designing general multipole shift equations, the other reason for presenting the multipole series expansion as an application of the lead field theory is to gain insight into the conditions for convergence of the multipole series expansion (3.5); we can expect that (3.5) converges if the current sources I are fully contained in the region of convergence of (3.3) (if a function series expansion 6 L f converges pointwise over a certain domain to the limit function g, then the integration of the series over the domain results in a series which converges to the integral of g). The latter region consists of the largest sphere centred around O which does not intersect the conductivity inhomogeneities of the body conductor (Kellogg, 1929; Martensen, 1968), in other words, this is the largest sphere in which 9 is harmonic (the conductivity inhomogeneities are the interfaces between homogeneous subregions of the conductor: these include the interface of the body surface and the internal interfaces). We therefore suggest the following criterion for the convergence of the multipole series expansion (3.5): The multipole series expansion converges, if we can find a separating surface S = S R which has the shape of a spherical surface centred around O.

(3.7)

Here, we use the concept of a so-called separating surface S, which is defined as any closed surface that separates the current sources I (inside S) from the conductivity inhomogeneities (outside S). In Fig. 1, the criterion (3.7) is illustrated for the case of a piecewise homogeneous conductor: it is fulfilled for case b, but not for case a.

Multipole Approach to the Inverse Problem in Electrocardiology σ=0

(a) S

σ=0

(b) S = S R St

St +O σ1

551

+O

σ2

σ σ1 2

Figure 1. Inhomogeneous volume conductor consisting of two homogeneous subregions with conductivities σ1 and σ2 . The location of the cardiac current sources I is indicated by the gray region.

R Although the proof of (3.7) can be obtained from the I d V integration of (3.3), this approach is not satisfactory because we do not want to involve the entire spatial distribution of I in the proof, but only those aspects of I which determine the external potential. These aspects are exactly represented by the multipole components qlm , and also by the infinite medium potential of I . The information content of qlm and 8∞ is identical and the relationship between the two is given by the multipole series expansion of 8∞ : 8∞ =

+∞ X l X

qlm · φlm ,

(3.8)

l=0 m=−l

where φlm is the infinite medium potential of the unit multipole current sources located at O. The series (3.8) converges outside the sphere S R that is the smallest sphere centred around O and containing all the current sources I (Kellogg, 1929; Martensen, 1968). Equation (3.8) results from applying (3.5) for the case of an infinite homogeneous medium, i.e., for 9 P equal to the Green function G ∞ (r−r0 ) = (4π σ |r−r0 |)−1 (σ is the conductivity value, r and r0 are, respectively, the position vector of P and the variable position vector with respect to O). The application of (3.6) yields multipole infinite medium potentials agreeing with the normalization convention of Jackson (1975): φlm (r) = (−1)l ∂lm G ∞ (r) =

1 Ylm (θ, φ) · . (2l + 1)σ r l+1

(3.9)

To prove (3.7), we consider the following corollary of the reciprocity principle for the measurement µ of an electrocardiographic lead: µ=σ

Z  S

 ∂9 ∂8 8 − 9 d S. ∂n ∂n

(3.10)

For our application we identify µ = 8(P) and 9 = 9 P , but the proof below holds for the multipole series expansion of µ in general [i.e., equation (3.5) with 8(P)

552

P. De Muynck et al.

and φlm (P) replaced by µ and (∂lm 9)(O), respectively]. In (3.10), S denotes the separating surface (see above) and σ is the conductivity value within S. The righthand side of (3.10) is a symmetrical bilinear operator E operating on 8 and 9; µ = E(8, 9). Furthermore, it is known that 8 = 8∞ +8sec , where 8sec is the part of 8 caused by the secondary sources located at the conductivity inhomogeneities of the conductor. Considering that both 9 and 8sec are harmonic inside S, one has E(8sec , 9) = 0, so that: µ = E(8∞ , 9). (3.11) The dependence of µ on I is now exclusively represented via 8∞ , whereas 9 expresses the dependence on the body conductor and the lead. Moreover, the multipole series for µ can be related to the series expansions (3.3) for 9 and (3.8) for 8∞ by applying E to the respective series partial sums 6 L 9 and 6 L 8∞ . Via the Green theorem it is easy to prove the identity E(φlm , yl 0 m 0 ) = δll 0 δmm 0 , so that one obtains E 6 L 8∞ , 6 L 9

!

= 6 L µ.

The notation 6 L µ denotes the partial sum of the series (3.5). The residual error of this series can be related to the residual errors δ L 9 and δ L 8∞ on S as ! δ L µ = E(8∞ , 9) − E 6 L 8∞ , 6 L 9

!

= E 6 L 8∞ + δ L 8∞ , 6 L 9 + δ L 9 − E 6 L 8∞ , 6 L 9 = E(8∞ , δ L 9) + E(δ L 8∞ , 9) − E(δ L 8∞ , δ L 9). More explicitly, one has δL µ = σ

Z  S

8∞

∂ ∂8∞ (δ L 9) − δL 9 ∂n ∂n

+δ L 8∞

∂9 ∂ − (δ L 8∞ )9 ∂n ∂n

 ∂ ∂ −δ L 8∞ (δ L 9) + (δ L 8∞ )δ L 9 d S. ∂n ∂n

!

Multipole Approach to the Inverse Problem in Electrocardiology

553

From the Harnack theorem (Kellogg, 1929), it follows that the application of the normal derivative operator ∂/∂n on (3.3) yields a series expansion of ∂9/∂n which converges in the same region as the series for 9 itself (De Muynck, 1994). Furthermore, we can identify the residual error of this series as   ∂9 ∂ δL = (δ L 9). ∂n ∂n The same considerations apply for the normal derivative of the series expansion of 8∞ , so that finally the residual error δ L µ can be written as   Z  δL µ ∂9 ∂8∞ = 8∞ δ L − δL 9 σ ∂n ∂n S   ∂9 ∂8∞ +δ L 8∞ − δL 9 ∂n ∂n      ∂9 ∂8∞ + δL δ L 9 d S. (3.12) −δ L 8∞ δ L ∂n ∂n It is clear that, if both the series for 8∞ and 9 converge on S (and hence also the derived series for ∂9/∂n and ∂8∞ /∂n converge), then the presence of a factor δ L (. . .) in each term of (3.12) will make δ L µ converge to zero. Considering the regions of convergence of the above-mentioned series, we conclude that the surface S must lie between two concentric spheres centred around O; an inner sphere which contains all the current sources I and an outer sphere which is fully contained in one homogeneous subregion of the body conductor. Obviously, this is equivalent to the condition of the convergence criterion (3.7), which is hence proved. An upper bound for the residual error δ L µ can be obtained by the Schwarz inequality. Since the series expansions for ϕ = 9, 8∞ , ∂9/∂n, and ∂8∞ /∂n are uniformly convergent over S, they are also convergent ‘in the mean’ which means that the RMS error sZ 1 L (ϕ) = (δ L (ϕ))2 d S S

will converge to zero as L → +∞. We can now apply the Schwarz inequality to the first term in the right-hand side of (3.12): "Z

8∞ δ L S



 #2    Z ∂9 ∂9 2 d S ≤ 1L (8∞ )2 d S ∂r ∂r S 

≤ 1L



∂9 ∂r

2

|8∞ |2RMS . Surface(S)

Since 8∞ is a continuous function over S, this function will have a finite root mean square norm |8∞ |RMS . Thus, it appears that we have derived an upper bound

554

P. De Muynck et al.

which is convergent to zero as L → +∞. Similar upper bounds can be derived for the other terms of (3.12); the upper bound of δ L µ is obtained by adding all the individual upper bounds together. The series expansions for 9 and ∂9/∂n are only dependent on the behaviour of 9 on S R , while the series expansions for 8∞ and ∂8∞ /∂n are only dependent on the behaviour of 8∞ on S R (De Muynck, 1994). Thus, the rate of convergence of (3.5) can be exclusively related to the behaviour of 9 and 8∞ on S R , and, apparently, an increase on S R of the amplitude or the normal gradient magnitude of either of the potentials 9 and 8∞ will have a negative effect on the convergence. The factor 8∞ represents the effect of the cardiac current sources I , and the factor 9 represents the effects of the body conductor and of the measurement point position P. If S R comes close to the generating current sources of either 8∞ (the actual cardiac sources I ) or 9 (the secondary sources induced at the conductivity inhomogeneities by the lead in P), then the rate of convergence will decrease. Therefore, a good convergence can be expected if we are able to find a separating surface S R that is reasonably far from both the actual cardiac sources and the torso conductivity inhomogeneities. 3.3. Multipole shift equations. Since the introduction of the multipole theory in electrocardiology (Geselowitz, 1960), considerable attention has been given to the properties of the multipole components under a translation of the multipole origin. More particularly, one wants to obtain the linear equations which express the set of multipole components ql(1) (for the multipole origin O1 ) in terms of the set of 1m1 multipole components ql(2) (for an alternative multipole origin O2 ). Geselowitz 2m2 (1960) and Brody et al. (1973) derived these equations for a maximum order of l = 2 and l = 4, respectively. In Morse and Feshbach (1953), the general set of equations for arbitrary values of l is given for the special case of a translation in the direction of the z-axis. Below we derive the set of equations for arbitrary values of l and for an arbitrary translation direction. An approach to the derivation of the shift equations is based on the composition rule for the derivative operators ∂lm , which expresses the operator ∂l 0 m 0 ∂lm as a linear combination of operators ∂l ∗ m ∗ ∂l 0 m 0 ∂lm = α

Fl 0 m 0 Flm Fl 0 m 0 Flm ∂l 0 +l,m α + β ∂l 0 +l,m β + R(∇)∇ 2 , Fl 0 +l,m α Fl 0 +l,m β

(3.13)

where Flm =

r

4π 1 √ 2l + 1 (l − m)!(l + m)!

and R is a homogeneous polynomial of degree l − 2, which is of no further interest here because we will use (3.13) only for harmonic arguments. The quantities α, β, m α and m β are functions of m 0 and m only; they are derived in Appendix A.

Multipole Approach to the Inverse Problem in Electrocardiology

555

z1

z2 O1

y1 a21

r1

x1

P

y2

O2

r2 x2

Figure 2. Multipole origins O1 ad O2 , with their respective frames of reference.

For the derivation of the shift equations, we assume that the frames of reference O1 x1 y1 z 1 and O2 x2 y2 z 2 are oriented in an identical way; this is shown in Fig. 2. From (2.2), it is evident that the translation properties of qlm and ylm are identical; (1) i.e., we need to express yl(1) in terms of the functions yl(2) with ylm ≡ ylm (r1 ) and 1m1 2m2 (2) ylm ≡ ylm (r2 ) and given r1 = r2 + a21 . To this end, we consider the Taylor series expansion (3.3) of yl(1) around O2 : 1m1 yl(1) 1m1

=

l1 l2 X X

∂l2 m 2 yl1 m 1 (a21 )yl(2) . 2m2

(3.14)

l2 =0 m 2 =−l2

It is straightforward to verify that the superposition coefficients in (3.14) are indeed identical to (∂l2 m 2 yl(1) )(O2 ). We see that the multipole shift equations make 1m1 up a set of linear equations and that the derivatives in (3.14) are the equation coefficients. To further elaborate these coefficients, we consider the multipole series for the infinite medium Green function (Jackson, 1975): +∞ X l X

G ∞ (r − r0 ) =

ylm (r0 )φlm (r).

(3.15)

l=0 m=−l

There are two ways to look at the series (3.15): it is the Taylor series of the potential 9 P (r0 ) of a monopole source at r, whereas it can also be considered as the multipole series expansion of the potential 8∞ (r) of a monopole source at r0 (both 8∞ and 9 P are infinite medium potentials). The series is convergent for r 0 < r . In accordance with (3.9), the application of the derivative operator (−1)l ∂l2 m 2 (r) to the left-hand side of (3.15) yields φl2 m 2 (r − r0 ). Furthermore, we can replace the derivative operator (−1)l ∂l2 m 2 (r) by ∂l2 m 2 (r0 ), because of the dependency in r − r0 . Finally, we obtain φl2 m 2 (r − r ) = ∂l2 m 2 (r )G ∞ (r − r ) = 0

0

0

l1 +∞ X X

l1 =0 m 1 =−l1

∂l2 m 2 yl1 m 1 (r0 )φl1 m 1 (r).

(3.16)

556

P. De Muynck et al.

Obviously, the superposition coefficients of the function φl1 m 1 in (3.16) are exactly the coefficients of the shift equations that we are looking for. To identify these factors, we apply (−1)l ∂l2 m 2 (r) to the right-hand side of (3.15), so that the following alternative expression for the series (3.16) is obtained: φl2 m 2 (r − r0 ) = (−1)l2 ∂l2 m 2 (r)G ∞ (r − r0 ) =

+∞ X I X

ylm (r0 )(−1)l2 ∂l2 m 2 φlm (r)

+∞ X l X

ylm (r0 )(−1)l+l2 ∂l2 m 2 ∂lm G ∞ (r).

l=0 m=−l

=

(3.17)

l=0 m=−l

We can use (3.13) and (3.9) to express the terms of (3.17) in terms of the functions φlm as follows:  α l+l2 l+l2 (−1) ∂l2 m 2 ∂lm G ∞ (r) = (−1) Fl2 m 2 Flm ∂l +l,m α G ∞ (r) Fl2 +l,m α 2  β + ∂l +l,m β G ∞ (r) Fl2 +l,m β 2   α β = Fl2 m 2 Flm φl +l,m α (r) + φl +l,m β (r) . Fl2 +l,m α 2 Fl2 +l,m β 2 The coefficients ∂l2 m 2 yl1 m 1 (r0 ) now follow from the identification of (3.17) with (3.16): ∂l2 m 2 yl1 m 1 =

(lX 1 −l2 )

δm 1 m α α

m=−(l1 −l2 )

+δm 1 m β β

Fl2 m 2 Fl1 −l2 ,m yl1 −l2 ,m Fl1 m 1

Fl2 m 2 Fll −l2 ,m yl1 −l2 ,m . Fl1 m 1

(3.18)

In conclusion, the shift equations for the multipole components are given by: ql(1) 1m1

=

l1 l2 X X

∂l2 m 2 yl1 m 1 (a21 ) · ql(2) . 2m2

(3.19)

l2 =0 m 2 =−l2

The coefficients of (3.19) are expressed in terms of the harmonic polynomials yl(1) (a21 ) via (3.18), so that they are suitable for numerical computation. 1m1 A computer algorithm for equation (3.19) has been implemented and verified. For the verification, we considered a number of random unit multipole sources l1 m 1 , . . . , ln m n located at the origins O1 , . . . , On . The computer algorithm was

Multipole Approach to the Inverse Problem in Electrocardiology

557

1 2 4 1 4

3 2 5 6 left

right

back

front

Figure 3. The body conductor model is shown via a horizontal cross-section at the left and a vertical cross-section with sagittal view at the right (the position of the intersecting planes is indicated). The numbers i = 1, . . . , 6 indicate the boundary surfaces si which divide the model into homogeneous subregions: torso surface (1), muscle layer (2), lungs (3 and 4), and intracavitary bloodmasses (5 and 6). The original torso surface (before it is extended outside the physical boundaries) is shown by a dashed line.

used to compute the multipole components of the sources with respect to a single multipole origin O. Using the obtained multipole components, we consequently computed the infinite medium potential of the multipole sources at various test points using the multipole series expansion (3.8). The resulting potentials of the multipole series expansion were compared with the infinite medium potentials obtained by the direct computation. The comparison showed that the multipole expansions converge to the directly computed values, at all the test points that were chosen.

4.

M ODELS AND M ETHODS

4.1. Modelling of the body conductor and cardiac current sources. In our approach, existing models are taken as a starting point, and modifications have been made in order to take into account all aspects relating to the convergence of the multipole equivalent generator. The original model data were acquired in the Medical Physics & Biophysics Department of the University of Nijmegen from magnetic resonance images of a human subject and from the consequent computation of the isochrones of ventricular depolarization on the basis of the measured body surface potential (Cuppen and van Oosterom, 1984). The isochrones constitute the model data for the cardiac current sources. The model for the body conductor is a piecewise homogeneous medium, which consists of the following homogeneous subregions: the lungs, the intracavitary

558

P. De Muynck et al.

bloodmasses of the heart, and the remaining part of the torso which includes the heart muscle (Fig. 3). To these subregions the following conductivity values are assigned: 0.044, 0.66, and 0.22 −1 · m−1 , respectively (Cuppen and van Oosterom, 1984; Rush et al., 1963). This model was modified by the addition of a skeletal muscle layer. We take this additional inhomogeneity into account, because the layer is close to the heart current sources; it obviously plays a role in the multipole convergence criterion, and, probably, it affects the rate of convergence of the multipole generator. The anisotropic conductivity properties of the muscle layer are taken into account by dilating the layer with a factor 1.94 (McFee and Rush, 1968; Cornelis, 1980; Gulrajani and Mailloux, 1983). This has two consequences: (1) the muscle layer is now considered to have an isotropic conductivity of 0.155 −1 · m−1 and (2) the boundaries of the conductor model extend outside the physical boundaries of the torso. The operation performed here is different from Gulrajani and Mailloux (1983) in two ways. First, whereas in Gulrajani and Mailloux (1983) the thickness of the layer is constant (1 cm), we use a variable thickness between 0.8 and 3 cm. The thickness values are derived on the basis of anatomical observations (e.g., the layer is much thicker at the back of the torso than at the front). Secondly, we use the more recent conductivity values of Epstein and Foster (1983) for the muscle layer, namely 0.52 and 0.08 −1 · m−1 for the tangential and normal conductivity, respectively. This results in different values for the isotropic conductivity and for the dilation factor [in Gulrajani and Mailloux (1983) these are 0.125 −1 · m−1 and 3, respectively]. For the cardiac current sources the uniform double-layer model is used. The model data comprises the geometry of the endocardial and epicardial surfaces and the distribution of the depolarization times on these surfaces. The latter were computed in the University of Nijmegen, using the procedure described by Cuppen and van Oosterom (1984); the procedure is based on the measurement of the body surface potentials in a set of 64 electrodes throughout the QRS complex. Evidence for the physiological correctness of the model data obtained in this way is given by Cuppen and van Oosterom (1984). It is evident that the presence of the intracavitary bloodmasses is in total contradiction with the conditions for the multipole convergence criterion. In fact, the bloodmasses are prevalent inhomogeneities within the body conductor (Brody, 1956), but at the same time it is impossible to separate them from the current sources. Therefore we are forced to consider the cardiac current sources as a combination of the uniform double layer itself and of secondary current sources which are located on the boundary surfaces of the bloodmasses. The study of the multipole convergence is consequently carried out on the modified conductor model, from which the bloodmass inhomogeneities are omitted. The secondary current sources (due to the conductivity differences between bloodmass regions and heart muscle) are described by their infinite medium potential which can be expressed as

Multipole Approach to the Inverse Problem in Electrocardiology

follows in a given point r outside the bloodmasses: 6 Z X 1 (σ− − σ+ )8(r0 )dr (r0 ). 4π σ (r) i=5 Si

559

(4.1)

The symbols σ− and σ+ denote the conductivity values inside and outside the boundary surfaces si , and dr (r0 ) is an elementary solid angle around r0 on the surface si seen from an observer at r. Furthermore, σ (r) is the conductivity value at the point r (which is equal to the mean of σ− and σ+ if r is on the bounding surface). Since (4.1) is a function of 8 on the boundary surfaces s5 and s6 of the bloodmasses, it is evident that our cardiac current sources are no longer solely dependent on the electric activity of the heart, but also on the conductivity properties of the conductor. The same conclusion applies to the multipole components. However, it appears that this dependency is almost exclusively defined in the geometry of the bloodmasses (this will be explained below), and therefore we can still assume that the multipole components are only dependent on the ‘heart’ itself, that is, on its electric activity and its internal geometry, but not on the surrounding body conductor. The finite medium potentials on the body conductor were computed via the boundary element method (BEM) (Barr et al., 1966; Lynn and Timlake, 1968a,b; Barr et al., 1977; Meijs et al., 1989), which is based on the discretization of the following integral equation (Barnard et al., 1967): p Z X 1 8(r) = 8∞ (r) + (σ− − σ+ )8(r0 )dr (r0 ). (4.2) 4π σ (r) i=1 Si The BEM computes 8 as a function of the infinite medium potential 8∞ (the uniform double-layer potential), which has to be specified in the nodes of the triangulated conductor boundary surfaces si , i = 1, . . . , p. Thus, 8 was computed for eight different time instants of the QRS complex (26, 38, 44, 48, 52, 60, 70, and 78 ms), starting from the computed values of 8∞ in the nodes. It should be noted that, while for the solution of the integral equations (4.2) only the finite medium potentials 8 on the six boundary surfaces are involved, the identity (4.2) furthermore constitutes a valid expression for 8(r) in any arbitrary point r outside the primary source region. Moreover, the expression for 8(r) is decomposed into the contribution of the primary current sources and of the secondary current sources due to the six boundary surfaces, as if they were all located in one infinite homogeneous medium. The secondary current sources are expressed as non-uniform double-layers on the boundary surfaces, where the double layer amplitude varies proportionally to 8 over the boundary surface. It is obvious now that the addition to the primary current sources of the secondary current sources due to the bloodmasses is accomplished by adding the two surface integrals over s5 and s6 of (4.2) to the infinite medium potential of the uniform double layer, i.e., the new infinite

560

P. De Muynck et al.

medium potential is the old infinite medium potential plus (4.1). From this point on, the notation 8∞ will refer to the new infinite medium potential, and qlm will refer to the new cardiac multipole components. Correspondingly, we will henceforth work with the modified conductor model without the bloodmasses, noting that we needed the original model only to compute 8 on the boundary surfaces and subsequently derive the new ‘integrated’ 8∞ , with two important remarks to be made: — The value of 8∞ can be computed in any arbitrary point r because 8∞ is the infinite medium potential of known double-layer sources. Hence it is very easy to compute the values of the multipole components qlm by computing 8∞ on a large enough sphere around the multipole origin, and consequently performing a regression with the functions φlm [see (3.8)]. — It is clear that the dependence of 8∞ and hence qlm on the conductivity properties of the conductor is expressed in the values of 8(r) over the surfaces s5 and s6 [via (4.1)], whereas the influence of the remaining (non-bloodmass) inhomogeneities on 8(r) can be estimated in first order via (4.2), namely by considering the relative contribution of the combined integrals p = 1, . . . , 4 to the total expression (4.2). It has appeared from our computations, however, that for r ∈ s5 and r ∈ s6 the relative contribution was consistently very small, hence our assertion that the dependency of qlm on the torso conductivity properties is almost exclusively defined in the geometry of the bloodmasses. Finally, the finite medium potentials 8lm are computed, also via the BEM [we repeat that the surfaces s5 and s6 have been removed from the conductor now; i.e., p = 4 in (4.2)]. The problem with high values of l is that the spatial frequency of the potentials φlm on the surfaces si can become quite elevated, especially for those parts of the surfaces which are closest to the multipole origin (in our case, the maximum value of l is 10). Special measures have been taken to guarantee that the computation of 8lm remains accurate when the order l increases. Therefore, we used an implementation of the BEM which includes the ‘auto solid angle enhancement’ described in (Meijs et al., 1989). In order to further increase the accuracy, we have applied an iterative refinement of the triangulation of the surfaces si . The goal of having smaller triangles is to guarantee that the potentials show a more or less linear behaviour within the triangles [the underlying assumption for the validity of BEM (Barr et al., 1977; Meijs et al., 1989)]. The triangles are subdivided into smaller triangles (2, 3, or 4) until the sides of the triangles are less than a specified critical size. Moreover, this critical size is not a constant over the surfaces, but is defined so that a selective refinement is obtained, the triangles becoming smaller when one approaches the multipole origin. In the refined triangulation, a total of 12978 triangles and 6497 nodes were used for the combined surfaces of the lungs, the skeletal muscle, and the torso boundary.

Multipole Approach to the Inverse Problem in Electrocardiology g

f (ii)

(iii)

561

k (ii)

(iii)

(ii)

(iii)

(i)

(ii)

(iii)

Figure 4. Position of the sphere S R with a radius of 9 cm, for three different origins f , g, and k. The position of S R is depicted in three cross-sections (i), (ii), and (iii). The crosssections (i) are shown in the top row and correspond to a frontal view. The planes of (ii) and (iii) are perpendicular to the plane of (i); their intersections with (i) are shown in the top row. The second and the third row show the cross-sections (ii) and (iii), respectively.

4.2. Application of the convergence criterion. For the verification of the convergence criterion, we tried to find a suitable separating surface S R , i.e., a spherical surface which fits in between the current sources and the body inhomogeneities. If such a surface S R exists, then its centre will be a valid choice for the position of the multipole origin. In our attempts to find S R , we tried various values for the radius R of the sphere and for the position of its centre. From our attempts, it became clear that the proposed separating surface S R does not exist for our body conductor model. To illustrate this conclusion, we depict in Fig. 4 three different attempts for positioning S R . Here, the current source region comprises the two intracavitary bloodmasses in conjunction with the myocardium (indicated in dashed lines). The atrial part of the heart muscle is omitted in the model because, during QRS, depolarization only takes place in the ventricular part. The inhomogeneities closest to the sources are the lungs and the precordial part of the skeletal muscle layer. Figure 4 clearly shows that it is impossible to keep these inhomogeneities out of S R , if S R is at the same time assumed to contain the entire current source region. Therefore, we conclude that the separating surface S R does not exist, and that the condition of the convergence criterion is not fulfilled.

562

P. De Muynck et al.

This observation does not give a clear indication of whether the multipole series is totally divergent or there are merely problems with multipole convergence— which may only be of minor importance for the practical application of the multipole model. An indication in support of the latter case is that the condition of the criterion (3.7) is only slightly violated: for the attempts f and k, the sphere S R has all the sources inside, while only a small fraction of the inhomogeneities are not outside S R . Also, the previous numerical simulations of the multipole approach on an inhomogeneous body conductor model by Dub´e et al. (1985) did not indicate major convergence problems (note, however, that no muscle layer was included in the model used by Dub´e et al.) 4.3. Choice of the multiple origin. The multipole convergence criterion failing, we shall turn in Section 5.1 to numerical simulations to find out how the multipole equivalent generator potential behaves on the body surface. For this purpose we first have to specify the location of the multipole origin; we will choose O to be optimal for the multipole approach to the inverse electrocardiological problem (in which the potential is estimated on a spherical quasi-epicardial surface S R , see Section 1). Hence, we now turn to the problem of finding the optimal surface S R , whereas the resulting choice for O is nothing else than the centre of S R . We consider the Lth order approximation of 8∞ on a sphere S R with radius R: 8∞;e ≈

L X l X Q lm Ylm (θ, ϕ), R l+1 i=0 m=−1

(4.3)

where 8∞;e denotes the infinite medium potential on the quasi-epicardial surface, i.e., 8∞ on S R (the index e refers to ‘epicardial’), and Q lm =

1 qlm . (2l + 1)σ

The right-hand side of (4.3) is the multipole series expansion (3.8) of 8∞ , which has been rewritten in terms of the surface spherical harmonics Ylm on S R . In the multipole approach, the estimated multipole components are used to obtain an estimate (4.3) of the quasi-epicardial potential, which is more meaningful diagnostically than the values of the multipole components by themselves (Damen, 1980; Tanaka et al., 1980; Pilkington and Morrow, 1982; Tanaka et al., 1982). An interesting aspect of the multipole approach presented here is that the estimated quasiepicardial potential is intrinsically only dependent on the electric activity of the heart itself because (4.3) depends on the qlm only. Our choice of the quasi-epicardium S R will be based on two requirements. The first one is that the multipole series expansion (3.8) must quickly converge on S R , so that (4.3) yields a good approximation for moderate values of L. The second requirement is that the quasi-epicardium should be as close to the heart as is

Multipole Approach to the Inverse Problem in Electrocardiology 1

1 t = 26 t = 38 t = 44 t = 48 t = 52 t = 60 t = 70 t = 78

O=f

0.9 0.8 0.7 0.6

0.8 0.7 0.6 0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1 0

2

4

6

t = 26 t = 38 t = 44 t = 48 t = 52 t = 60 t = 70 t = 78

O=g

0.9

0.5

0

563

8

10

0

0

2

4

6

8

10

1 t = 26 t = 38 t = 44 t = 48 t = 52 t = 60 t = 70 t = 78

O=k

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

2

4

6

8

10

Figure 5. Relative RMS residual error of the multipole series expansion of 8∞ over S R for the eight selected time instants as a function of L.

possible, so that the obtained potential is diagnostically relevant. To find an appropriate quasi-epicardium, we computed the residual error of (4.3) for various radii and positions of S R with L = 10. We found that R = 9 cm was the smallest radius for which small enough residuals could be obtained. The radius of the quasi-epicardium was hence chosen as R = 9 cm, and its position was selected on the basis of the smallest residuals. In Fig. 5, we show the computed residuals for three different positions of S R . The three spheres S R coincide with the tentative separating spheres which were previously shown in Fig. 4 and they are respectively labelled f , g, and k. The actual choice that we made for the quasi-epicardium is the sphere k, which is situated in between the spheres f and g. The multipole origin is consequently at the centre of the sphere k. In Fig. 5, large residual errors show up when S R is close to the regions of the heart which are electrically active at the designated time. In the case of S R = f , for instance, S R is close to the apex, and large residual errors are observed for the times t = 44, 48, 52. These are exactly the time instants when the depolarization wave front is close to the apex. The errors disappear for S R = g, because S R is moved away from the apex; on the other hand, large residuals now show up for other time instants (mainly at the end of QRS). This must be ascribed

564

P. De Muynck et al.

to the fact that S R is now closer to the intracavitary bloodmasses, which are located in the heart at the opposite side of the apex. The sphere k, finally, is positioned in between f and g; with this choice a quasi-epicardium is obtained which completely contains all the heart current sources and gives optimal convergence over the eight selected times instants for the multipole series expansion of 8∞ . From the case k of Fig. 5 it is clear that for L = 10 the multipole series expansion (4.3) yields a good approximation of 8∞ , with residuals of up to 10% at the most. This means that (4.3) is a valid way to represent the quasi-epicardial potential, using the multipole components. In the inverse problem of electrocardiology, the problem is not with the representation of 8∞ , but with the fact that only the lower order multipole components can be estimated in a reliable way (these orders are actually much lower than 10). In order to use (4.3) for the estimation of 8∞ , the series consequently has to be truncated at a lower value of L. This truncation of the higher order terms corresponds to omitting components of high spatial frequency in the reconstruction of 8∞ [the suppression of high spatial frequency components is typical of various approaches to estimating epicardial potentials, see Rudy and Messinger-Rapport (1988)].

5.

C ONVERGENCE R ATE OF THE M ULTIPOLE E QUIVALENT G ENERATOR

5.1. Numerical simulations. To investigate the convergence of the multipole equivalent generator for the proposed realistic models of the cardiac electric generator and the body conductor (i.e., a model taking into account the lung and muscle layer inhomogeneities), the torso residual errors of the generator potential were computed numerically. The results are shown in Fig. 6. We note that the multipole origin is chosen at the point k (see Section 4.3). The RMS errors of Fig. 6 are computed over the nodes of the triangulated torso surface, so that the approximation errors in the precordial area are especially stressed (the distribution of nodes is more dense in the precordial area than in the other areas of the body surface). This makes sense, because in this area the signal content of the electrocardiographic potential is highest and the majority of the electrocardiographic leads are placed here. The multipole equivalent generator only really converges for two of the eight time instants in the QRS interval: for the other six time instants we see a stagnation of the residual error between L = 5 and L = 8, whereas for the higher values of L the error again increases. The non-convergence is not surprising because the multipole convergence criterion is not fulfilled. However, the residual errors for L = 5 are sufficiently small, so that we can consider the multipole equivalent generator to be convergent for practical applications such as the regression method. Indeed, the behaviour of the residual error for L > 5 is not really relevant for the regression method, because due to the ill-posedness of the inverse problem it does not make sense to consider values of L which are higher than 5 anyway. It is more important to have low residual errors for low values of L than to have

Multipole Approach to the Inverse Problem in Electrocardiology

565

0.9 t =26 t =38 t =44 t =48 t =52 t =60 t =70 t =78

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

1

2

3

4

5

6

7

8

9

10

Figure 6. Relative RMS residual error of the multipole equivalent generator potential |δ L 8|/|8| on the torso model surface for the eight selected time instants as a function of the order L of the multipole approximation.

a multipole series which is convergent in the strict sense, but which nevertheless has large residuals for low values of L. In addition, it is evident that at a certain point, increasing the value of L will decrease the accuracy of the regression method instead of increasing it, even in the case that the multipole equivalent generator converges monotonically δ L 8 → 0. The reason for this is that at a certain point the regression linear equations [i.e., 6 rL 8(Pi ) = 80 (Pi ), see Section 2] become nearly singular when 8lm of higher order are added (this is implied by the ill-posedness of the problem). Even if the series residual δ L further decreases for increasing L, the presence of the noise term |ν| in (2.4) implies large errors for the estimated qlm as soon as high orders of L become involved. The spatial distribution of the residual error δ L 8 has been inspected, and it is found that the error is mainly concentrated in the precordial area of the torso, especially for higher values of L. Apart from this, it is difficult to recognize any other noteworthy features of the spatial distribution of the error. For higher values of L, we observe that the error consists of patterns in the precordial area; without having a clear explanation for this, we observe that these patterns can be either of a low or a high spatial frequency. For the time instants t = 48 and t = 60, we also inspected the distribution of the residual error on the internal surfaces of the torso model, i.e., the surfaces of the lungs and muscle layers. The results are presented in Fig. 7. Due to the proximity of the left lung to the current sources and multipole origin (see Fig. 4), it is not surprising that the multipole equivalent generator is completely divergent on the left lung surface. On the other surfaces, we see that the residual error initially decreases as a function of L, but after L = 5 there is first a stagnation of the error and then an increase of the error (the fact that this is more

566

P. De Muynck et al.

0.9

0.9

0.8

t =48

0.7

left lung right lung muscle layer body surface

0.8 0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

1

2

3

4

5

6

7

8

t =60

9

10

0

1

2

3

4

5

6

left lung right lung muscle layer body surface

7

8

9

10

Figure 7. Relative RMS residual error of the multipole equivalent generator potential |δ L 8|/|8| on the internal surfaces of the torso model for the time instants t = 48 and t = 60 as a function of L.

evident for t = 48 than for t = 60 is in agreement with the graphs in Fig. 6 for these two time instants). As to the spatial distribution of the residual error on the internal surfaces, it was observed that the error was most prominent in the regions that are closest to the cardiac current sources. This observation also applies for the potential on the body surface, i.e., there is a concentration of the error in the precordial area. The significance of the convergence criterion (3.7) of Section 3.2 is further illustrated by an additional numerical simulation, where the residual error was computed for a model with the internal inhomogeneities removed from the torso. The cardiac current sources remain the same (including the bloodmass secondary current sources), but we now use a homogeneous torso conductor with the same external shape as the inhomogeneous torso of Fig. 3. The results obtained for the errors are shown in Fig. 8. Obviously, the conditions of the convergence criterion are now satisfied because of the omission of the internal surfaces. We see, accordingly, that the residual errors now decrease monotonically; from the graphs in logarithmic scale (not included here), there is evidence that the errors decrease exponentially as a function of L. In contrast to the inhomogeneous case, the multipole equivalent generator is now fully convergent for all the time instants during QRS. Up to L = 4, the magnitude of the errors is comparable with the inhomogeneous case; however, now the error does not stagnate but continues to decrease after L = 5. We therefore conclude that for our realistic torso model the problems with the equivalent generator convergence are due to the presence of the internal inhomogeneities. The significance of the convergence criterion lies in the fact that it gives an a priori indication that there is convergence for the homogeneous case, whereas there might be complications for the inhomogeneous case. In conclusion, the numerical simulations that were carried out show that the convergence problems are not so severe that we must reject the multipole approach to the inverse problem a priori. Furthermore, some additional numerical simulations

Multipole Approach to the Inverse Problem in Electrocardiology

567

0.7 t =26 t =38 t =44 t =48 t =52 t =60 t =70 t =78

0.6 0.5 0.4 0.3 0.2 0.1 0

1

2

3

4

5

6

7

8

9

10

Figure 8. Relative RMS residual error of the multipole equivalent generator potential |6 L 8 − 8|/|8| on the homogeneous torso model surface for the eight selected time instants as a function of L.

were also carried out in order to study the accuracy of the regression method for the estimation of multipole components, and we encountered no problems related to the convergence. In order to do this, three different lead systems were defined, one called ‘U64’ with 64 electrodes which are spread in a uniform way over the torso surface [this is similar to the lead system used by Dub´e et al. (1985)], and two lead systems ‘N64’ and ‘N32’ (De Muynck, 1994) with 64, respectively 32 electrodes, which are more concentrated in the precordial area (i.e., similar to what has been conventionally used in electrocardiology for the purpose of body surface potential mapping). In the simulations, a different noise pattern was added to each of the eight computed potentials 8 (t = 26, . . . , 78), where the patterns were generated randomly as uncorrelated white noise with a standard deviation of 20 µV. For each different degree l, the RMS error on the estimated multipole moments was computed (the mean is taken over the eight selected times and over the 2l + 1 different values of m), and the result is divided by the root mean square value of the multipole moments of that degree, in order to obtain relative errors: v u P78 Pl r 2 u m=−1 (qlm − qlm ) R E(l) = t t=26 . P78 Pl 2 t=26 m=−1 qlm In Fig. 9, the estimation errors R E(l) of the regression method are plotted as a function of L, i.e., the maximal degree of the fitting function 6 L 8 as defined in Section 2. The optimal point L min (l), i.e., the value of L where the error becomes minimum, is obviously dependent on l. Comparison of the curves U64 and N64 in Fig. 9 points out that a uniform electrode placement leads to better results, while the comparison of N32 and N64 indicates that the addition of the 32 extra elec-

568

P. De Muynck et al.

1

2 l =1

l =2

0.5

1

0.2

0.5

0.1

0.2

0.05

N32 N64 U64

0.02

N32 N64 U64

0.1 0.05

1

2

3 L

4

5

2

3

4 L

5

6

10 l =3 5 2 1 0.5

N32 N64 U64

0.2 3

4

5

6

L Figure 9. Relative root mean square error R E(l) on estimated multipole moments for l = 1, 2, 3, using the lead systems N64, N32, and U64.

trodes improves accuracy. For the best case U64 and L = L min (l), the regression method estimates the multipole moments qlm with an average relative error of 4%, 15%, 45% for a degree of l = 1, 2, 3, respectively (note that the value of L min has been determined a posteriori here; in practice, the value of L min will not be available). For l > 3 the estimates become inaccurate and unstable. However, the observation which is important in the context of the convergence considerations is the following: the estimation errors decrease in function of L at first and attain a minimum in L min , then, beyond L = L min , they immediately start to increase at a fast rate. This could be attributed to two reasons: (1) the multipole series residual δ L 8 increases after L = 6 (see Fig. 6) and (2) the addition of higher degree multipole terms increasingly reveals the ill-posedness of the electrocardiological inverse problem. However, since L min is typically much smaller than 6, we can safely assume that (2) is the real reason that explains this observation, whereas the non-ideal convergence properties of the multipole equivalent generator are not really significant in the context of regression. Our final conclusion, therefore, is that the multipole equivalent generator can be considered convergent in the context of practical applications.

Multipole Approach to the Inverse Problem in Electrocardiology

569

5.2. Further theoretical considerations. Whereas in the previous paragraph the simulations with the realistic models have reassured us on the use of regression for practical applications, this paragraph describes a theoretical construction showing that, in the general case, situations can occur where small residuals are not necessarily a reliable indication for the accuracy of the estimated multipole components. We note, however, that in our simulations with realistic torso models these situations did not occur, so the conclusions of this paragraph do not contradict the practical results of Section 5.1. The use of the multipole potentials as basis functions for the interpolation of the body surface potential is discussed by Pilkington and Morrow (1982). Such basis functions will be efficient if the multipole convergence is sufficiently fast. The motivation for choosing the multipole potentials, rather than empirically obtained sets of basis functions, is that the sought for parameters pretend to have a physical meaning: if we interpolate on the basis of a minimum mean square residual, then the interpolation parameters will be the multipole components—as estimated by the regression method. However, from our discussion in Section 2, it is clear that the physical meaning of these parameters is lost if the multipole equivalent generator fails to converge (even when the interpolation yields a good approximation to the body surface potential). A strong indication that this might actually happen is that in certain situations, where the multipole equivalent generator fails to converge, it may still be possible to approximate the electrocardiological potential 8 via a superposition of multipoles up to an arbitrarily small error ε. Indeed, a ‘multipole approximation criterion’ can be derived which guarantees the existence of such approximations under conditions which are less restrictive than the conditions of the multipole convergence criterion. If the conditions of this criterion are fulfilled, there are two possible cases: — The multipole equivalent generator converges, hence a valid choice for multipole approximations are obviously the partial sums 6 L 8 of the multipole equivalent generator: the functions 6 L 8 for L = 1, 2, . . . form a sequence of functions, which converges to 8 for L → +∞. — The multipole equivalent generator does not converge, however multipole approximations can still be found, of the following form: 6 Le 8 =

L X l X

e qlm (L)8lm .

(5.1)

l=1 m=−l

We refer to the approximations 6 Le 8 as ‘generalized multipole approximations’, where the index e is used to distinguish them from the multipole equivalent gene erator approximations 6 L 8 (similarly qlm is introduced as a notation). It is exe plicitly indicated that the superposition coefficients qlm are now dependent on L. Although it is still possible in the second case to construct a sequence of functions

570

P. De Muynck et al.

e 6 Le 8 which converges to 8 for L → +∞, it will require that the values of qlm e change in between consecutive approximations. Having qlm independent of L is only possible if the multipole equivalent generator converges [and in that case one e will necessarily have qlm = qlm ; see De Muynck (1994)]. The multipole approximation criterion is declared as follows:

If the cardiac current sources and the multipole origin are situated in one and the same connected homogeneous subregion V of the torso conductor, then for any arbitrarily small ε > 0 a superposition 6 Le 8 can be found which approximates the electrocardiological potential 8 with an error that is less than ε everywhere outside V

(5.2)

The proof of (5.2) is given in Appendix B. If the conditions of the criterion are satisfied, then a sequence 6 Le 8 (generalized multipole approximation) which converges to 8 for L → +∞ can be constructed. Although (5.2) can somehow be considered as an extension of the multipole convergence criterion, it is not a positive indication with regard to the estimation of the multipole components. In the e case that 6 Le 8 converges to 8 and 6 L 8 does not, the values of qlm and qlm will differ (from the considerations in Appendix B, this appears to be the case, except for the orders l which are relatively low compared with L). Furthermore, it is guaranteed on the basis of the convergence of 6 Le 8 that the interpolation of the body surface potential with multipole potentials works and that small residuals can be obtained. Based on a similar reasoning as in Section 2, where 6 Le 8 now takes over the role of 6 L 8, we may assume that the estimated multipole components e r qlm will approach qlm as L increases (no noise is assumed). Thus we find that the estimated multipole components do not approach the right values, because of the e above-mentioned discrepancies between qlm and qlm . We conclude, therefore, that under the described conditions [(5.2) is satisfied, but (3.7) is not] the regression method may yield rather small residuals, whereas the estimation of the multipole components is inaccurate. This brings evidence to our assertion in Section 2 that the value of the residual Res is not a reliable indication for the accuracy of the r estimated multipole components qlm . As to the use of the multipole functions for interpolation, it appears to us that if the multipole convergence is not assured, then there is no particular reason to prefer these functions over empirically obtained interpolation functions such as those mentioned by Pilkington and Morrow (1982).

Multipole Approach to the Inverse Problem in Electrocardiology

6.

571

S UMMARY AND C ONCLUSIONS

Extensions to the existing multipole theory are needed to understand the convergence properties of the multipole generators in finite media. We elaborate the multipole theory specifically for the case of the inhomogeneous body conductor, by involving the concept of the lead field potential and the related lead tensors and components. The primary result of our approach is the convergence criterion for the multipole equivalent generator in the inhomogeneous body conductor. On the basis of this theoretical result, it appears that the proximity of the internal inhomogeneities (lungs and skeletal muscle layer) to the cardiac electric current sources may pose a problem for the convergence of the multipole equivalent generator. In particular, it has been verified that the conditions of the criterion are not satisfied in the realistic model that we use. At this stage, we conclude that the basic validity of the multipole approach in electrocardiology has to be reconsidered. Our starting point is the fact that it is mandatory to use an inhomogeneous body conductor model in the multipole approach: for both the integration method and regression method, it has been observed that considerable errors may result if an originally inhomogeneous model is replaced by a homogeneous model (Dub´e et al., 1985). However, using such an inhomogeneous model, one should be seriously concerned with the multipole convergence problem. To obtain more clarity in this matter, we have investigated the multipole convergence rate in a quantitative way by means of the numerical simulations described in Section 5.1. Admittedly, the obtained multipole convergence rates are specific for the particular models that we used for the cardiac current sources and body conductor. It is our hope, however, that the models used represent rather well all the aspects that are relevant for the multipole convergence. For instance, we have included the skeletal muscle layer in the conductor model, and care has been taken to model the spatial distribution of the cardiac current sources as accurately as possible (as pointed out before, the spatial extent of the current sources plays a major role in the convergence). Our conclusion as to the convergence is dual. On the one hand, the doubts raised by the theoretical considerations are confirmed, and it is shown that there is no convergence in the strict mathematical sense of the word. On the other hand, the torso residual errors at L = 5 are sufficiently small, so that the multipole equivalent generator can be considered convergent in the context of practical applications, such as the regression method.

ACKNOWLEDGEMENTS We kindly thank Thom Oostendorp and Adriaan van Oosterom for making the model data and the BEM software described in Section 4.1 available to us.

572

P. De Muynck et al.

A PPENDIX A:

M ULTIPOLE S HIFT E QUATIONS : DERIVATION OF E QUATIONS (3.2) AND (3.13)

Although we have provided analytic expressions for the ∂lm and ylm in (3.2) and (3.4) respectively, neither the derivation of the series expansion of the multipole equivalent generator in Section 3.1 nor the proof of the multipole convergence criterion in Section 3.2 depend on the particular choice of ∂lm and ylm . On the other hand, a particular choice for the ylm needs to be made in order to formulate the multipole shift equations in Section 3.3. Having made this choice in (3.4), it remains to be shown that the corresponding ∂lm are given by (3.2). For this, we consider the following series expansion for the infinite medium Green function G ∞ (r − r0 ) = (4π σ |r − r0 |)−1 (Morse and Feshbach, 1953, p. 1274): +∞

l

X X 2 − δm0 (l − m)! 1 r 01 m m 0 0 = P (cos θ)P (cos θ ) cos(m(φ−φ )) . l 4π σ |r − r0 | 4π σ (l + m)! l r l+1 l=0 m=0 After the substitution of cos m(φ −φ 0 ) by cos(mφ) cos(mφ 0 )+sin(mφ) sin(mφ 0 ), it is possible to express the dependence in r0 of the series expansion via the functions ylm (r0 ): +∞ l 1 1 XX = ylm (r0 )wlm (r), 4π σ |r − r0 | 4π σ l=0 m=−l where the coefficients wlm (r) are given by: m ≥0:

p wlm (r) = (−1)m 2 − δm0

r

4π 2l + 1

s

(l − m)! Plm (cos θ) cos(mφ) (l + m)! r l+1 (A1)

m >0:

p wl,−m (r) = (−1) 2 − δm0 m

r

4π 2l + 1

s

(l − m)! Plm (cos θ) sin(mφ) . (l + m)! r l+1 (A2)

It is obvious that the described series expansion is just the multipole expansion (3.5) for the potential 8 in a point r, where 8 is the potential generated by a unit monopole current source in r0 whose multipole components are given by qlm = ylm (r0 ). The coefficients wlm (r) hence correspond to (3.6) (times 4π σ ), and after some manipulations of the derivatives it can be shown that wlm (r) = 4π σ (−1)l ∂lm G ∞ (r) = ∂lm [(−1)l r −1 ]. To find the ∂lm we need to find derivative operators which generate expressions (A1)–(A2) when operating on (−1)l r −1 , and for this purpose we consider the following identity (Morse and Feshbach, 1953, p. 1281):   Plm (cos θ)eimφ ∂ ∂ m ∂ l−m 1 l +i = (−1) (l − m)! . (A3) ∂x ∂y ∂z l−m r r l+1

Multipole Approach to the Inverse Problem in Electrocardiology

573

Obviously the ∂lm are given by the real and imaginary parts of the operator in the left-hand side of (A3), provided with appropriate multiplication factors. The elaboration of this finally leads to (3.2). It should be noted that (3.2) is uniquely defined up to a term R(∇)∇ 2 , where R is an arbitrary homogeneous polynomial of degree l − 2: this makes sense because the ∂lm only operate on harmonic arguments in the current context. Likewise, in the derivation of the composition rule (3.13), derivative operators are considered to be equal if they differ only by a term R(∇)∇ 2 . For the derivation of (3.13) we first proceed to the ‘normalized’ derivative operators defined by ∂ lm =

∂lm . Flm

It follows that the proof (3.13) is equivalent to the proof of the following composition rule: ∂ l 0 m 0 ∂ lm = α∂ l 0 +l,m α + β∂ l 0 +l,mβ .

(A4)

c In order to prove (A4) we make the link to complex derivative operators ∂ lm :

c ∂ lm

c ∂ l,−m

 ∂ ∂ m ∂ l−m = (−1) +i (m ≥ 0) ∂x ∂y ∂z l−m   ∂ m ∂ l−m ∂ m c∗ −i (m > 0). = (−1) ∂lm = ∂x ∂y ∂z l−m m



(A5) (A6)

c It can easily be verified that the relation between the ∂ lm and the ∂ lm is given by (m > 0):

c ∂ l0 = ∂ l0

and



∂ lm ∂ l,−m



=

√

2



  c 1 1 Re(∂ lm ) = √ 1 c Im(∂ lm ) 2 i

  c  1 ∂ · lm c∗ . (A7) i ∂ lm

c The composition rule for the ∂ lm is as follows: c ∂ lc0 m 0 ∂ lm = ∂ lc0 +l,m 0 +m .

(A8)

The proof of (A8) can be obtained on a case-by-case basis where the cases correspond to the different combinations for the signs of m 0 and m. If m 0 and m are zero or have the same sign, (A8) is proven in a trivial way, using (A5)–(A6). If m 0 and m have different signs, the proof additionally requires substitutions of the type (valid for harmonic arguments): 

∂ ∂ +i ∂x ∂y



∂ ∂ −i ∂x ∂y



→−

∂2 . ∂z 2

574

P. De Muynck et al.

As an example we consider the following case with m ≥ m 0 > 0:   0  0 0 ∂ ∂ m ∂ ∂ m ∂ l −m ∂ l−m c ∂ lc0 ,−m 0 ∂ lm = (−1)m −i +i ∂x ∂y ∂x ∂y ∂z l 0 −m 0 ∂z l−m  m 0   0  ∂ ∂ ∂ ∂ −m +m ∂ ∂ m −i +i +i = (−1) ∂x ∂y ∂x ∂y ∂x ∂y 0

0

∂ l +l−m −m × l 0 +l−m 0 −m ∂z   0 0 0 ∂ ∂ −m +m ∂ l +l+m −m −m 0 +m +i = ∂ lc0 +l,−m 0 +m . = (−1) ∂x ∂y ∂z l 0 +l+m 0 −m The other cases are similar. The composition rule for the ∂ lm can now be derived based on (A8) and the relations (A7). Again a case-by-case approach is followed and the nine combinations for the signs of m and m 0 are covered in four matrices D 1 , . . . , D 4 (M > 0 and M 0 > 0):   ∂ L0 M0 1 D = · [∂ L0 ] ∂ L 0 ,−M 0    c    c    c  1 1 1 1 1 1 ∂ L 0 M 0 ∂ cL0 ∂ L0 M0 · c∗ =√ · c∗ · ∂ L0 = √ 1 1 ∂ L 0 M 0 ∂ cL0 ∂ L0 M0 2 i i 2 i i     c   1 1 1 ∂ L 0 +L ,M 0 ∂ L 0 +L ,M 0 =√ · , = 1 ∂ c∗ ∂ L 0 +L ,−M 0 2 i i L 0 +L ,M 0   D 2 = ∂ L 0 0 · [ ∂ L M ∂ L ,−M ]   1  c  1 1i c c∗ = √ ∂ L 0 0 · [ ∂ L M ∂ L ,M ] · 1 i 2 = [ ∂ L 0 +L ,M ∂ L 0 +L ,−M ] ,     D 3 = ∂ L 0 0 · ∂ L0 = [∂ cL 0 0 ] · [∂ cL0 ] = [∂ cL 0 +L ,0 ] = [∂ L 0 +L ,0 ]   ∂ L0 M0 4 · [ ∂ L M ∂ L ,−M ] D = ∂ L 0 ,−M 0    c    1 1 1 ∂ L0 M0 1 1i c c∗ = · c∗ · [ ∂ L M ∂ L ,M ] · 1 i ∂ L0 M0 2 1i i    c    1 1 1 ∂ L 0 M 0 ∂ cL M ∂ cL 0 M 0 ∂ c∗ 1 1i L M = · c∗ c · c∗ 1 i ∂ L 0 M 0 ∂ L M ∂ c∗ 2 1i i L0 M0 ∂ L M    c  1 c 1 1 1 ∂ L 0 M 0 ∂ cL M + ∂ cL 0 M 0 ∂ c∗ (∂ L 0 M 0 ∂ cL M − ∂ cL 0 M 0 ∂ c∗ ) L M L M i = · c∗ c 1 c∗ c∗ c∗ (∂ L 0 M 0 ∂ cL M − ∂ c∗ ∂ L 0 M 0 ∂ L M + ∂ c∗ 2 1i i L0 M0 ∂ L M L0 M0 ∂ L M ) i

Multipole Approach to the Inverse Problem in Electrocardiology

575

Re(∂ cL 0 M 0 ∂ cL M ) + Re(∂ cL 0 M 0 ∂ c∗ LM) c c c c∗ Im(∂ L 0 M 0 ∂ L M ) + Im(∂ L 0 M 0 ∂ L M )



Im(∂ cL 0 M 0 ∂ cL M ) − Im(∂ cL 0 M 0 ∂ c∗ LM) c c c −Re(∂ L 0 M 0 ∂ L M ) + Re(∂ L 0 M 0 ∂ c∗ LM)   Re(∂ cL 0 +L ,M 0 +M ) Im(∂ cL 0 +L ,M 0 +M ) = Im(∂ cL 0 +L ,M 0 +M ) −Re(∂ cL 0 +L ,M 0 +M )   c c M Re(∂ L 0 +L ,M 0 −M ) −Im(∂ L 0 +L ,M 0 −M ) +(−1) Im(∂ cL 0 +L ,M 0 −M ) Re(∂ cL 0 +L ,M 0 −M )

=



= D 4a + D 4b , with D

 1 ∂ L 0 +L ,M 0 +M =√ 2 ∂ L 0 +L ,−M 0 −M

4a

 ∂ L 0 +L ,−M 0 −M . −∂ L 0 +L ,M 0 +M

The conversion to operators ∂ lm in D 4b is a bit more complex because it depends on the sign of M 0 − M; one has (m > 0): c Re(∂ lm ) = √12 ∂ lm c Re(∂ l0 ) = ∂ l0 c c∗ Re(∂ l,−m ) = (−1)m Re(∂ lm )= c Im(∂ lm ) = √12 ∂ l,−m c Im(∂ l0 )=0 c c∗ Im(∂ l,−m ) = (−1)m Im(∂ lm )=

m (−1) √ ∂ lm 2

m+1 (−1) √ ∂ l,−m . 2

After synthesis of the three different cases for the sign of M 0 − M we obtain: (−1)min(M ,M) √ 2  √ 1 + δ M 0 M ∂ L 0 +L ,|M 0 −M| × sign(M 0 − M)∂ L 0 +L ,−|M 0 −M| 0

D 4b =

−sign(M − M 0 )∂ L 0 +L ,−|M 0 −M| √ 1 + δ M 0 M ∂ L 0 +L ,|M 0 −M|



with the value of sign, respectively, given by 0, −1 and 1 for zero, negative and positive arguments. The nine cases for the signs of m and m 0 are now combined in one matrix D, which relates as follows to D 1 , . . . , D 4 :  ∂ L 0 M 0 ∂ L0 ∂ L 0 M 0 ∂ L ,−M ∂ L0 M0 ∂ L M D =  ∂ L00∂ L M ∂ L 0 0 ∂ L0 ∂ L 0 0 ∂ L ,−M  ∂ L 0 ,−M 0 ∂ L M ∂ L 0 ,−M 0 ∂ L0 ∂ L 0 ,−M 0 ∂ L ,−M  4b   4a  4b 4a 1 D11 0 D12 D11 D11 D12 3 2 2   0  D11 0 0  D11 D12 = + . 4b 4b 4a 4a 1 D21 0 D22 D21 D21 D22 (D α ) (D β ) 

576

P. De Muynck et al.

Comparing this expression with (A4), we see that D contains all individual cases for ∂ l 0 m 0 ∂ lm : D α represents the first α-dependent term of the right-hand side of (A4), and D β represents the second term. This means that we have finally proven the composition rule (A4), and consequently also (3.13), with values for α, β, m α and m β that are given in the tables below (M > 0 and M 0 > 0): α

m=M

m=0

m = −M

mα

m=M

m=0

m = −M

m0 = M 0

√1 2

1

√1 2

m0 = M 0

M0 + M

M0

−M 0 − M

m0 = 0

1

1

1

m0 = 0

M

0

−M

m 0 = −M 0

√1 2

1

−1 √ 2

m 0 = −M 0

−M 0 − M

−M 0

M0 + M

β m0 = M 0

m=M q 0 (−1)min(M ,M) 1+δ M 0 M √ 2

m0 = 0

0

0 0 (−1)min(M ,M) √ sign(M −M) m 0 = −M 0 2

A PPENDIX B:

m=0

m = −M

mβ

m=M

m=0

m = −M

0

0 0 (−1)min(M ,M) √ sign(M−M ) − 2

m0 = M 0

|M 0 − M|

-

−|M 0 − M|

m0 = 0

-

-

-

-

|M 0 − M|

0 0

0 q 0 (−1)min(M ,M) 1+δ M 0 M √ 2

m 0 = −M 0 −|M 0 − M|

G ENERALIZED M ULTIPOLE A PPROXIMATION : P ROOF OF S TATEMENT (5.2)

The proof of (5.2) is given in two steps. In the first step, the criterion is proved for the particular case that the current source region can be contained within a homogeneous sphere. We denote this sphere and its centre by B1 and O1 , respectively. Hence, a homogeneous region C can be constructed, as the union of a sequence of spheres B1 , . . . , Bn which satisfy the following requirements: (1) the first sphere (B1 ) contains the cardiac current sources; (2) each subsequent sphere contains the centre of the preceding sphere, i.e., ri ∈ Bi+1 , i = 1, . . . , n − 1; (3) the centre of the last sphere must coincide with the multipole origin, i.e., On = O; and (4) each sphere lies in V , i.e., C is a subregion of V , where V is the connected homogeneous subregion of the conductor mentioned in (5.2). An example of such a sequence of spheres is shown in Fig. 10. For the proof of the criterion (5.2), we first consider the multipole series expansion of 8 at the multipole origin O1 , which then, in turn, is approximated by a multipole series expansion at O2 , etc. Partial sums of these multipole series expansions at the consecutive multipole origins O1 , . . . , On are respectively denoted as 81 , . . . , 8n (the (i) multipole coefficients of 8i are denoted as qlm and the maximum multipole order of 8i is denoted as L i ). This means that 8i , i = 1, . . . , n, is an approximation for 8i−1 (80 is used here as an alternative notation for 8). Moreover, from the

Multipole Approach to the Inverse Problem in Electrocardiology

577

2

Figure 10. Sequence of spheres B1 , . . . , Bn with their respective centers O1 , . . . , On .

requirements (1) and (2) it follows that the multipole convergence criterion is satisfied for all i = 1, . . . , n (with Bi as the separating sphere), so that the respective approximation errors |8i − 8i−1 | can be made arbitrarily small over Bi (the overbar designates the complementary region). Considering the described sequence of approximations and the fact that On = O, we suggest 8n as the multipole approximation 6 Le 8 for 8. Hence, we must prove that for any arbitrary small approximation error ε, one can find a 8n such that: sup|8n − 80 |{V } < ε,

(B1)

where V is the region in which the generalized multipole approximation is being considered [see criterion (5.2)]. In the notation used, the supremum symbol is followed by the function and, further, by the domain over which the supremum is taken (between brackets). We will actually prove that one can find a 8n for which: sup|8n − 80 |{C} < ε,

(B2)

n where C = ∪i=1 Bi . Since V is a subregion of C, it is evident that (B1) follows from (B2). The starting point of the proof of (B2) is the following triangular inequality:

sup|8n − 80 |{C} ≤

n X

sup|8i − 8i−1 |{C}.

(B3)

i=1

For each i = 1, . . . , n, the relationship Bi ⊆ C implies that C ⊆ Bi , so that sup|8i − 8i−1 |{C} ≤ sup|8i − 8i−1 |{B i }.

(B4)

578

P. De Muynck et al.

From (B3) and (B4), it follows that sup|8n − 80 |{C} ≤

n X

sup|8i − 8i−1 |{Bi }.

(B5)

i=1

The inequality (B5) shows that (B2) can be satisfied if the approximation error |8i − 8i−1 | over Bi can be made smaller than ε/n for each i = 1, . . . , n. Our previous considerations pointed out that this can be done, so that the first step of the criterion is proved. It is understood that in the construction of 8n , each 8i depends on 8i−1 and on (i) the desired value of ε/n. The multipole coefficients qlm of 8i are a function of the multipole coefficients of 8i−1 , namely via the multipole shift equations (3.19). Likewise, the number of terms that are used in 8i (i.e., the maximum multipole order L i ) also depends on the multipole coefficients of 8i−1 , as well as on ε/n. For the multipole coefficients of the multipole approximation 6 Le 8 = 8n , one has L n−1

qlen m n (L)

=

qµ(n)n

=

X

ln−1 X

ln−1 =0 m n−1 =−ln−1

qµ(n,n−1) n µn−1

...

L1 l1 X X

qµ(2,1) · qµ1 1 . (B6) 2 µ1

l1 =0 m 1 =−l1

Here, the factors qµ(i,i−1) are the coefficients of the multipole shift equations (3.19) i ,µi−1 for the origins Oi−1 and Oi (note that µ = l(l + 1) + m is a shorthand notation (1) for the combined indices l and m). By definition, the multipole coefficients q1m are the cardiac multipole components with respect to the multipole origin O1 . The (2) multipole coefficients qlm are equal to the cardiac multipole components at O2 for l ≤ L 1 , but for L 1 < l ≤ L 2 they will be different [see (3.19), and the fact that the multipole components of 81 are set to zero for l > L 1 ]. Similar considerations are true for the multipole series expansions at O3 , . . . , On . It consequently appears that only the lower order multipole coefficients of 6 Le 8 will coincide with the cardiac multipole components; for orders higher than L 0 = min(L 1 , . . . , L n−1 ) they will differ. We conclude that the obtained multipole approximation 6 Le 8 coincides with the multipole equivalent generator for the lower order terms l < L 0 , and that a discrepancy will occur for the higher order terms l > L 0 . It is possible to increase the number of the coinciding terms L 0 . On the other hand, the discrepancy for the higher order terms remains, because L = L n will inevitably increase as well. In regard of this, one can consider the procedure of the multipole approximation 6 Le 8 as a method to make the series 6 L 8 convergent; the convergence is achieved by a regrouping of the higher order terms of 6 L 8 (which are no longer added to the approximation in a strictly sequential way, as is the case of the multipole series expansion). In the second step of the proof of (5.2), i.e., for the case that the current source region cannot be contained within a homogeneous sphere, the current source region is partitioned into subregions such that each individual subregion can be contained

Multipole Approach to the Inverse Problem in Electrocardiology

579

2

Figure 11. Partitioning of the current source region into B11 and B12 and construction of the corresponding sequences of spheres.

in a homogeneous sphere. Consequently, each of these spheres is used as the starting sphere B1 of a sequence of spheres such as those of Fig. 10. This is illustrated in Fig. 11. The current source I are partitioned according to the partitioning of the source region, I = I1 + · · · + I N. The electrocardiological potential can accordingly be decomposed as 8 = 81 + · · · + 8 N .

(B7)

The multipole approximation 6 Le 8 is defined as the sum of the approximations for the individual potentials 8 j , 6 Le 8 =

N X

8nj j .

(B8)

j=1

j

Here, n j denotes the number of spheres that is used to go from the first sphere B1 j (the sphere containing I j ) to the last sphere Bn j , which is centered at the multipole j origin. On the basis of the first part of the proof, we know that 8n j approximates 8 j up to an arbitrarily small error outside the region C j , which is the region of the jth row of spheres, nj j C j = ∪i=1 Bi . We will now prove (5.2), showing that a multipole approximation 6 Le 8 of the type (B8) can be constructed, so that sup|6 Le 8 − 8|{V } < ε.

(B9)

580

P. De Muynck et al.

To prove the inequality (B9), we impose a maximum approximation error of ε/N j for each 8n j , ε sup|8nj j − 8 j |{C j } < . (B10) N We now consider the joint region of all the spheres C = ∪ Nj=1 C j . On the basis of (B10), it can be shown that ε is an upper limit for the error |6 Le 8 − 8| over C, the complementary region of C, sup|6 Le 8

N N X X j j − 8|{C} = sup 8n j − 8 {C} j=1

≤ sup

N X

j=1

|8nj j − 8 j |{C}

j=1

≤

N X

sup|8nj j − 8 j |{C}

N X

sup|8nj j − 8 j |{C j } < ε.

j=1

≤

j=1

j

Since all the spheres Bi lie within V , the region V is a subregion of C, and it follows that sup|6 Le 8 − 8|{V } ≤ sup|6 Le 8 − 8|{C} < ε. The criterion (5.2) is therefore generally proven for the case that the cardiac current sources and the multipole origin are in one and the same homogeneous subregion of the body conductor.

R EFERENCES Arthur, R. M. (1999). Spherical-harmonic approximation to the forward problem of electrocardiology. J. Electrocardiol. 32, 103–114. Arthur, R. M., D. B. Geselowitz, S. A. Briller and R. F. Trost (1972). Quadrupole components of the human surface electrocardiogram. Am. Heart J. 83, 663–677. Arzbaecher, R. C. and D. A. Brody (1976). The lead field: vector and tensor properties, in The Theoretical Basis of Electrocardiography, Nelson and Geselowitz (Eds), Oxford: Clarendon Press, pp. 175–201.

Multipole Approach to the Inverse Problem in Electrocardiology

581

Barnard, A. C. L., I. M. Duck, M. S. Lynn and W. P. Timlake (1967). The application of electromagnetic theory to electrocardiology—derivation of integral equations. Biophys. J. 7, 443–462. Barr, R. C., T. C. Pilkington, J. P. Boineau and M. S. Spach (1966). Determining surface potentials from current dipoles, with application to electrocardiography. IEEE Trans. Biomed. Eng. 13, 88–92. Barr, R. C., M. Ramsey and M. S. Spach (1977). Relating epicardial to body surface potential distributions by means of transfer coefficients based on geometry measurements. IEEE Trans. Biomed. Eng. 24, 1–11. Brody, D. A. (1956). A theoretical analysis of intracavitary blood mass influence on the heart lead relationship. Circ. Res. 4, 731–738. Brody, D. A., J. C. Bradshaw and J. W. Evans (1961). The elements of an electrocardiographic tensor theory. Bull. Math. Biophys. 23, 31–42. Brody, D. A., J. W. Cox and L. G. Horan (1973). Hexadecapolar shift equations applicable to equivalent cardiac generators of lower degree. Ann. Biomed. Eng. 1, 481–488. Cornelis, J. (1980). Studie van de verwerking van elektrocardiologische potentialen, PhD dissertation, Vrije Universiteit Brussel, Belgium (in Dutch). Cuppen, J. J. M. and A. van Oosterom (1984). Model studies with the inversely calculated isochrones of ventricular depolarization. IEEE Trans. Biomed. Eng. 31, 652–659. Damen, A. (1980). On the observability of electrical cardiac sources, PhD dissertation, Technische Hogeschool Eindhoven, The Netherlands. de Guise, J., R. M. Gulrajani, P. Savard, R. Guardo and F. A. Roberge (1985). Inverse recovery of two moving dipoles from simulated surface potential distributions on a realistic human torso model. IEEE Trans. Biomed. Eng. 32, 126–135. De Muynck, P. (1994). Multipooltheorie in de Elektrocardiologie, PhD dissertation, Vrije Universiteit Brussel, Belgium (in Dutch). De Muynck, P. and J. Cornelis (1996). Multipole approach to the inverse electrocardiological problem for the realistic inhomogeneous body conductor, in Computers in Cardiology, Vienna 1995, Murray and Arzbaecher (Eds), IEEE Computer Society Press, pp. 91–94. De Muynck, P., E. Nyssen, J. Cornelis, P. Block and O. Steenhaut (1986). Exploring the limits of ECG preprocessing methods, insensitive to interindividual differences in thoracic properties, in Computers in Cardiology, Boston, 1986, IEEE Computer Society Press, pp. 219–222. Dub´e, B., P. Savard, R. Guardo, R. M. Gulrajani and J. P. Drouhard (1985). Surface integration and least-squares procedures for the inverse recovery of cardiac multipole components. Ann. Biomed. Eng. 13, 43–58. Epstein, B. R. and K. R. Foster (1983). Anisotropies in the dielectric properties of skeletal muscle. Med. Biol. Eng. Comput. 21, 51–55. Geselowitz, D. B. (1960). Multipole representation for an equivalent cardiac generator. Proc. IRE 48, 75–79. Geselowitz, D. B. (1997). The theoretical basis of bioelectromagnetic fields. Invited lecture in Biomedizinische Technik 42 (Erg. 1).

582

P. De Muynck et al.

Gulrajani, R. M. and G. E. Mailloux (1983). A simulation of the effects of torso inhomogeneities on electrocardiographic potentials, using realistic heart and torso models. Circ. Res. 52, 45–56. Gulrajani, R. M., P. Savard and F. A. Roberge (1988). The inverse problem of electrocardiography: solution in terms of equivalent sources. CRC Crit. Rev. Biomed. Eng. 16, 171–214. ¨ Helmholtz, H. (1853). Uber einige Gesetze der Vertheilung elektrischer Str¨ome in k¨orperlichen Leitern mit Anwendung auf die thierisch-elektrischen Versuche. Pogg. Annln. Physik Chemie 89, 211–233 and 353–377. Jackson, J. D. (1975). Classical Electrodynamics, John Wiley & Sons Inc. Kellogg, O. D. (1929). Foundations of Potential Theory, New York: Dover Publications. Kneppo, P. and L. I. Titomir (1979). Integral characteristics of the human cardiac electrical generator from electric field measurements by means of an automatic cylindrical coordinator. IEEE Trans. Biomed. Eng. 26, 21–28. Lynn, M. S. and W. P. Timlake (1968a). The numerical solution of singular integral equations of potential theory. Numer. Math. 11, 77–98. Lynn, M. S. and W. P. Timlake (1968b). The use of multiple deflations in the numerical solution of singular systems of equations, with applications to potential theory. SIAM J. Numer. Anal. 5, 303–322. Malmivuo, J. and R. Plonsey (1995). Bioelectromagnetism. Principles and Applications of Bioelectric and Biomagnetic Fields, Oxford: Oxford University Press. Martensen, E. (1968). Potentialtheorie, Stuttgart: B.G.Teubner (in German). McFee, R. and S. Rush (1968). Qualitative effects of thoracic resistivity variations on the interpretation of electrocardiograms: the low resistance surface layer. Am. Heart J. 76, 48–61. Meijs, J. W. H., O. W. Weier, M. J. Peters and A. van Oosterom (1989). On the numerical accuracy of the boundary element method. IEEE Trans. Biomed. Eng. 36, 1038–1049. Morse, P. M. and H. Feshbach (1953). Methods of Theoretical Physics, McGraw-Hill Book Company. Nolte, G. and G. Curio (1997). On the calculation of magnetic fields based on multipole modeling of focal biological current sources. Biophys. J. 73, 1253–1262. Nyssen, E., J. Cornelis, Y. Christophe and P. Block (1986). ECG Preprocessing based on a physical model, Proc. IV Mediterranean Conference on Medical and Biological Engineering—MECOMBE’86, Sevilla-Spain: pp. 200–203. Pilkington, T. C. and M. N. Morrow (1982). The usefulness of multipoles in electrocardiography. CRC Crit. Rev. Biomed. Eng 175–192. Plonsey, R. (1966). Limitations on the equivalent cardiac generator. Biophys. J. 6, 163– 173. Plonsey, R. and D. B. Heppner (1967). Considerations of quasi-stationarity in electrophysiological systems. Bull. Math. Biophys. 29, 657–664. Rudy, Y. and B. J. Messinger-Rapport (1988). The inverse problem in electrocardiography: solutions in terms of epicardial potentials. CRC Crit. Rev. Biomed. Eng. 16, 215–268. Rush, S., J. Abildskov and R. Mc Fee (1963). Resistivity of body tissues at low frequencies. Circ. Res. 12, 40–50.

Multipole Approach to the Inverse Problem in Electrocardiology

583

Takano, N., R. Lawson, J. Hyttinen and J. Malmivuo (1996). Multipole like term estimation from the standard ECG records. Med. Biol. Eng. Comput. 34 (Suppl. 1), 89–90. Tanaka, H., K. Yajima, T. Ihara and T. Furukawa (1980). An inverse solution in electrocardiology—estimation of epicardial potential mapping using spheroidal harmonics, in MEDINFO 80, Lindberg and Kaihara (Eds), North Holland Publishing Company, pp. 254–258. Tanaka, H., K. Yayima, T. Ihara and T. Furukawa (1982). An inverse solution of electrocardiology—epicardial potential estimation by using the orthogonal expansion method. Japanese Heart Journal 23 (Suppl.), 352–354. Titomir, L. I. and P. Kneppo (1985). Simultaneous analysis of the cardiac electric and magnetic fields using the scalar multipole expansion. Bulletin of Mathematical Biology 47, 123–143. Titomir, L. I. and P. Kneppo (1994). Bioelectric and Biomagnetic Fields. Theory and Applications in Electrocardiology, Boca Raton: CRC Press. Received 5 May 1998 and accepted 16 December 1999