Localization of magnetocardiographic sources for myocardial infarction cases using deterministic and Bayesian approaches

In this paper, the inverse problems of cardiac sources using analytical and probabilistic methods are solved and discussed. The standard Tikhonov regularization technique is solved initially to estimate the under-determined heart surface potentials from Magnetocardiographic (MCG) signals. The results of the deterministic method subjected to noise in the measurements are discussed and compared with the probabilistic models. Hierarchical Bayesian modeling with fixed Gaussian prior is employed to quantify the uncertainties in source reconstructions. A novel application of Variational Bayesian inference approach has been presented to estimate the heart sources. The reconstruction results of Variational Bayesian model with non-stationary priors are compared with solutions of simplistic Bayesian approach; and the performances are evaluated using Root Mean Square Error (RMSE) and correlation co-efficient metrics. The Bayesian solutions in the study are also extended to localize the MCG sources for two types of Myocardial infarction cases.

www.nature.com/scientificreports/ One of the objectives of this paper is to reconstruct epicardial sources using Bayesian frameworks to address the uncertainties that occur in the deterministic approach. The other objective is to apply Hierarchical models using the simplistic Bayesian and the novel (in this domain) Variational Bayesian methods in solving the inverse problem of epicardial activities. Further, the proposed inverse models have been analyzed based on noise-free and noisy MCG signals of normal and Myocardial Infarction cases.

