Non-Invasive Electroanatomical Mapping: A State-Space Approach for Myocardial Current Density Estimation

Electroanatomical mapping is a method for creating a model of the electrophysiology of the human heart. Medical professionals routinely locate and ablate the site of origin of cardiac arrhythmias with invasive catheterization. Non-invasive localization takes the form of electrocardiographic (ECG) or magnetocardiographic (MCG) imaging, where the goal is to reconstruct the electrical activity of the human heart. Non-invasive alternatives to catheter electroanatomical mapping would reduce patients’ risks and open new venues for treatment planning and prevention. This work introduces a new system state-based method for estimating the electrical activity of the human heart from MCG measurements. Our model enables arbitrary propagation paths and velocities. A Kalman filter optimally estimates the current densities under the given measurements and model parameters. In an outer optimization loop, these model parameters are then optimized via gradient descent. This paper aims to establish the foundation for future research by providing a detailed mathematical explanation of the algorithm. We demonstrate the feasibility of our method through a simplified one-layer simulation. Our results show that the algorithm can learn the propagation paths from the magnetic measurements. A threshold-based segmentation into healthy and pathological tissue yields a DICE score of 0.84, a recall of 0.77, and a precision of 0.93.


Introduction 1.Motivation
Over the past century, our understanding of the human heart has increased dramatically.Today, we can diagnose and treat many more diseases than a hundred years ago.Despite this, heart disease remains the leading cause of death [1].Cardiac arrhythmias may occur as a consequence of underlying heart disease such as cardiomyopathy or be a primary phenomenon such as in long QT syndrome, an inherited ion channel disorder.Accurately diagnosing arrhythmias often requires using invasive electroanatomical mapping techniques, thereby putting the patient at risk for complications such as bleeding or infection.
Non-invasive methods that provide the same result would spare them the risks associated with invasive measurements [2,3].In combination with Stereotactic Body Radiation Therapy(SBRT) [4], this could potentially allow for completely non-invasive treatment of cardiac arrhythmias [5,6].These methods would also allow measurements in cases where they are not medically indicated, providing an excellent opportunity for research and prevention.

State of the Art
The gold standard for electroanatomical mapping is invasive catheterization with systems such as EAMS, EnSite Velocity ® , and CARTO ® 3.By tracking the position of the catheter and measuring the potential of the myocardium with the catheter head, a computer program reconstructs a map of the local activation time.Medical doctors can then use this information to select the region for ablation.To minimize risk and stress for the patients, mapping and ablation are usually performed during the same procedure.Therefore, an extensive analysis of the data prior to ablation is not possible.The invasive nature of this method also prevents its use as a diagnostic measurement.The reported accuracies of commercial systems are between 1 mm and 2 mm, while the reported threshold for clinical relevance is 3 mm [7].
Non-invasive measurements aim to reconstruct electrical activity inside the myocardium indirectly from data such as electrocardiography (ECG), magnetocardiography (MCG) and magnetic resonance imaging (MRI).While the MRI data is used to construct a spatial model of the heart, the MCG or ECG data is used to reconstruct the electrical activity of the heart.The use of a patient-specific torso model is necessary to accurately calculate body surface potentials from myocardial currents [8][9][10].Since the relative permeability of the human torso does not vary much, not including a torso model incurs fewer errors for forward magnetocardiographic models [11].For this reason, we will focus on magnetic measurements [12], although an extension to include electric measurements might prove beneficial [11,13,14].
There are different approaches to solving the problem of non-invasive electroanatomical mapping.One way to circumvent the prohibitively large number of unknown parameters is to compute the current density distribution from a low-dimensional set of parameters.Gillette et al. [15] use a sophisticated electrophysiological model to compute the distribution of cardiac sources and ECG potentials over time.They then use parameter sampling to match the simulated ECG signals to the measurements.A major advantage of this digital twin approach is that the optimized parameters have well-understood electrophysiological meanings.The main limitation of this approach is the poor scaling of parameter sampling to higher-dimensional parameter vectors.Gillette et al. acknowledge that accurate modeling of substrate-based diseases would require regions of deviating parameters.At the opposite end of the spectrum from digital twins, there are supervised machine-learning approaches.Instead of fitting meaningful electrophysiological parameters to the measurements, models, usually deep neural networks, are trained to predict cardiac sources directly from non-invasive measurements.For example, Chen et al. [9] trained a convolutional neural network on pairs of heart and body surface potentials from five pigs.This approach is not limited to specific types of diseases, such as the cardiac digital twin approach, since no explicit assumptions about propagation pathways or mechanisms are made.The main challenge is acquiring suitable training data.

