A method for magnetocardiography functional localization based on boundary element method and Nelder-Mead simplex algorithm.

BACKGROUND: The magnetocardiography (MCG) functional localization can transfer the biomagnetic signal to the electrical activity information inside the heart. The electrical activity is directly related to the physiological function of the heart.
METHODS: This study proposes a practical method for MCG functional localization based on the boundary element method (BEM) and the Nelder-Mead (NM) simplex algorithm. Single equivalent moving current dipole (SEMCD) is served as the equivalent cardiac source. The parameters of SEMCD are adapted using the NM simplex algorithm by fitting the measured MCG with the calculated MCG obtained based on BEM. The SEMCD parameters are solved in the sense that the difference between measured and calculated MCG is minimized.
RESULTS: The factors affecting the localization accuracy of this BEM-NM method were first explored with synthetic signals. Then, the results with real MCG signals show a good agreement between the SEMCD location and the region where ventricle depolarization starts, demonstrating the feasibility of this idea.
CONCLUSIONS: This is the first three-dimensional localization of the onset of ventricular depolarization with the BEM-NM method. The method is promising in the noninvasive localization of lesions for heart diseases.

INTRODUCTION

Better diagnostic tools for heart health need to be developed urgently since heart disease is the leading cause of death worldwide, causing one in four deaths in the United States (Murphy et al., 2018). Magnetocardiography (MCG) is the magnetic field measured near the thorax produced by the electrical activity of the human heart. MCG is considered a complementary means of detecting the physiological condition of the heart to the electrocardiogram (ECG).

Similar to ECG, the MCG signal can be used for the diagnosis of different heart diseases. The typical clinical applications of MCG include arrhythmogenic risk assessment, cardiac source location, and fetal cardiac health assessment (Fenici et al., 2005). Although the MCG and ECG share the same underlying source, the ionic currents in the heart, MCG still has advantages in specific application scenarios. MCG measurement avoids electrode contact and the associated noise. More importantly, MCG enables more precise localization of heart currents, which makes it promising for cardiac functional localization.

To achieve MCG functional localization, modeling the process of MCG generation with the help of computers is necessary. The modeling of MCG consists of solving two problems: the forward problem and the inverse problem.

The solution to the forward problem is to obtain the MCG with a known source model and proper modeling of the human torso (Malmivuo & Plonsey, 1995). On the one hand, the heart currents are the primary source. Several types of source templates are introduced to model the heart currents, including the equivalent current dipole, equivalent magnetic dipole, and current dipole distribution (Lu, 2010). On the other hand, heart currents form an electric field in the torso, resulting in volume currents contributing to the MCG. The finite element method (FEM) and the boundary element method (BEM) are used to calculate the volume currents distribution in the torso. The MCG signals are calculated with both types of sources considered.

The objective of the inverse problem is reconstructing the underlying electrical activity of the heart from the MCG, including cardiac source localization. The inverse problem in MCG was elaborated, and methods to solve the inverse problem were reviewed in Nenonen (1994). Recently, a method to solve the MCG inverse problem using the Bayesian approach was proposed in Bhat and Anitha (2020). However, there are still many unresolved difficulties for the MCG inverse problem. One of the most critical difficulties is the ill‐posed nature of the inverse problem, which means that with the measured MCG data, and there are kinds of solutions for the source. To guarantee a unique and consistent solution, the restrictions of source configuration are essential. The heart source templates can be applied to limit the number of source parameters.

Moreover, an effective solution to the inverse problem is generally based on the proper solving of the forward problems. The choice of either cardiac source or torso modeling method in the forward problem profoundly affects the inverse problem solution. The single equivalent moving current dipole (SEMCD) is the most commonly used source model since it is simple and suitable for describing the source when the current is constrained in a relatively small region. The SEMCD is applied in localizing focal activity for WPW (Wolff–Parkinson–White) syndrome, arrhythmogenic centers, and sites of origin of ventricular electrical activation. More studies based on SEMCD were reported on the study of ECG (Armoundas, Feldman, Mukkamala & Cohen, 2003; Armoundas, Feldman, Mukkamala, Mullen, et al., 2003; Fukuoka et al., 2006; Tysler & Svehlikova, 2013).