Methodology
In this section, the forward model of MCG has been described which discusses the source models, the explanation of the Lead Field/Transfer Matrix followed by the mathematical depiction of the forward model. In the next subsection, the inverse problem is described and solved using the standard Tikhonov method, the Hierarchical Bayesian method and the Variational Bayesian Linear Regression method. Two different heart models have been used, one for the forward problem and the other to solve the inverse problem in order to avoid the 'inverse crime' .
Forward problem. Source model. The source model is designed based on the prior knowledge of the electrophysiological nature of the heart. The source parameters assumed in the current study consisted of a ventricular surface model extracted from ECGsim 25 , covered by Q = 257 nodal locations with each discretized node assigned to epicardial potentials. The thorax model of size (25,45,45) cm consisting of 300 nodes with 596 tetrahedral meshes is considered in the study. The heart mesh was placed inside the thorax model at a location (0.3, 0.3,0) cm that lay in the realistic anatomical position between the lungs and behind the sternum 13 (geometrical models were assembled in SCIRun software 26,27 ). The volume of the thorax were filled with assumed conductivity σ values of 0.6, 0.04 and 0.2 S/m for the ventricles, the lungs and the torso, respectively 28 .
The discretized model of the heart is defined as the vector of q = (q 1 , q 2 , . . . q Q ) sources. The ventricular surface nodes at q points oriented in â q directions have been assigned with epicardial potentials. In this study, two types of Myocardial infarction (MI) cases are simulated using ECGsim software 25 . The activations are created at a node with 1.5 cm potential spread (shown as patch in the inset of Fig. 1a) near the anterior left ventricular wall by varying the transmembrane potentials waves: (1) ST elevated MI is simulated by decreasing the magnitudes of the transmembrane potentials wave by 48% 29 , and (2) the time instants of the transmembrane potentials wave was shortened to get increased T peak MI case. The amplitudes of the abnormal epicardial potentials are shown in Fig. 1c.
The current inverse study is tested on another heart model (all simulations were designed and executed in SCIRUN/BioPSE software 26 ) and the models were imported from ECGSim software 25 as source with the same diseased cases as shown in Fig. 1b. The model consisted of 585 discretized nodes with 1156 triangular surface mesh elements. The source sampling in this structure is more than the previous heart model. The reason for performing the inverse study on this test model is explained in the succeeding sections.
Transfer matrix. After modeling the structural and functional prior assumptions of the heart aligned with the torso, the spatial detector vector m = (m 1 , m 2 , . . . m M ) are placed in parallel to the thorax. The MCG detectors are collected from a set of 9 × 9 observation points ( M = 81 ) with 3 cm intervals in between (available in ECGsim 25 ). Now, a transfer matrix is designed to construct the spatial sensitivity between Q sources and M detectors. The generic field is recorded by placing unit current dipole vectors â q at the nodes of the myocardium  www.nature.com/scientificreports/ with the help of Biot-Savarts law. The magnetic field observed at sensor array (position vector in three dimensional Cartesian coordinates denoted by r m = (r 1 , r 2 , r 3 , . . . r M ) ) due to unit dipoles at r ′ q = (r ′ 1 , r ′ 2 , r ′ 3 , . . . r ′ Q ) heart points are appended for all the sources to construct a lead field spatial tensor of dimension M × Q × 3: The tensor L(r m , �r ′ q ,â q �) is expressed as: The elements of the constructed spatial transfer tensor L(r m , has the dimension R M×Q×3 ; includes the locations and orientations of the sources, denoted as L in further discussions. Forward linear model. The forward magnetic fields (MCG) with dimensions R M×t×3 are defined as an array of M detector coils, B f (r m , t) = [B fx (r m , t), B fy (r m , t), B fz (r m , t)] which are computed from the lead field matrix L and prior scalar epicardial potential distributions s(�r ′ q ,â q �, t) ∈ R Q×t . In further discussions throughout the article, B f (r m , t) and s(�r ′ q ,â q �, t) are denoted as B f and s, respectively. The following equation represents the linear model to construct the forward magnetic field: The MCG waves are simulated for normal and abnormal cases and considered as true observed signals in the inverse problem. But, the usage of the same model to obtain the inverse solution leads to a problem called 'inverse crime' i.e., the results may yield best estimates which could be equivalent to the truth. One can overcome this problem by adding noise to the linear forward problem 30 . However, the actual modeling error of using the same models (unrealistic scenarios) in the inverse estimations is still a 'crime' and the results may end up providing over-estimated solutions 31,32 .
In our study, it has been attempted to solve the 'inverse crime' . First the forward calculated magnetic field intensity (subjecting to uncertainties) is distorted by Gaussian noise with substantial amount of Signal-to-Noise Ratio (SNR) levels (ranging from 6 to 16 dB) and then the hidden potentials are determined. However as explained before, this does not solve the inverse crime so we further utilized the forward model of the different heart structure and the forward calculated magnetic field intensity obtained from the first source to estimate the unknowns. Figure 2a,b represent the MCG maps at mid ST instants of normal and ST elevated MI case from the first heart model, respectively. Similarly, MCG fields at increased T peak maps of normal and abnormal of model 1 are illustrated in Fig. 2d,e, respectively. The forward calculated magnetic field intensity from the second source model (shown in Fig. 1b) are mapped in Fig. 2c,f depicting the ST elevated and increased T instants, respectively. The variations in the amplitudes of the field intensities can be visualized in the maps.
The linear forward field ( B f ) representing the cardiac cycle from t = 1 s to t = T s, subjected to noise n(r m , t) denoted as n is formulated to get B fn and is expressed as: The noise n of time frame Ts is derived from the SNR dB with signal power P s = 1 T T t=1 |B f | 2 dt and noise power P n = 1 T T t=1 |n| 2 dt is given by: Since the noise power is not known, the parameter of a specific noise amplitude is formed using the above equation that reduces to for a range of SNR from 6 to 16 dB. The obtained noise power amplitude P n is then multiplied with the random noise sequence of the length similar to the signal's time points.
(1) L(r m , �r ′ q ,â q �) = µ 0 4π Deterministic approach: Tikhonov regularization. The unknown source activations are estimated with the help of transfer matrix L and from known observations B(r m , t) denoted as B using squared error function 14 . The general cost function to be minimized is computed using: Since the system is under-determined ( M < Q ), it causes over-fitting in the function. In order to overcome this problem, a penalty term called solution norm (L2 norm) is introduced that contains squared magnitude of the epicardial source weights. The standard Tikhonov (L2 norm) regularization implemented in the study, is expressed as: where R is the regularization parameter which controls and helps in fitting the sources by minimizing the cost function.
Hierarchical Bayesian framework. In the Bayesian estimation, the magnitudes of unknown source activities are estimated, which are priorly assumed to be random with known prior distributions. The joint probability density function (pdf) of the magnetic field data B and source activities s are assumed to be normally distributed. In this, p(s) is the prior distribution of s and p(B | s) is the likelihood function of MCG conditioned on the sources s. The posterior distribution to estimate the activities is written as: where p(B) in the denominator is the normalizing factor of observed signals.
The conditional probability of the complete MCG cardiac cycle from ith detector B it = B(r m = r i , t) with time points ( t = 1, . . . , T ) given the whole latents of the epicardial potentials s qt = s(�r ′ q ,â q �, t) considering the lead field ( L ) is given by, We assume the noise n (Eq. 3) is Gaussian with zero mean and is i.i.d across time expressed as www.nature.com/scientificreports/ Here, is a diagonal precision matrix in which the diagonal entries are equal to the inverse of the noise variances for the corresponding observed MCG data.
The forward likelihood over the MCG at M detectors in an instant is defined as: Similarly, the prior distribution of the epicardials s jt for a jth source is assumed to be Gaussian and i.i.d across time and the model for the entire time series of s jt ( t = 1, . . . T ) is given by Here, is the precision or covariance matrix which contains the inverse of the source variances ( σ 2 s ). The explicit form of the Gaussian distribution of the epicardial priors at Q sources in a specific instant is given by Two scalar parameters α and β called hyperparameters are introduced to control the distribution of the parameter s. The hyperparameters supporting the precision matrices are assumed as � = αI and � = βI for prior source and noise covariances, respectively.
These precision hyperparameters decide the variance estimates of the posterior distribution. The posterior distribution of the epicardial weights (Gaussian) for any source at r ′ q is formed as where the mean and variance are given by The values of the hyperparameters are unknown and are derived from the data. This can be done by introducing prior distributions over the hyperparameters α and β , and predicting the posteriors by marginalizing with respect to these hyperparameters and epicardial weights s.
We employed the method of evidence approximation, since the complete marginalization over s, α , and β by integrating analytically is not tractable 21,30 .
If the hyperpriors are modeled over α and β , the predictive distribution by marginalizing over s, α , and β gives If the hyperposteriors p(α, β | B) are most probable around the values α and β , the posterior distribution of epicardials is computed by marginalizing over s with α and β . From Bayes' theorem, hyperposteriors are expressed as Thus, by maximizing the p(B | α, β) provides the evidence for α and β. The marginal likelihood of the MCG data p(B | α, β) is obtained by integrating over the sources s, From Eqs. (9) and (12), the evidence function is written in the form where (10) n ∼ N (n | 0, � −1 ) www.nature.com/scientificreports/ The expression exp{−E(s)}ds cannot be obtained analytically and practical solution of the same over s is larger. To simplify this, Taylor expansion of E(s) is considered around the minimum value by retaining up to the second order 23 . The expansion of E(s) around its minimum value is given by where s MP denotes the most probable solution (epicardial potentials). The second term in the expansion is the minimum value and can be discarded or equated to zero. The second order derivative function is and By differentiating the log evidence of Eq. (18) with respect to α , we get, are the eigen values of βL ⊤ L. Similarly, by maximizing the log evidence with respect to β , results in The Algorithm 1 explains the source estimation procedure by updating the hyper-parameters in iterative manner until it reaches to a convergence criteria.
Variational Bayesian linear regression. The Variational Bayesian linear regression method is a probabilistic algorithm previously used in solving other applications as seen in 33 . So in this paper, this method has been explored for the first time to localize the sources from MCG signals.
Bayesian model. The linear model explained in the forward model (Eq. 2) that assumes linear relation between Q dimensional s sources and M dimensional MCG observations with independent noise at the detectors is considered in the inverse procedure. The magnetic field intensity likelihood of the observations B (similar to the previous approach) is assumed with constant-variance Gaussian noise distribution � = βI 20 , defined as: The graphical representation of Bayesian models is illustrated in Fig. 3a,b to understand the relationships between variables.
Priors and hyper-priors. In the previous section (simplistic Bayesian approach), the prior source weights in the linear model were assumed to be Gaussian distribution. Due to this stationary prior, it is possible to estimate www.nature.com/scientificreports/ marginal likelihood and obtain posteriors 19 . It is also important to extract the sparser characteristics of the unknown signals such as region of abnormal spread estimations in the epicardial potentials that appear to be smoother solutions in stationary prior. To overcome this problem, a non-stationary conjugate normal inversegamma distribution is assumed on the prior sources s and variance � = βI , parametrized by hyper prior � = αI graphically represented in Fig. 3. The prior source distribution is modeled as: In this prior, β acts as inverse variance on s with zero-mean. Due to imposition of Gamma on β , β −1 is defined as the inverse-gamma function with shape a 0 and scale b 0 and is expressed as: where η = β 2 (αs ⊤ s) + 2b 0 The convergence criteria of the algorithm 1 is evaluated by checking the change in updates of hyperparameters. This is done by setting a stop threshold ǫ = 10 −3 value to track the difference between the hyperparameter updates. Then the algorithm is constrained to stop at iteration j if P(j) < ǫ where: Variational inference. The final step is to estimate the cardiac sources in terms of variational posteriors termed as p(s, β, α | B) . Since it is difficult to obtain the closed form of posteriors, a distribution g(s, β, α) is introduced in the model to solve it as variational approximation. With the help of Kullback-Leibler (KL) divergence, the difference between g(s, β, α) and posterior p(s, β, α | B) is minimized, defined as: where g = g(s, β, α) and p = p(s, β, α | B) . The distribution g(·) is computed by minimizing the KL divergence: The main idea is to find the distribution g that is closer to the distribution p. The divergence in Eq. (27) reduces to: Now, with the help of known ln p(B) , the variational posteriors g(s, β) and g(α) are estimated by maximizing the variational lower bound L(g) 23 (equivalent to minimizing the KL divergence function in Eq. (28)).
The variational posteriors for s and β , with fixed g(α) is given by: The posterior on sources with variance g(s, β) reduces to: The hyperparameter update is terminated when the variational lower bound L(g) remains unchanged for more than 0.001% between the two successive iterations.

Results and discussion
In this section, the performance of algorithms 1 and 2 for epicardial source reconstructions are compared with that of deterministic method in noise-free and noisy conditions. The inverse algorithms are computed in MAT-LAB and the results are visualized in SCIRun software 25,30 . The true epicardial maps were generated using the geometry of model 1 (mappings at t= 85ms and t= 125 ms shown in Fig. 4A.i,B.i, respectively) while the transfer matrix generated from the geometry of model 2 was used to solve the inverse problem. Due to this reason the inverse crime was avoided however it led to some significant modeling errors. The results obtained from the Tikhonov approach in noiseless condition ranged from − 1.4 to 2.9 mV; where the region of spread showed a similar kind of ruptured information as that of true potentials (shown in Fig. 4A.a). The magnitudes of the spread regions of the Bayesian inference ( Fig. 4A.b) yielded a reasonable amount of improvement in the results than the Tikhonov method and the solutions ranged from − 1.4 to 3 mV in noise-free case. The main study was to visualize and evaluate the results for uncertainty conditions. This is performed by adding a noise of 8 dB SNR to the forward magnetic field intensity and solving the inverse problems. It can be observed from the Fig. 4A.d that the Tikhonov method in the noisy cases produced the contorted reconstruction spread at the diseased node. Figure 4A.e represents the Bayesian estimation results where the distributed area surrounding the diseased node was noticeable with magnitude ranges [− 3.8 3.2] mV. The VBLR approach performed well in both noise-free and noisy conditions evident from the maps shown in Fig. 4A.c,A.f, respectively. Similarly, Fig. 4B.a-c describe the estimations solved by Tikhonov, Bayesian inference, (31) ln g(s, β) = ln N (s |ŝ, β −1Ĉ ) · Gam(β |â,b) (32) ln g(α) = E s,β ln p(s, β | α) + ln p(α) + constant www.nature.com/scientificreports/ and VBLR methods, respectively at t = 125 ms. Also, the inverse solutions obtained at this instant from Tikhonov, Bayesian inference, and VBLR techniques due to the impact of noises in the forward field are shown in Fig. 4B.d-f, respectively. Figure 5 depicts the time series reconstructions of the heart surface potentials at the desired node index. The potentials activated in the model 1 (shown in Fig. 5A.i,B.i) are considered as the reference truth values for ST elevated and increased T cases, respectively. As compared to the true values, the amplitudes reconstructed by the VBLR were less than the Bayesian solutions (visually-for elevated ST case shown in Fig. 5A.iii,iv); but the temporal RMSEs of the VBLR method outperformed the Bayesian method in both noisy and noise-free conditions which are later discussed in Table 2.
It can be observed that the Tikhonov approach attempted to produce good solution as shown in Fig. 5A.v,B.v visually, but the amplitudes are not satisfactory which still showed the elements of noise present in the estimated results. This shows that uncertainties in the forward magnetic field intensity are less handled by deterministic methods. The solution of the Bayesian inference (Fig. 5A.vi,B.vi) traced out better estimates than the Tikhonov method. The temporal reconstruction of the potentials showed remarkably better results in VBLR method even under 8 dB SNR (Fig. 5A.vii,B.vii) than the other methods.  Fig. 6. As the regularization parameter is varied, the graphical tool displays the trade-off occurring between the regularized solution and its fit to the observed data. In the work, the regularization parameter was varied from = 1E − 10 to = 10 with 100 linearly spaced vectors in between them. The optimal regularizer can be found at a corner of the curve isolating the horizontal and vertical lines in the log − log scale. The vertical part of the curve in Fig. 6c corresponds to the solutions with small values. As the values are tuned to the small values, the solution norm starts increasing, thereby decreasing the norm of the residuals as shown in the Fig. 6a,b (for = 0.01 , the values of �s� = 401.582 , and the residual norm found was �B − (Ls + n)� = 4.8466).
The solution norm decreases slightly for large values causing the solution to over regularize (and increases the residual norm); corresponds to the horizontal part of the L-curve. The corner in Fig. 6c is considered as an optimum value = 1.6681 that yields a good solution fit for the latent t = 85 ms (mid ST elevated point). www.nature.com/scientificreports/ The inverse algorithms are tested by applying different levels of noise (SNR ranging from 6 to 16 dB) to the MCG system and Root Mean Square Error (RMSE) between the true (assumed) and estimated potentials are computed for quantitative evaluations and is defined as: where s q denotes the assumed epicardial potentials for two different heart models with geometrical nodes Q 1 and Q 2 , and ŝ q represents the estimated potentials at the same spatial points Q 1 and Q 2 .  It can be seen that both the Bayesian and Tikhonov methods take time to adjust themselves from noisy to less noisy conditions however VBLR always maintains similar values and hence it can be said that it works best even in noisy conditions. The VBLR method thus outperforms the other methods by showing good sRMSE s of 0.5428 mV even in worst cases from 6 to 8 dB and reached to a constant sRMSE s of 0.538 mV after 8 dB.
Similarly, the VBLR showed good estimations in the second case at t = 125 ms, with sRMSE s of 0.347 mV for all the noise levels whereas the hierarchical Bayesian inference provided negligibly better solutions than the Tikhonov approach from 6 to 8 dB and reached constant sRMSE s of 0.3672 mV after 8 dB in both the methods. Thus we can conclude that the VBLR method provides far better estimates of sRMSE values than Tikhonov and Bayesian methods in noisy cases.
Another metric called Correlation Co-efficient (CC) used in the study is tabulated in Table 1. The results of CC obtained for ST elevated case reached to 0.14 which is a very low value and it got improved by 1.7 times in the Bayesian approaches under noisy conditions. Along with this, the temporal RMSEs are evaluated to check the performances of the estimated epicardials with respect time formulated as:   Fig. 8e). The results of the epicardial estimations for different conditions are discussed in Fig. 4. Here, we discuss the MCG forward reconstructions from the estimated heart activities (model 2) and the results are mapped on to the detector planes for field visualization. It was found that there were only small deviations in magnitudes of reconstructions from deterministic and Bayesian methods in noise-free conditions as shown in Fig. 8b-d. Figure 8f illustrates the reconstruction results of Tikhonov regularization (8 dB SNR) which showed reduction in amplitude range [− 84 27] pT. The MCG waves reconstructed from Bayesian model handled uncertainties and amplitudes appeared near to the range of true maps as shown in Fig. 8g. However, the VBLR reconstructions (Fig. 8h) showed an improvement in the map than Bayesian method and the distribution range found was [− 51 22] pT.
The optimal values of hyperparameters are updated in the simplistic Bayesian modeling by maximizing the measurement likelihood and ratio corresponding to the regularization parameter. The final RMSE measure of Bayesian method in elevated ST case converged to 0.5463 mV, for hyperparameters α = 6.5676 and β = 1.907E15 in noise-free condition. For 6dB, the RMSE reached to 2.5081 mV for α = 0.0232 and β = 5.38E12. But in VBLR method, the distribution function q(α) known as hyperpriors is approximated using KL divergence. The hyperposterior mean ( E α (α) =ˆcd ) obtained was 0.0023 and ( E β,s (βs ⊤ s) = 780.3) at 6 dB (lowest SNR) in ST elevated case. To evaluate the hypers, the optimal ratios of α β are captured at different SNRs along with the converged RMSEs in the hierarchical Bayesian method as tabulated in Table 2. Similarly the means of hyperposteriors are obtained in VBLR technique to study the behaviour at uncertainty conditions. It was observed that there was slight improvement in the RMSE of VBLR method than the Bayesian approach. The nature of variations of α β (hierarchical Bayesian) and E α (α) (VBLR) with respect to iterations showed similar activities; where the values of both ratio and hyperparameters increase in noisy signals, whereas decrease in noiseless measurements.