Contributions
In this paper, we present the mathematical foundation of a novel algorithm to estimate the current density distribution in myocardial tissue from MCG measurements.We demonstrate the capabilities and limitations using a simplified simulation.The medical problem we address primarily is the localization of arrhythmogenic tissue.
Similar to the cardiac digital twinapproach, the parameters of our model have an electrophysiological meaning.Unlike the cardiac digital twin approach, these parameters can be different for each voxel, allowing for inhomogeneities in the myocardial tissue.This is possible because we do not have to sample the parameter space.Intead, we optimize the model parameters using gradient descent.In contrast to supervised machine learning approaches, no access to heart surface potentials is required.

Methods
We divide the myocardium into N v equally sized voxels.For each voxel, the threedimensional current density j , where N m is the number of considered time steps, and v ∈ [0 . . .N v − 1] is the voxel index.The vector j(n) = [j 0 (n) ⊺ , . . ., j N v −1 (n) ⊺ ] ⊺ is the concatination of these current densities, and the matrix J = [j(0), . . ., j(N m − 1)] ⊺ is the vector of current densities over time.Similarly Z = [z(0), . . ., z(N m − 1)] ⊺ denotes all measurements over time.The measurement vector for one time step is written as z(n) = [z 0 (n), . . ., z N s −1 (n)] ⊺ , where N s is the number of individual sensors.The measurements can be written as z(n) = d(n) + s(n), where d(n) is the desired signal and s(n) is the noise signal.
A forward model maps the estimated current densities ĵ into the measurement space through a measurement matrix H, resulting in the estimated measurements ẑ: Obtaining the estimated current densities by computing the pseudoinverse of the measurement matrix H and multiplying it by the measurements z is not feasible: The reason for this is the ill-posed nature of the problem.Instead, the problem of non-invasive measurements is stated as finding the, in a sense, best estimate of J for some given measured signals Z:

Overview
The proposed current density estimation algorithm consists of three main parts: the Model Initialization, the State Estimation , and the Model Refinement (cf. Figure 1).
The Model Initialization uses image data I, and knowledge of the electrophysiology of a healthy heart to construct the initial state space model M(0) .In the State Estimation part, a sparse Kalman filter estimates the current densities Ĵ(i) from the model M(i) and the measurements Z, where the superscript i indicates the epoch.Using the estimated current densities Ĵ(i) and the measurements Z, the new model M(i+1) is calculated in the Model Refinement block.The main loop consists of the steps State Estimation and Model Refinement.These two steps execute alternately.After termination, the final estimated current densities Ĵ(∞) and the final model M(∞) can be used for further analysis.

Forward Model
The forward model describes the mapping of the current densities j into the measurements z.The measurement vector consists of N s magnetic measurements at positions We calculate the magnetic fields using a superposition of Biot-Savart's law for all voxels: where N v is the number of voxels p v,k is the position of the k-th voxel, and V is the volume of the voxels.Using Equation (4) we calculate measurement matrix entries during the Model Initialization.In the State Estimation and Model Refinement step, the measurement matrix is used to calculate the estimated measurements by applying Equation (1).

