J. Magn. Soc. Jpn., 41, 75-80 (2017)<Paper>

Simulation of Extended Source Localization using sLORETA Method forMagnetocardiography
W. Sun and K. Kobayashi
Faculty of Engineering, Iwate University, 4-3-5 Ueda, Morioka, Iwate 020-8551, Japan

In this study, cardiac source localization was simulated using the spatial filter method. Three types of spatial
filters were obtained using the standardized low-resolution brain electromagnetic tomography (sLORETA) method,
based on different examination procedures. In Type A filter, the examination was conducted at the front of the torso.
In both Type B and Type C filters, the examinations were conducted at the front and back of the torso; however, thedistance from the frontal observation plane to the center of the heart model was different for each type. In thesimulation experiments, first the goodness of fit (GOF) value was introduced to determine the proper threshold foreach spatial filter. Then, single and multiple dipole sources were simulated at different depths with and without noise.
The extension of the solutions computed using these spatial filters was investigated. Finally, the performances ofthese spatial filters, with the corresponding averaged thresholds, were evaluated using the GOF. Type B and Type C
spatial filters demonstrated reduction in the extension of source dependency on source depth and improvement in theaccuracy of source localization with noisy data.
Key words: Magnetocardiography, sLORETA method, source localization, spatial filter
1. Introduction
Magnetocardiography (MCG) is a noninvasive
technique that detects the magnetic field generated by
the electrical activity in the heart1). Recently, MCG has
been attracting a lot of attention in relation to the early
detection of heart diseases. As MCG is a multichannel
measurement technique and as the signals are notaffected by the shape of the lungs and torso2), it has highpotential for clinical applications. MCG aims at
obtaining functional mapping, starting from thevisualization of signal sources. Because the MCG
inverse problem is ill posed and does not have a uniquesolution, the estimation of MCG signal sources is thebiggest challenge for its application3). MCG signals are
generated by currents flowing within the myocardial
fibers during cardiac activity4). Therefore, we canconsider that an extended source model will be suitablefor estimating MCG sources.
The spatial filter methods5, 6) were developed for
MCG source localization, and could be used to obtain areliable three-dimensional (3D) outline of the heart.
Among the many spatial filter methods, standardized
low-resolution brain electromagnetic tomography(sLORETA)7) is widely used for visualizing brain activity,
as it has zero localization error for a single dipole withnoiseless data. However, it has some limitations,especially under noisy conditions and in case of deep
sources8, 9). For conventional method, MCG examination
is just conducted in the front of the torso, in order toobtain strong signals. In a previous study10), the spatial
filter always favors the sources that are close to themagnetic field sensors and leads to reconstruction of alarge extension for deep sources. Namely, the extensionof the source depends on the source depth. MCG candetect heart signals not only at the front but also at theback of the torso, even though the signals detected at
the back of the torso are very weak11). The assessmentsof the spatial filter based on the different examination
method have not been studied. In order to assess thesource localization abilities of the spatial filters, andinvestigate the dependency of source extension on thesource depth, the simulation study is needed.
In this study, we make sure the distances from thefrontal and back observation planes to the heart are thesame, to reduce the dependency of the source extension
on the source depth, and obtain a spatial filter that canestimate a cardiac source well with noisy data.
2. Simulation ModelIn this study, we simulated MCG source localization
using the sLORETA method based on a conventional
64-channel MCG system. Superconducting quantum
interference device (SQUID) modules were aligned in an
8 × 8 configuration to create a square observation area
sized 175 mm × 175 mm12).

Fig. 1 Frontal and back observation planes and theposition of heart in torso.

Journal of the Magnetics Society of Japan Vol.41, No.4, 201775

Simulations were carried out for the 64-channel
MCG system within a 168 mm × 168 mm × 120 mmvolume. Six thousand six hundred and fifteen voxels
were divided with the size of each voxel being 8 mm × 8
mm × 8 mm. Fig. 1 shows three types examination
procedures. Type A is the conventional examination
method, which is only conducted at the front of the torso.
In both Type B and Type C, the examinations are
conducted not only at the front but also at the back of
the torso; however, the distance from the frontal
observation plane to the center of the heart model isdifferent for each type. For Type C, in order tocompensate for the varying sensitivities of the sensors
to the current sources at different depths, we make surethe distances from the frontal and back observation
planes to the heart are the same. The distance from thefrontal observation plane to the center of the heartmodel was set as 80 mm for Type A and Type B and as
160 mm for Type C. The distance from the backobservation plane to the center of the heart model was
set as 160 mm for Type B and Type C. For the forward
problem, we calculated the lead-field matrix K using the
Biot–Savart Law.
3. Methods
sLORETA is a spatial filter method used to solve theill posed problem. It can obtain a blurred source image
by using a linear inverse operator7).

