Physics of the Heart: From the macroscopic to the microscopic

Physics of the Heart: From the macroscopic to the microscopic Xianfeng Song Advisor: SimaSetayeshgar April 17, 2007

Outline • Part I: Transport Through the Myocardium of Pharmocokinetic Agents Placed in the Pericardial Sac: Insights From Physical Modeling • Part II: Electrical Wave Propagation in a Minimally Realistic Fiber Architecture Model of the Left Ventricle: Dynamics of Phase Singularities • Part III: Calcium Dynamics in the Myocyte

Xianfeng Song, Department of Physics, Indiana University Keith L. March, IUPUI Medical School Sima Setayeshgar, Department of Physics, Indiana University Part I:Transport Through the Myocardium of Pharmocokinetic Agents Placed in the Pericardial Sac: Insights From Physical Modeling

Motivation: Diffusion in Biological Processes Diffusion is the dominant transport mechanism in biology, operative on many scales: Intracellular[1] • The rate of protein diffusion in the cytoplasm constrains a variety of cellular functions and limit the rates and accuracy of biochemical signaling in vivo. Multicellular [2] • Diffusion plays an important role during the early embryonic pattern formation in establishing and constraining accuracy of morphogen prepatterns. Tissue-level [3] • Diffusion controls delivery of glucose and oxygen from the vascular system to tissue cells and also governs movement of signaling molecules between cells. Need for careful characterization of diffusion constants governing various biophysical processes. • [1] Elowitz, M. B., M. G. Surette, et al. (1999). J. Bact. 181(1): 197-203. • [2] Gregor, T., W. Bialek, R. de Ruyter van Steveninck, et al. (2005). PNAS 102(51). • [3] Nicholson, C. (2001), Rep. Prog. Phys. 64, 815-884.

Background: Pericardial Delivery • Thepericardial sacis a fluid-filled self-contained space surrounding the heart. As such, it can be potentially used therapeutically as a “drug reservoir.” • Delivery of anti-arrhythmic, gene therapeutic agents to • Coronary vasculature • Myocardium via diffusion. • Recentexperimental feasibilityof pericardial access [1], [2] Vperi (human) =10ml – 50ml • [1] Verrier VL, et al., “Transatrial access to the normal pericardial space: a novel approach for diagnostic sampling, • pericardiocentesis and therapeutic interventions,” Circulation (1998) 98:2331-2333. • [2] Stoll HP, et al., “Pharmacokinetic and consistency of pericardial delivery directed to coronary arteries: direct comparison • with endoluminal delivery,” Clin Cardiol (1999) 22(Suppl-I): I-10-I-16.

Part 1: Outline • Experiments • Mathematical modeling • Comparison with data • Conclusions

Experiments • Experimental subjects: juvenile farm pigs • Radiotracer method to determine the spatial concentration profile from gamma radiation rate, using radio-iodinated test agents • Insulin-like Growth Factor (125I-IGF, MW: 7734 Da) • Basic Fibroblast Growth Factor (125I-bFGF, MW: 18000 Da) • Initial concentration delivered to the pericardial sac at t=0 • 200 or 2000 mg in 10 ml of injectate • Harvesting at t=1h or 24h after delivery

ExperimentalProcedure • At t = T (1h or 24h), sac fluid is distilled: CP(T) • Tissue strips are submerged in liquid nitrogen to fix concentration. • Cylindrical transmyocardial specimens are sectioned into slices: CiT(x,T) x denotes i CT(x,T) = Si CiT(x,T) x: depth in tissue

Mathematical Modeling • Goals • Determine key physical processes, and extract governing parameters • Assess the efficacy of agent penetration in the myocardium using this mode of delivery • Key physical processes • Substrate transport across boundary layer between pericardial sac and myocardium: • Substrate diffusion in myocardium:DT • Substrate washout in myocardium (through the intramural vascular and lymphatic capillaries):k

Idealized Spherical Geometry Pericardial sac: R2 – R3 Myocardium: R1 – R2 Chamber: 0 – R1 R1 = 2.5cm R2 = 3.5cm Vperi= 10ml - 40ml

Governing Equations and Boundary Conditions • Governing equation in myocardium: diffusion + washout CT: concentration of agent in tissue DT: effective diffusion constant in tissue k: washout rate • Pericardial sac as a drug reservoir(well-mixed and no washout): drug number conservation • Boundary condition: drug current at peri/epicardial boundary

Example of Numerical Fits to Experiments Agent Concentration Error surface

Fit Results Numerical values for DT, k,  consistent for IGF, bFGF