System Model
The system model describes the prediction of system states j(n|n − 1).The system model needs to facilitate arbitrary propagation paths and propagation velocities to accurately model the propagation of the action impulse through the myocardium [16][17][18].Modeling the interactions between all voxels becomes intractable as the number of voxels increases.Therefore, we use First-order Thiran all-pass filters [19] to model the interaction between neighboring voxels.
These infinite impulse response (IIR) filters have an amplitude response of one over the whole frequency range and a constant group delay for 0 ≤ f ≲ f s /10.Within these limits, arbitrary group delays τ ≥ 1 sample can be achieved by applying the following equations: Equation ( 5) is written for a single all-pass output y.The parameter k is the integer part of the desired group delay τ, while the parameter a determines the fractional part of the group delay [19]: We calculate the predicted current densities by applying: Here, ĵv (n) is the three-dimensional current density in x-, y-, and z-direction in one voxel v.The vector y v ∈ R 243 contains the all-pass outputs that contribute to j v .The sparse matrix C v ∈ R 3×243 contains the 243 gains c, that map the all-pass outputs into the current densities.The influence of the control function u(n) on the current density ĵv (n) is modelled by the corresponding entry in the control vector b v .
The control function determines the shape of the action potential.We use Myokit [20] to calculate an action potential based on the dynamic O'Hara-Rudy model [21].The control function u(n) is calculated by differentiating the action potential and scaling it to a maximum value of 1 A/mm 2 .

Model Initialization
In the Model Initialization, the initial values for the delays τ (cf.Equations ( 5) and ( 6)), the gains c (cf. Equation ( 7)), and the control vector b are computed based on the voxel positions We consider six voxel types ζ in this study: sinoatrial node, atrium, atrioventricular node, His-Purkinje system, ventricle, and pathological.These six types are the minimum number needed to model the electrophysiology of the heart sensibly.They allow for a localized excitation of the heart, a localized connection between the atria and ventricles, a coordinated excitation of the ventricles via the His-Purkinje system, and a disturbed function via the pathological voxels.The assumed connectivities between the voxel types [15,18,22,23] are shown in Figure 2. indicate allowed connections between neighboring voxels.For example, a direct connection from atrium to ventricleis not allowed.Instead, the propagation has to pass through the atrioventricular node and the HIS-Purkinje system before reaching the ventricles.Voxels of type pathological are not connected to other voxels in this model and, therefore, never excited.(Figure adapted from [24]).
First, we calculate the delays τ (0) for all pairs of neighboring voxels based on their positions p and the propagation velocity ϕ of the input voxel type ζ i .For all voxels v i , i ∈ [0. . . ., N v − 1] and their up to 26 neighbors v o , the delay τ (0) i,o is calculated according to: The variable τ (0) i,o describes the delay of the propagation from voxel v o to voxel v i in samples.
Next, the gains C are calculated using an iterative algorithm.We use the symbol C i,o ∈ R 3×3 to denote the nine gain entries that describe the propagation from voxel v o to voxel v i .It is necessary to keep track of the activation times of each voxel.We denote these as σ ∈ R N v .Furthermore, it is necessary to keep track of the current directions in each voxel.We denote these as δ ∈ R N v ×3 .We initialize the activation time of the sinoatrial node with σ sinoatrial = 0 ms and calculate the others during the iterative procedure detailed below.This ensures that the activation in the initial model starts in the sinoatrial node.We arbitrarily choose the current direction of the sinoatrial node to be δ sinoatrial = [1, 0, 0] ⊺ .The entries of the gain matrix C are calculated by the following procedure.

1.
The iteration time t i is set to 0 ms.

2.
For all voxels v o with an activation time σ o equal to the iteration time t i the following steps are carried out.
(a) Identify all neighboring voxels v i , that have not been connected yet.(b) Discard all neighbors with incompatible voxel types (cf. Figure 2).(c) Calculate the current direction in v i according to the following equation, where p i is the position of the voxel v i and p o it the position of the voxel v o : Calculate the gains C i,o between two voxels v i and v o according to the following equation: The indices i and o are dropped for better readability and the indices j, k ∈ {0, 1, 2} are used to index the three by three matrix Calculate the activation time of the voxel v i as Set the iteration time t i to the smallest activation time σ that is bigger than the current iteration time t i .
We repeat steps two and three until we do not find a new activation time in step three.This procedure ensures that fast propagation paths are preferred, no voxel is activated twice, and all voxels that can be connected are connected.
Since the control function is scalar, the control matrix B simplifies to a vector b ∈ R 3N v .We set the entry corresponding to the x component of the current density of the sinoatrial node to one and all others to zero.

