A Probabilistic Model of Cardiac Physiology and Electrocardiograms

An electrocardiogram (EKG) is a common, non-invasive test that measures the electrical activity of a patient's heart. EKGs contain useful diagnostic information about patient health that may be absent from other electronic health record (EHR) data. As multi-dimensional waveforms, they could be modeled using generic machine learning tools, such as a linear factor model or a variational autoencoder. We take a different approach:~we specify a model that directly represents the underlying electrophysiology of the heart and the EKG measurement process. We apply our model to two datasets, including a sample of emergency department EKG reports with missing data. We show that our model can more accurately reconstruct missing data (measured by test reconstruction error) than a standard baseline when there is significant missing data. More broadly, this physiological representation of heart function may be useful in a variety of settings, including prediction, causal analysis, and discovery.


Introduction
Heart disease annually claims the lives of over 600,000 people in the United States and over 17 million people worldwide [1,2]. Early detection of cardiovascular disease is critical, and accurate characterization of patient risk can improve care. Consequently, a commonly used diagnostic tool is the electrocardiogram (EKG), which uses electrodes placed on the patient's body to measure the heart's electrical activity. EKGs are widely administered to diagnose a variety of cardiac abnormalities, including ischemia, heart attack, and arrhythmias. EKGs (and other waveform and image data) likely carry additional information about the overall health of a patient, beyond what is already in electronic health records (e.g. past diagnoses, medications, and lab tests). As we begin to use EHR data to inform health care decisions, EKG features can be used alongside patient EHR data to construct more accurate predictors, similar to unstructured physician notes [3] or the dynamics of ICU vital signs [4].
However, unlike discrete records, the EKG is a high-dimensional, dynamic object that is difficult to summarize. While an EKG could be summarized with standard ML tools, a generic approach would ignore the underlying physiology that drives EKG observations: the waveforms are a measurement of an electrical signal originating in the heart and propagating through the body.
In this work, we pursue a physiological representation of an EKG by modeling the location, orientation, and strength of the electric potential induced by cardiac activity. Such a model has many benefits. First, because our model represents physical electrode locations, our approach can flexibly and coherently incorporate data from systems with more (or fewer) electrodes than standard 12-lead EKGs,  Figure 1: (a) Graphical depiction of the latent moving dipole model. The instantaneous electrical activity of the heart is summarized by a latent moving dipole, with location s t and moment p t . This object induces a potential field (depicted by dotted lines), which is measured by the electrodes located at r e according to Equation 2 -we depict the right arm (e = ra), left leg (e = ll) and two precordial electrodes, e = v1 and e = v2. The electrode measurements are linearly combined to produce the EKG lead observations (mentioned in Section 2 and detailed in Appendix A). (b) Prior distribution over the 3-d spatial location of the precordial electrodes, v1 -v6.
including higher-dimensional (e.g.256-electrode body surface EKGs [5]) and lower-dimensional (e.g. 5-electrode Holter monitors) electrocardiographic signals. Second, we show that our model is robust to large amounts of missing observations. Finally, it produces a physically interpretable representation of the EKG that can be used in a variety of applications, such as prediction or causal analysis -for example, measuring the effect of a therapy on some aspect of cardiac function. We evaluate our model with a lead reconstruction task -given an EKG where data from some leads are missing, how accurately can we reconstruct the held out data? We find that our model can more accurately reconstruct missing EKG leads than non-physical factor models when we have a lot of missing data (which is common in our clinical emergency department setting). We conclude with a discussion of statistical inference challenges and model extensions.

