Modeling The Invers Electrophysiology of the Heart Activation As Anisotropic Volume Source Inside Inhomogeneous Volume Conductor

1. Introduction

There are two problems in modeling the electrophysiology of the heart: the Forward Problem and the Inverse Problem [1]. The heart Forward Problem involves designing a model that is capable of determining the field on the surface of known body (conductor) that is generated by electrical-sources inside the heart. Solving the Forward Problem requires the development of electrical models which are capable of describing the bioelectrical criteria of the heart and the body (Figure 1).

The Inverse Problem of the heart electrophysiology can be summarized as how to find detailed information about the distribution of the bioelectric sources in the heart from the electric potentials on the body surface (Body Surface Potential Maps BSPM) or the magnetic fields (MCG) outside the body.

To find the source, given the measured field, there are several approaches that have been used [1]. These include:

1-

An empirical approach based on the recognition of typical signal patterns that are known to be associated with certain source configurations.

2-

Imposition of physiological constraints based on the information available about the anatomy and physiology of the active tissue. This imposes strong limitations on the number of available solutions.

3-

Modeling the source and the volume conductor using simplified models. The source is characterized by only a few degrees of freedom (for instance a single dipole which can be completely determined by three independent measurements).

4-

Examining the lead field pattern, from which the sensitivity distribution of the lead and therefore the statistically most probable source configuration can be estimated.

Where it is difficult to compare one heart model against another [2], validation of the Inverse Problem solution is usually done by comparing results of the system to other real measurements. Real measurements are usually obtained from measuring the volume source potential (usually the Epicardium wall or the Endocardium wall potential) using invasive techniques [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]. Other methods tend to use simulators of the Forward Problem in the process and then compare the results obtained from the Inverse Model with the ones generated by the Forward Model [18,19,20,21,22,23,24,25].

2. Methods

There is a linear relation between the source and the body surface potential. This relation can be introduced as

$$b=Ax$$

(1)

where b is the body surface potential vector (dependent variables vector [Mx1]), x is the source vector (independent variables vector [Nx1]), and A is the lead field transfer matrix (MxN) that maps x to b. Thus, solving the Inverse Problem of the heart requires solving the equation

$$x=Tb$$

(2)

where T is the inverse transfer matrix. Different forms of the transfer matrix A has been introduced in terms of Transmembrance Potential (TMP) instead of Current Destiny vector [26] and it is apparently singular matrix. Thus, finding T cannot be accomplished by regular linear system solution such as Gaussian Elimination.

2.1. Regularization

In general, the Inverse Problem of the heart’s electrophysiology is an ill-posed problem because it breaks one or more of Hadamard conditions of well-posed problem [27,28,29]. As a consequence of the solution not being unique, the error of body surface potential reading being presented in all configurations of the Inverse Problem, and in the multiple pole distribution configuration, the problem becomes strongly underdetermined [22,25] as the number of body surface reading is much less than the number of sources.

Solving the ill-posed Inverse Problem requires reformulating the problem by including some assumptions such as smoothness of solution, which is called the regularization [27,28,29]. Tikhonov regularization [30] is the most common technique used to solve ill-posed problems such as this Inverse Problem is solved by minimizing both the residual error norm and the regularized solution norm. Using this technique the objective function can be written as

$$\begin{array}{l}\stackrel{\Lambda }{x}={\min }_{\arg (x)}\left\{{‖b-Ax‖}^2+\alpha .{‖\mathrm{\Gamma }x‖}^2\right\}\\ ={\min }_{\arg (x)}\left\{{‖b-Ax‖}^2+\alpha .x^T{\mathrm{\Gamma }}^T\mathrm{\Gamma }x\right\}\end{array}$$

(3)

where $\stackrel{\Lambda }{x}$ is an estimate of the solution, $\alpha$ is the regularization parameter and it is usually determined by L-curve [30,31,32,33] or U-curve [34] and $\mathrm{\Gamma }$ is the regularization matrix.

2.2. The Minimum Norm (MN) Solution

Three methods for solving this function are described here. The Minimum Norm (MN) solution of the Tikhonov regularization is provided by assigning the regularization matrix to the identity matrix. The objective function then becomes:

$$\stackrel{\Lambda }{x}=\min \left\{{‖b-Ax‖}^2+\alpha .{‖x‖}^2\right\}$$

(4)

and the solution of this objective function will be:

$$\begin{array}{l}\stackrel{\Lambda }{x}=A^T{(A{A}^T+\alpha .I)}^{-1}b\\ ={(A^{T}A+\alpha .I)}^{-1}{A}^{T}b\end{array}$$

(5)

2.3. The Weighted Minimum Norm (WMN) Solution

