## Info

Fig. 7.7. The results obtained by applying the second-order Tikhonov scheme to the problem discussed in Section 7.2.5. In the second row noise has been added to the observation data (still recorded along the body surface dT, cf. Figure 7.4). a) Estimates obtained by using the regularization parameter e = 10~8. b) The body surface potentials associated with the epicardial potentials shown in Figure a). c) In this case, 1% noise has been added to the observation data, and e = 1. The figure shows the epicardial potential computed by the second order scheme along with the correct potential. d) The noisy observation data and the simulated body surface potential associated with the epicardial potentials depicted in Figure c).

For both the simulation results shown in Figure 7.8, the associated body surface potentials were reproduced accurately. In fact, they were very similar to the potential depicted in Figure 7.7d).

7.3. Computing the Location and Orientation of a Dipole 253 Potentials at the heart surface Potentials at the heart surface

100 200 300 400 500 600 100 200 300 400 500 600

Epicardial node number Epicardial node number a) b)

Fig. 7.8. a) Temporal regularization with e = 0.01 in equation (7.60). b) "Hybrid" regular-ization with e = 0.01 and ¡3 = 1 in equation (7.67).

7.3 Identifying a Simplified Source - Computing the Location and Orientation of a Dipole

The electrical signal that can be measured on the surface of the human body is generated by the heart. The heart can thus be viewed as a current source positioned in a volume conductor. In Chapter 2 a simple model of this phenomenon was introduced,:

J n cf. the discussion preceding equations (2.5)-(2.6). Here u, M, f and n denote the electrical potential, the conductivity of the involved tissues, the current source(s) and the region occupied by the human body, respectively. Throughout this section we will focus on this simple elliptic boundary value problem. Note that n = H U T (see Figure 7.4), and thus the solution domain of (7.68) consists of both the heart H and its surrounding body T.

In these kinds of models, the function f represents the electrical activity of the heart. As mentioned in Chapter 1, this activity is frequently modelled by a so-called dipole; see Figure 3.4. This means that f is, at least approximately, defined in terms of two opposite point charges15.

15 In many cases more than two point charges are used to represent the electrical activity of the heart. This leads to models of the form (7.68)-(7.70), involving several dipoles. In the present text, we will only consider the case of a single dipole.

Models of this kind, and the heart vector interpretation of them, go back as far as Einthoven's work in the early 20th century. They have always been a central part of ECG analysis; see, e.g. Keener and Sneyd [72]. The purpose of the present section is not to discuss the physiological or biological reasons for studying models of the form (7.68)-(7.70) in detail. Rather, we will focus on various mathematical and computational properties of these equations.

Assume that the characteristics of M, Q and f are known. Then we can use the Finite Element Method (FEM), discussed in Chapter 3, to solve (7.68)-(7.70) numerically for u. Consequently, we can simulate ECG recordings along the surface dQ of the body, and study how various changes in the dipole (heart) influence the electrical signal. In the present situation, this is the direct problem: we compute the effect (ECG) of a dipole (the cause) on the electrical potential.

In this section, our ultimate goal is not to solve this direct problem, but to design methods for computing the characteristics of the current source f from body surface measurements of the electrical potential u. That is, we will try to use observations of u\sn to compute the properties of f. From a medical point of view, solving this inverse problem can hopefully improve our qualitative and quantitative understanding of the heart; see, e.g. Gulrajani [54], Pullan, Paterson and Greensite [116], Yamashita and Geselowitz [147], and references therein.

