Investigations of sensitivity and resolution of ECG and MCG in a realistically shaped thorax model

Solving the inverse problem of electrocardiography (ECG) and magnetocardiography (MCG) is often referred to as cardiac source imaging. Spatial properties of ECG and MCG as imaging systems are, however, not well known. In this modelling study, we investigate the sensitivity and point-spread function (PSF) of ECG, MCG, and combined ECG+MCG as a function of source position and orientation, globally around the ventricles: signal topographies are modelled using a realistically-shaped volume conductor model, and the inverse problem is solved using a distributed source model and linear source estimation with minimal use of prior information. The results show that the sensitivity depends not only on the modality but also on the location and orientation of the source and that the sensitivity distribution is clearly reflected in the PSF. MCG can better characterize tangential anterior sources (with respect to the heart surface), while ECG excels with normally-oriented and posterior sources. Compared to either modality used alone, the sensitivity of combined ECG+MCG is less dependent on source orientation per source location, leading to better source estimates. Thus, for maximal sensitivity and optimal source estimation, the electric and magnetic measurements should be combined.


Introduction
Electrochemical activity of cardiac muscle cells gives rise to an electromagnetic field. To study the electrical activity of the heart, this field is measured and analyzed: Electrocardiography (ECG) measures differences of electric potential on the body surface, and magnetocardiography (MCG) measures the magnetic field outside the body. In clinical use, the most common electrocardiographic application is the 12-lead ECG, where the electrical activity of the heart is interpreted directly from signal tracings. The inverse problem of ECG/MCG aims at characterizing the cardiac activity by combining multichannel measurements with a computational model of the origin and transfer of signal (Nenonen 1994, Gulrajani 1998, MacLeod and Brooks 1998: First, the forward problem is solved, i.e., signal topographies are modelled using a template source (distribution) and a volume conductor model of the thorax. Then, sources underlying the measured signals are estimated using this model and prior information on the source. Forward and inverse modelling thus yields information that may help the interpretation of measured ECG and MCG signals.
In clinical use, invasive study of the cardiac electrical activity is considered the gold standard. Treatment of several disorders relies on this kind of functional information. For example, catheter ablation procedures have shown major progress over the past decade; emerging techniques enable electrophysiological study and ablation even without x-ray fluoroscopy on which these procedures traditionally depend (Anselmino et al 2013). Most often the invasive diagnosis and treatment are performed during the same operation (see e.g. Berger et al 2006) and the decision between operative and other treatment options has to be made based on prior knowledge and other, usually noninvasive, diagnostic methods. Noninvasive cardiac source imaging with ECG or MCG could provide information that facilitates allocation, planning and execution of this kind of invasive diagnosis and treatment, leading to reduced operational and fluoroscopy time (Berger et al 2006, Kim et al 2007. Noninvasive cardiac source imaging can perhaps also be used to monitor the treatment outcome and it can also be applied in general population. In addition to its potential direct clinical benefit, it thus has the potential to increase our understanding of the physiology and pathophysiology of the heart. The inverse problem of electrocardiography is ill-posed and does not have a unique solution. To render the solution unique and stable, the source estimate needs to be constrained using prior information and possibly additional regularization. The prior information typically affects both the source space and source estimates. The inverse problem of ECG has been reviewed, e.g., in MacLeod and Brooks (1998), Dössel (2000) and Pullan et al (2010); here we briefly introduce the main concepts. The most restricted inverse solution is obtained with moving-single-dipole modelling (Gulrajani et al 1984), in which the whole ECG/MCG signal at one time-instant is assumed to come from one focal source. The single-dipole model can be successfully used for e.g., localizing the earliest ventricular activation or a small ischemic lesion, but in connection to such a restrictive model, it is essential to have independent information that warrants the use of the model (Gulrajani et al 1984). If there is no prior information that would support a model with only one or a couple of focal sources, a distributed source model that allows sources everywhere in the heart (Han et al 2012) or on a cardiac surface , Ghanem et al 2005, Berger et al 2006, van Dam et al 2009 is typically applied. The distributed surface source can be implemented as a primary current distribution (Gulrajani 1998, transmembrane voltage distribution (Berger et al 2006) or, in the case of ECG, as an epicardial potential distribution (MacLeod et al 1995, Ghanem et al 2005, Stenroos and Haueisen 2008. To make a distributed-source inverse problem well-posed, further restrictions on the source estimate are needed. For example, if focal sources are assumed but also extended sources may exist, an approach that minimized the L1 norm of the solution (e.g. Ghosh and Rudy 2009) may be optimal. Loosening the constraint, we may choose to promote overall smoothness (the 2nd order Tikhonov regularization that minimizes the L2 norm of the Laplacian of the source; Horácek and Clements 1997, Ghanem et al 2005, Stenroos and Haueisen 2008 or the overall L2 norm of the solution (minimum-norm estimate, also called the 0th order Tikhonov regularization). Also in the case of distributed solutions, the constraints need to be posed a priori based on independent information; if no specific prior information is available, the minimum-norm estimator provides the least-constrained and thus the most general solution. Regarding model optimization, a major benefit for L2-minimizing methods is that they can be expressed in closed form as linear estimators, independent of the actual measurement data, while L1 methods depend on data in a non-linear way and are thus considerably more difficult to assess.
Most earlier ECG or MCG inverse studies have concentrated on solving a specific problem, thus enabling the use of specific constraints based on assumptions. Typical examples of these assumptions include the focality of ectopic activity (Lai et al 2013), propagation of depolarization wavefront Pesola 2001, van Dam et al 2009) and static phase of activity , Konttila et al 2014. The use of strong constraints and prior assumptions may, however, complicate generalization of the methods and results.
The properties of source estimates depend on the amount of information brought into the process by three main parts: measurement, model of signal transfer, and prior knowledge about the sources. Several studies have suggested that MCG and ECG may provide complementary information (Lant et al 1990, Brockmeier et al 1997, Tenner et al 1999, Oostendorp and Pesola 2001. It has been interpreted that the sensitivity of MCG and ECG is different for particular source configurations: MCG has been considered more sensitive to 'tangential' , Tenner et al 1999 and vortex (Brockmeier et al 1997, Liehr et al 2005 currents whereas ECG has been considered more sensitive to 'radial' , Tenner et al 1999 currents. However, it is not always clear whether the term 'current' refers to cardiac activity (primary currents) or to the total current close to the body surface and sensors. In addition, 'radial/tangential' sometimes refers to orientation of this current with respect to the heart, to the ECG/MCG sensors, or to the body surface. To our knowledge, the difference of ECG/MCG sensitivity and its possible effect on solving the inverse problem has not been studied quantitatively and systematically as a function of source location and orientation in realistic geometry.
Solving the inverse problem of ECG and MCG is often referred to as cardiac source imaging. However, the spatial properties of ECG and MCG as imaging systems, especially when combined, are not well known. The spatial properties of an ECG/MCG measurement can be studied in terms of sensitivity of sensor arrays to elementary sources at different positions. For linear estimators, the point-spread function (PSF) describes the response of an imaging system to a point-like elementary source; the response to any source can be written as a linear combination of PSFs, and the location and spread of each PSF can be characterized with simple metrics. This kind of analysis has recently been applied to electro-and magnetoencephalography (see, e.g., Molins et al 2008, Hauk et al 2011, Stenroos and Hauk 2013, but, to our knowledge, it has not been applied to ECG or MCG.
In this study, we investigate spatial properties of ECG, MCG, and combined MCG+ECG measurements using sensitivity and PSF analysis. Signal topographies are modelled using a realistically-shaped volume conductor model, and the inverse problem is solved using a distributed source model and linear source estimation with minimal use of prior information. To our knowledge, this is the first time that the ECG and MCG are assessed together as imaging systems and their spatial properties are assessed globally around the ventricles. This analysis brings out spatial features and biases of modalities that need to be understood when interpreting the results of source estimation and that facilitate interpretation of any multichannel ECG/MCG measurement.

Models and methods
The MCG and ECG signal topographies were simulated in sensor geometries of measurement devices in BioMag Laboratory, Helsinki, Finland, illustrated in figure 1. For MCG, the 99-channel Neuromag MCG device was used. The measurement head of the device was positioned over the left side of the anterior chest and tilted slightly to follow the curvature of the thorax, like we do in our experimental MCG studies. The measurement head holds a slightly curved array of thin-film sensor triplets in 33 locations over the area of about 30 cm in diameter; in each location there is a magnetometer detecting B z and two orthogonal planar gradi- ; here x, y, and z refer to local coordinates of the sensor, z being the normal vector of the sensor unit. The spatial extent of the MCG sensors was taken into account by using numerical quadrature with magnetic field computed at four points per sensor. For ECG, we used the 120-lead Helsinki body-surface potential mapping (BSPM) electrode layout covering the whole thorax; the layout is based on using 18 vertical strips of electrodes placed according to anatomical landmarks, with 5 cm interelectrode spacing; for details and figure, see e.g. Takala et al (2001).
The anatomical model was based on the Dalhousie standard thorax model (Hren et al 1998), illustrated in figure 1. To ensure adequate numerical accuracy we refined 3 the triangulated surface meshes of the model compartments. For volume conductor model, the compartments for thorax (with 1402 final vertices), lungs (656), and intra-cardiac blood (688) were assumed homogeneous with electrical conductivities of [1 1/4 4] × 0.23 S/m, respectively. For source space, the epicardial surface of the ventricles (802) was slightly modified: The part of the mesh that intersected the blood meshes and the vertices that were too close to the blood for stable numerical computation were excluded (leaving 740 vertices for source space).

Forward computation
Electrochemical activity of myocardial cells gives rise to an electromagnetic field. The source of the fields can be expressed in macroscopic scale as primary current density J p (Tripp 1983, Sarvas 1987, Gulrajani 1998; the tissues of the thorax are assumed resistive and thus the problem obeys the quasistatic Maxwell equations (Sarvas 1987, Gulrajani 1998). The primary current density gives rise to a primary charge distribution and electric field Figure 1. The model geometry. On the left, MCG sensors and anterior BSPM electrodes are shown with blue squares and black dots, respectively; in the middle, the anterior view of the anatomical model is shown; on the right, the discretized source locations are shown with black dots on the surface mesh of the ventricles along with example source triplets with red arrow corresponding to the normal (n), blue to the first tangential (t 1 ), and green to the second tangential (t 2 ) basis vector orientation.
where ϕ is the electric potential and σ conductivity. The total current is thus = + J J J p V . Due to the conservation of charge we then get the Poisson equation for potential the differences of potential on the body surface are measured as ECG. After solving ϕ, the magnetic field B can be integrated from the Biot-Savart law where r is the observation point and the integral is calculated with respect to primed coordinate ′ r over the whole thorax volume ′ V . As all currents give rise to magnetic field, the field measured as MCG consists of two components: the field generated directly by the primary current density J p and the field due to volume current density J V (see, e.g. van Oosterom et al 1990, Nenonen 1994. In this work, the source was modelled as a freely-oriented primary current density inside the myocardium , Han et al 2012, Goernig et al 2013, Lai et al 2013. For numerical computations, the source was discretized into a set of unit-amplitude dipoles that can be interpreted as source atoms (the ECG/MCG topography of any source distribution can be described with a linear combination of topographies of such elementary sources): Three orthogonal dipoles were placed in every vertice of the source mesh (see figure 1), effectively lumping together transmural segments of primary current. Similar surface-constrained approaches have been used earlier to describe cardiac source activity (for example, Nenonen et al 2001, Oostendorp and Pesola 2001, van Dam et al 2009, Goernig et al 2013. The unit-dipole orientations were set according to local geometry of the heart: For each vertex of the source mesh, an orthonormal local basis consisting of two tangential unit vectors (tangential sources, denoted as t 1 and t 2 ) and a unit normal vector (normal source, n) was defined. Orientation of the first tangential basis vector was set roughly towards the apex of the heart, and the second tangential basis vector was then acquired through cross product.
The forward problem was solved using the boundary element method (BEM), formulated using the linear Galerkin approach (Stenroos et al 2007, Stenroos andHaueisen 2008). Signal topographies generated by dipole triplets in each vertex of the source mesh were computed and collected into (columns of) a forward operator L (often called the lead-field matrix). The zero level of the potential was set to the mean value of all electrodes. The measurement/forward model for signal d is then where vector s contains the amplitudes of the discretized primary current density and n is the measurement noise.

Source estimation
We estimate the sources using linear inverse operator (i.e. estimator) G: This problem is ill-posed and the solution needs to be constrained. The well-known minimum-L2-norm (MN) estimator , Lin et al 2006 provides such a solution with minimal constraints: The MN estimator minimizes the L2-norm of the estimate while balancing between the reconstruction of noise and measured signal. Assuming all source candidates are treated equally, the MN estimate can be written (van Oosterom 1999, Lin et al 2006 where T denotes matrix transpose, λ 2 is the regularization parameter, and C is the estimate of noise covariance matrix. This MN estimate may also be interpreted as a Bayesian maximum-a-posteriori estimate (van Oosterom 1999). Using the noise covariance matrix C as the regularization operator is equivalent to transforming the signal and the forward model into signal-to-noise basis through whitening (Lin et al 2006); the whitening also allows combining measurements with different sensor types, such as magnetometers and gradiometers or MCG and ECG, in a single estimate. For whitened data = ∼ − d C d 1 2 and lead-fields = ∼ − L C L 1 2 , the above estimate can be written as, which is the classical MN estimate. In this work, we assumed the same average signal-tonoise ratio (SNR) for all sensor types. This was implemented by scaling the topographies and signals (Fuchs et al 1998, Oostendorp and based on the average signal power in each sensor type, i.e., the mean signal power over the whole source space was computed for each sensor type and the same value was set in all corresponding diagonal elements of the (diagonal) scaling matrix replacing C in the above equations. The regularization parameter λ 2 was set as suggested by Lin et al (2006), assuming the mean power-SNR-ratio of 10.

Sensitivity and resolution metrics
Spatial properties of MCG/BSPM were assessed using two types of metrics that characterize the signal generation and source estimation: The gross signal strength over all sensors was quantified for each elementary source separately; this sensitivity was defined as the L2-norm of the whitened (scaled) topography: where l i is the i th column of the lead-field matrix. Both forward and inverse model were used to assess the source estimation performance with the help of resolution analysis that has previously been used in geophysics (Backus andGilbert 1968, Tarantola 2005) and electro-and magnetoencephalography (Lin et al 2006, Molins et al 2008, Hauk et al 2011, Stenroos and Hauk 2013. The resolution operator is defined (Tarantola 2005) as R = GL. Columns of R are called point-spread functions (PSF) whereas rows are called cross-talk functions (CTF): Each column/PSF shows how the estimate for the corresponding elementary source spreads to the source space. Each row/CTF shows how elementary sources over the whole source space would contribute, or leak, to the estimate at target vertex. Thus the performance of any linear estimator and measurement modality can be characterized with the PSFs and CTFs that are independent of the actual measured signals; the response to any source distribution represented by the source space can be expressed as linear combination of PSFs. As the PSF and CTF for each elementary source are identical for the MN estimator used here (Liu et al 2002), we present the results for the PSFs only; the results apply for CTFs as well.
PSF is a general description of an imaging system. Because it would be impractical if not impossible to examine every PSF visually (in our case 2220 distributions on the heart surface), specific features of PSF (distribution) are quantified with simple metrics. In brain studies, PSF/CTF analysis has been applied to scalar-valued (cortically-constrained) source model (Molins et al 2008, Hauk et al 2011, Stenroos and Hauk 2013. In this work, we extend the analysis to vector-valued sources and compute three resolution metrics: localization error, spatial deviation, and orientation error. The localization performance is assessed by computing the distance between the center-of-mass of the estimate and the true source. For a source at r true , the localization error LE is (Lin et al 2006, Stenroos andHauk 2013) where ŝ i is the estimated (discretized three-component) source at i th vertex of the source mesh. The spread of the estimate is assessed by computing the spatial deviation (SD) (Stenroos and Hauk 2013) The performance of the estimator in reconstructing the source orientation is assessed using a novel metric, the orientation error (OE): For a (unit) source s true , true whose value of 0 corresponds to parallel and 1 to orthogonal orientation of the PSF/estimate components in comparison to the true source orientation. All metrics were focused on the area where PSF amplitude exceeds half-maximum value ('hot spot'): in the above the index i runs through all locations in which ŝ i has values above relative threshold t, ≥ŝ s t , i max with = t 0.5

Computations and results
Forward and inverse operators were built for MCG and BSPM. The sensitivity and PSF as well as the associated metrics were then computed separately for each source location and orientation; the metrics as function of source position were plotted as distributions on the epicardial surface. In addition, the computations were repeated for a combined MCG+ECG estimator. All computations and visualizations were carried out using Matlab and in-house BEM tools (Stenroos and Haueisen 2008) further developed from Helsinki BEM library (Stenroos et al 2007).

Example: signal topographies and point-spread functions (PSFs)
In figure 2 we show an example of MCG/BSPM signal topographies and point-spread functions for one source location on the anterior wall. This example shows differences in sensitivity and source-estimation performance between modalities that are typical for a large part of anterior myocardium. The signal topographies for three source orientations are presented on the upper panel of the figure. In MCG, topographies show higher amplitude for the tangential source components than for the normal component, whereas BSPM shows relatively higher amplitude for the normal component. For this source location, MCG is thus relatively more sensitive to tangential sources whereas BSPM is relatively more sensitive to the normal source, as indicated also in the corresponding sensitivity values shown below the topographies.
The point-spread functions for the three source orientations are presented on the lower panel of figure 2. In MCG, PSFs of tangential sources are peaking closer to the source location compared to the PSF of the normal source, whereas in BSPM the opposite holds. This relative difference is also indicated in the corresponding values of LE. In MCG, the PSFs of the tangential sources are slightly less spread around the source location than the PSF for the normal source, as indicated by the corresponding SD values. Again, the opposite holds for BSPM.
The orientation distribution of the three-component PSF is visualized along with the true source orientation in the close-ups on the lower panel of figure 2. This relationship is characterized with the OE metric: In MCG, the tangential source orientations are better reconstructed than the normal source orientation, as indicated in the corresponding values of OE. Signal topographies, PSFs, and corresponding metrics for one source location and three orientations. Signal topographies for tangential unit sources (t 1 , t 2 ) and a normal source (n) are shown in the top row; the corresponding PSFs are shown in the middle row, with close-up of the source region in the lower row. In each PSF plot, the source location and orientation are indicated by a black dot and line, respectively. The signal topographies (for 1 Am unit sources) per modality have the same symmetric color scaling about zero. Each PSF plot shows the color-coded amplitude of the estimate, normalized with respect to its maximum; the orientations of the estimated source distribution are indicated with arrows on top of each PSF plot (note that the arrows pointing into the heart, as in BSPM for tangential sources (t 1 , t 2 ), are not visible through the colormap). Below the plots, values for sensitivity, PSF maxima, and, metrics LE (cm), SD (cm), and OE (%) are shown.
In BSPM, the PSF amplitude for the tangential sources peaks on both sides of the source location, roughly on the axis along the source orientation; a closer inspection reveals that the strongest PSF elements are roughly parallel to the source-surface normal and of opposite polarity within each peak. Thus the strongest PSF elements are almost orthogonal to the tangential (normal) source orientation in BSPM (MCG); this is indicated also in the relatively high value of the corresponding OE values. Figure 3 presents the sensitivity of MCG and BSPM as a function of source position and orientation, plotted as distributions on the source surface. The results show that MCG is more sensitive to tangential sources than to normal sources, whereas BSPM is more sensitive to normal sources than to tangential sources. This difference is most pronounced in the anterior myocardium where the tangential sources dominate in MCG and normal sources in BSPM. Sensitivity to posterior sources is clearly lower compared to anterior sources, especially in MCG.

Resolution metrics
Localization errors (LE) are presented in figure 4. In MCG, tangential sources (t 1 , t 2 ) on the anterior wall are better localized compared to normal sources (n), whereas in BSPM the normal sources are relatively better localized. In comparison, LE is clearly higher for posterior sources with similar relative difference between tangential and normal sources in BSPM; in MCG the localization of posterior sources fails, apart from the area near the apex of the heart. A comparison between the modalities shows that MCG tends to localize the anterior tangential sources better whereas BSPM outperforms MCG in localization of anterior normal sources and posterior sources.
Spatial deviations (SD) are presented in figure 5. The results for MCG show that the PSFs for tangential sources are less spread with slightly lower SD compared to the normal sources whereas in BSPM the SD is relatively lower for normal sources. Again, for both modalities the error is larger on the posterior wall. A comparison between modalities shows that in MCG the PSFs for anterior tangential sources are slightly less spread compared to BSPM; in BSPM the SD is lower for normal and for posterior sources. Orientation errors (OE) are presented in figure 6. On the anterior wall, the OE shows similar behavior to the LE and SD above: The results for MCG show that the orientation of the tangential sources is better reconstructed compared to the normal sources, whereas in BSPM the opposite holds. Contrary to LE and SD results above, in BSPM the OE is lower for most posterior sources in comparison to anterior sources; in MCG this holds for normal sources. A comparison between modalities shows that the orientation of the anterior tangential sources is better reconstructed in MCG, whereas BSPM outperforms MCG for anterior normal sources and for posterior sources.

Combined MCG and BSPM
The resolution metrics presented above show clearly that MCG and BSPM perform in a complementary way. Thus, to study the potential improvement in performance, the analysis was repeated for a combined MCG and BSPM. The combining of the modalities was done in the same way as the combination of gradiometers and magnetometers in MCG, using sensor weighting based on mean signal energies of the modalities. The results are presented in figure 7: The sensitivity distributions for different source orientations are clearly more similar  with locally roughly equal sensitivity to tangential and normal sources (compare to figure 3). The estimates of the combined measurement have clearly lower errors LE, SD, and OE in areas that are problematic with either modality alone (compare to figures 4-6).

Discussion
In this study, we assessed spatial properties of ECG and MCG in terms of sensitivity of sensor arrays to elementary sources distributed around ventricles and in terms of source estimation performance using linear inversion and point-spread function (PSF). The results show that the modalities are different and complementary in nature: For sources on the anterior wall of the ventricles, the sensitivity of MCG is higher for tangential than for normal sources whereas BSPM is relatively more sensitive to normal sources; the sensitivity of the combined MCG+BSPM is clearly less dependent on source orientation. In comparison, the sensitivity to posterior sources is clearly lower, especially in MCG. The PSF analysis shows that MCG performs better in estimation of anterior tangential sources whereas BSPM is superior in estimating anterior normal sources and posterior sources. With combined MCG+ECG the localization error, spread, and orientation error of source estimates are clearly smaller, especially in areas that are problematic with either modality used alone; with the added information in the measurement and model, we can build an estimator that better distinguishes between signal topographies of different sources.
The results show that the sensitivity profile of the modality is clearly reflected in source estimates so that the estimates are biased towards strong elementary sources (those with high sensitivity). This may lead to source estimates with erroneous location and orientation; our results reveal this risk for MCG with normal sources and BSPM with tangential sources on the anterior wall of the ventricles which is traditionally considered as the region easiest to characterize in electrocardiography (for example, see figure 2). In addition, the spread of the estimate is higher in areas with low sensitivity such as the posterior wall. These uncertainties should be taken into account when interpreting source estimates. Compared to either modality used alone, the sensitivity of the combined MCG+ECG is clearly less biased to different source orientations per source position. This complementarity leads to superior source estimates and, in addition, is likely exploitable in most signal analysis approaches also in sensor domain. Thus, for maximal sensitivity, optimal source estimation and largest information content, the electric and magnetic measurements should be combined.

Models
Primary current distribution can be regarded as the most general macroscopic source model for describing the bioelectromagnetic field on and outside the body surface, i.e., in the farfield. The freely-oriented primary current distribution is also capable of expressing the socalled electrically silent sources like a current loop (magnetic dipole). Volume current and ECG arise from sources and sinks, i.e., the divergence of the primary current as in (1), whereas the magnetic field depends in addition on the path of the primary current as in (2). The elementary source used for forward simulations in this study is a current dipole, i.e., a short, straight segment of primary current. Special cases of curved sources (Brockmeier et al 1997, Figure 7. Sensitivity and resolution metrics for combined MCG and BSPM. For further explanation, see the captions of figures 3-6. Liehr et al 2005), associated with electrically silent magnetic source components, are out of scope of this work and should be investigated later; the far field of such sources, like any primary currents, can be modelled with a couple of current dipoles. The applied source model is thus capable of explaining such electrically silent source components in source estimation, if the measured data supports them. The volume conductor model used in this work is realistically shaped and it contains the major conductivity inhomogeneities, i.e., the lungs and intra-cardiac blood. This corresponds to the level of detail in experimental inverse ECG/MCG studies. We used our BSPM and MCG sensor setups that are similar to most existing measurement setups: typically a BSPM setup covers both anterior and posterior thorax, whereas MCG is measured over the anterior chest only; our omission of posterior MCG sensors or inclusion of posterior ECG electrodes does thus not introduce any comparison bias with respect to the real-world measurements. The forward problem was solved using state-of-the-art linear Galerkin formulation which, with the used models, guarantees numerically stable results.
To do PSF analysis, a linear inverse model is needed. From the family of linear models, MN estimator is the most simple and general one: it poses minimal restrictions to the solution and thus the results reflect the actual contents of the measured topographies better than in other approaches. Other approaches with strong priors, e.g., nonlinear models, likely give much better results in their specific applications (Pullan et al 2010) but it is difficult to assess the relative contribution of prior information and actual measurement in the solution.
Finally, as this is to our knowledge the first study that deals with sensitivity and pointspread of ECG and MCG in realistically-shaped geometry in shared framework, we have focused on basic spatial properties of MCG and ECG and have omitted some issues that should be studied in the future: First, we have not studied the effect of noise or volumeconductor modelling errors. The latter can be studied using PSF/CTF analysis; for examples in brain studies, see e.g. Irimia et al (2013) and Stenroos and Hauk (2013). Geometric model is an essential part of inverse modelling (Pullan et al 2010); for example, Cheng et al (2003) showed that errors in relative size and positions of the heart and torso dominate over other sources of error. Nevertheless, there are still several open questions like the required level of model resolution and fidelity, including the role of anisotropic properties of myocardium in volume conduction. Second, in the case of MCG, the myocardial anisotropy has also another role: The fiber structure of the cardiac muscle gives rise to situations, where near-field measurement of magnetic field reflects a current pattern that may differ considerably from that observed with electric measurement (Holzer et al 2004, McBride et al 2010. In terms of bidomain simulations, this has been connected to an equivalent current for local MCG that flows along the activation wavefront, while the local ECG is associated with a source that is normal to the wavefront Woods 1999, Dos Santos andKoch 2005). The simulations have, however, been carried out in heavily simplified (2D, infinite extent, constant fiber direction) slabs that do not produce realistic volume currents and that do not take into account the effects of twisting fiber-directions on summation/cancellation of far-field MCG signals. To assess the role of cardiac anisotropy in MCG signal generation and MCG and ECG source imaging in general, further measurements and simulations in realistically-shaped anisotropic 3D heart and finite thorax models at realistic source-sensor distances are needed.

PSF analysis
For our analysis, the source discretization must be dense enough to characterize the shape, location, and spread of the PSFs adequately. Judging by the smoothness as well as by the level of detail present in the visualized PSFs, the discretization density along the surface is adequate. Moreover, according to our results, the spatial extent of the PSF is generally much larger than the ventricular wall thickness. This suggests that, with purely spatial approaches, also the transmural depth resolution of ECG/MCG imaging is limited and that constraining the source space on a surface instead of using the myocardial volume is likely an adequate descriptor of measured ECG and MCG; we are currently investigating this.
The differences between MCG and ECG are often discussed in terms of sensitivity to so-called tangential and radial currents (Wikswo Jr 1983, Kwong et al 2013. Apparently, the term 'current' has sometimes been used to refer to the source and sometimes to the volume or total current; this ambiguity is removed in the present study with the clear separation of primary current, i.e., the source, and volume current, enabling quantitative study. Moreover, the exact orientation that is referred to seems to have been rather ambiguous, especially in realistic geometry (where the volume conductor clearly contributes to the signal; see, e.g. van Oosterom et al 1990): it is not clear whether 'tangential/radial' has referred to orientation with respect to sensors, to thorax surface, or perhaps to the heart surface. To clarify this, we defined the source orientations according to the local geometry of the epicardial (source) surface. Thus, our elementary sources are also roughly tangential/normal with respect to the well-conducting blood inside the heart. Our sensitivity results agree with the earlier interpretation that ECG is less sensitive to tangential sources compared normal ones (Wikswo Jr 1983); this is often referred to as the Brody effect. Moreover, the normallyoriented elementary sources correspond to the (epicardial) double-layer source elements used in several studies (van Oosterom et al 1990, Oostendorp and Pesola 2001, Oostendorp et al 2002, van Dam et al 2009, Konttila et al 2014. Our results apply also for (elements of) distributed sources; because of the linear inverse, the contribution of a source element on the estimate is the same regardless of other active sources. The metrics are, however, usually non-linear and cannot be combined linearly like PSFs. Thus, for extended sources, specific metrics should be constructed (see e.g. Konttila et al 2014). In addition to PSF, crosstalk function (CTF) is useful in context of distributed sources. As the PSF and CTF for each elementary source are identical for the MN estimator used here (Liu et al 2002), the results for PSF metrics can be interpreted in terms of crosstalk as well.

Complementarity of MCG and ECG
It could be argued that the superior source estimation performance of the combined MCG+ECG as predicted by our results is explained by the larger number of channels-this is, however, not the case: We recomputed the resolution metrics also using only 75 anterior ECG electrodes and using 352 ECG electrodes positioned at the vertices of the original thorax mesh, and compared them to the results (obtained with 120 electrodes) presented in this paper. While the source estimation performance improved with the increasing number of ECG channels in both scenarios, i.e., in ECG and in MCG+ECG, this improvement was fairly small and the error metric distributions were very similar to the ones presented in the results section (the omission of posterior ECG electrodes resulted in worse posterior estimates). Furthermore, in agreement with the presented results, adding MCG to each electrode set improved the estimates of tangentially-oriented sources clearly but resulted in only a small improvement in the estimates of normally-oriented sources; for tangential sources, even with 75-lead ECG combined with MCG (a total of 174 sensors), the average LE, SD, and especially OE were smaller than in the 352-lead ECG used alone. In addition, we tested adding MCG sensors on the back for a total of 198 MCG sensors and 120 ECG electrodes; this further improved the estimates of posterior tangentially-oriented sources (results not shown); in practice, MCG measurement covering the whole thorax may be possible in the future with, e.g., room-temperature atomic magnetometers (see e.g. Bison et al 2009). The above findings are in line with the results by Oostendorp and Pesola (2001) who showed that adding more ECG leads after a certain number does not yield more information about the source, whereas adding MCG channels will increase the information content.
Several experimental studies have suggested that MCG and ECG contain complementary information (Lant et al 1990, Brockmeier et al 1997, Oostendorp and Pesola 2001: MCG has been shown to be more sensitive than ECG (body-surface potential mapping) to electrophysiological changes related to physical  and pharmacological (Brockmeier et al 1997) stress and prior myocardial infarction (Lant et al 1990); for a recent review on the diagnostic value of MCG, see Kwong et al (2013). It has been interpreted that this difference between the modalities could be due to unique sensitivity of MCG to vortex currents (Brockmeier et al 1997, Liehr et al 2005 or due to relatively higher sensitivity to tangential sources (Kwong et al 2013, Wikswo Jr 1983. While we did not study vortex sources and can thus not comment on their role, our results support at least the latter interpretation. The question is where the tangential sources could arise from. A commonly-used approach is to express the ECG (and MCG) based on uniform double-layer (UDL) hypothesis (Wikswo 1983, van Oosterom et al 1990, Oostendorp and Pesola 2001, van Dam et al 2009. In the UDL model, the source is assumed a double-layer with uniform strength located at the activation wavefront, and this double-layer is replaced by an equivalent double-layer at the surface of the myocardium; the equivalent source is thus normally-oriented with respect to the heart surface. Furthermore, it is assumed that the whole myocardium is homogeneous and activated at some time during cardiac cycle. This kind of equivalent source formulation has been shown to work in activation time imaging where prior temporal information can be effectively utilized (Oostendorp and Pesola 2001). Using the UDL model, van Oosterom et al (1990) argued that the MCG of normal (healthy) ventricular depolarization can be derived from the ECG/BSPM, disputing the complementarity of the modalities. The results of van Oosterom et al (1990) do, however, not contradict with the present work: In the case of normal activation, the electrical activity seems to be adequately explained by normally-oriented equivalent sources (UDL model) Pesola 2001, van Dam et al 2009). Our results show that adding MCG channels to an ECG leadset clearly improves the estimation performance for tangential sources but only slightly for normally-oriented sources. Thus, if the heart is healthy and the UDL model is used in the inverse problem, the benefit of adding MCG may indeed be small. However, the UDL model is clearly not valid in case of ischemia or infarction or any other condition, where the local action potential amplitude is decreased; the presence of (necrotic) scar tissue due to prior infarction may result in severely compromised source estimates (Oostendorp et al 2002). This is because the equivalent double-layer on the heart surface may not be capable of representing the border of the healthy (assumed homogeneous) myocardium around the necrotic area; to correct for this, the source surface would need to be altered to exclude the necrotic tissue (Oostendorp et al 2002). Given that the double-layer is by definition normally-oriented, any alteration to this source surface, such as described above, will lead to tangential (equivalent) source components not included in the source space. Another formulation for surface double-layer source model was presented by Geselowitz (1989): He showed that, with certain assumptions, the double-layer amplitude depends on the local transmembrane voltage, i.e., the assumption of uniformity (as in the UDL model) can be removed. This model allows also the inclusion of repolarization phase in activation time imaging (van Dam et al 2009). Although this kind of surface double-layer (with variable amplitude) is capable of explaining also the tangential sources, the interpretation of the results may be difficult as our example illustrates (see figure 2).
The freely-oriented primary current density source model is capable of expressing different source orientations, including the tangential ones, explicitly. However, interpretation of the source estimates may be challenging , Han et al 2012, Goernig et al 2013, even with the use of prior information such as the spatio-temporal continuity of activity (Han et al 2012). To alleviate this, more information is needed in the process of solving the inverse problem; in this, enhancing the measurement by combining ECG and MCG is a potential candidate. While the inverse of problem of ECG has been in the focus of research in recent decade, the use of MCG is marginal and thus its potential, especially when combined to ECG, is largely unknown. Our results indicate that combining MCG and ECG will improve source estimates significantly. Specifically, our results show that the combined MCG+ECG is superior to ECG in probing tangential sources.

Conclusions
In this work, we studied the spatial properties of MCG, ECG, and combined MCG+ECG in terms of sensitivity and point-spread functions. These properties were analyzed, for the first time, systematically as function of source location and orientation. MCG excels with anterior tangential sources, while ECG sees better the normally-oriented and posterior sources; the sensitivity of combined ECG+MCG is less dependent on source orientation per source location, leading to clearly improved source estimates. For maximal sensitivity, optimal source estimation and largest information content, the electric and magnetic measurements should be combined.