In MCG functional localization, the SEMCD was used to localize the pre‐excitation size of WPW syndrome patients (Bruder et al., 1994; Nenonen et al., 1991). Although the volume conductor was modeled as a realistic torso in both studies, the inhomogeneous conductivity was not considered. In another study, SEMCD location was inferred based on Levenberg–Marquardt (LM) algorithm (Mariyappa et al., 2012). However, the volume conductor effects were not considered in detail. The inhomogeneous conductivity of the volume conductor was partially taken into account (Chen et al., 2014), but only the boundaries of the heart and torso were considered.

In this article, a practical idea to achieve MCG functional localization is proposed with the cardiac source modeled as SEMCD and the volume conductor effect solved by the BEM. Estimation of the parameters of SEMCD is obtained by minimizing the difference between the calculated MCG solved by the forward problem and the measured one. The nonlinear least square optimization method used here is the NM simplex algorithm. The inhomogeneous conductivity is considered more detailed compared to previous studies. This approach is called the BEM–NM method in the following.

METHODS

A realistic forward model to solve the forward problem is necessary. To compromise between computational consumption and the accuracy of results, the BEM is applied. Meanwhile, the inverse problem is solved using the NM simplex algorithm. The source parameters are obtained by fitting the measured MCG with the calculated MCG solved by the forward model.

The forward model

In this part, the forward model, including the cardiac source model and volume conductor model, is presented in detail.

The cardiac source is modeled as SEMCD, that is, a single current dipole with varying magnitude, orientation, and position. Since the current dipole is thought of as an infinitesimally short length of wire segment that carries a current, the solution of SEMCD location is considered to be the center of the concentrated cardiac current activity (Nenonen, 1994). The SEMCD has six independent parameters. The three coordinates describe the position of SEMCD, while the three dipole moment components describe the magnitude and direction.

The influence of the volume currents on MCG is calculated numerically using the BEM. The torso is considered as a piecewise homogeneous volume conductor and divided into several regions. The conductivity inside each region is constant. The realization is based on a MATLAB library proposed in Stenroos et al. (2008).

For an observation point with a position vector , the magnetic field is

is the magnetic field due to the primary current density assumed to be SEMCD here and is the magnetic field brought by the volume currents. Only the volume currents through the boundary surface are considered to contribute to the MCG in the BEM. Thus, the is calculated as where is the position vector of the source point at the surface S and is the electrical potential of the source point, and are conductivities of the internal and external region of the surface, N is the number of surfaces.

As shown in (2), the MCG components are decided by the potential distribution at each surface. The calculation of the potential distribution needs to specify the basis function and weight function and minimize the residual after discretization. For the details of calculating the potential distribution, please refer to Stenroos et al. (2008). In this study, the linear basis functions are chosen, and the Dirac δ weight function is applied.

Six boundary surfaces with realistic shapes are considered. For numerical solutions, the surfaces are discretized to nodes and triangular elements. The data of the surfaces are downloaded from a website: www.ecgsim.org. The material is described in van Oosterom and Oostendorp (2004). The six surfaces include interfaces of the torso, ventricle, lungs, and ventricular cavities. The volume conductor model is illustrated in Figure 1.

FIGURE 1

The torso volume conductor model and measurement points of MCG

The measurement points of MCG form a 9 × 9 array on a plane 30mm away from the front plane of the torso, also shown in Figure 1. The distance between adjacent nodes is also 30mm. To ensure consistency, the forward model calculates the MCG at the same locations. All the MCG signals considered in this article are the components perpendicular to this measurement plane. The conductivity values of different regions are listed in Table 1, mainly referring to the data in Ramon et al. (1998).

TABLE 1

The conductivity values of different regions in the forward model

Region	Ventricle	Lungs	Torso	Ventricle cavities	Air outside body
Conductivity (S/m)	0.239	0.067	0.033	0.649	0

Solving the inverse problem

This part aims to estimate the six parameters of SEMCD at a given moment with the measured MCG signal.

The typical approaches to estimate the source parameters include the LM algorithm (Marquardt, 1963), Nelder–Mead (NM) simplex method (Lagarias et al., 1998). LM algorithm is more commonly used. However, the work using the LM algorithm, like Nenonen et al. (1991) and Mariyappa et al. (2012), considers at most one boundary surface. As the number of boundary surfaces grows, the relationship between the dipole parameters and the magnetic field becomes more complex. In addition, a derivative is required for the LM algorithm, which can be numerically approximated. Therefore, the implementation of the LM algorithm will be more difficult due to derivative operations. The NM simplex algorithm is chosen here since it only requires calculating the function value in iteration.