Data Description and Model
EKG Background An EKG lead measures the electric potential difference between two (or more) electrodes. An electrode is a conductive pad placed on a part of the body (e.g. left leg or chest) that measures instantaneous electric potential. A lead records the difference between a pair (or a specific combination) of electrode readings (typically in millivolts); cardiologists visually inspect the lead outputs for abnormalities. The standard twelve leads are commonly referred to as L = {I, II, III, aV R, aV L, aV F, V 1, . . . , V 6}. 2 An EKG tracing of length T consists of |L| = 12-dimensional observations at T time steps, We model lead EKG observations by describing the physical location, orientation, and strength of the electric potential induced by cardiac activity. We introduce model components that (i) describe the underlying state of electrical activity; (ii) map latent electrical activity to electrode potential values; and (iii) map electrode values to lead observations. To address point (i) we model electrical activity with a single moving dipole, following [6]. We address point (ii) by modeling the 3-d location of electrodes and using a simplified torso conductance model. Point (iii) is straightforward and can be modeled by a fixed linear map -we describe this step in detail in Appendix A.
Moving Dipole Generative Model Each observation x t is generated by the electric potential induced by a single latent dipole with state z t (s t , p t ), where s t ∈ R 3 is the spatial location of the dipole and p t ∈ R 3 is its moment. A dipole is an idealized point whose location and moment induce a spatial electrostatic potential that is measured by each electrode. At time t, the latent dipole induces a potential at each of the |E| electrodes that are spatially located at points r e where e = 1, . . . , |E| and r e ∈ R 3 . Assuming a uniform conductance torso model, a single dipole with state z t = (s t , p t ) will have a potential measured at electrode e, located at r e , given bỹ where κ is the (constant) electrical conductivity of the torso [6] (which we set to .2, following [7]).
To get from electrode to lead, there is a fixed linear function which produces the twelve standard leads, namely Priors To complete the model, we specify a prior distribution over the latent dipole state, s t and p t , and the location of each electrode r e . We put zero-centered spherical Gaussian priors over s t and p t and weakly informative priors over electrode locations, p e (r), described in Appendix B.1. The data generating process is x t ∼ N (Ox t , Iσ 2 ) wherex t,e = g(s t , p t , r e , κ) In this current iteration we ignore temporal correlation in s t and p t , leaving this for future work.
Inference The above generative model will admit tractable (though approximate) variational inference [8]. As a first step, however, we use MAP inference to show the utility of this modeling approach, solving the optimization problem ln p e (r e ) .
To solve this optimization problem, we use L-BFGS [9] with gradients computed by automatic differentiation [10]. See Section D for further discussion of inference.

Experiments
We assess model performance by measuring held-out lead reconstruction error. For each EKG, we hold out portions of the waveform during training (depicted in Figure 2). We then use the model to impute the held-out values and compute the root mean squared error (RMSE). A model that better captures the underlying phenomenon will more accurately impute missing lead values.
We apply our model to a small sample of 10-second EKG segments from two datasets: • ptb: the PTB Diagnostic ECG Database [11], a database of 549 full twelve-lead EKG records. • ed: a subset of 997 EKG recordings from emergency department patients from a large metropolitan hospital, each exhibiting "normal sinus rhythm." These EKGs are derived from a traditional report that saves the full time series for leads II, V1, and V5, and short (approximately 2.5 second) segments of the other 9 leads. This missingness can frustrate analysis of the full EKG object, which we address with both our model and baseline approaches. Figure 2b shows a typical ed record.
Baseline linear factor model We compare the latent dipole model to probabilistic PCA (with K = 3 and K = 6 latent dimensions), which describes lead observations as Here F are the K latent factors, and z t is the K-dimensional latent state that describes observation x t (analogous to our dipole location and moment). We use an implementation of PCA designed to handle missing entries [12].
Empirical Results For each dataset, we fit the dipole model and two PCA models (K = 3 and K = 6) to each EKG. We compare the median RMSE over all patients in each dataset. Figure 3 depicts the median RMSE for each model for the ed and ptb. We observe that the dipole model more reliably reconstructs out of sample lead observations in the ed setting, where the short leads are missing most entries (see Figure 2). We observe similar performance for all three models in the ptb dataset. The dipole model is constrained to produce a spatially consistent reconstruction of EKG leads, which may help in the absence of many lead observations.
We produce a full imputation of an ed EKG in Figure 5, which shows that the dipole model can have more stable imputations. To contrast this, the ptb EKG imputation in Figure 6 shows that both models have an easier time reconstructing out of sample leads when most of the tracing is present. Figure 4a shows the dipole model reconstruction on an ed EKG. Figure 4b compares the residuals for the dipole model and two PCA models on the same EKG. We see that the dipole model accurately reconstructs leads II and aVF (almost no structure in the residual), but has a difficult time with the precordial leads V1 and V5, which may be due to inaccurate reconstruction of electrode location. The PCA models have structured residuals for all leads around regions of high voltage.

