^{1}, Elizabeth M. Cherry

^{2}, Alain Karma

^{3}and Wouter-Jan Rappel

^{4}

### Abstract

We present a novel algorithm for modeling electrical wave propagation in anatomical models of the heart. The algorithm uses a phase-field approach that represents the boundaries between the heart muscle and the surrounding medium as a spatially diffuse interface of finite thickness. The chief advantage of this method is to automatically handle the boundary conditions of the voltage in complex geometries without the need to track the location of these boundaries explicitly. The algorithm is shown to converge accurately in nontrivial test geometries with no-flux (zero normal current) boundary conditions as the width of the diffuse interface becomes small compared to the width of the cardiac action potential wavefront. Moreover, the method is illustrated for anatomically realistic models of isolated rabbit and canine ventricles as well as human atria.

The mechanisms for the generation and maintenance of cardiac arrhythmias, the leading cause of death in the industrialized world, are poorly understood. Computer modeling can furnish valuable information, provided that the models used are sufficiently realistic. Modeling electrical wave propagation in models of cardiac geometries requires not only an accurate description of the electro-physiological properties of the heart but also a precise implementation of the geometrical structure and fiber orientation within the tissue. This paper addresses these geometrical requirements and presents a new algorithm to model wave propagation in anatomical models of the heart. The algorithm is based on the diffuse interface phase-field approach that has been used in a wide range of contexts. Simulation results are presented that quantify the accuracy of the method and illustrate its application to realistic heart geometries ranging from rabbit and canine ventricles to the human atria.

We gratefully acknowledge the hospitality of the Aspen Center for Physics where part of this research was performed. This research was facilitated through an allocation of advanced computing resources by the National Computational Science Alliance, through the support of the National Science Foundation, and the computations were performed in part on the National Science Foundation Terascale Computing System at the Pittsburgh Supercomputing Center. We acknowledge the National Biomedical Computation Resource (National Institutes of Health Grant No. P41 RR-08605). This research was supported in part by National Science Foundation Grants No. PHY99-07949 and No. 0320865 (F.H.F. and E.M.C.), by NIH Grant No. 5 F32 HL073604-02 (E.M.C.), by the NSF-sponsored Center for Theoretical Biological Physics (Grants No. PHY-0216576 and No. 0225630) (W.-J.R.), and by NIH Grant No. HL075515-01 (W.-J.R.). A.K. acknowledges support of NIH Grant No. P50-HL52319.

I. INTRODUCTION

II. METHODS

III. APPLICATION TO SAMPLE GEOMETRIES

A. One-dimensional cable

B. Two-dimensional irregular tissue geometries

IV. APPLICATION TO ANATOMICAL MODELS IN THREE DIMENSIONS

V. CONCLUSION

### Key Topics

- Heart
- 35.0
- Boundary value problems
- 22.0
- Anisotropy
- 12.0
- Finite difference methods
- 11.0
- Wave propagation
- 11.0

## Figures

(Color online). (a) Example of the phase field along a one-dimensional section from the slice of rabbit ventricles shown in (b). The solid line indicates the original values of 1 or 0 assigned to and the open circles represent the final steady-state values of after solving Eq. (3) using . (c) Different solutions for using , 0.05, and corresponding to diamonds, circles and stars, respectively. Note that for the choice of in Eq. (3) , matches , as shown by the long dashed curve when using and the short dashed curve when using . The width of the interface is approximately .

(Color online). (a) Example of the phase field along a one-dimensional section from the slice of rabbit ventricles shown in (b). The solid line indicates the original values of 1 or 0 assigned to and the open circles represent the final steady-state values of after solving Eq. (3) using . (c) Different solutions for using , 0.05, and corresponding to diamonds, circles and stars, respectively. Note that for the choice of in Eq. (3) , matches , as shown by the long dashed curve when using and the short dashed curve when using . The width of the interface is approximately .

(Color online). Fiber orientation information for the anatomical model of rabbit ventricles of Ref. ^{ 20 } . (a) and (b) Slices (a) and (b) from the apex showing the projections of the three-dimensional fibers in the plane. (c) The three components of the fiber orientation vectors along the line in (b) for points inside the original domain (white area, length between 0.57 and 0.9 cm) and points in the enlarged domain (gray areas) whose fiber values originally were undefined. Continuity is present in all components. See text for details of the iterative process used to determine values for previously undefined fiber orientations.