Conclusion
In this paper, we presented inverse problems of MCG using deterministic and probabilistic methods employed on two types of MI cases. The results are estimated in terms of epicardial maps from the observed MCG signals in noisy and noiseless conditions. In competition to the simplistic Bayesian technique, a new method called Variational Bayesian Linear approach is applied which contained non-stationary normal inverse-gamma priors that used KL divergence to approximate the posteriors. The other drawback of selecting the optimal parameter is addressed in Tikhonov regularization. It is shown that the problem is overcome in Bayesian models that learn the posteriors automatically and estimate the unknowns from prior knowledge and observations. The performance of VBLR method is shown to be better than the simplistic Bayesian approach since it uses varying distribution of priors that estimate the posteriors of heart activities which fit well using KL divergence. The current work on inverse algorithms has been tested on different heart shapes. Further, the method assumed zero effect of the conductivity profile to the magnetic measurements in the forward simulation steps, due to which, the transfer matrix was built to map between different heart shapes and the detectors. The complicated electrodynamic design of heart and torso with respiratory movements, skin conductivity profiles, bio-electric and bio-magnetic fields has to be modelled in future studies to mimic a realistic environment for cardiac source localization.
To avoid the inverse crime, two different geometrical models were used. The forward measurements were simulated from the potentials of the first heart model. The observed MCG data was tested on the different heart geometry (model 2) by inverting the forward matrix constructed from the first model in order to determine the source activities.
The current outcome of this research limits its study to the simulations and solving the inverse problems of MCG on different realistic heart models. The future study of this research involves localization of abnormal cardiac sources due to the effect of other artefacts such as realistic noise levels and breathing movements with physically acquired MCG measurements. The results of this paper promise the clinicians an advantage of diagnosing  www.nature.com/scientificreports/ the infarction region or any other cardiac-related diseases as the reconstructed heart surface potential maps on the myocardium. The study provides the advantage of source localization with high accuracy. This serves as an inspiration to design an online system for non-invasive monitoring of inverse solutions during MCG recording where the clinician will be able to view the unknown functional activities on the screen. However, the models depicted in the paper are limited only to generic models. There is a motive to obtain good quality localization for modeling the geometries of sources acquired from the individual Magnetic Resonance Imaging (MRIs) systems. Further, the inverse problems of the cardiac activations presented in this paper from the MCG signals are to be thoroughly investigated for the subject-specific assessments and could be validated with the recorded epicardial potentials of the structurally healthy hearts.

Data availability
The datasets generated during the current study are available in the ECGSIM software https:// www. ecgsim. org/ index. php and the computations were simulated and analysed using SCIRUN https:// www. sci. utah. edu/ cibcsoftw are/ scirun. html.