State Estimation
The State Estimation uses a sparse Kalman filter to estimate the current densities Ĵ(i) using the current model M(i) and the measurements Z.For each timestep, we carry out the following steps.First, we predict the current densities ĵ(n|n − 1) by applying Equations ( 5) and (7).Next, we predict the state covariance matrix P(n|n − 1) by applying: Since the gains c are defined only for neighboring voxels, the summations include only said neighbors.We denote this with N .Then, the Kalman gain K(n) is computed by applying: We use s −1 m,o to denote the entry (m, o) of the inverse S −1 .After that, we correct the current densities ĵ(n|n) using the measurements z(n) by applying: Finally, we correct the state covariance matrix P(n|n) using the measurements z(n) by applying:

Model Refinement
The Model Refinement uses the estimations of the current densities Ĵ(i) and the measurements Z to adjust the model M (i) , yielding the refined model M (i+1) .To this end, we calculate a loss function L(n) for each time step.Based on this loss we calculate the gradients ∂L ∂c and ∂L ∂a .We then adjust the gains c and the all-pass coefficients a using these gradients.The goal is to minimize the difference between the actual current densities J and the estimated current densities Ĵ.Since we do not have access to the actual current densities J we minimize the difference between ground truth measurements Z and predicted measurements Ẑ instead.We denote this part of the loss as L m : We place an electrophysiologically motivated constraint on the estimated current densities to restain the possible solutions.The absolute current density in a voxel may be lower than the one imposed by the control function u(n) due to the voxel not being fully filled with myocardial tissue or changes of the myocardial tissue filling the voxel [16,25,26].There is, however, no reason for the current density to exceed the nominal maximum value.Therefore, we add a term to the loss that increases once the absolute current density in a voxel || ĵv || 1 exceeds the normalized value of 1.01, and is zero otherwise.We denote this part of the loss as L v : The final loss is the sum of both parts weighted by the regularization strength γ: In order to update the gains c and delays τ, we have to calculate the respective partial derivatives of the loss function.Starting with ∂L m ∂c and applying the chain rule: we arrive at: where y is the all-pass output corresponding to the gain c and h is the column of the measurement matrix H corresponding to the current density j influenced by the gain c (cf. Equation ( 7)).Next, ∂L v ∂c evaluates to: where y is the all pass output corresponding to c, ĵ is the estimated current density corresponding to c, and ĵv are the three current densities, including ĵ corresponding to the voxel influenced by c.For ∂L m ∂a we can again apply the chain rule: arriving at: The summation over the corresponding all-pass outputs y is because the all-pass coefficients are shared across nine all-pass filters that connect two voxels.Calculating ∂y ∂a is done by applying the following iterative equations: where ∂y ∂a = ∂y ∂a IIR + ∂y ∂a FIR ([27], Chapter 15).Lastly, since the parameters τ only influence the propagation velocity, the regularization loss Lv is not used in the calculation of ∂L ∂a .In other words, ∂L v ∂a is set to zero.Using these gradients we then update the gains c according to: and the all-pass coefficients according to where η is the learning rate.The all-pass coefficient a is only valid for values between zero and one.If, after applying Equation ( 27), the coefficient a is less than zero, it is increased by one, and the integer part of the delay (cf.Equation ( 5)) is also increased by one.If, on the other hand, the coefficient exceeds one, we have to consider two cases.If the integer part of the delay is at least one, we reduce the coefficient and the integer part by one.Otherwise, we set the coefficient is set to one.

