Transcript “Distance” between two tips

Physics of the Heart: From the macroscopic
to the microscopic
Xianfeng Song
Advisor: Sima Setayeshgar
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
Part I:
Transport Through the Myocardium of
Pharmocokinetic Agents Placed in the
Pericardial Sac:
Insights From Physical Modeling
Xianfeng Song, Department of Physics, Indiana University
Keith L. March, IUPUI Medical School
Sima Setayeshgar, Department of Physics, Indiana University
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

The pericardial sac is 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.

Recent experimental feasibility of 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
Experimental Procedure
 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  = 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, = 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 !!
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.
Epi
Endo
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:  = (4.6±3.2) x 10-6 cm s-1, bFGF:  =(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
W.F. Witkowksi, et al., Nature 392, 78 (1998)

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 selfsustained 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.
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
J.P. Keener, et al., in Cardiac Electrophysiology, eds.
D. P. Zipes et al. (1995)
Anatomical canine ventricular model
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
Fibers on a nested pair of
ventricles, revealing the nested
surfaces in the LV, from
ventricular fiber surfaces, from
C. E. Thomas, Am. J. Anatomy (1957).
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 inversesine 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
2
L   f ( ,
1
d
,  )d
d
f
d  f 



 d   ' 
z
subject to:
 0  0

Fiber trajectory:
L  0 

  1
 1 
  1  a 2 sec 1 
Fiber trajectories on nested pair of conical surfaces:
inner surface
outer surface
Electrophysiology: Governing Equations
 Transmembrane potential propagation
Cm
u
   ( Du )  I m
t
Cm: capacitance per unit area of membrane
D: conductivity tensor
u: transmembrane potential
Im: transmembrane current
 Transmembrane current, Im, described by simplified
FitzHugh-Nagumo type dynamics [1]
I m  ku(u  a)(u  1)  uv
v 
mv 
    1  v  ku(u  a  1
t 
m2  u 
v: gate variable
Parameters: a=0.1, m1=0.07, m2=0.3,
[1] R. R. Aliev and A. V. Panfilov, Chaos Solitons Fractals 7, 293 (1996)
k=8, =0.01, Cm=1
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
Dlocal
 D//

 0
 0

0
D p1
0
0 

0 
D p 2 
Lab Coordinate
Dlab  R 1 Dlocal R
Parallelization
 The communication can be minimized when parallelized along azimuthal
direction.
 Computational results show the model has a very good scalability.
CPUs
Speed up
2
1.42 ± 0.10
4
3.58 ± 0.16
8
7.61 ±0.46
16
14.95 ±0.46
32
28.04 ± 0.85
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
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
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
t=2
t = 999
FHN Model:
Numerical Convergence
The results for filament length agree to
within error bars for three different mesh
sizes.
The results for filament number agree to
within error bars for dr=0.7 and dr=0.5. The
result for dr=1.1 is slightly off, which could
be due to the filament finding algorithm.
Filament Number and Filament Length
versus Heart size
The computation time for dr=0.7 for one
wave period in a normal heart size is less
than 1 hour of CPU time using FHN-like
electrophysiological model. Fully realistic
model requires several days per heart cycle
on a high-performance machine [1]
[1] Hunter, P. J., A. J. Pullan, et al. (2003), Annual Review of
Biomedical Engineering 5(1): 147-177.
Scaling of Ventricular Turbulence
Log(total filament length) and Log(filament number)
The average filament length, normalized by
Both filament length
versus Log(heart size)
average heart thickness, versus heart size
These results are in agreement with those obtained with the fully realistic
canine anatomical model, using the same electrophysiology. [1]
[1] A. V. Panfilov, Phys. Rev. E 59, R6251 (1999)
Conclusions so far…
 We have constructed and implemented a minimally realistic fiber
architecture model of the left ventricle for studying electrical wave
propagation in the three dimensional myocardium.
 Our model adequately addresses the geometry and fiber architecture of the
LV, as indicated by the agreement of filament dynamics with that from
fully realistic geometrical models.
 Our model is computationally more tractable, allowing reliable numerical
studies. It is easily parallelizable and has good scalability.
 As such, it is more feasible for incorporating
 Realistic electrophysiology
 Bidomain description of tissue
 Electromechanical coupling
Work in Progress
 Computational:
Investigate role of geometry and fiber architecture on scroll
wave stability
(Preliminary results indicate filament instability is suppressed
in minimally realistic model versus rectangular slab!)
 Analytical:
Extend perturbation analysis of scroll waves in the presence of
rotating anisotropy [1] to include filament motion
[1] Setayeshgar, S. and Bernoff, A. J. PRL (2002).
Rotating anisotropy
Coordinate System
Governing Equations
Perturbation Analysis
Scroll Twist Solutions
Twist
Twist
Scroll Twist, Fz
Rotating anisotropy generated scroll twist,
either at the boundaries or in the bulk.
Significance?
In isotropic excitable media ( = 1), for twist > twistcritical,
straight filament undergoes buckling (“sproing”) instability [1]
What happens in the presence of rotating anisotropy ( > 1)??
Henzi, Lugosi and Winfree, Can. J. Phys. (1990).
Filament Motion
Filament motion (cont’d)
Filament Tension
Destabilizing or restabilizing role of rotating anisotropy!!
Part III:
Calcium Dynamics in the Myocyte
Xianfeng Song, Department of Physics, Indiana University
Sima Setayeshgar, Department of Physics, Indiana University
Part III: Outline
Importance and background on calcium
signaling in myocytes
Future directions
Overview of Calcium Signaling
From Berridge, M. J., M. D. Bootman, et al. (1998). "Calcium - a life and death signal." Nature 395(6703): 645-648.
Elementary events (red) result from the entry of external Ca2+ across the plasma membrane or
release from internal stores in the endolasmic or sarcoplasmic reticulum (ER/SR).
Global Ca2+signals are produced by coordinating the activity of elementary events to produce a
Ca2+ wave that spreads throughout the cell.
The activity of neighboring cells within a tissue can be coordinated by an intercellular wave that
spreads from one cell o the next.
Fundamental Elements of Ca2+
Signaling Machinery
Borisyuk, A. (2005). Tutorials in mathematical biosciences. Berlin, Springer
Calcium stores: External and internal stores, i.e. Endoplasmic Reticulum (ER), Sarcoplasmic Reticulum (SR),
Mitochondria
Calcium buffers: Calcium is heavily buffered in all cells, with at least 99% of the available Ca 2+ bound to
large Ca2+-binding proteins., such as Calmodulin, Calsequestrin.
Calcium pumps: Ca2+ is moved to Calcium stores by varies pumps.
Calcium channels: Ca2+ can enter the cytoplasm from calcium stores via varies channels, i.e. ryanodine
receptors (RyR) and inositol trisphosphate receptors (IP3R).
Ventricular Myocyte
Some facts about myocytes
From Borisyuk, A. (2005). Tutorials in mathematical biosciences. Berlin, Springer
The typical cardiac myocyte is a cylindrical cell
approximately 100 mm in length by 10mm in
diameter
Three physical compartments: the cytoplasm, the
sarcoplasmic reticulum (SR) and the mitochondria.
The junctional cleft is a very narrow space
between the SL and the SR membrane.
 Calcium Induced Cacium Release (CICR)
Ventricular Myocyte Structure
A small amount of Ca2+ goes into the junctional
cleft thus induce large scale of Ca2+ release from
calcium stores (mainly SR).
Excitation-Contraction Coupling (ECC)
The depolarization of the membrane initial a small
amount of Ca2+, thus induce CICR and initiate
contraction.
Calcium induced Calcium release
Future Directions
What is the role of receptor clustering on
calcium signaling?
What is the role of the buffer Calsequestrin
in facilitating calcium release?
Receptor Clustering
RyR and IP3R channels are spatially organized in clusters, with the
distance between clusters are approximate two order magnitude larger
than the distances between channels within one cluster.
High resolution image showing a Ca2+ puff evoked by
photoreleased InsP3 which demonstrate an IP3R cluster (From Yao,
Y. etc, Journal of Physiology 482: 533-553.)
Analogy with chemotaxis receptor clustering in E.coli shown to be
important in [1]
 Signal amplification
 Noise reduction
[1] Skoge, M. L., R. G. Endres, et al. (2006), Biophys. J. 90(12): 4317.
The Role of Calsequestrin


Calsequestrin is the buffer inside SR, most of which are located close to
RyRs.
Borisyuk, A. (2005). Tutorials in mathematical biosciences. Berlin, Springer
Calsequestrin plays an important role during CICR.
(A) The channel opens, Ca2+ adsorbed to linear CSQ polymers feeds rapid release. (B) The polymers are depleted Ca2+ thus
disassemble. (C) Depletion becomes deeper as Ca2+ replenishes the proximate store and the CSQ polymers reassembles
(from Launikonis et al, PNAS 103(8) 2982-7 (2005))

Role of Calsequestrin polymerization/depolymerization on its diffusive
uptake of Ca2+ as a store?
Thanks!!