Time Course from Simulation Parameters: DT = 7×10-6cm2s-1 k = 5×10-4s-1 a = 3.2×10-6cm2s2

Effective Diffusion, D*, in Tortuous Media • Stokes-Einstein relation D: diffusion constant R: hydrodynamic radius : viscosity T: temperature • Diffusion in tortuous medium D*: effective diffusion constant D: diffusion constant in fluid : tortuosity For myocardium, l= 2.11. (from M. Suenson, D.R. Richmond, J.B. Bassingthwaighte, “Diffusion of sucrose, sodium, and water in ventricular myocardium, American Joural of Physiology,” 227(5), 1974 ) • Numerical estimates for diffusion constants • IGF : D ~ 4 x 10-7 cm2s-1 • bFGF: D ~ 3 x 10-7 cm2s-1 Our fitted values are in order of 10-6 - 10-5 cm2sec-1, 10 to 50 times larger !!

Epi Endo Transport via Intramural Vasculature Drug permeates into vasculature from extracellular space at high concentration and permeates out of the vasculature into the extracellular space at low concentration, thereby increasing the effective diffusion constant in the tissue.

Diffusion in Active Viscoelastic Media Heart tissue is a porous medium consisting of extracellular space and muscle fibers. The extracellular space consists of an incompressible fluid (mostly water) and collagen. Expansion and contraction of the fiber bundles and sheets leads to changes in pore size at the tissue level and therefore mixing of the extracellular volume. This effective "stirring" [1] results in larger diffusion constants. [1] T. Gregor, W. Bialek, R. R. de Ruyter, van Steveninck, et al., PNAS 102, 18403 (2005).

Part I: Conclusions • Model accounting for effective diffusion and washout is consistent with experiments despite its simplicity. • Quantitative determination of numerical values for physical parameters • Effective diffusion constant IGF: DT = (1.7±1.5) x 10-5 cm2s-1, bFGF: DT = (2.4±2.9) x 10-5 cm2s-1 • Washout rate IGF: k = (1.4±0.8) x 10-3 s-1, bFGF: k = (2.1±2.2) x 10-3 s-1 • Peri-epicardial boundary permeability  IGF: a = (4.6±3.2) x 10-6 cm s-1, bFGF: a =(11.9±10.1) x 10-6 cm s-1 • Enhanced effective diffusion, allowing for improved transport • Feasibility of computational studies of amount and time course of pericardial drug delivery to cardiac tissue, using experimentally derived values for physical parameters.

Part II:Electrical Wave Propagation in a Minimally Realistic Fiber Architecture Model of the Left Ventricle: Dynamics of Phase Singularies Xianfeng Song, Department of Physics, Indiana University Sima Setayeshgar, Department of Physics, Indiana University

Part II: Outline • Motivation • Model Construction • Numerical Results • Conclusions and Future Work

The Heart as a Physical System

Motivation • Ventricular fibrillation (VF) is the main cause of sudden cardiac death in industrialized nations, accounting for 1 out of 10 deaths. • Strong experimental evidence suggests that self-sustained waves of electrical wave activity in cardiac tissue are related to fatal arrhythmias. • Goal is to use analytical and numerical tools to study the dynamics of reentrant waves in the heart on physiologically realistic domains. • And … the heart is an interesting arena for applying the ideas of pattern formation. W.F. Witkowksi, et al., Nature 392, 78 (1998) Patch size: 5 cm x 5 cm Time spacing: 5 msec

Big Picture What are the mechanisms underlying the transition from ventricular tachychardia to fibrillation? How can we control it? Tachychardia Fibrillation (Courtesty of Sasha Panfilov, University of Utrecht) Paradigm: Breakdown of a single spiral (scroll) wave into disordered state, resulting from various mechanisms of spiral wave instability

Focus of Our Work Distinguish the role in the generation of electrical wave instabilities of the “passive” properties of cardiac tissue as a conducting medium • geometrical factors (aspect ratio and curvature) • rotating anisotropy (rotation of mean fiber direction through heart wall) • bidomain description (intra- and extra-cellular spaces treated separately)* from its “active” properties, determined by cardiac cell electrophysiology. *Jianfeng Lv: Analytical and computational studies of the bidomain model of cardiac tissue as a conducting medium

Motivated by … “Numerical experiments”: Winfree, A. T. in Progress in Biophysics and Molecular Biology (1997)… Panfilov, A. V. and Keener, J. P. Physica D (1995): Scroll wave breakup due to rotating anisotropy Fenton, F. and Karma, A. Chaos (1998): Rotating anisotropy leads to “twistons”, eventually destabilizing scroll filament Analytical work: In isotropic excitable media Keener, J. P. Physica D (1988) … Biktashev, V. N. and Holden, A. V. Physica D (1994) … In anisotropic excitable media Setayeshgar, S. and Bernoff, A. J. PRL (2002)