Discussion
We developed a model of EKG observations with a physiological inductive bias. We compared our model to non-physical latent factor models, showing its potential to more accurately characterize patient physiology. There are many directions of future work. We will model temporal dynamics in the latent variables z t , incorporating local temporal smoothness and quasi-periodicity into the prior. Such dynamics may better characterize arrhythmias or other abnormalities not present in normal sinus rhythm EKGs. The non-conjugacy of the dipole-electrode model complicates inference; a focus will be on efficient approximate inference methods for this model, following [13] and [14]. We also aim to use model output in downstream statistical analyses -exploring the structure of patients with adverse cardiac events, building risk predictors alongside historical patient medical record that can aid physician decision-making, and measuring the causal effect of therapies on aspects of cardiac function.

A Electrode to Lead transformation
Electrocardiogram data are presented as leads, which are differences in potential values between electrode measurements (or averages of electrode measurements). This is a somewhat straightforward transformation that we will account for in the likelihood. The output potential difference for each lead are simple functions of the electrode readings: We can express these relationships as a simple linear transformation of the electrode valuesx into lead values x where the transformation matrix O has a fixed block structure, where the blocks are defined

B.1 Electrode location prior
Fitting data with the above appearance model will require inferring the electrode locations R. We can use the general guidelines of where to attach electrodes to place a prior on their relative location.
We use a simple elliptical model of the thorax to set prior values on the location of V1 -V6. We set the chest width to be 25 centimeters and the major-minor axis ratio to be 2.75. The V1 -V6 locations are at even angular locations from 260 degrees to 360 degrees (all the way to the left side of the chest). Figure 1b depicts the prior (with 2 standard deviations for the precordial (chest) electrodes (v1 -v6).
The location of the right arm, left arm, and left leg electrodes will be given diffuse priors flanking the modeled torso depicted in Figure 1b.

B.2 Conductance Parameter Priors
Average thoracic tissue impedance can be approximated as resistive with conductivity values around .2 S/m [15].  Figure 4b also depicts the residuals for the dipole model and both PCA models. Figure 5 displays the full missing data imputation for an emergency department patient. Similarly, Figure 6 shows the missing data imputation for the ptb example, which has no missing data.

D Inference
We currently fit the model with MAP inference over both the latent z t and global r e parameters. This joint maximization, however, may lead to overfitting -we would prefer to maximize the following marginal likelihood objective {r} = arg max r ln p(x|r)p(r) = arg max r ln p(x|z, r)p(z)p(r)dz .
As an example, this is the objective that PCA essentially fits -the local weights are integrated out (in the E-step) and the global factors are optimized (in the M-step).
Our dipole is challenging, however, because the latent factors z and global variables r enter into the likelihood term, p(x|z, r) with a complex function, given by Equation 5.
The posterior distribution over z given the observations x and electrode locations r does not admit a closed form p(z t |x t , r) ∝ p(x t |z, r)p(z t ) where we define g to be the |E|-length electrode observations g(z t , r) = g(s t , p t , r 1 , κ), . . . , g(s t , p t , r |E| , κ) , recalling that the latent state z t (s t , p t ) is divided into a location and moment component. Because g(s t , p t , r e , κ) is nonlinear in s t , the posterior distribution over s t becomes intractable. We do note that g(·) is linear in p t when you fix s t , r e , and κ. If we place a Gaussian prior over p t , this conditional posterior should be tractable. We plan to use this tractable conditional structure in future work to scale inference in this model.