Simulations
The question we want to answer with the following simulations is if our state-space current density estimation approach is able to localize regions of reduced electrical activity, which would indicate the presence of arrhythmogenic tissue [28].While our approach can learn propagation paths and velocities from magnetic measurements, we focus on the propagation paths in these simulations.The comparison with other sophisticated approaches, as introduced in Section 1.2, remains challenging due to the lack of available source code.Therefore, we compare our results to the pseudoinverse solution (cf.Equation ( 2)) instead.The current density estimation should yield results that facilitate differentiation between healthy and pathological tissue.Therefore, we employ a simple threshold-based segmentation.

Setup
To answer this question, we use an abstracted one-layer heart simulation as depicted in Figure 3.The heart model consists 25 × 37 × 1 (x, y, z) voxels of size 2.5 mm, resulting in a total heart size of 65 mm × 92.5 mm × 2.5 mm (cf. Figure 3a).This model represents the unrolled surface of the heart.The voxel types are assigned based on the approximate shape of the regions in the real heart to generate a morphologically sound sinus rhythm.We differentiate between the ground truth model M, used to generate the data, and the initial model M(0) .While the initial model assumes a healthy heart without any pathological regions, the ground truth model includes a region without any conduction in the right ventricle.We show the difference in voxel types between ground truth and the initial model in Figure 3c.The assumed propagation velocities are 1.1 m/s for the sinoatrial node, atrium, and ventricle, 0.012 m/s for the atrioventricular node, and 4.5 m/s for the HIS-Purkinje system.While these values are simplified, we based them on values reported in literature [15,18,22,23].We run both models for a simulation time of 1 s with a simulation frequency of 2000 Hz.  (d-f) During the time of the maximal difference between ground truth and initial model (t = 0.3675 s).(g-i) The maximum of this L1-norm over the complete cardiac cycle.(j-l) The simulated magnetic fields z around the time of the QRS-complex.The black dashed line corresponds to the time of the maximum difference between the ground truth and the initial model (t = 0.3675 s).

Setup
Figure 3d,e depict the L1-norm of the simulated current densities in each voxel for t = 0.3675 s.At this time, the propagation already passed through the HIS-Purkinje system and is spreading outwards through the ventricles.The propagation is symmetrical for the initial model, while the ground truth has no electrical activity in the pathological region.Figure 3f shows the difference in these L1-norms.The current density in the initial model is higher than in the ground truth model in the region of the pathology.Since we cannot show all time steps in this fashion, Figure 3g-i show the maximum over these L1-norms during the complete cardiac cycle.While this maximum is constant over all voxels in the initial model, the ground truth model has no electrical activity in the pathological region for the complete cardiac cycle.
We use a 4 × 4 × 3 array of ideal magnetic sensors.The sensors are placed equally across a 250 mm × 250 mm × 100 mm region centered 200 mm above the heart.Since we simulate one layer of voxels, the measurements have to be scaled up to match reasonable field amplitudes as reported in [29][30][31].We scale the measurements by a fixed factor of 70 to achieve a maximum R-peak amplitude of 100 pT for the ground truth model.We then superimpose the simulated signals with white noise of an equivalent noise spectral density of 40 fT/ √ Hz, which is in the order of low temperature superconducting quantum interference devices (SQUIDs) [8]. Figure 3j-l depict the magnetic measurements around the time of the QRS-complex.The dashed black line indicates the time of the maximum difference (22.5 pT) between the ground truth and the initial model at time t = 0.3675 s.
The number of voxels is N v = 962, resulting in 2886 system states.The number of sensors is N s = 144, which is on the high end in the number of sensors of available MCG systems [29].

Pseudoinverse
Figure 4 summarizes the results of the pseudoinverse solution.Since the number of unknown values greatly exceeds the number of knowns, the differences between the predicted measurements and the ground truth are vanishingly small.The same, however, can not be said for the current densities.Here, no discernable structure is present.This solution favors using the outer voxels to reconstruct the magnetic measurement, but no segmentation into healthy and pathological tissue is possible based on these results.The difference between the predicted measurements and the ground truth measurements (cf. Figure 3j).(c) Maximum of the L1-norm of the current densities in each voxel over the complete cardiac cycle.(d) The difference between these maxima and those of the ground truth (cf. Figure 3g).