(Color online). Fiber orientation information for the anatomical model of rabbit ventricles of Ref. ^{ 20 } . (a) and (b) Slices (a) and (b) from the apex showing the projections of the three-dimensional fibers in the plane. (c) The three components of the fiber orientation vectors along the line in (b) for points inside the original domain (white area, length between 0.57 and 0.9 cm) and points in the enlarged domain (gray areas) whose fiber values originally were undefined. Continuity is present in all components. See text for details of the iterative process used to determine values for previously undefined fiber orientations.

Representative action potentials as measured in uniform cables from (a) the Nygren *et al.* human atrial model (Ref. ^{ 38 } ), (b) the Fox *et al.* canine ventricular model (Ref. ^{ 39 } ), and (c) the phenomenological ionic model in Ref. ^{ 12 } with parameters to reproduce canine epicardial cells. Note that the Fox *et al.* model has a faster upstroke than the phenomenological model (values of are 270 and , respectively). (d) The upstrokes corresponding to the models shown in (a) (long dashes), (b) (solid), and (c) (dashes). Models with different values of are used to show its effect on the accuracy of the solutions obtained using the phase-field method (see Fig. 5 ).

Representative action potentials as measured in uniform cables from (a) the Nygren *et al.* human atrial model (Ref. ^{ 38 } ), (b) the Fox *et al.* canine ventricular model (Ref. ^{ 39 } ), and (c) the phenomenological ionic model in Ref. ^{ 12 } with parameters to reproduce canine epicardial cells. Note that the Fox *et al.* model has a faster upstroke than the phenomenological model (values of are 270 and , respectively). (d) The upstrokes corresponding to the models shown in (a) (long dashes), (b) (solid), and (c) (dashes). Models with different values of are used to show its effect on the accuracy of the solutions obtained using the phase-field method (see Fig. 5 ).

Membrane potential distribution along a -long cable at different times using the phase-field method with , , and (symbols) and using a standard zero-flux finite-difference code with the same discretization. The initial condition is a brief excitation at the center of the cable at . This produces a symmetric excitation that propagates to the edges. (a) Voltage distribution during depolarization (all voltage values are ). Initial time is , final time is , and the voltage distribution is plotted every . (b) Voltage distribution during repolarization (all voltage values are ). Initial time is , final time is , and the voltage profile is plotted every . (c) Comparison of action potentials at the boundary using finite differences (solid) and the phase-field method (dotted), with the two upstrokes highlighted in the inset. The phenomenological model is used.

Membrane potential distribution along a -long cable at different times using the phase-field method with , , and (symbols) and using a standard zero-flux finite-difference code with the same discretization. The initial condition is a brief excitation at the center of the cable at . This produces a symmetric excitation that propagates to the edges. (a) Voltage distribution during depolarization (all voltage values are ). Initial time is , final time is , and the voltage distribution is plotted every . (b) Voltage distribution during repolarization (all voltage values are ). Initial time is , final time is , and the voltage profile is plotted every . (c) Comparison of action potentials at the boundary using finite differences (solid) and the phase-field method (dotted), with the two upstrokes highlighted in the inset. The phenomenological model is used.

(a) Relative error in the maximum upstroke velocity as a function of the phase field width for one-dimensional cable simulations as shown in Fig. 4 for the Fox *et al.* (short dashes), Nygren *et al.* (long dashes), and phenomenological (solid) models. Note that for the values of used throughout this paper, the relative error is less than 10%. (b) Cumulative error in action potential for the same cases. The cumulative error is obtained by computing the absolute error in voltage over the time course of one action potential and then computing the ratio of the area under that curve to the area under the curve of the action potential calculated with the finite difference code. (c) Relative error in action potential duration for the same cases. In all cases, and . Slight increases in error can be observed for the smallest values of because for these values the interface is very steep and is not adequately resolved by the fixed used.

(a) Relative error in the maximum upstroke velocity as a function of the phase field width for one-dimensional cable simulations as shown in Fig. 4 for the Fox *et al.* (short dashes), Nygren *et al.* (long dashes), and phenomenological (solid) models. Note that for the values of used throughout this paper, the relative error is less than 10%. (b) Cumulative error in action potential for the same cases. The cumulative error is obtained by computing the absolute error in voltage over the time course of one action potential and then computing the ratio of the area under that curve to the area under the curve of the action potential calculated with the finite difference code. (c) Relative error in action potential duration for the same cases. In all cases, and . Slight increases in error can be observed for the smallest values of because for these values the interface is very steep and is not adequately resolved by the fixed used.