(1)Φ  ΚJ  nwhere Ф is the detected signals, K is the lead-field
matrix, J is the original current sources within theheart, and n is the noise at each sensor. Using thewell-known minimum norm estimate (MNE) method,the estimated current sources Ĵ can be written as:Ĵ  GΦ  K T KK T  2 I Φ(2)MNE solution using its variance in the following way:Ĵ sLORETA ,i  Sii

1/2

Ĵ i

(3)

where Sii is the diagonal element of the covariancematrix S withS  GK(4)

In this study, we first calculate the lead-field matrix K
that describes the sensitivity of the sensors to eachcurrent source location. For each current source location,
K includes two columns that represent the X component
and Y component, respectively. Then, Type A, Type B,
and Type C spatial filters are obtained using thesLORETA method.
In the coordinate system, the center point of the top
voxels layer is (0, 0, 20) and the unit is mm. A single
dipole source (0, 0, 20) in the X direction is simulated bychanging the depth in 8 mm increments.
Since the real cardiac sources are considered to beextended sources, simulations with multiple dipoles arealso performed. A normal MCG wave is obtained using
the Iwate create project 64-channel MCG system11). An
isomagnetic field map of the MCG wave at the R peak isshown in Fig. 2(a). The generated magnetic field of themultiple sources is shown in Fig. 2(b). Multiple
sources—(−16, 8, 20) in the X direction, (−8, 8, 20) in the
X direction, and (0, 0, 20) in the Y direction—are
simulated by changing the depth in 8 mm increments.
4. Simulation Results and Discussions
4.1 Simulation without Noise
Fig. 3 shows the simulation results of a single dipole
at different depths (depth = 28 mm, 76 mm, and 124mm). The relative strength was normalized by the
maximum of the solutions. The white rectangle denotes
the maximum of the solutions. All the maximum

where G is a spatial filter matrix, T denotes the matrix
transpose, λ2 is the Tikhonov regularization parameter,and I represents an identity matrix. The sLORETA
solution can be obtained from the standardization of theFig. 2 (a) Isomagnetic field map of MCG wave at R
peak. (b) Isomagnetic field map of multiple dipole
sources. GOF of two patterns is 0.93. One step of theisomagnetic field contour line is 10 pT.

76

Fig. 3 Left column is the true source in the X direction
placed at different depths (from top to bottom depth = 28mm, 76 mm, and 124 mm). The three right columns arecross sections of simulation results using Type A, TypeB, and Type C spatial filters. The white rectangle
denotes the maximum of the solutions.

Journal of the Magnetics Society of Japan Vol.41, No.4, 2017

solutions are located at real source positions. It meansthat there is no localization error without noise.
However, the solutions have large spatial extensions.
Obviously, the extensions of the solutions computed by
the Type A spatial filter have larger dependencies on thedepths of the true source than that computed by theType B and Type C spatial filters. Because a source thatis far from the sensor contributes less to the measured
data than one that is near the sensor, the solution needs
large extensions to fit the measured data. By contrast,
Type B and Type C spatial filters produce solutions withextensions that have fewer dependencies on the depthsof true sources.
For an extended source, a threshold of relative
strength has to be determined. The full width at half
maximum (FWHM) is used as the threshold13, 14). It
means that only the excess 50% of the maximum
sources are analyzed to avoid influence from the weakinsignificant source estimates. In our study, thegoodness of fit (GOF) value is introduced to determine
the proper threshold. GOF describes the similarities in
the measured and estimated signals at all N sensors15).
GOF is defined as follows:
N

N

i 1

i 1

2
GOF  1   Bmea,i  Bcal,i 2 /  Bcal,
iFig. 4 Left column is the true source in the X
direction placed at different depths (from top tobottom depth = 28 mm, 76 mm, and 124 mm). The
three right columns are cross sections of simulation
results using Type A, Type B, and Type C spatial
filters (threshold = 0.91, 0.90, and 0.93, respectively).
The white rectangle denotes the maximum of thesolutions.

(5)

where Bmea,i denotes the measured magnetic field andTable 1 Threshold of relative strength simulated withthe true source at different depths when thecorresponding GOF is close to 1. (All of the GOF is 0.99
in this study).
depth (mm)Type A

Type B

Type C

20

0.847

0.892

0.940

28

0.830

0.856

0.946

36

0.850

