International audienceSophisticated models for the electrocardiographic inverse problem are available, but their reliance on imaging data and large numbers of electrodes limit their use. Simple models such as the equivalent dipole model (EDM) therefore remain relevant. We developed a probabilistic approach to the equivalent unbounded uniform single dipole problem and developed a natural extension to the bounded nonuniform case that relies on a patientspecific statistical inference of the propagation mechanism between the location of the dipole and the electrode locations. The two models were tested on data simulated with a detailed heart-torso model with four different activation sequences and three different sets of tissue characteristics. We observed a throughout enhancement of the ability to reconstruct the ECG of the patient-specific model when compared to the uniform unbounded dipole model

Arrieula, Andony

Cardoso, Gabriel

Dubois, Rémi

Haïssaguerre, Michel

Moulines, Eric

Potse, Mark

Robin, Geneviève

Robin, Geneviève

INRIA a CCSD electronic archive server

HAL Id: hal-03936912https://hal.inria.fr/hal-03936912Submitted on 12 Jan 2023HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoirespublics ou privés.A Patient-Specific Equivalent Dipole ModelGabriel Cardoso, Geneviève Robin, Andony Arrieula, Mark Potse, MichelHaïssaguerre, Eric Moulines, Rémi DuboisTo cite this version:Gabriel Cardoso, Geneviève Robin, Andony Arrieula, Mark Potse, Michel Haïssaguerre, et al.. APatient-Specific Equivalent Dipole Model. Computing in Cardiology, Sep 2022, Tampere, Finland.￿hal-03936912￿A Patient-Specific Equivalent Dipole ModelGabriel Cardoso1,4,6, Geneviève Robin2, Andony Arrieula3,4,5, Mark Potse5,4,3,Michel Haı̈ssaguerre4,6,7, Eric Moulines1, Rémi Dubois4,61 CMAP Ecole Polytechnique, Palaiseau, France2 Université d’Evry, Evry, France3 CARMEN Research Team, Inria Bordeaux – Sud-Ouest, Talence, France4 IHU Liryc, fondation Bordeaux Université, Pessac, France5 Univ Bordeaux, IMB, UMR 5251, Talence, France6 Univ Bordeaux, CRCTB U4045, INSERM, Bordeaux, France7 Bordeaux University Hospital (CHU), Electrophysiology and Ablation Unit, Pessac, FranceAbstractSophisticated models for the electrocardiographic in-verse problem are available, but their reliance on imag-ing data and large numbers of electrodes limit their use.Simple models such as the equivalent dipole model (EDM)therefore remain relevant. We developed a probabilis-tic approach to the equivalent unbounded uniform sin-gle dipole problem and developed a natural extension tothe bounded nonuniform case that relies on a patient-specific statistical inference of the propagation mechanismbetween the location of the dipole and the electrode loca-tions. The two models were tested on data simulated witha detailed heart-torso model with four different activationsequences and three different sets of tissue characteristics.We observed a throughout enhancement of the ability to re-construct the ECG of the patient-specific model when com-pared to the uniform unbounded dipole model.1. IntroductionModeling the electrical activity of the heart as a sin-gle, time-varying electric current dipole dates back toEinthoven et al.[1]. Einthoven’s triangle [2, volume 1section 10.3.6] assumes a single fixed dipole in an un-bounded uniform body. Information from leads I, II, IIIof the standard ECG allows measurement of the individualcomponents of such a dipole. However, the assumption ofan unconfined uniform body significantly limits the inter-pretability of the resulting dipole parameters. A currentdipole with magnitude pt ∈ R3 and position xp ∈ R3 in aconductor with an arbitrary bounded volume, generates apotentialyit = c(xi,xp) · pt, c : R3 × R3 → R3 (1)at position xi ∈ R3, where i ∈ {1, ...,K} is the index ofthe measurements [2, Chapter 10]. This is the main groundfor the development of vectorcardiography and Frank’slead system, whose goal is to measure each of the com-ponents of pt in an orthogonal basis in the equivalent elec-trical space. Frank [3] established a method for estimatingan orthogonal lead systems, by inducing a current dipolein three directions in an artificial torso and measuring theresulting potential at the surface. This led to the creationof a 3d lead system called the Frank’s lead system definedby Fx, Fy, Fz leads that represents an orthogonal basis ofthe equivalent electrical space.The vectorcardiogram obtained using this lead systemhas been successfully applied to a number of problems [4].An important limitation of this approach is that it is notpatient-specific, being experimentally derived from an arti-ficial torso. Therefore we developed a new patient-specificEDM by using Bayesian statistics to infer the function yitfor each patient. We tested the method on simulated ECGsproduced with a detailed human heart and torso model,with different activation origins and different conductionproperties.2. Equivalent Unbounded dipole modelOur goal is to estimate the posterior distribution of theparameters given the observed data.The following notation will be used throughout therest of this paper: If ait is a scalar variable indexed byi ∈ {1, ...,K} and t ∈ {1, ..., T}, we write A =(a11, a21, ...aKT ). Assuming propagation through an un-bounded homogeneous medium, the torso potentials areyit|pt∼N(cHOM(xi) · pt, σ2i,t), cHOM(xi) = xi/‖xi‖3(2)where pt represents the dipole vector of a current dipoleplaced at 0.Under this model, the log-likelihood of the observationsis`(Y |P ) ∝ −T∑t=1K∑i=1(yit − pt · cHOM(xi))2σ2i,t(3)If the measurements (yit)i=1,...,Kt=1,...,T start at the beginningof the QRS complex of the ECG then it is reasonable toexpect pt to be small at the beginning, since it correspondsto the onset of the activation of the ventricles. There is noa priori reason to expect that the components of the dipoleare positive or negative, so a Gaussian distribution centeredon 0 is a natural candidate for the prior distribution p0.Most of the common sampling rates in the ECG domain(sampling rate > 250 Hz) are higher than the frequency ofthe dipole curve pt. Therefore, we expect a small deviationbetween pt and pt+1, i.e., we expect ‖pt+1 − pt‖ to beclose to 0. This leads us to define the prior distribution aspt+1|pt = N (pt, σ2pI), p0 = N (0, I). (4)The log probability of the posterior can be derivedthrough Bayes’ equation and it is negative if the observa-tion Y is given, up to an arbitrary normalizing constant,by:T∑t=1[K∑i=1(yit − pt · cHOM(xi))2σ2i,t+‖pt − pt−1‖2σ2p]+ ‖p0‖2(5)This setting is the well known case of Bayesian linear re-gression [5] with Gaussian measurements noise and statedisturbances. It can be shown thatP |Y ∼ N (µP = ΣPMTΣ−1m Y ,ΣP = (Rreg +MTΣ−1m M)−1)(6)when ΣP is invertible, which is the case in our settings.3. Our modelTo make the EDL model patient-specific we considerc as a random function that must be estimated for eachpatient. To develop a Bayesian approach, we need to definea prior over a function space and a likelihood. The mainidea is to use the model developed in section 2 to definea Gaussian process as a prior distribution over a functionspace. For f : R3 → R3 :yit|f ,pt,xi∼N(cINH(xi) · pt, σ2i,t)cINH(xi;f) = {xi + f(xi)}/‖xi‖3(7)The prior over the function f will be defined as, follow-ing Bishop [5] section 6.4:f ∼ GP(0, kθ(x,x′)) (8)Variable Distribution Interpretationlk LogNormal(µ∗ = l0(X), σ∗2 = 10) Expect values to be in the in-terval [ l0(X)10 , 10l0(X)]νk LogNormal(µ∗ = σ̂2kernel, σ∗2 = 1.01) Expect values to in the inter-val [ σ̂2kernel1.01 , 1.01σ̂2kernel]p0 N (0, I) Same as in Section 2pt N (0, σ2pσ̂2kernelI) Same as in Section 2 but wefound numerically that it wasbest to scale it by variance asthe Gaussian process kernelsTable 1. Description of the priors used in the numericalexperiments.The prior and posterior are thus determined by the choiceof the kernel function k. We use a component-wise inde-pendent scaled radial basis function (RBF) kernel definedaskν,l(x, x′) = νexp(−(x− x′)2/2l2)kθ(x,x′) = diag(kνi,li(x[i],x′[i]), i = 1, 2, 3).(9)The parameter ν determines the variance scale of the co-variance matrix and in our case controls the magnitude off , i.e., how far we are from the uniform dipole model.The parameter l indicates the rate at which the covariancebetween two points decreases. A small l would result inmore variation being allowed for f between two neighbor-ing points x, x′. In what follows, θ = (ν1, l1, ν2, l2, ν3, l3).Since the model errors are Gaussian, we end up with anexact posterior Gaussian over a functional space [5, section6]. The likelihood of model (7) is:P(Y |θ,P ) = Ef(X)|θ [P(Y |f , θ,P )] (10)From (8), f(X)|θ is distributed according to a mul-tivariate Gaussian distribution and the same is true forY |f , θ,P . Using formula (2.115) from Bishop [5] thelikelihood has the distributionY |θ,P ∼ N (0,Σ(θ,P ,X)) (11)where we refer to appendix B for the definition ofΣ(θ,P ,X). We proceed in the same manner as in Sec-tion 2 by defining a prior over the parameters in order tocompute the posterior distribution. In Table 3 we use themultiplicative mean and multiplicative standard deviationparametrization of the log normal distribution. The typicalgeometry scale l0(X), σ̂2kernel as well as the procedure toestimate σp are defined in the appendix.In this setting, our goal is to estimate the parameter θ̂ bymaximizing the posterior log-likelihood:θ̂ = arg maxθlogP(Y |θ,P ) + logP(θ) + logP(P ) (12)4. Numerical EvaluationTo solve (12), we used the Adam gradient descent al-gorithm [6]. For computation efficiency, a batch of elec-trodes was sampled for each optimization step. Other pa-rameters used in the optimization are given in appendix C.We initialized the dipole parameter P using µP definedin (6). The metric used to evaluate goodness of fit wasthe coefficient of determination, known as R2, defined asR2(y, z) = 1−∑Tt=1(yt−zt)2∑Tt=1(yt−ȳ)2where ȳ = 1T∑Tt=1 yt.We used a detailed anatomical model of the heart andtorso created from computed tomography data. Simula-tions were performded on a bi-ventricular mesh with a uni-form edge length of 200 µm. A human ventricular ionicmodel [7] was used with a monodomain reaction-diffusionequation. The potentials at 252 electrode locations on thetorso were computed using lead fields [8]. Three anatomi-cal configurations were created: one with a full healthy tis-sue, and two with low-conductivity zones with radii of 12and 18 mm. In these zones the conductivity was 10 timessmaller than in the healthy tissue. For each configuration,3 paced beats from different origins (A, B, C) and one nor-mal ventricular beat were simulated. The pacing locationsand the low-conductivity zones are shown in Fig. 3. In to-tal, 12 sets of signals were simulated. For computationalefficiency, the signals were sampled at 500 Hz.We placed the dipole at the barycenter of the three stim-ulation points A, B, and C. We use the following nota-tion for the figures: INH-DIP corresponds to the solutionof(12), and HOM-DIP corresponds to the solution of (6).We compared the median per patient R2 score throughthe different beat origins and different sizes of the low-conduction zones in figure 1, where we see that the INH-DIP model performed considerably better than the HOM-DIP model. The corresponding reconstructed electrogramsare shown in figure 2.Figure 1. Comparison between the inter-patient mean ofthe per simulation median of the R2 scores for each lead.We can see an enhancement of the R2 scores throughout.5. DiscussionThe model developed in section 2 is simplistic and sev-eral approaches have been developed to reduce the modelerror, e.g. by Frank [3]. However, to the best of our knowl-edge, none of them is based on probabilistic modeling ofthe function c in (1). We have presented a probabilisticframework for solving the EDM problem that allows us toFigure 2. The means of the Likelihoods Y |P , θ,X foreach model as well as the confidence intervals (shaded)corresponding to µ± 3σ in two sample ECG leads.Figure 3. Mesh of the modeled heart used for simula-tions, with the location of the pacing sites (A, B, C) andthe low-conduction zone with a radius of 12 mm in blueand a radius of 18 mm in red.derive a natural extension to a patient-specific model thatoutperforms the equivalent unbounded dipole model in allsituations tested. There are several ways to further developand validate this model. One of them is the development ofthe kernel function, for example, to handle the correlationbetween the different components of the function f . Vali-dation with clinical data is also a major challenge, as is theability of the model to handle even more uncertain data (noknowledge of heart placement, varying geometries).References[1] Einthoven W, Fahr G, de Waart A. On the direc-tion and manifest size of the variations of potential inthe human heart and on the influence of the position ofthe heart on the form of the electrocardiogram. Amer-ican Heart Journal 1950;40(2):163–211. ISSN 0002-8703. URL https://www.sciencedirect.com/science/article/pii/0002870350901657.[2] Comprehensive electrocardiology, 2010. URL http://dx.doi.org/10.1007/978-1-84882-046-3.[3] Frank E. An accurate, clinically practical system for spatialvectorcardiography, May 1956. URL http://dx.doi.org/10.1161/01.CIR.13.5.737.[4] Armoundas A, Feldman A, Mukkamala R, Cohen R. A sin-gle equivalent moving dipole model: An efficient approachfor localizing sites of origin of ventricular electrical activa-tion. Annals of biomedical engineering 06 2003;31:564–76.[5] M. C. Pattern Recognition and Machine Learning. Springer,2016. ISBN 9780387310732.[6] Kingma DP, Ba J. Adam: A method for stochastic optimiza-tion, 2014. URL https://arxiv.org/abs/1412.6980.[7] ten Tusscher KHWJ, Noble D, Noble PJ, Panfilov AV.A model for human ventricular tissue. American Jour-nal of Physiology Heart and Circulatory Physiology 2004;286(4):H1573–H1589. URL https://doi.org/10.1152/ajpheart.00794.2003. PMID: 14656705.[8] Potse M. Scalable and accurate ECG simulation for reaction-diffusion models of the human heart. Front Physiol 2018;9:370.AppendixA. Matrix formulation equivalent uni-form single dipoleY = N (M(X)P ,Σm), M(X)i,j = xji/‖xi‖3M = blockdiag(M, · · · ,M︸ ︷︷ ︸T times) Σtm = diag(σ21,t, ..., σ2K,t)Σm = blockdiag(Σ1m, · · · ,ΣTm) Σp =σ20 0 0 0 . . . 00 σ20 0 0 . . . 00 0 σ20 0 . . . 00 0 0 σ2 . . . 0. . . . . . . . . . . . . . . 00 0 0 0 . . . σ2Ri,j =1, if i = j−1, if i = (j − 3); i ≥ 30, otherwise, Rreg = RTΣ−1p R(13)B. Matrix formulation Gaussian processCk(P ) =pk0IK,K· · ·pkT IK,K ,diag(‖x1‖3, ..., ‖xK‖3)−1Σ(θ,P , X) = Σm +3∑k=1Ck(P )Kνk,lk(X,X)CTk (P )(14)C. Optimization GPFor the ensuring l > 0 and ν2 > 0 of (12) foreach kernel we optimize the log of those quantities. Thetypical scale length for a given geometry X is l0 =∑Ki=1 min{‖xi − xj‖|j 6= i}/K The other parametersare specified below:Variable Value Variable ValueBatch size 32 Learning rate dipole 0.5Learning rate logl 0.01 Learning rate logν 0.01Number of iterations 10000 β (from Adam) (0.9, 0.999)ε (from Adam) 10−8 σkernel 0.1Tr(Σm)Tr(S)P̄σp4‖Y maxt‖T‖M(X)(1,1,1)‖where Y maxt is the K-dimensional vector having eachcomponent being the maximum over time of yit.AcknowledgementsThis work was supported by the Région Nouvelle-Aquitaine, grant nr. 2017 – 1R50109 – 00013434; the Eu-ropean Research Council (H2020 grant agreement number715093, ECSTATIC); and the French National ResearchAgency, grant references ANR-10-IAHU04-LIRYC andANR-11-EQPX-0030. This work was granted access tothe HPC resources of CINES under the allocation 2021-A0110307379 made by GENCI.

A Patient-Specific Equivalent Dipole Model

https://hal.inria.fr/hal-03936912v1/file/paper.pdf

A Patient-Specific Equivalent Dipole Model

Abstract

Similar works

Full text

Available Versions

INRIA a CCSD electronic archive server