From Idealized to Fully Realistic Geometrical Modeling Rectangular slab Anatomical canine ventricular model J.P. Keener, et al., in Cardiac Electrophysiology, eds. D. P. Zipes et al. (1995) Courtesy of A. V. Panfilov, in Physics Today, Part 1, August 1996 Minimally realistic model of LV for studying electrical wave propagation in three dimensional anisotropic myocardium that adequately addresses the role of geometry and fiber architecture and is: • Simpler and computationally more tractable than fully realistic models • Easily parallelizable and with good scalability • More feasible for incorporating realistic electrophysiology, electromechanical coupling, • bidomain description

LV Fiber Architecture Early dissection results revealed nested ventricular fiber surfaces, with fibers given approximately by geodesics on these surfaces. 3d conduction pathway with uniaxial anisotropy: Enhanced conduction along fiber directions. From Textbook of Medical Physiology, Guyton and Hall. cpar = 0.5 m/sec cperp = 0.17 m/sec Anterior view of the fibers on hog ventricles, revealing the nested ventricular fiber surfaces, from C. E. Thomas, Am. J. Anatomy (1957). Fibers on a nested pair of surfaces in the LV, from C. E. Thomas, Am. J. Anatomy (1957).

Peskin Asymptotic Analysis of the Fiber Architecture of the LV: Principles and Assumptions • The fiber structure has axial symmetry • The fiber structure of the left ventricle is in near-equilibrium with the pressure gradient in the wall • The state of stress in the ventricular wall is the sum of a hydrostatic pressure and a fiber stress • The cross-sectional area of a fiber tube does not vary along its length • The thickness of the fiber structure is considerably smaller than its other dimensions.

Peskin Asymptotic Model: Results • The fibers run on a nested family of toroidal surfaces which are centered on a degenerate torus which is a circular fiber in the equatorial plane of the ventricle • The fiber are approximate geodesics on fiber surfaces, and the fiber tension is approximately constant on each surface • The fiber-angle distribution through the thickness of the wall follows an inverse-sine relationship Cross-section of the predicted middle surface (red line) and fiber surfaces (solid lines) in the r, z-plane. Fiber angle profile through LV thickness: Comparison of Peskin asymptotic model and dissection results

Model Construction Nested cone geometry and fiber surfaces • Fiber paths • Geodesics on fiber surfaces • Circumferential at midwall Fiber trajectories on nested pair of conical surfaces: subject to: Fiber trajectory: inner surface outer surface

Electrophysiology: Governing Equations • Transmembrane potential propagation • Transmembrane current, Im, described by simplified FitzHugh-Nagumo type dynamics [1] • Cm: capacitance per unit area of membrane • D: conductivity tensor • u: transmembrane potential • Im: transmembrane current v: gate variable Parameters: a=0.1, m1=0.07, m2=0.3, k=8, e=0.01, Cm=1 [1] R. R. Aliev and A. V. Panfilov, Chaos Solitons Fractals 7, 293 (1996)

Numerical Implementation Working in spherical coordinates, with the boundaries of the computational domain described by two nested cones, is equivalent to computing in a box. Standard centered finite difference scheme is used to treat the spatial derivatives, along with first-order explicit Euler time-stepping.

Conductivity Tensor Transformation matrix R Local Coordinate Lab Coordinate

Parallelization • The communication can be minimized when parallelized along azimuthal direction. • Computational results show the model has a very good scalability.

Phase Singularities Tips and filaments are phase singularities that act as organizing centers for spiral (2D) and scroll (3D) dynamics, respectively, offering a way to quantify and simplify the full spatiotemporal dynamics. Color denotes the transmembrane potential. Movie shows the spread of excitation for 0 < t < 30, characterized by a single filament.

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Find all tips

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Random choose a tip

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Search for the closest tip

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Make connection

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Continue doing search

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Continue

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface The closest tip is too far

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Reverse the search direction

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Continue

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Complete the filament

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Start a new filament

Filament-finding Algorithm “Distance” between two tips: If two tips are not on a same fiber surface or on adjacent surfaces, the distance is defined to be infinity. Otherwise, the distance is the distance along the fiber surface Repeat until all tips are consumed

Filament-finding result FHN Model: t = 2 t = 999

Physics of the Heart: From the macroscopic to the microscopic