0.860

0.950

44

0.866

0.870

0.940

52

0.902

0.880

0.925

60

0.910

0.883

0.917

68

0.924

0.885

0.911

76

0.934

0.886

0.900

84

0.936

0.910

0.910

92

0.942

0.925

0.920

100

0.942

0.929

0.931

108

0.938

0.938

0.940

116

0.937

0.950

0.952

124

0.927

0.950

0.951

132

0.923

0.947

0.947

average

0.91

0.90

0.93

Fig. 5
Relationship of volume percentage and depthof true source.

Fig. 6 GOFs of solutions computed using Type A,
Type B, and Type C spatial filters when the true
source is placed at different depths.

Journal of the Magnetics Society of Japan Vol.41, No.4, 201777

Bcal,i denotes the calculated magnetic field at the i thsensor. A high GOF indicates that the estimated sources
can explain the measured data well. The thresholdvalues with true sources at different depths are shownin Table 1, when the corresponding GOF is close to 1.
The average value is chosen for computing the solutions.
The relative strength, which does not exceed theaverage value of thresholds, is considered to be 0.
In the simulation, a single dipole source in the X
direction is placed at different depths (depth = 28 mm,
76 mm, and 124 mm). Fig. 4 shows the cross sections of
the solutions computed by Type A, Type B, and Type C
spatial filters. The extension of the solution has a large
dependency on the depth of the true source for the TypeA spatial filter because of the varying solid angles
between the dipole at different depths and theobservation plane. In addition, the solutions computedby the Type B and Type C spatial filters show smaller
extensions in the Z direction than that computed by theType A spatial filter, owning to the data measured from

Fig. 7 Left column is the true sources at different
depths (from top to bottom depth = 28 mm, 76 mm,
and 124 mm). The three right columns are crosssections of simulation results using Type A, Type B,
and Type C spatial filters (threshold = 0.91, 0.90, and0.93, respectively). The white rectangle denotes the
maximum of the solutions.

Fig. 8 Relationship of volume percentage and depthsof true sources.

78

Fig. 9 GOFs of solutions computed using Type A,
Type B, and Type C spatial filters when the true
sources are placed at different depths.
the back of torso. Because the distances from the frontal
and back observation planes to the heart model are thesame, the extensions show less dependency on the depthof the true source for the Type C spatial filter.
By computing the volumes of the voxels in which thesolutions are localized, the relationship between thevolume percentage and the true source depth can beobtained as shown in Fig. 5. The extensions of solutions
show less dependency on the depth of the true source for
the Type B and Type C spatial filters. The GOFs of TypeA, Type B, and Type C spatial filters with the proper
thresholds are shown in Fig. 6. The GOF values of theType A and Type B spatial filters are far away from 1
and even drop to negative values for deep sources.
Therefore, the Type C spatial filter shows better
performance.
Fig. 7 shows the cross sections of the solutions
computed using the Type A, Type B, and Type C spatial
filters. By computing the volumes of the voxels in which
the solutions are localized, the relationship between thevolume percentage and the true source depth can beobtained as shown in Fig. 8. Obviously, the Type B andType C spatial filters show extensions that have lessdependency on the depth of the true sources, comparedto the Type A spatial filter. Fig. 9 shows the GOFs of theType A, Type B, and Type C spatial filters simulatedusing sources at different depths. The GOFs of Type A
and Type B spatial filters are far away from 1 and evendrop to negative values for deep sources. It means thatthese two spatial filters cannot fit the measured datawell for deep sources. By contrast, the GOF of the TypeC spatial filter has a good performance despite theslight drop in the case of deep sources.
4.2 Simulation with Noisy Data
Since the MCG inverse problem is ill posed, weaknoise signals can cause large changes in the solution. In
this study, Gaussian white noise is added to thesimulated magnetic field data and the inverse solutions

Journal of the Magnetics Society of Japan Vol.41, No.4, 2017

are computed using the above three spatial filters. The
SNR ranges from 0 to 40 in 5 dB increments. The SNR is
defined as:SNR  10log varBexact  / σ 2(6)where Bexact is the variance of the simulated noise-free
measurements and σ2 is the variance of the added
Gaussian white noise. In this study, for a single dipole
source, the localization error is the distance from the
maximum estimated source to the true source. For
multiple sources, the localization error is the distancefrom the maximum estimated source to the center pointof the true sources. The mean value is obtained from100 calculations. In order to obtain the samemeasurement conditions, the same σ value of noise is set
up for these three spatial filters at a certain noise level.
Since the magnitudes of the signals detected at thefrontal and back observation planes are different, the
SNR is different. In this study, the x-axis values of Fig.
10 and Fig. 11 are calculated only using themeasurement data that obtained based on Type A
examination procedure.