The flowchart for applying the NM algorithm to solve the inverse problem is shown in Figure 2a. The NM algorithm is to find the optimal solution from one initial guess to minimize the function value. The function value corresponds to the sum of squares of the difference between calculated and measured MCG in this work. The process of function value calculation is illustrated in Figure 2b. Moreover, the parameters of SEMCD can be regarded as coordinates of a point in a six‐dimensional space. A simplex with seven points is first generated and then modified in the following iterations. Please refer to Lagarias et al. (1998) for the detailed process of NM algorithm implementation.

FIGURE 2

(a) The flowchart of solving the parameters of SEMCD using NM algorithm, (b) The flowchart of calculating the function value, noted as f(x) which is the sum of squares of the difference between measured and calculated MCG

RESULTS

In this section, the influence of the initial guess and noise on localization accuracy will be given based on synthetic MCG signal, and the validity of the BEM–NM method is verified based on measured MCG signal.

Synthetic MCG signal

The synthetic MCG signal is calculated by the forward model in part 2.1, which is also used in the inverse problem. In other words, the model is perfect in this part. The SEMCD with known parameters serves as a source, and the BEM–NM method is applied to search it. The objective was to evaluate the probability of finding an accurate solution using the NM algorithm in the ideal modeling state and how the initial guess and noise affect it.

The location of SEMCD is generated within the region of the ventricular. Three points are chosen, noted as Point 1, Point 2, and Point 3, which also refer to SEMCD in the following. Point 1 is located in the ventricular septum. Points 2 and 3 are located on the left and right free wall of the ventricle, respectively. The dipole moment components are generated randomly in the range of [−10,10] µAm to ensure the synthetic MCG is relatively close to the actual MCG in terms of magnitude. The specific parameters of the three SEMCDs are shown in Table 2.

TABLE 2

The parameters of the three SEMCDs used in this article

	x/mm	y/mm	z/mm	px /µAm	py /µAm	pz /µAm
Point 1	49.9	36.6	−15.5	−9.3111	−1.2251	−2.3688
Point 2	39.5	75.1	5.3	−2.1555	3.1096	−6.5763
Point 3	83.2	24.3	−10.3	8.2675	2.6472	−8.0492

Abbreviation: SEMCD, single equivalent moving current dipole.

Then, the synthetic MCG data of the three SEMCDs are calculated using the forward model. The MCG will be added noise when the influence of noise is considered. After that, the synthetic MCG data are treated as the “measured MCG.” The distance d between the position found by the BEM–NM method and the actual position is obtained to evaluate localization accuracy.

Range of the initial guess

The initial guess is essential for the NM simplex algorithm. The effect of the initial guessing range on the localization accuracy is considered here. A parameter e is presented to describe the range of the initial guess. Specifically, when e is identified, the positions of all initial guesses are located in the area with the coordinates of x, y, z range from [x 0 − e, x 0 + e], [y 0 − e, y 0 + e], [z 0 − e, z 0 + e], respectively, where (x 0, y 0, z 0) is the coordinates of the true SEMCD. For each SEMCD and e, one hundred independent trials were conducted for evaluation. The position of the initial guess is generated randomly around the true SEMCD, and the dipole moment parameters are also obtained randomly in the same way as the counterpart of the true SEMCD.

For each trial, if d is smaller than 10 mm, the SEMCD is thought to be found roughly. If d is smaller than 0.1 mm, the SEMCD is thought to be found exactly. The rate of finding is obtained by calculating the percentage of eligible trials. The rate of rough finding and the rate of exact finding in the case of different SEMCD and e are shown in Figure 3a. The rate of rough finding decreases as e increases, indicating a larger initial guess range reduces the possibility of rough finding. Here, the probability of finding roughly is greatest when e is set to be 2 or 3 cm.

FIGURE 3

(a) The rate of rough finding in the case of different SEMCD and e, (b) The rate of exact finding in the case of different SEMCD and e. e is the description of the initial guess range

Although the rate of exact finding still decreases with increasing e, the rate of exact finding and the rate of rough finding are not positively correlated. The SEMCD with a more significant rate of rough finding does not have a larger rate of exact finding. Besides, the rate of exact finding is relatively low, no more than 25% under all experimental conditions in this study. Even for perfect modeling, the probability of exactly find the actual location of SEMCD is still low. In the absence of information about SEMCD, it is necessary to ensure that the initial guess is close to the actual value and to increase the number of searches.

Noise in MCG signal