(Color online). Propagation of a point stimulus applied off-center in a quarter-annulus. (a)–(d) Electrical excitation (orange) propagating into quiescent tissue (black) at times , 20, 35, and , respectively. (e) Comparison of wave front contours at intervals for the solution using the phase-field method in Cartesian coordinates (red, symbols) and the reference solution using polar coordinates (black, solid). Grid spacings are and or 0.45°, with .

(Color online). Propagation of a point stimulus applied off-center in a quarter-annulus. (a)–(d) Electrical excitation (orange) propagating into quiescent tissue (black) at times , 20, 35, and , respectively. (e) Comparison of wave front contours at intervals for the solution using the phase-field method in Cartesian coordinates (red, symbols) and the reference solution using polar coordinates (black, solid). Grid spacings are and or 0.45°, with .

(Color online). Propagation of a point stimulus applied (a) to the lower right corner of a quarter-annulus, (b) to the lower right corner of a quarter-annulus with a hole, and (c) near the right edge of an anisotropic square domain with a hole. Wave front contours are shown at intervals for the solution using the phase-field method (red, symbols) and for the reference solution using finite differences (black, solid). Reference solutions are obtained using polar coordinates in (a) and (b) and using standard finite differences in (c). Note that the contours are normal to all the boundaries for both solutions in (a) and (b). The ratio of diffusion constants parallel and perpendicular to the fibers in (c) is and the fiber angle is . Grid spacings are and or 0.45° with for (a) and (b), while for (c) the grid spacing is with and .

(Color online). Propagation of a point stimulus applied (a) to the lower right corner of a quarter-annulus, (b) to the lower right corner of a quarter-annulus with a hole, and (c) near the right edge of an anisotropic square domain with a hole. Wave front contours are shown at intervals for the solution using the phase-field method (red, symbols) and for the reference solution using finite differences (black, solid). Reference solutions are obtained using polar coordinates in (a) and (b) and using standard finite differences in (c). Note that the contours are normal to all the boundaries for both solutions in (a) and (b). The ratio of diffusion constants parallel and perpendicular to the fibers in (c) is and the fiber angle is . Grid spacings are and or 0.45° with for (a) and (b), while for (c) the grid spacing is with and .

Maximum relative error in propagation velocity in the domain shown in Fig. 7(c) as a function of the , the width of the phase field. For values of below , the error is less than 10%. The inclusion of anisotropy increases the error by 2%–4% compared to an isotropic simulation for nearly all values of tested.

Maximum relative error in propagation velocity in the domain shown in Fig. 7(c) as a function of the , the width of the phase field. For values of below , the error is less than 10%. The inclusion of anisotropy increases the error by 2%–4% compared to an isotropic simulation for nearly all values of tested.

(Color). Example simulations using the phase-field method in complex cardiac geometries. (a) Single scroll wave in the rabbit ventricular model. The left and right images show posterior and anterior views, respectively. (b) Slabs of the rabbit ventricles during scroll wave propagation (posterior view). The slabs are perpendicular to the apex-base axis and proceed toward the apex. (c) Propagation of an electrical wave in the canine ventricles after a stimulus along a simulated Purkinje network. The left image shows an anterior view of the ventricles with a small portion cut out to allow the endocardium to be seen. The cut-out view on the right shows the anterior endocardium. (d) Two views of a spiral wave in the anatomical model of the human atria. The Nygren *et al.* model of human atrial cells is used. Electrical potential is color-coded with red corresponding to strongly depolarized tissue and blue corresponding to repolarized tissue. In all cases, grid spacing is , and the phase-field control parameter is . Time steps are for (a) and (b), for (c), and for (d).

(Color). Example simulations using the phase-field method in complex cardiac geometries. (a) Single scroll wave in the rabbit ventricular model. The left and right images show posterior and anterior views, respectively. (b) Slabs of the rabbit ventricles during scroll wave propagation (posterior view). The slabs are perpendicular to the apex-base axis and proceed toward the apex. (c) Propagation of an electrical wave in the canine ventricles after a stimulus along a simulated Purkinje network. The left image shows an anterior view of the ventricles with a small portion cut out to allow the endocardium to be seen. The cut-out view on the right shows the anterior endocardium. (d) Two views of a spiral wave in the anatomical model of the human atria. The Nygren *et al.* model of human atrial cells is used. Electrical potential is color-coded with red corresponding to strongly depolarized tissue and blue corresponding to repolarized tissue. In all cases, grid spacing is , and the phase-field control parameter is . Time steps are for (a) and (b), for (c), and for (d).

Article metrics loading...

Full text loading...

Commenting has been disabled for this content