The Weighted MN (WMN) solution is by considering the regularization matrix as a weighting matrix. The objective function then becomes:

$$\stackrel{\Lambda }{x}=\min \left\{{‖b-Ax‖}^2+\alpha .{‖Wx‖}^2\right\}$$

(6)

and the solution of this objective function will be:

$$\begin{array}{l}\stackrel{\Lambda }{x}={(W^{T}W)}^{-1}{A}^T{[A{(W^{T}W)}^{-1}{A}^T+\alpha .I]}^{-1}b\\ ={[A^{T}A+\alpha .(W^{T}W)]}^{-1}{A}^{T}b\end{array}$$

(7)

where W is the weighting matrix. The weighting matrix can be smoothed by any smoothing functions such as Laplacian smoothing operator.

2.4. The Low Resolution Brain Electromagnetic Tomography (LORETA) solution

In addition to the MN and WMN solution there is the Low Resolution Brain Electromagnetic Tomography (LORETA) solution [35,36,37,38,39,40] which was developed by Pascual-Marqui. The LORETA methods are linear imaging methods, used mainly to localize the electric neuronal activity in the brain. There are three versions of LORETA methods. The first is the original LORETA, that has good accuracy in localizing test sources even when they are deep and this makes the overall average localization error to be smaller than one voxel. The comparisons between LORETA method and other methods have proven that LORETA has the smallest localization error [36,37,38,41]. The second version is the standardized LORETA (sLORETA) which has zero localization error under ideal conditions of no noise. The last version is the exact LORETA (eLORETA), which is the most advanced version of these models and was derived to include the reading noise error in account making the solution more realistic. The objective matrix form of eLORETA is

$$\stackrel{\Lambda }{x}=W^{-1}{A}^T{(A{W}^{-1}{A}^T+\alpha H)}^{-1}b$$

(8)

where

$$H=I-{\displaystyle \frac{1^{T}1}{\begin{array}{cc}1& 1^T\end{array}}}=\left[\begin{array}{cccc}1& 0& \mathrm{..}& 0\\ 0& 1& \mathrm{..}& 0\\ :& :& :& :\\ 0& 0& \mathrm{..}& 1\end{array}\right]-\left[\begin{array}{cccc}{\displaystyle \frac{1}{N}}& {\displaystyle \frac{1}{N}}& \mathrm{..}& {\displaystyle \frac{1}{N}}\\ {\displaystyle \frac{1}{N}}& {\displaystyle \frac{1}{N}}& \mathrm{..}& {\displaystyle \frac{1}{N}}\\ :& :& :& :\\ {\displaystyle \frac{1}{N}}& {\displaystyle \frac{1}{N}}& \mathrm{..}& {\displaystyle \frac{1}{N}}\end{array}\right]$$

(9)

such that

$$W_j^2=A_j^T{(A{W}^{-1}{A}^T+\alpha H)}^{-1}{A}_j$$

(10)

where H is the average centering operator, and the suffix j means the column number in the matrix as the matrix W is the diagonal matrix, and I is the unity vector.

There are two important matters regarding the Inverse Problem solution. The regularization algorithms are efficient in localizing single primary source [42], but more sources produce inaccurate results. The other matter is that there are some difficulties in localizing sources that are in some zones such as boundaries and Septum of the heart [24].

2.5. Data Generation

The implementation of the Inverse Problem solution is similar to the Forward Problem solution and accomplished using the C/C++ programming language in addition to the OpenGL library to visualize Inverse Model results. Matrix operations in the implementation are performed by some of matrix operations libraries that are found in the ALGLIB site [43]. These libraries have been verified to give correct solution compared with MATLAB.

Localizing the distribution of transmembrance potential is the problem in concern to be solved. Solving the Inverse Problem in terms of transmembrance potentials (TMP) reduces the size of the problem to the third of the current density solution as stated earlier [26]. A low resolution heart geometry that is employed in the solution process by resizing the resolution one [45]. The conduction network of the heart [46,47] could be used in a later stage for modeling the propagation wavefront [48,49]. The number of voxels of the low resolution version of the heart is 960 elements (Figure 2) giving 960 independent variables in the case of TMP and 2880 (960 x 3) independent variables in the current density formulation.

The numbers of electrodes (dependent variables vector) that are used to measure the body surface potential are less than this number [50,51]. The 200, 100, 64 and 32 electrodes configurations were tested for the Inverse Problem solution (Figure 3). The 12-Lead ECG readings are not included because they were not measured with respect to the same reference, and the six electrograms which were measured with respect to WCT are too few to generate a good solution (the system will be extremely underdetermined).