Here, the influence of noise on the possibility of rough finding is explored by adding Gaussian white noise to synthetic MCG signals. The inverse problem is solved based on noisy MCG signals. Figure 4 shows the rate of rough finding with additive noise of different amplitudes. The signal‐to‐noise ratio (SNR) of the noisy MCG is given by where x(n) is the clean synthetic MCG signal calculated in various measurement points and u(n) is the Gaussian white noise, N is the number of measurement points. Each rate is also obtained after 100 trials. For all SEMCD, the adding noise does not significantly reduce the likelihood of finding the correct solution when the SNR is larger than 15DB. However, the possibility of rough finding shows a clear downward trend when the SNR decreases below 15DB.

FIGURE 4

The rate of rough finding in the case of different SEMCD and noise amplitude

Real MCG signal

The real MCG signals of a normal male are downloaded from www.ecgsim.org, measured at 9 × 9 points on the plane parallel to the torso front surface. The signal's sampling frequency is 1000Hz, and the measuring duration is 0.1 s of the QRS interval.

The SEMCD and BEM model used is simplified compared with the actual situation. To ensure the applicability of SEMCD, only the signals at the beginning of the QRS interval are used. Since ventricular depolarization has just begun, the activation region is small, which may satisfy the assumption of SEMCD.

First, the SNR of the MCG needs to be assessed to choose the moment of the QRS interval beginning. The P(t) is a description of the measured MCG amplitude at time t, calculated by . v(n) is the measured MCG of measurement point with the index n at time t and N is the total number of measurement points. The P(t) of t = 1, 2, …, 9 msec is listed in Table 3.

TABLE 3

Measured MCG energy within milliseconds of the start of the QRS interval

t/msec	1	2	3	4	5	6	7	8	9
P(t)/pT 2	0.7554	0.9014	1.422	2.48	4.3717	7.4905	12.8749	20.3378	29.4503

Abbreviations: MCG, magnetocardiography; msec, milliseconds.

When t = 1 or 2 msec, P(t) is mainly brought by background noise. Consequently, the components of the background noise in P(t) are around 0.8 pT 2. Then, the amplitude of the MCG signal increases and the value of P(t) keeps growing. Considering the SNR should be as large as possible to make sure localization accuracy. If the noise energy is considered not to change rapidly, 7 msec is the first instant meeting the minimum requirement of SNR, which is 10 DB. To compromise with the requirement to be as early as possible to ensure the applicability of SEMCD, the beginning of the QRS interval is finally chosen to be 7 msec.

Then, a series of initial guesses need to be identified. Since there is no prior information for the origin of ventricular activation here, the SEMCD is searched throughout the ventricle. The 3D space of ventricular is divided into 100 cubes with 2 cm sides, and 30 points are chosen randomly inside each cube as the locations of initial guesses. All the dipole moment components are generated randomly in the range of [−5,5] µAm. Totally, 3000 sets of initial guesses are obtained in this way.

Finally, the program completed 3000 runs with different initial guesses. The smallest function value is 1.1898 pT 2, and the parameters of the SEMCD are (34.91, 34.10, −18.02) mm for the location coordinate and (−7.13, 6.30, 5.80) µAm for dipole moment. The position of this SEMCD in the ventricle is shown in Figure 5. The SEMCD is located near the left surface center of the interventricular septum, closer to the base and the posterior wall. The position may correspond to the region where activation of the ventricle starts.

FIGURE 5

The position of the found SEMCD based on measured MCG signal, corresponding to the blue‐violet dot in the graph

Another research is referred to verify the result. The propagation of excitation in normal human hearts was studied in (Durrer et al., 1970), based on the measurement of isolated human hearts. The SEMCD solved by our BEM–NM method is in reasonable agreement with one location of earliest activation in the ventricular septum. The fit of the position of the normal ventricular excitation origin obtained by two completely different approaches confirms the feasibility of the BEM–NM method for noninvasive localization based on measured MCG signals.

DISCUSSIONS

The results based on real MCG signals show that the BEM–NM method can obtain the heart currents information with only the MCG signal and the internal structure of the torso, indicating the feasibility of noninvasive cardiac function monitoring. The reasonable results derive from a more realistic forward model and extensive trials of initial parameters. Cardiac current activity at the onset of the QRS interval partially in line with the assumption of SEMCD is also quite important.

However, there are still many limitations of using the BEM–NM method for MCG functional localization. The primary limitations come from the inaccuracy of the forward modeling, the distortion of the measured signal, and the calculating consumption.

