Impulse data models for the inverse problem of electrocardiography
11 Basis Function Based Data Driven Learning for theInverse Problem of Electrocardiography
Tommy Peng, Avinash Malik, Laura Bear, and Mark L. Trew
Abstract — Objective:
This paper proposes an neural networkapproach for predicting heart surface potentials (HSPs) frombody surface potentials (BSPs), which re-frames the traditionalinverse problem of electrocardiography into a regression problemthrough the use of Gaussian 3D (G3D) basis function decompo-sition.
Methods:
HSPs were generated using G3D basis functionsand passed through a boundary element forward model toobtain corresponding BSPs. The generated BSPs (input) andHSPs (output) were used to train a neural network, which wasthen used to predict a variety of synthesized and decomposedreal-world HSPs.
Results:
Fitted G3D basis function parameterscan accurately reconstruct the real-world left ventricular pacedrecording with percent root mean squared error (RMSE) of . ± . %. The basis data trained neural network was ableto predict G3D basis function synthesized data with RMSEof . ± . %, and G3D representation of real-world datawith RMSE of . ± . %. Activation map produced fromthe predicted time series had a RMSE of 17.0% and meanabsolute difference of . ± . ms when compared to thatproduced from the actual left ventricular paced recording. Conclusion:
A Gaussian basis function based data driven modelfor re-framing the inverse problem of electrocardiography as aregression problem is successful and produces promising timeseries and activation map predictions of real-world recordingseven when only trained using Guassian data.
Significance:
TheHSPs predicted by the neural network can be used to createactivation maps to identify cardiac dysfunctions during clinicalassessment.
Index Terms —Electrocardiogram decomposition, Pro-arrhythmic risk, Gaussian Mesa Functions, Electrocardiogramprediction.
I. I
NTRODUCTION
Predicting heart surface potentials (HSPs) from body surfacepotentials (BSPs) is termed as the inverse problem of electro-cardiography [1]. Achieving good predictions would allow fornon-invasive cardiac surface assessments which can assist insource localization for cardiac procedures such as ablation [2].The forward relationship between HSP and BSP is describedin Eq 1, where A characterizes the geometric relationshipbetween the body and heart surface potential at any instancein time. The problem of estimating A − given instances ofconcurrent BSPs and HSPs is known to be ill-posed and Tommy Peng is with the Department of Electrical and Computer Engineer-ing, University of Auckland. E-mail: [email protected] Malik is with the Department of Electrical and Computer Engi-neering, University of Auckland.Laura Bear is with IHU-LIRYC.Mark L. Trew is with the Auckland Bioengineering Institute, University ofAuckland.Manuscript submitted November 10, 2020. This work was funded by theNational Science Challenges, Science for Technological Innovation, NZ, grant3713917. Mark L. Trew was supported by a grant from the Fondation Leducq. therefore difficult to solve [1], [3]. Current approaches suchas electrocardiographic imaging (ECGI) [4], [5] estimates theHSP ( ˆ HSP λ ) by minimizing the error or L2-norm residualwith an additional a penalty term λR as shown in Eq 2, where λ is the regularization parameter and R is the regularizor. Reg-ularization parameters must be chosen carefully, as too largeof a regularization parameter results in overly smooth HSPpredictions, and too small of a regularization parameter resultsin oscillating HSP predictions [1], [6]. ECGI resolves thisissue by using a number of electrophysiological constraintsand special smoothing operators [5]. While methods exist forestimation of the quasi-optimal regularization parameter, thereis no single estimation technique which performs best for allgeometries and signal-to-noise ratios (SNRs) [7], [8], [3]. BSP = A × HSP (1) ˆ HSP λ = || BSP − A × HSP || + λR (2)Ultimately, predicting HSPs from BSPs can be framed as aregression problem [9], where we find some approximation for A − given simultaneously observed BSP-HSP pairs. Here, we hypothesize that the prediction of HSPs from BSPs (tradition-ally, the inverse problem) can be framed as a forward problemand solved by using a data driven model. Neural networks havebeen shown to be effective at solving data driven regressionproblems [10]. During training, neural networks implicitlylearn the linear or non-linear relationship between HSP andBSP. This data driven approach for solving the inverse problemof electrocardiography is successful in recovering HSPs fromBSPs [6], [11], [9].However, a common problem encountered by the neuralnetwork approach is the richness of heart states found in thetraining data set [11], [2], [12]. Neural networks can onlypredict data similar to that found in the training set, thereforetraining a neural network to solve the cardiac inverse problemrequire a large amount of BSP and HSP recordings acrossdifferent heart states. Unfortunately, detailed simultaneousHSP and BSP recordings from real hearts in different heartstates are limited. As a solution, we propose to decompose theHSPs into a basis function across heart conditions and trainthe neural network with examples from that common space.Decomposition is important because while a particular heartstate HSP may not be in the training set, the decomposed HSPbasis function will likely have been. To that end, Gaussianbasis functions have been shown to be effective at capturingcardiac signals [13], [14]. These basis functions can be usedto model and generate real-world-like cardiac signals underdifferent drug and disease states [15], [16]. Furthermore, a r X i v : . [ q - b i o . Q M ] F e b boundary element forward models using recording electrodegeometries have been relatively successful in capturing theHSP to BSP relationship (Laura/Mark’s forward model). Thisrelationship allows us to generate training data pairs for theneural network.We aim to generalize this solution for heart surface poten-tials using signal decomposition, signal projection, and neuralnetworks. Our contributions are: 1. we show that HSPs canbe effectively modelled and generated using Gaussian 3D basisfunctions; 2. we demonstrate that a single neural networktrained on Gaussian basis functions can be used to predictnever seen before physiologically relevant synthesized HSPsfrom BSPs; 3. we extend the prediction of HSPs from signalsrecorded in the real-world. We show that a data drive modeltrained using a suitable basis function (HSPs) and forwardmodel (BSPs) can be effective for prediction of real-worldHSPs from BSPs. In essence, our approach is predicting theoutput of a system, by learning the system characteristicsthrough basis function driven responses.II. M
ETHODS
Figure 1 is a visual representation of the Gaussian basisneural network pipeline. Note here that the forward model isonly being used for data generation purposes. The training andtesting data can ultimately be replaced by recordings from realworld experimental recordings.
A. Data Context
The original heart-in-tank data was recorded by Bear et.al. at IHU-LIRYC [17]. In the experiment, an excised pigheart was perfused in Langendorff mode, with an epicardialelectrode sock of 108 electrodes attached to the ventricles. Theheart was placed inside a human-shaped torso tank with 128body surface electrodes. Simultaneous body and heart signalswere recorded. The geometry of the epicardium and bodysurface was extracted from CT full form images.Heart and body surfaces are abstracted to regular 2D matrixof points. These provide a standard context for training andapplying neural network data models. They are also consistentwith the topology of electrode matrices used for electricalmapping on both the heart and body surfaces [17]. Conse-quently, the heart surface is represented by a grid of 9-by-12points and the body surface a grid of 16-by-16 points. Thebody surface grid is arbitrary and has more points than typ-ical numbers of body surface electrodes. See SupplementaryMaterials for more details on the 2D abstractions.
B. Neural Network Data Model1) Neural Network Design and Implementation:
A neu-ral network with a dense fully connected one hidden layerarchitecture was implemented in Python using Keras witha Tensorflow back-end. This architecture is mathematicallyproven to be capable of approximating any smooth continuousfunction given enough neurons and appropriate activationfunctions [18], [19]. The network is designed to capture therelationship in the x - y plane between 2D Gaussian basis HSP and its simultaneous 2D BSP. Therefore, the input layerconsists of 256 neurons corresponding to the 16-by-16 bodysurface 2D abstraction points, and the output layer consistsof 108 neurons corresponding to the 9-by-12 heart surfaceelectrodes. Neuron bias was turned off, as we expect a zeroBSP input to produce a zero HSP output.
2) Training Data Generation:
The Gaussian basis func-tion [13], [15] and Gaussian mesa functions [14] have beenused in the past to successfully decompose cardiac signals.The decomposition process allows for cardiac signals fromdifferent cardiac states to be expressed in a common parameterspace, which allows for objective representation and quanti-tative comparison of signals. Here, we consider a Gaussian3D basis function (G3D) for describing 3D cardiac signalmatrices as it is smooth and continuous in both space andtime. We define the G3D to be Eq 3, where x and y are the2D abstraction dimensions, t is the the normalized heartbeattime in the recording. Normalized heartbeat time at a samplecan be calculated by dividing the current sample number bythe total number of samples in the recording of the heartbeat,and can vary between 0 and 1 for a heart beat. Parameter A isthe amplitude, µ i is the location of the Gaussian peak in the i dimension ( i ∈ x, y, t ), σ i is the standard deviation or width ofthe Gaussian in the i dimension ( i ∈ x, y, t ) expressed relativeto the size of that dimension. Here, we limit σ x to be equalto σ y to reduce the number of impulses needed to train theneural net. G D ( x, y, t ) = A × exp ( − ( ( x − µ x ) σ x + ( y − µ y ) σ y + ( t − µ t ) σ t )) (3)HSP training data was synthesized from G3D basis functionparameters. The G3D training data included components withpeaks ( µ ) at each combination of integer value of the x (i.e. 1-9) and integer value of the y (i.e. 1-12) dimension. The σ x and σ y at each peak location was allowed to vary in 0.02 intervalsbetween 0.01 and 0.5. µ t was set to 0.5 and σ t was set to 0.3for all basis components, this was done so that the trainingdata would always have a non-zero x - y plane 2D surface withthe peak located at the temporal center. Amplitude ( A ) wasallowed to be 1 or -1, due to the real-world HSPs exhibitingboth positive and negative potentials for a single heartbeatrecording. The training data generation used all combinationsof given µ , σ , and A values. This generated a total of 129600 x - y plane 2D HSPs.The calculation of heart surface potentials given bodysurface potentials uses a volume conductor torso model [20].The model is based on a heart in a torso tank experimental setup [17]. The specific model is not a critical component of thiswork and could be replaced by other models or experimentalheart and body surface recordings. The 129600 HSPs sampleswere passed through the volume conductor forward modelto produce 129600 BSP-HSP pairs in the training set, eachcontaining a Gaussian basis HSP and its simultaneous BSP.
3) Training Regime:
A 70-30 training-validation split wasapplied to the training data. A 0.1 dropout rate was im-
Fig. 1:
Summary neural network pipeline. plemented to reduce over-fitting. The hyperparameters forthe neural net were tuned through grid search of possiblehyperparameter pairings. The grid search was performed inPython with
GridSearchCV from scikit-learn. The hyperpa-rameter search space can be found in Table I, with besthyperparameters highlighted in green. Training was done usingthe Adaptive Momentum Estimation (ADAM) optimizer withroot mean squared loss function. Training was completed over100 epochs.
TABLE I:
Hyperparameter Search SpaceHidden Neurons [100, 150, 200, 250, 300, 350, 400]Activation Function [ReLU, Sigmoid, TanH]Regularization [None, L1, L2]Batch Size [16, 32, 64]Learning Rate [0.1, 0.01, 0.001]
C. Predictions using Parameterized Neural Networks
Prediction was done using the testing data sets, which were not presented to the neural network during training. Duringprediction, the input to the neural network is the 2D BSPsurface, and the output from the neural network is the 2DHSP surface. All predicted HSPs were evaluated using error,which was calculated per 2D potential surface ( x - y plane) byusing percent root mean squared error (RMSE) found in EQ 4,where ˆ p is the predicted or model generated 2D signal surfaceand p is the ground truth. RM SE (ˆ p, p ) = 100 (cid:112) mean ((ˆ p − p ) ) max ( p ) − min ( p ) (4)
1) Pure Moving Gaussians:
In order to create a testing setwith physiologically relevant generated Gaussian based timeseries, we synthesized 3D matrices of 9-by-12-by-100 sizein which the peak of the Gaussian in the x , y dimensionsmoves in a straight line between two points (p1 to p2) on theheart electrode grid. The combination of points can be foundin Table II, and 2D contour representation of the Gaussianbasis function at different points along the 100 samples inthe t dimension for path 1 can be seen in Fig 2. The σ x and σ y was set to 0.1768 for the synthetic HSPs for thismoving Gaussian test data set. The paths were designed toreflect different pacing locations commonly found in clinicalsettings such as sinus, right ventricle, and left ventricle. Thesynthesized HSPs in the testing set were passed through thevolume conductor forward model for their respective BSPs.This resulted in 600 x - y plane 2D surface BSP-HSP pairs inthe pure Gaussian testing set. These 2D surface BSPs were used as the input for the neural network to predict 2D surfaceHSPs. TABLE II: x , y Dimension Location of Test I PointsPath 1 2 3 4 5 6p1 [1, 1] [1, 6] [1, 12] [4.5, 1] [4.5, 12] [9, 1]p2 [9, 12] [9, 1] [1, 6] [9, 12] [9, 1] [9, 12]
Fig. 2:
Synthesized moving Gaussian data path 1 (from point [1, 1] to [9,12]) with contours of the moving Gaussian. a) first sample of the moving pathwith peak centered at point 1, 1; b) 50 samples into the moving path withpeak centered at point 5, 6.5; c) last sample of the moving path with peak at12, 12.
Here, we explore the effects of modifying the testing andtraining data on the predicted HSPs in the following ways:1)
Noise Robustness Test : This is a modification on thetest data. Signals recorded in the real world are oftennoisy, and that noise often contributes to poor predic-tions from neural networks. Similar to [12], forwardmodel BSPs from the set of moving pure Gaussiansignals were augmented with Gaussian white noise withMATLAB’s awgn at 2dB, 5dB, 10dB, 20dB, and 50dBsignal-to-noise-ratio (SNR), to create an additional ro-bustness test dataset to observe neural net sensitivityarising from errors in BSP measurement.2)
Geometry Robustness Test : This is a modification onthe test data. Often times during real-world experiments,there is noise in the geometry measurements of the heartand body. We simulate changes or measurement errorsin heart body geometry through: 1. shifting the heartlocation within the forward model towards the back ofthe body at 5mm increments between 0 and 40mms; 2.rotating the heart within the forward model at ◦ in-crements between − ◦ and ◦ . The rotation axis wasfound using principle component analysis on the heartelectrode locations in MATLAB. A positive rotation isclockwise along the axis, while a negative rotation iscounterclockwise along the axis. The subsequent alteredforward models were then used to produce BSPs from HSPs along all 6 paths of the moving Gaussian testingset.3)
Reduced Training Set Test : Here, explore how reducingthe training set would effect the prediction results. It iswell known that an effective training data set for neuralnetworks includes examples from all expected predictionoutcomes [21]. We perform 3 reduction schemes for the x - y dimension peak location as shown in subplots b, c,d of Fig 3. We split the original G3D training set into 25different σ reduction schemes consisting of training setswith a single x - y dimension σ at 0.02 intervals between0.01 and 0.5 for HSPs. Each of the reduced training setswere used to train a separate neural network. Fig. 3:
The locations of training set G3D peak ( µ ) locations in the x - y (9-by-12) plane of the HSP surface. Red dots indicate the points in the x - y planewith a training set G3D peak. a) original training set from Section II-B1; b)reduction scheme 1; c)reduction scheme 2; d) reduction scheme 3. The prediction results from the three data modification testswere compared against the pure moving Gaussian predictionsfrom the neural net trained on the full training set. The resultsfrom the geometry robustness and reduced training set testswere were evaluated by comparing the difference in signalpeak locations in the 2D x - y plane.
2) Experimental Recorded Data:
In this study, we makeuse of the left ventricular (LV) paced averaged beat recordingfrom the experimental data [17] which is 649 samples in lengthrecorded at a sampling rate of 2048Hz. This resulted in an LVHSP recording matrix of 9-by-12-by-649 after 2D abstractionat each sample in time. The 9-by-12-by-649 LV HSP matrixwas fitted to G3D basis functions through an extension of thegeneralized orthogonal forward regression (GOFR) technique[14], [16]. The potential values were normalized between -1and 1 while maintaining scale. To fit each G3D component,the current unexplained signal was first correlated against alibrary of G3D basis functions. A description of the librarycan be found in the Supplementary Materials. Then, thelibrary function with the largest absolute correlation valuewas set as the initial point for loss function minimizationusing MATLAB’s fmincon . Sum of squared differences (Eq 5),where v is the 1D form of the 3D matrix to be fitted and ˆ v is the 1D form of the current fit of the 3D matrix, was usedas the loss function. The bounds of optimization are shown inTable III. TABLE III:
Table of fmincon Optimization Bounds
Bound µ x µ y µ t A σ xy σ t Upper 9 12 1 Max pot. 0.5 0.5Lower 1 1 0 Min pot. 0 0 The bounds on µ limits the peak of the G3D components toalways be contained within the recording, the A bounds allowthe amplitude of the G3D peak to vary between the range ofpotentials within the recording, and the σ bounds allow fittedG3D components to have visible effect throughout the wholerecording. The bounds chosen here are similar to those foundin [16], which have been shown to be effective for describingelectrocardiographic signals. The fitted G3D component wasthen subtracted from all basis functions in the G3D library andthe unexplained signal. The process iterated until the numberof predefined G3D components were fit to the 3D recordingmatrix. The final G3D fit of the recording was evaluated per2D potential surface ( x - y plane) using percentage root meansquared error ( RM SE ) where the range ( max ( p ) − min ( p ) )is 2 due to the normalization of potential between -1 and 1. SSD ( v, ˆ v ) = (cid:88) ( v − ˆ v ) (5)The left ventricular averaged heartbeat recording HSP wasfitted with 100 G3D components. The 649 time instances of x - y plane 2D HSP surface for each of the 100 G3D compo-nent were individually passed through the volume conductorforward model to simulate their respective BSPs. This testingset created a testing set of 64900 x - y plane 2D surface BSP-HSP pairs. These 2D surface BSPs were used as the inputfor the neural network to predict 2D surface HSPs. Predictedreal-world HSPs were evaluated via activation time, which isdefined to be the time point for the most negative gradient(dV/dt). III. R ESULTS
A. Prediction of Moving Pure Gaussian HSPs
This section shows the prediction results for the movingpure Gaussian testing set. Fig 4 show the root mean squarederror (RMSE) between the predicted and generated puremoving Gaussian HSPs. Across the 6 paths from Table II,the mean RMSE is . ± . %. An example of time seriescomparison for path 3 can be seen in Fig 5. B. Robustness Tests
For the noise robustness test, the moving Gaussian signalaugmented with white noise (Section II-C1) is supplied as theinput to the basis function data trained neural network, withthe generated moving pure Gaussian HSP set as the expectedprediction outcome. The RMSE across all 6 paths at differentnoise levels can be found in subplot a of Fig 6. As expected,the error associated with the prediction decreases as the noiselevel decreases. Subplot b from Fig 6 shows an example of theHSP time series predicted based on BSPs with varying noiselevels.For the geometry robustness test, the prediction is carriedout using the neural network trained with the training setwith BSP-HSP pairs found with the forward model with noshift or rotation; however, the testing set is produced withforward models with rotations or shifts in the heart locationas described in Section II-C1. Subplot a of both Fig 7 andFig 8 show that the error of prediction increases as the
Fig. 4:
RMSE for the predicted 2D sample slices (along x-y axis) for the G3Dsynthesized HSPs moving from point to point. The path numbers correspondto those found in Table II. change in heart geometry increases. Subplot b of Fig 7 showsthat the predicted time series share similar morphology, butare increasingly different in amplitude for as shifts increasebetween 0mm and 40mm. Subplot b of Fig 8 shows that anegative rotation angle increases the amplitude and decreasesthe Gaussian width of the signal in the time domain, while apositive rotation angle decreases the amplitude and increasesthe Gaussian width of the signal in the time domain.
C. Effects of Reduced Training Sets
For peak location reduction schemes, as an intuitive guide,Fig 9 shows a subsection of the electrode grid formed by fourneighbouring electrodes (shown as red circles). The trainingset includes examples of G3D basis functions with peaks atthe electrodes. The application of the peak location reductionschemes shown in Fig 3 can then be thought of as increasingthe distance between neighbouring electrodes ( d x and d y ). InFig 9, if we want to predict a HSP Gaussian with peak thatis not at any of the electrodes (shown as the yellow star), themaximum Euclidean distance error in peak location that canbe made for the prediction is (cid:113) d x + d y , which is the largestEuclidean distance to any of the four neighbouring electrodesfrom any point within the grid.Subplot a of Fig 10 are the box plots for the Euclideandistance between predicted and expected peaks for the syn-thetic moving Gaussian testing data set with different reductionschemes. The proposed euclidean distance error limit for eachreduction scheme is shown as the blue line, note how none ofthe 600 total samples along the 6 paths exceed this limit inall x - y dimension peak reduction scheme. Subplot b of Fig 10shows that the RMSE increases as the number of examplesfound in the training set decreases.Individual neural networks were trained using the 25 dif-ferent σ reduction sets. Subplot a of Fig 11 shows the meanand standard deviations for the Euclidean distance betweenpredicted and expected Gaussian peak locations as a function of the σ in the training set for the 25 different neural networks.The σ for the moving Gaussian in the testing set is 0.1768(shown as red star). Note that the in both subplots of Fig 10prediction results become better as the σ for the reducedtraining set approaches 0.1768. D. G3D Decomposition of Recorded HSPs
The 3D matrix of size 9-by-12-by-649 for the left ventric-ular recording was fitted using 100 G3D components. G3Dreconstruction of the recording was done by summing the100 components. Figure 12 shows fitted G3D isosurfaces forthe recording, with further examples in the SupplementaryMaterials. Over the 649 HSP recording samples, the recon-structed signal from 100 G3D components deviates no morethan 5% from the normalized recording when comparing thepotential surfaces (9-by-12 x - y axis surface), with a meanof . ± . %. The isosurfaces of the first 3 fitted G3Dcomponents are shown in the bottom row of Figure 12. Notehow the 3 components are a part of the summation process tofrom the final 100 component G3D representation. E. Prediction of Recorded HSPs
Subplot a of Figure 13 summarizes the RMSE per 2Dabstraction slice ( x - y plane) between the neural network HSPprediction and expected G3D HSP fit. Here, we expect theG3D fitted HSPs to be the ideal outcome from the neural net,as the testing input to the neural net is the BSPs generatedby the forward model from the G3D components of the leftventricular HSP recording. The mean RMSE is . ± . %.Subplot b from Figure 13 contrasts the predicted HSP timeseries with its corresponding G3D fit, and normalized actualrecording. Figure 14 compares the activation times derivedfrom predicted HSPs and recorded HSPs. The mean abso-lute difference in activation times across all electrodes was . ± . ms. The RMSE in predicted activation times is17.0%. IV. D ISCUSSION
In this paper, we present a proof of concept for a basisfunction based neural network approach for solving the in-verse problem of electrocardiography. We have shown that aneural network can be used to modeling the data relationshipsbetween electrical signals on the body and heart surfaces,for a simplified torso model. We trained the neural networkwith static Gaussian basis and used the neural network topredict synthetic signals and recorded heart surface signals.In the context of our problem the method is robust to noiseand perturbed heart locations. This is an important first stepfor developing a new approach to the inverse problem ofelectrocardiography.We have shown that a G3D basis model can effective capturethe morphological behavior found within HSPs expressed as a3D matrix ( x , y , t ). The G3D basis model offers description ofthe HSP 3D matrix which improves with the number of G3Dcomponents fit (Supplementary Materials), this is a commonbehavior for basis function models [22]. Due to the smooth Fig. 5: a) Path 3 on the heart surface 9-by-12 grid, along with the time series correlation between predicted and synthetic data at each electrode. The electrodewith the worst prediction in correlation is marked with a red triangle, while the electrode with best prediction in correlation is marked with a red star. b)Time series prediction results for best prediction correlation at electrode (x = 3, y = 11) for path 3 in Table II. c) Time series prediction results for worstprediction correlation at electrode (x = 11, y = 11) for path 3 in Table II. nature of the Gaussian basis, its representation of a givenrecording is fair but not exact, as seen in Figure 12. This issimilar to the behavior noted in [16], where small componentselectrocardiographic signals were sometimes not captured bythe Gaussian model. The G3D basis function decompositiontechnique proposed here is not signal dependent, and cantheoretically decompose, and therefore express, HSPs fromvarying heart states in the same basis function space.The G3D basis function impulses on the heart surface alongwith calculated (through the forward model) simultaneousrecordings on the body surface allowed for training of a neuralnetwork which was able to predict both synthesized HSPs andreal-world decomposed HSPs from BSPs. The network per-forms well when predicting the synthesized moving Gaussiandata set described by Table II, with low mean RMSE value of8.46% across all 6 paths. As the neural net was only trainedwith basis functions having peaks at integer x and y dimen-sion locations, the prediction error and euclidean distance isexpected to fluctuate in Fig 4 along all 6 paths which includepeaks at non-integer x and y dimensions. Interestingly, subplota of Fig 5 suggests that the time series at electrode locationsthat are closer to the points of activation (i.e. along the path)achieve better prediction results when compared to electrodeswhich are far from the point of activation. This is most likelydue to the distant electrodes having a close to zero time series(subplot c of Fig 5), which is difficult to predict for an neuralnet largely trained on impulse or non zero time series data.The trained neural network is susceptible to noise in theinput BSPs as seen in Fig 6. Note that there is a large dropin predicted HSP quality when going from BSP inputs with20dB to 10dB signal-to-noise ratio in subplot a of Fig 6.Furthermore, subplot b of Fig 6 suggests that, as expected,noisier BSPs will produce noisier HSPs as neural networkoutputs. However, this is expected as network is trained purelyon G3D basis function data without noise. Methods to improverobustness of predictions at signal-to-noise ratios smaller than20dB will be explored in the future. In subplot a of Fig 7,it shows that the neural network prediction results becomeincreasingly erroneous when predicting testing set HSPs from hearts that are increasingly distant from the spatial locationof the heart used to generate the training set. Encouragingly,the predicted time series in subplot b of Fig 7 show thatthe shape of the predicted HSP signal is largely unchangedbetween 0mm and 40mm heart location shifts. Furthermore,in both Fig 7 and Fig 8, there is a non-random change inpredicted signal amplitude, peak time, and width as the heartshifts and rotates. Therefore, in the future, there is a possibilityof training a generic neural network for solving the cardiacinverse problem for different heart and body geometries byincorporating the geometry information as an input.Here, we also observe several important factors for con-structing an effective training set for prediction of Gaussianbasis HSPs using neural nets. We evaluated the effectivenessof training sets based their effect on the Euclidean distancebetween predicted and expected Gaussian peaks. This is animportant measure, as most interesting heart behaviors occurwhen the electrical potential on the heart surface is non-zero.We show that the maximum Euclidean distance error betweenpredicted and expected Gaussian peaks on the heart surfaceis related to the distance between Gaussian peaks found inthe training set. Specifically, if the training set has examplesof peaks at distance d x in the x dimension of the heartsurface, and d y in the y dimension of the heart surface, thenthe maximum error in Euclidean distance of the predictedpeak is (cid:113) d x + d y . This limit is made for the proposedmodel under the assumption that the training set HSP-BSPpairs at neighbouring electrodes are sufficiently different, andthe difference between training set HSP-BSP pairs betweenneighbours’ neighbours are larger than for neighbours alone.In Fig 11 we show that the neural network trained on thefull training set outperforms any of the neural nets trainedon reduced sets with G3D basis functions having a single σ .Therefore, for our proposed data driven model, an effectivetraining set must include a fair representation of data fromvarying σ values. This is similar to the guidelines for effectivetraining sets proposed by [21].The neural network prediction of G3D decomposed HSPleft ventricular recordings is especially encouraging, as the Fig. 6: a) Boxplot of root mean squared (RMSE) value across 2D HSP slices( x - y plane) between predicted and expected HSPs at different noise levels forthe input BSP. b) Time series example of predicted HSPs at different BSPnoise levels for path 1 from Table II. The expected prediction is the G3Dgenerated moving HSP (blue dotted line). neural network is trained on the G3D basis function HSPsand forward model BSPs solutions, and not via a data set ofsimultanteously recorded BSPs and HSPs. While the RMSEseen in subplot a of Fig 13 is larger than that found inpredictions of synthesized moving G3D HSPs in Fig 4, theyare similar to those found for HSPs predicted by other neuralnetworks trained on real-world HSP and BSP data [11], [12].The dynamic behavior of the recorded HSP is predicted well.This can be seen in the volatile activation within the first200 samples, the relatively flat period between 200 and 400samples, and the repolarization period after 400 samples inwhich the potential returns to baseline for both the recordedand predicted time series in subplot b of Fig 13. This is furthersupported by the mean absolute difference in predicted andexpected activation time being 10.7ms.Here, we present the prediction of real-world signals basedon neural networks trained with basis function data to serveas a case study for how such an approach can be usedto solve clinical problems. Our model can be replaced bymeasurements. When the developing the method a model reduces some uncertainties. Based on the viability of the datadriven method shown here, it is potentially feasible to bypassthe forward model for calculating BSPs and train the neuralnetwork based on recordings of real world basis function HSPsand their corresponding BSPs in a heart-in-tank experimentin the future. Here, we find that the quality of predictionswill increase as the number of HSP-BSP samples pairs in thetraining set increases. In particular, we show that there is amaximum error in Euclidean distance of predicted Gaussianpeaks. This limit can be used to calculate the appropriatedistance between stimulus electrodes for real-world HSP-BSPtraining set recordings to train the proposed model to produceHSP predictions at a desired error level. Furthermore, basedon results from the decomposition of HSPs, it seems feasiblethat BSPs from clinical recordings can be decomposed viafitting techniques into basis functions. These decomposed BSPcomponents can be useful for prediction of clinical HSPs.V. C ONCLUSION
In this work we propose a Gaussian basis function decom-position approach to bypass the rich training data problemof neural network approaches for the inverse problem ofelectrocardiography. We have shown that a network trainedon purely generated basis function HSPs and their respectivegenerated BSPs can be used to predict real world recordings ofHSPs decomposed into the same basis function set. Activationmaps of the HSPs, predicted from BSPs with a neural network,can be used to identify cardiac dysfunctions during clinicalassessments. In future works, this method would be tested onreal world HSP Gaussian recordings and their respective BSPswhen the data becomes available.R
EFERENCES[1] B. J. Messinger-Rapport and Y. Rudy, “Regularization of the inverseproblem in electrocardiography: A model study,”
Math. Biosci. , vol. 89,no. 1, pp. 79–118, 1988.[2] N. Zemzemi and R. Dubois, “A machine learning regularization ofthe inverse problem in electrocardiography imaging,”
Comput. Cardiol.Conf. , pp. 1135–1138, 2013.[3] A. J. Pullan et al. , “The Inverse Problem of Electrocardiography,”
Compr. Electrocardiol. , pp. 299–344, 2010.[4] S. Ghosh and Y. Rudy, “Application of L1-Norm Regularization to Epi-cardial Potential Solution of the Inverse Electrocardiography Problem,”
Annu. Biomed. Eng. , vol. 37, no. 5, pp. 902–912, 2009.[5] R. Yoram and J. E. Burnes, “Noninvasive electrocardiographic imaging,”
Annu. Noninvasive Electrocardiol. , vol. 4, no. 3, pp. 340–359, 1999.[6] A. Malik et al. , “A machine learning approach to reconstruction of heartsurface potentials from body surface potentials,” in , 2018.[7] P. Colli-Franzone et al. , “Finite Element Approximation of RegularizedSolutions of the Inverse Potential Problem of Electrocardiography andApplications to Experimental Data,”
Calcolo , vol. 22, no. 1, pp. 91–186,1985.[8] ——, “A mathematical procedure for solving the inverse potentialproblem of electrocardiography. analysis of the time-space accuracyfrom in vitro experimental data,”
Math. Biosci. , vol. 77, no. 1-2, pp.353–396, 1985.[9] N. Zemzemi et al. , “From body surface potential to activation mapson the atria: A machine learning technique,”
Comput. Cardiol. (2010). ,vol. 39, pp. 125–128, 2012.[10] T. Y. Kwok and D. Y. Yeung, “Constructive algorithms for structurelearning in feedforward neural networks for regression problems,”
IEEETrans. Neural Networks , vol. 8, no. 3, pp. 630–645, 1997.[11] K. Bujnarowski et al. , “CT-Scan Free Neural Network-Based Recon-struction of Heart Surface Potentials from ECG Recordings,” pp. 1–4,2020.
Fig. 7: a) Boxplot of Euclidean distance between predicted and expected peaks of HSPs at different spatial shifts for the heart geometry in the forward model.b) Time series example of predicted HSPs at different shifts in heart geometry for path 3 from Table II. The peaks of the expected outcome from the testingset and the predictions are labelled as red stars. The peak of the expected outcome is also marked as the red dotted vertical line for comparison.
Fig. 8: a) Boxplot of Euclidean distance between predicted and expected peaks of HSPs at different rotation angles for the heart geometry in the forwardmodel. b) Time series example of predicted HSPs at different shifts in heart geometry for path 3 from Table II. The peaks of the expected outcome from thetesting set and the predictions are labelled as red stars. The peak of the expected outcome is also marked as the red dotted vertical line for comparison.
Fig. 9:
A small section of the HSP electrode grid, showing 4 neighbouringelectrodes. Training set for the neural net contains HSP-BSP pairs with HSPshaving Gaussian basis function peaks at the the electrodes (shown as redcircles). The expected peak of the predicted HSP is somewhere within thegrid of four electrodes (shown as the yellow star). The distance betweenelectrodes in the x direction is d x and the distance between electrodes in they direction is d y .[12] A. Karoui et al. , “A Spatial Adaptation of the Time Delay NeuralNetwork for Solving ECGI Inverse Problem,” Lect. Notes Comput. Sci.(including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics) ,vol. 11504 LNCS, pp. 94–102, 2019. [13] P. E. McSharry et al. , “A dynamical model for generating syntheticelectrocardiogram signals,”
IEEE Trans. Biomed. Eng. , vol. 50, no. 3,pp. 289–294, 2003.[14] F. Badilini et al. , “Automatic analysis of cardiac repolarization morphol-ogy using Gaussian mesa function modeling,”
J. Electrocardiol. , vol. 41,no. 6, pp. 588–594, 2008.[15] T. Peng et al. , “Parametric Modeling of Electrocardiograms usingParticle Swarm Optimization,” in , 2018, pp. 2695–2698.[16] ——, “Predictive modeling of drug effects on electrocardiograms,”
Comput. Biol. Med. , vol. 108, pp. 332–344, 2019.[17] L. R. Bear et al. , “Cardiac electrical dyssynchrony is accurately detectedby noninvasive electrocardiographic imaging,”
Hear. Rhythm , vol. 15,no. 7, pp. 1058–1069, 2018.[18] G. Cybenko, “Approximation by Superpositions of a Sigmoidal Func-tion,”
Mathematics of Control, Signals, and Systems , vol. 2, pp. 303–14,1989.[19] K. I. Funahashi, “On the approximate realization of continuous mappingsby neural networks,”
Neural Networks , vol. 2, no. 3, pp. 183–192, 1989.[20] L. R. Bear et al. , “Forward Problem of Electrocardiography: Is ItSolved?”
Circulation: Arrhythmia and Electrophysiology , vol. 8, no. 3,pp. 677–684, 2015.[21] R. Thirumalainambi and J. Bardina, “Training data requirement for aneural network to predict aerodynamic coefficients,”
Indep. Compon.Anal. Wavelets, Neural Networks , vol. 5102, p. 92, 2003.
Fig. 10:
Box plot representation of prediction error using neural networks trained with training sets reduced with schemes as shown in Fig 3. The left-mostbox plot with no training data reduction scheme corresponds to error found from all 6 paths found using the full training set. a) Euclidean distance in the x - y plane between predicted and expected Gaussian peak locations for the synthetic moving Gaussian testing data set. The proposed limits on maximumEuclidean distance between predicted and expected Gaussian peaks is shown as the blue line for each box plot. b) RMSE per 2D slice. Fig. 11:
The mean of the errors when predicted using neural networks trained with σ reduction scheme sets is shown as the blue solid line, the 95% confidenceintervals are marked by the black dotted lines. The σ of the testing set Gaussian is shown as the red star. For comparison, the mean of the errors when thefull training set is used is the red dotted line. a) Euclidean distance in the x - y plane between predicted and expected Gaussian peak locations for the syntheticmoving Gaussian testing data set; b) RMSE per 2D slice.[22] E. K. Roonizi and R. Sameni, “Morphological modeling of cardiacsignals based on signal decomposition,” Comput. Biol. Med. , vol. 43,no. 10, pp. 1453–1461, 2013.[23] L. S. Lee et al. , “Department of Electrical and Computer EngineeringPart IV Research Project Project Report Project Number : 10 DeepLearning Neural Nets for Detecting Heart Activity,” Tech. Rep., 2018. Fig. 12:
Top row from left to right shows the isosurface for the full 9-by-12-by-649 matrix for the left ventricular recording, the isosurface of the full G3Drepresentation using 100 components, and the RMSE per 2D slice ( x - y plane) for the fitted signal. Bottom row shows the isosurfaces for the first 3 G3Dcomponents fitted. The isosurfaces are drawn at values of -0.1, -0.05, 0.05, and 0.1 to show the various potentials found within the 3d signal matrices. Fig. 13: a) RMSE values per HSP 2D slice between the neural networkpredicted and expected G3D fit potentials for each temporal sample of theleft ventricular averaged beat recording. b) Time series examples of predictedHSPs for the left ventricular average beat recording at electrode (x = 7, y =11). The G3D fit HSP is the ideal prediction from the neural network (trainedbased on G3D basis functions only). The normalized recorded HSP is alsopresented as a reference. Fig. 14:
3D activation map mesh comparison of the activation maps produced from predicted HSP (Left) and real world recording (Right). The electrodeclosest to the point of activation (left ventricle) is marked as the red star. Supplementary Materials
A. Overall Methodology
An overview of the methodology is shown in Figure S1.
Fig. S1:
Summary of methodology.
B. Body and Heart Surface 2D Projections
We performed 2D geometry projections for the heart surfaceelectrodes and body surface electrodes separately for eachinstance in time. At every instance in time, there are 108HSP and 128 BSP electrode potential values. The heart surfaceelectrodes were projected onto an evenly spaced 9-by-12 meshas shown in Fig S2. Using this heart 2D projection, the LVpaced averaged beat HSP was projected into a 3D array ofsize 9-by-12-by-649. The body mesh from the experimentaldata was used to project the body surface electrodes ontoa cylinder approximation of the body using pointCloud inMATLAB [23], as shown in Fig S3. The cylinder was un-wrapped via a cut along the left anterior descending arteryline identified in the experimental data. Bad leads identifiedvia the experimental data were removed. The 2D unwrappedcylinder surface was re-sampled onto a 16-by-16 mesh vialinear interpolation using griddata in MATLAB.
Fig. S2:
Summary of heart electrode 2D projection, where electrodes in amesh were projected into 9-by-12 matrix positions.
C. G3D Library
A combination of different G3D parameters were used tocreate a G3D library which can describe a wide range of signalshapes and sizes. Here, the G3D basis functions within thelibrary all have an
Amp of 1. The µ x was allowed to varyin 8 even locations between 1 and 9 (9 is the size of the xdimension for signal to be fitted). The µ y was allowed to vary Fig. S3:
Summary of body electrode 2D projection. Electrodes in 3D spacewere first projected onto a cylinder. The cylinder was then unwrapped alongthe left anterior descending artery line. The unwrapped 2D surface was thenre-sampled into 16-by-16 matrix positions. in 8 even locations between 1 and 12 (12 is the size of they dimension for signal to be fitted). The µ t was allowed tovary evenly in 2, 4, 8, 16, 32, 64 splits between 0 and 1 ( t hasunits in heart beat time, varies between 0 and 1). The σ xy wasallowed to vary evenly in 2, 4, 8, 16, 32, 64 splits between 0and 0.5. The σ t was allowed to vary evenly in 2, 4, 8, 16, 32,64 splits between 0 and 1. The final G3D library containedG3D basis functions generated from all given combinations ofthe 6 parameters. D. Synthetic Moving G3D Paths
Figure S5 summarizes the synthetic moving Gaussian pathsalong with their respective contours at 3 points in time foreach of the 6 paths found in Table II. The 2D representationof these paths on the 9-by-12 heart projection surface can beseen in Fig S4.
Fig. S4:
Synthesized moving Gaussian data paths on the HSP surface (9-by-12). The synthesized data have Gaussians with peaks along these indicatedpaths, smoothly translating from the start point to end point in 100 samples.
E. G3D Decomposition Representation
Fig S6 shows that error of the summed basis functionrepresentation of the HSP signal decreases as the number ofG3D’s components fit to the signal increase. This is expectedas each new basis function was only fit to the unexplained partof the signal in an iterative process.Fig S7 shows the fitted G3D component isosurfaces for the9-by-12-by-649 HSP matrix. In the subplots, the number of Fig. S5:
Synthetic moving Gaussian paths (100 samples in total length each) along with contours at 3 points in time. The subplot labels are the path numberfollowed by a) 1st sample, b) 50 samples in, c) last sample. summed components progressively increase. Note how the iso-surface begins to resemble that from the actual recording fromFigure 12 as the number of summed components increase,which indicates better representation of the original recording.Subplot a from Fig S8 shows the similarities in morphologyand scale between G3D reconstructed and recorded HSP timeseries for the electrode at x=7, y=11. Subplot b from Fig S8presents a subset of 5 G3D waves at electrode x=7, y=11which were used to reconstruct the G3D fit time series shownin subplot a.
F. Neural Network Training and Validation Loss
Fig S9 shows the mean squared error loss for both thetraining and validation sets over the 100 epochs. Fig. S6:
RMSE of whole left ventricular average beat recording 3D Arrayper G3D components fit. Fig. S7:
Summed isosurface of the first 10 G3D components. Each subplot presents the partial sum of the G3D representation of the HSP matrix with thegiven components in the title.
Fig. S8: a) Time series comparison between G3D fit and left ventricular averaged beat HSP at electrode position x = 7, y = 11 b) Five largest fitted G3Dcomponents in amplitude at electrode x = 7, y = 11 as time series. c) RMSE for 2D sample slices (along x-y axis) for the 100 G3D component fit of leftventricular averaged beat HSP.