Typically, the position and strength of the point charges in the dipole f will depend on time; they will move around during a heart beat. Consequently, f should ideally be a function of both space x e Q and time t e [0,t*j, where [0,t*j represents the time interval of interest. This means that the inverse problem described above must be solved for every t e [0,t*j. Assuming that we have a method for solving this inverse problem at a fixed time, this can be easily be accomplished by applying similar techniques to those introduced in Section 7.2.4. However, to keep things simple, we will throughout this section focus on the stationary version of this task.

Note that, if u satisfies (7.68)-(7.69), the function u + c, where c is an arbitrary constant, also satisfies these two equations. Thus, condition (7.70) has been added to ensure the uniqueness of the solution. The role of the constant c in the present framework will be studied in detail in Section 7.3.5 below.

A Remark. In Chapter 2 we introduced the highly complex model (2.65)-(2.72) for the electrical activity in the human body. Now we want to study the heart by considering the simple system of equations (7.68)-(7.70). Is it possible to justify this approach? Comparing these two models in a precise mathematical sense is, of course, a very challenging scientific problem. Nevertheless, it is actually quite simple to establish a rough link between (2.65)-(2.72) and (7.68)-(7.70). This can be accomplished as follows.

From equation (2.67) in Chapter 2 it follows that

Moreover, in the tissue surrounding the heart, the propagation of the electrical potential is governed by a model of the form

cf. (2.65)-(2.72) and Figure 7.4. Thus, by simply putting

u0 in T, we see that u, as defined in (7.73), satisfies an equation on the form (7.68).

A further investigation of the link/relationship between these two models would require a deep and lengthly analysis of the involved PDEs. That topic certainly lies beyond the scope of this text.

### 7.3.1 Preliminaries

Prior to presenting the mathematics of the inverse problem described above, we will derive a property that the source term f in (7.68)-(7.70) must satisfy; cf. also Hackbusch [56] or Marti [99].

Assume that u solves (7.68)-(7.70), and let fa be a smooth function defined on Q. From (7.68) it follows that ffadx = V ■ (MVu)fa dx, Ja

J n J n and by Gauss' divergence theorem, and equation (7.69), we conclude that

Jn Jn Jen Jn

Since fa was arbitrary, we may choose fa(x) = 1 for all x e Q, and then (7.74) implies that i fdx = 0. (7.75)

Consequently, if (7.68)-(7.70) has a solution, f must satisfy (7.75). That is, this boundary value problem cannot have a solution unless f fulfills (7.75). Thus, we will throughout this section assume that the involved source term belongs to the following subspace V of the classical L2(Q) space,

7.3.2 The Inverse Problem

Our goal is to derive a mathematical formulation of the inverse problem described above. Clearly, the solution u of (7.68)-(7.70) depends on the source term f, i.e.

Having computed u(f) we can read off the trace (restriction)

u(f) a, of this function along the boundary dQ of Q. Thus, we may define the operator

an, where V is the space of functions defined in (7.76). In the present context, the task of computing Q(f), for a given f, is the direct problem. In other words, the direct problem consists of solving (7.68)-(7.70), for a given source f, and thereafter recording the simulated electrical potential u along the boundary dQ of the torso Q.

Assume that we want to apply noninvasive methods for characterizing the current source. Only observations of the potential along the surface dQ of Q are then available to us as data. That is, u(f) (7.79)

an is the available data and f represents the unknown source. So, as mentioned above, our challenge can be described as follows: Can we use data in the form (7.79) and the model (7.68)-(7.70) to determine f ?

In mathematical terms this problem takes the form: Let d be an observation of the potential along dQ, solve the equation

for f. Note that, if Q is invertible, then f = Q-1(d).

Throughout this section we will assume that f depends on the spatial variable x, and that this function is parameterized by a set of parameters p1,... ,pn, i.e.

Our task is to determine p1,...,pn such that (7.80) is (approximately) satisfied. Depending on the sort of parameterization that is used, this will lead to either a linear system of algebraic equations or a nonlinear minimization problem defined in terms of a suitable cost-functional.

an an

### 7.3.3 Parameterizations Leading to Linear Problems

We will now show that the operator Q, defined in (7.77), is a linear mapping. Thereafter, this property will be used to design a discrete approximation of (7.80), which in turn leads to a linear system of algebraic equations.

Assume that f e V, and let u = u(f) denote the solution of (7.68)-(7.70). For an arbitrary constant c e R, consider the function v = cu. Clearly

Jn Jn

This means that v = cu satisfies the following set of equations

Thus, since

Next, consider the boundary value problems

/ u2 dx = 0, n where fi and f2 are functions in V; see (7.76). Clearly,

V • (MV(u1 + u2)) = V • (MVu1) + V • (MVu2) = f1 + f2 in tt, MV(u1 + u2) • n = MVu1 • n + MVu2 • n = 0 along dtt,

Jn Jn Jn and hence it follows that v = u1 + u2 solves the problem

V • (MVv) = f1 + f2 in ii, MVv • n = 0 along dQ,

This means that

0 0