Firstly, the forward modeling error comprises two parts, the cardiac model error and the volume conductor model error.

Although the SEMCD with a simple structure can describe a single and concentrated current activity in a small spatial area, it is not suitable as an equivalence source when the range of current activity is extended and spread over multiple spatial locations. As a result, the application of SEMCD lies in locating the current activity in a single concentration area, such as the current focal activity of WPW syndrome and arrhythmogenic. Due to the singular nature of the dipole, numerical solution difficulties arise when the dipole is very close to the potential observation point at the boundaries, which also needs to be handled appropriately.

Moreover, the torso is considered a piecewise homogeneous volume conductor in this study. However, such an assumption does not correspond to the actual situation. For example, the myocardium is anisotropic. As for the spatial information of the boundary surface, an individualized volume conductor model can be constructed from MR‐scans. However, the boundary surface is dynamically changing since breathing and heartbeat. All these factors need to be considered to refine the forward model. The good news is that the MCG is less sensitive to the conductivity properties of the torso comparing with ECG (Mäkijärvi et al., 2010), implying that the volume modeling can be relatively simplified.

Secondly, the distortion of the signal is considered here. The MCG is weak compared with ambient magnetic noise. A high‐quality signal is necessary for achieving MCG functional localization. Figure 4 shows that once the SNR falls below 15DB, the likelihood of finding the SEMCD position decreases. Proper handling of residual noise in the measured signal is essential.

Thirdly, the number of initial guesses is large, thus increasing the computational effort. The implementation of the NM simplex algorithm makes more complex forward models possible. However, the probability of finding the correct solution using the NM simplex algorithm is relatively low, which still leaves room for improvement.

Despite several limitations exist, the value of the BEM–NM method is still emphasized. The BEM–NM method can link the physiological functional information with the structural one and convert information from MCG signals into the spatial location of intracardiac currents. This property makes it promising for pre‐surgical planning of lesion localization.

Besides, some limitations can be overcome while maintaining the main idea of this method. The specific models, such as the volume conductor model, can be further improved. A complete heart boundary surface, including atria and other interfaces in the torso, can be introduced to enable more accurate forward models. The conductivity of each region can be adjusted to a more accurate value simply by modifying the corresponding parameters. The high‐quality MCG signal measured at more positions, for instance, adding a measurement plane behind the body, will also help obtain more accurate results.

As a generalizable localization method for physiologically focal electrical activity, the BEM–NM method itself is not limited to MCG localization and has applicability for magnetoencephalography source localization after adjusting the forward model.

The single dipole model and the BEM method were already applied in the ECG inverse solution (Tysler & Svehlikova, 2013). The differences between the study by Tysler et al. and this work are listed here: (i) The dipole position was predefined in the study of Tysler et al. and only three dipole moment components were computed. However, the SEMCD position is unknown and to be sought for the BEM–NM method, which leaves more freedom for the solutions. (ii) The difference of integral body surface potential maps served as the input data of the inverse localization proposed by Tysler et al. while the difference of the MCG maps of a particular moment is the equivalent in this study. Thus, two kinds of input data correspond to different time ranges. (iii) The NM simplex algorithm is first applied to search the inverse solution in the BEM–NM method.

CONCLUSIONS

A method to achieve MCG functional localization by fitting the measured MCG with the MCG calculated by the forward model is proposed in this article. The cardiac source is modeled as the SEMCD, and the volume effect is calculated using the BEM. The parameters of the SEMCD are modified based on the NM simplex algorithm to fit the measured MCG. Taking advantage of the NM simplex algorithm, a more realistic forward model is introduced in this work.

The forward model calculates the synthetic MCG signal with known SEMCDs. The parameters of SEMCDs are solved with synthetic MCG to explore the influence of the initial guess and the noise on the localization accuracy. Then, the BEM–NM method is applied to search the currents of the normal human heart with the measured MCG signal. One of the original excitation areas of the ventricle is found in the early time of the QRS complex.

It is the first time to localize the onset of ventricular excitation based on the BEM and NM algorithm combination, showing the reliability of this method to localize focal currents of the heart in a noninvasive way. Thus, the method is promising in localizing lesions of WPW syndrome and arrhythmogenic.

CONFLICTS OF INTEREST

The authors have stated that there are no conflicts of interest in connection with this article.

ETHICS

The study does not directly involve clinical trials. The data of the human participant used in this study are publicly available from www.ecgsim.org.