Fig. 10 Mean localization error of single dipole
source (depth = 124 mm) with different noise levels.

Fig. 10 and Fig. 11 show the simulation results of asingle dipole source and multiple sources, respectively.
The single dipole source and the multiple sources areplaced at relatively deep positions with a depth of 124
mm. For single dipole source estimation, very weaknoise led to Type A spatial filter generating a large
localization error. By contrast, the localization errors ofType B and Type C spatial filters were stable and about
0 at 20 ~ 40 dB. For multiple sources estimation, thelocalization errors of Type B and Type C spatial filters
were stable at 0 ~ 40 dB. Apparently, the Type B andType C spatial filters have improved source localization
accuracies. It may be owing to the data obtained fromthe back of the torso. In addition, Type B and Type C
spatial filters have good performances for multiple
source localizations even at high noise levels. Hence,
they are fit for estimating multiple sources.
5. Conclusion
In this study, Type A, Type B, and Type C spatial
filters were obtained using the sLORETA method basedon three types of MCG examination procedures. The
GOF was introduced to determine the proper thresholdfor each spatial filter. The extensions of the solutions
computed using these spatial filters with thecorresponding thresholds were investigated. In thesimulation of a single dipole and multiple dipoles placedat different depths, the Type C spatial filter
demonstrated an extension that had less dependency on
the depth of the true source. In addition, the GOF of theType C spatial filter was close to 1 even for deep sources.
It meant that the Type C spatial filter had the best
performance for extended source estimation, compared
to Type A and Type B spatial filters. From the results ofthe simulations with noisy data, we found that the TypeB and Type C spatial filters had improved source
localization accuracies and that they were fit for
estimating extended sources. The reduction of thesource extension for deep sources and the improvementof the source localization accuracies owned to theexamination at the back of the torso. The same distancefrom the frontal and back observation planes led to theless dependency of the source extension on the source
depth. In the future, we hope to detect the real MCG
signals and estimate the cardiac source using Type B
and Type C spatial filters.
Acknowledgements This work was supported by theJapan Society for the Promotion of Science KAKENHI(Grants-in-Aid for Scientific Research (C)) Grant
Number JP26350535.
References

Fig. 11 Mean localization error of multiple sources(depth = 124 mm) with different noise levels.

1) J. S. W. Kwong, B. Leithäuser, J.-W. Park, and C.-M.Yu: Int.
J. Cardiol., 167, 1835 ( 2013).
2) K. Nakai et al.: Int. J. Cardiovasc. Imaging, 21, 555 (2005).
3) V. Mäntynen, T. Konttila, and M. Stenroos: Phys. Med. Biol.,59, 7141 (2014).

Journal of the Magnetics Society of Japan Vol.41, No.4, 2017793) M. De Melis et al.: J. Magn. Soc. Jpn., 34, 453 (2010).
4) I. Tavarozzi et al.: Ital. Heart J., 3, 75 (2002).
5) K. Nakai et al.: Int. J. Cardiovasc. Imaging, 22, 573 (2006).
6) G. Shetty and M K: Int. J. Emerg. Technol. Innov. Res., 2,1156 (2015).
7) R. D. Pascual-Marqui: Methods Find. Exp. Clin. Pharmacol.,24,5 (2002).
8) E. Palmero-Soler, K. Dolan, V. Hadamschek and P. A. Tass:
Phys. Med. Biol., 52, 1783 (2007).
9) R. Terusawa et al.: J. Magn. Soc. Jpn., 36, 272 (2012).

10) M. Yoshizawa, K. Nakai, K. Kobayashi, Y. Nakamura, Y.
Uchikawa: J. Jpn. Biomagn. Bioelectromagn. Soc., 15, 109(2002).
11) K. Yoshida et al.: J. Arrhythm., 31, 345 (2015).
12) M. Yoshizawa et al.: Phys. C: Supercond. Appl., 426, 1572
(2005).
13) F.-H. Lin et al.: NeuroImage, 31, 160 (2006).
14) M. Fuchs, M. Wagner, T. Kohler, and H.-A. Wischmann: J.
Clin. Neurophysiol. 16, 267 (1999).
15) K. Pesola et al.: Med. Biol. Eng. Comput., 37, 2 (1999).

Received Jan. 19, 2017; Accepted May 15, 2017.

80

Journal of the Magnetics Society of Japan Vol.41, No.4, 2017