State-Space Approach
Figure 5 summarizes the results of the state-space approach.We optimized the model for 2000 epochs with a learning rate of η = 200 and a regularization strength of γ = 1. Figure 5a shows the loss function for the first 500 epochs.Although the loss decreases slightly until the last epoch to a final value of 5.8 × 10 −5 , the maximum decrease happens in the earlier epochs.The maximum of the measurement differences decreases significantly from around 20 pT to 1 pT.The differences in current densities to the ground truth in pathological tissue regions are significantly reduced, while differences in healthy regions increase slightly.

Current Densities
Figure 5f,g show the results of the simple threshold-based segmentation.The segmentation compares the maximum L1-norm in each voxel to a threshold.We classify the voxel as pathological if the L1-norm is below a threshold.Otherwise, we classify it as healthy.Based on this segmentation, a DICE score is calculated.Figure 5f shows this DICE score for all thresholds from 0.6 to 1.0.A threshold of 0.88 achieves the maximum DICE score of 0.84. Figure 5g shows the results of a segmentation based on this optimal threshold.Most pathological voxels are correctly classified.We achieve a recall of 0.77.Not many healthy voxels are incorrectly classified as pathological, achieving a precision of 0.93.All false positives (FP) are adjacent to the actual pathological region.
Executing one epoch of this algorithm takes around 1 s on an M2 MAX processor, using a single core.The algorithm is mostly converged after 500 epochs or 8 min.All 2000 epochs take around 30 min.

Discussion
We simulated a single layer of tissue leading to a relatively small number of voxels.At the same time, the number of sensors is on the high end for MCG sensor systems.Even in these favorable circumstances, the naive solution of using the pseudoinverse is not feasible.More sophisticated approaches are needed to extract useful information about the current density distribution in the human heart from non-invasive measurements.
While our state-space approach is not able to perfectly reconstruct the ground-truth current densities, the results show a clear structure.Throughout the optimization, the difference between the estimated magnetic measurements and the ground-truth measurements is reduced.The same holds for the current densities.The optimized model shows a reduced maximum current density in the pathological region of the ground truth model.For a perfect reconstruction, the maximum current density in this region would have to be zero.The optimized model only reduces the maximum to around 0.5 A/mm 2 in this region.At the same time, the maximum current density is also slightly reduced in regions that do not correspond to the pathological region in the ground truth model.The most obvious reason for this is that while our state-space approach applies a strong regularization and bias on the estimated current densities, the ill-posedness of the problem still leads to non-unique solutions.An extensive hyperparameter optimization over the learning rate, regularization strength, and Kalman parameters is expected to improve these results.With a DICE score of 0.84, segmentation solely based on the maximum current densities in each voxel is possible.These findings prove, that it is possible to differentiate between pathological and healthy tissue for this simplified model.This does suggest that regions with reduced electrical activity could also be found in real hearts, given sufficient measurements.
There are currently several limitations to our approach.First, we do not explicitly model any secondary currents caused by return paths.The update step of the State Estimation block implicitly calculates the return currents that flow through the myocardium.Our model cannot describe currents outside these voxels.Second, magnetic sensors do not measure the magnetic field perfectly.They only have a limited linear region and a non-flat frequency response [29].Furthermore, they do not measure the magnetic field at one point but accumulate it over their sensing area.Third, there is a gradient in the refractory period in the myocardium [32].Our model cannot capture this aspect of the electrophysiology of the heart.The control function encodes the de-and repolarization phases of the action potential.Since both are propagated using the same parameters, the refractory time cannot be varied depending on the location of the voxel.Since the signal generated by the depolarization phase is dominant over the one generated by the repolarization phase, the expected impact on the localization of arrhythmogenic tissue is low.
In order to apply this algorithm to patient data the following steps have to be taken.First, generating the spatial model of the patient's heart requires a segmented MRI scan.Manual segmentation is feasible for proof of concept studies.For widespread application, automatic segmentation [33] is needed to reduce the time effort.Second, the initial propagation velocities can be estimated using characteristic times of the measured MCG.For example, researchers can use the P-wave's duration and the spatial extent of the atrium to calculate the average propagation velocity in this area.After this, one can apply the algorithm as outlined in this paper.The resulting estimated current densities can then inform the decisions of medical professionals.With an execution time of 30 min, the algorithm is viable in clinical practice.The execution time will increase when applied to a whole heart instead of just one layer.There are, however, still opportunities to speed up the algorithm by using multithreading or fine-tuned code optimization.Another practical concern is that body movements (mainly respiration-induced) cause a subtle time-varying position change of the sensor array towards the heart.However, magnetic localization approaches (as conceptualized in [34]) can provide the position data required to compensate for such artifacts.While most clinics have access to MRI systems, they rarely have access to MCG systems [35].An extension of our algorithm to also work for the more readily available ECG data would necessitate the extension of the forward model, including the consideration of the torso geometry and permittivities [8][9][10].