As the number of the independent variables exceeds the number of the dependent variables, this problem is highly underdetermined and due to that, three regularization methods are used for localizing the primary sources locations inside the ventricles, namely, the Minimum Norm (MN), the Weighted Minimum Norm, (WMN) and the eLORETA. The matrix multiplication which is of O(N³) is the operation of the highest computing cost among other matrix operations where N is the number of elements in the matrix. This makes TMP’s formulation take the advantage over the current density formulation in reducing the computation cost.

The MN method has less matrix operations among all other methods including the WMN and the eLORETA, requiring less processing time. The eLORETA method requires the largest time as it contains several matrix multiplications. Another factor that affects the processing time is the number of electrodes readings (the dependent variables), and for example, with the MN method, 32 electrodes require a few minutes for solving the problem, while 200 electrodes it requires double that time.

2.6. Localization Problem of Boundary Points

Localizing the boundary points (including both external heart boundary, and blood volume boundary) that act as primary sources points is not accurate in most of cases, especially the ones located in the Endocardium of the ventricles. A novel approach is introduced to solve this problem, where small leakage currents between boundary points of the heart and neighboring non-excitable points are assumed. The solution then includes more points, which make the heart’s boundary points becomes sub-boundary points and so the problem is transferred to these added points which will never be activated. This provided a solution to the problem of localizing the activation points in some parts in the heart [24].

3. Results

For evaluating the correctness of the inverse solution, thousands of runs for each regularization method was performed. It was reported that the regularization methods are successful in determining single pacing point [42]. However, multiple pacing points were also included to evaluate the power of each of these regularization methods in determining such cases. The number of body surface electrodes is other factor that is used to in this evaluation as well. Some initial activation points are assigned to the heart model and the BSPM is being generated. The resultant BSPM is used as an input to recover the location of the activation points.

Selection of the regularization parameter for all of the used methods is accomplished by testing some small values and choice of the one that provide the best results. The weighting matrix of the eLORETA converges after few numbers of iterations (5 to 7). The evaluation of the efficiency of the each of these regularization methods is accomplished by comparing the recovered activation points from each of them with the original assigned activation points.

When an activation point (red pixel) is assigned as a volume source inside the heart Myocardium (Original), the Forward Model is used to generate the BSPM that corresponds to this activation. The generated BSPM is then used as the field input to the regularization method. It is clear from Figure 4 that these regularization methods are capable to determine an approximate location of the original activation point. The red color introduces positive component, the blue color introduces negative component and the green color introduces the zero component. The eLORETA is found to be the best method among these three methods that has excellent coloration with the original source points.

Two thousands different settings for each method were used to calculate the Inverse Problem solution. Each setting includes primary sources with random location ranging from one to five sources and body surface electrodes of 200, 100, 64 and 32 electrodes. Each 100 of these settings have similar number of primary sources and body surface electrodes, but with random locations for the primary sources. The MN results required almost 15 hours on a Pentium 4 – 1.7 GHz processor with 2GB RAM, while the WMN results took about 30 hours and the eLORETA method about 80 hours.

To measure the compatibility with the original pacing points (primary sources), the CC and the RE were calculated together with and the Mean and the Standard Deviation for the CC and the RE. Table 1 shows the mean and the standard deviation of each group of settings. It was found that the best localization results are produced when the number of primary sources equal to one and the more the number of electrodes the more the accurate the results (Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 andFigure 12). The eLORETA method was found to be the best method among the three methods for localizing the original source.

Finally, regularization methods are successful in determining a single primary source, but have less ability to localize multiple primary sources. The number of body surface readings affect the correctness of the regularization method, as the greater number or body surface electrodes produces better results. Identifying boundary pacing sites is done accurately by assuming a small leakage current to the surrounding surface and includes these points in the regularization process. The eLORETA method produces better results in the Inverse Problem solution than the MN or the WMN methods.

4. Conclusion

Solving the Inverse Problem of cardiac electrophysiology requires handling its inherently ill-posed nature through regularization techniques. This study has demonstrated that the MN, WMN, and eLORETA methods can successfully localize single internal sources based on body surface potentials, with eLORETA showing the highest accuracy and robustness across different test cases. While MN is computationally efficient, its accuracy is lower compared to eLORETA, especially in multi-source scenarios. The number of body surface electrodes significantly influences solution accuracy, with more electrodes yielding better localization results. Additionally, a novel boundary localization strategy—assuming small leakage currents to neighboring non-excitable points—significantly improves the accuracy in localizing activation points near the heart's boundary and septal regions. Despite limitations in multi-source localization, eLORETA offers a promising and more realistic approach to solving the cardiac Inverse Problem.