Conclusions
This paper introduces a novel state-space algorithm for the non-invasive estimation of cardiac current densities from MCG signals, allowing for arbitrary propagation paths and velocities.We demonstrated the algorithm's ability to differentiate between healthy and pathological tissue in a simulated environment.Segmentation with the optimal threshold achieves a DICE score of 0.84.The presented results prove the ability to learn propagation paths from magnetic measurements.An extension of the algorithm to deal with the current limitation and apply it to patient data is feasible.Potentially, the algorithm can non-invasively localize arrhythmogenic tissue in clinical settings after further development.

Figure 1 .
Figure 1.Schematic overview of the current density estimation algorithm.

Figure 2 .
Figure 2. Block diagram of the assumed connectivities between the six considered voxel types and the input function u.The input function u is connected to the sinoatrial node.The other arrows indicate allowed connections between neighboring voxels.For example, a direct connection from atrium to ventricleis not allowed.Instead, the propagation has to pass through the atrioventricular node and the HIS-Purkinje system before reaching the ventricles.Voxels of type pathological are not connected to other voxels in this model and, therefore, never excited.(Figure adapted from[24]).

Figure 3 .
Figure 3. Overview of the used simulation setup.(a) Voxel types of the ground truth model: The colors of the voxels correspond to their type as depicted in Figure 2. A region of pathological tissue is present in the right ventricle.(b) Voxel types of the initial model: The initial model assumes a healthy heart without any pathological regions.(c) The difference in voxel types: red regions stand for differences between the ground truth and initial model.Green regions stand for equalities.(d-i) The L1-norm of the current densities in each voxel.(d-f)During the time of the maximal difference between ground truth and initial model (t = 0.3675 s).(g-i) The maximum of this L1-norm over the complete cardiac cycle.(j-l) The simulated magnetic fields z around the time of the QRS-complex.The black dashed line corresponds to the time of the maximum difference between the ground truth and the initial model (t = 0.3675 s).

Figure 4 .
Figure 4.The results of applying the pseudoinverse solution to the data generated by the ground truth model introduced in Figure 3. (a) The estimated measurements in pT.(b)The difference between the predicted measurements and the ground truth measurements (cf.Figure3j).(c) Maximum of the L1-norm of the current densities in each voxel over the complete cardiac cycle.(d) The difference between these maxima and those of the ground truth (cf.Figure3g).

Figure 5 .
Figure 5.The results of the nested state-space optimization procedure after 2000 epochs.(a) The loss over the first 500 epochs.(b) The estimated magnetic measurements in pT.(c) The difference between these estimates and the ground truth (cf.Figure3j).(d) Maximum of the L1-norm of the current densities in each voxel throughout the cardiac cycle.(e) The difference between these maxima and those of the ground truth (cf.Figure3g).(f) The DICE score for different segmentation thresholds.A threshold of 0.88 achieves the maximum DICE score of 0.84.(g) The predicted voxel types for a threshold of 0.88.The color-coded regions stand for: green-true positive (TP), blue-true negatives (TN), orange-false positive (FP), and red-false negatives (FN).