Cardiac Dynamics Exercise
Background
Your heartbeat is a wave that passes across your heart muscle. The heart tissue
is an excitable medium: a small triggering impulse can lead to a large response
(an electrical discharge across the cell membranes, together with a
contraction of the heart muscle). A piece of heart tissue can be triggered by
the excitation of a neighboring piece of tissue, which is the basis for
the wave action.
In the normal heart, a pulse propagates outward from the sinoatrial node
(your natural pacemaker, which excites itself). The animation at right
shows the heart undergoing spiral waves, which occur when your heart
is malfunctioning (ventricular tachycardia).
This exercise first studies a simple model for an excitable medium,
the FitzHugh-Nagumo equations, using nullcline analysis to develop
insight into the phase portrait of the action potential. We then
couple these equations into a simple partial differential equation
representing the heart tissue. We learn techniques for solving
dynamical partial differential equations. We develop an interactive
graphical model for the heart allowing the student to apply localized
electrical "shocks" to create and break up spiral waves.
Learning Goals
Science: You will be introduced to the problem of pattern
formation generally, and excitable media specifically, in the context
of a simple model of electrical activity in the heart.
Pattern formation arises broadly in a number of fields including
biology, materials science, fluids, ecology, etc. The method of
nullclines will be first used to characterize the dynamical behavior
of individual cardiac cells in the model, and then you will simulate a
partial differential equation model to explore the
dynamics of spiral wave defects in cardiac tissue.
Computation: You will use root finding methods find the fixed
point of the dynamics of a single cardiac cell. Finite-difference methods
will then be developed for the simulation of the cardiac tissue PDE.
You will develop computationally efficient stencil schemes for computing
spatial derivatives, and integrate visualization tools to allow for
visualization and steering of the running simulation.
Procedure
Consult the following material for orientation and background:
Single-cell analysis of the FitzHugh-Nagumo model
- Download the file FitzNagHints.py (Hints for FitzHugh-Nagumo ODE model), and rename it as "FitzNag.py".
- Consider the FitzHugh-Nagumo model for the electrical activity of
a single cardiac cell (i.e., with no coupling of the transmembrane
potential between cells):
- As described in Exercise N.9(a) in Cardiac Dynamics Exercise, find and plot the nullclines of the single-cell
FitzHugh-Nagumo equations for the default values of the system
parameters: epsilon=0.2, gamma=0.8, and beta=0.7.
This will involve writing a routine
PlotNullclines(gamma, beta) to plot the nullcline curves for a
given values of gamma and beta.
- The resting state of the cell (V*, W*) is the fixed point of the
dynamics. Write a function
FindFixedPoint(gamma, beta) that, for given values of gamma and
beta, will find the values V* and W* such that dV/dt = 0 and dW/dt = 0,
using the root-finding routine scipy.optimize.brentq as
described in the Hints file.
- Write a function RunFitzNag to integrate the model equations
from a prescribed initial condition and given values of model parameters
epsilon, gamma and beta. This should return the full dynamical trajectory
produced by the integrator.
- Test you fixed point solution (V*,W*) by integrating using that as
an initial condition. Does the trajectory deviate from the initial condition?
- Write a function Stimulate to give a "kick" to the system
in its resting state (V*, W*). The kick should be an increase in the
transmembrane potential V by a prescribed amount delta.
- Integrate the
equations of motion for stimuli delta of various sizes. Plot the
trajectories V(t) and W(t) as a function of time, as well as in
the phase plane (the V-W plane) along with nullcline curves. What
happens for small stimuli? What happens for large stimuli? Does the
transition between different behaviors appear sharp?
(See this link for a discussion of this point.)
-
Change the parameter beta to beta=0.4, and integrate the system
dynamics from the resting state. (Run for a reasonably long time, e.g.,
t=100 time units.) Plot V(t) and W(t).
Does the system spontaneously oscillate? (Does the plot V(t) look
a little bit like an electrocardiogram?)
The fixed point (V*,W*)
becomes unstable - giving rise to a Hopf bifurcation -
for beta < beta_c = 7/15.
- There is a nice discussion of the dynamics of the FitzHugh-Nagumo
equations at Scholarpedia. If you're so inclined, you
could try making a movie reminiscent of that shown in Fig. 4.
(Hint, you can use the pylab.savefig command to save successive
frames, and the ImageMagick convert command from the command line
to save a series of gif images to an animated gif, or the
ImageMagick animate command to display an animation directly
to your screen (in X-Windows).
FitzHugh-Nagumo PDE model of cardiac tissue dynamics
- Download the file FitzNag2DHints.py (Hints
for FitzHugh-Nagumo PDE model), and rename it as "FitzNag2D.py".
- Now consider the extension to the FitzHugh-Nagumo model
which now couples the transmembrane potential among neighboring cells
:
- Having introduced a spatial derivative, we will approximate it
via finite differences:
or, changing notation to indicate array indices via subscripts:
- This finite-difference approximation shown above can be summarized
schematically via a stencil that indicates the action of the
differential operator at every site in the array. This is sometimes
called a 5-point stencil since it involves accessing values of the
array on 5 different lattice sites.
- See the notes below for discussion of
stencil representations of finite-difference operators.
Write a function del2_5(A,dx) that will compute the
laplacian of an array A with lattice spacing dx, for the 5-point
stencil.
- Rather than explicitly looping over the indices of the array
A to compute the derivative, use array slicing to create
"shifted" copies of A and array addition to combine those copies.
- Write a alternative function del2_9(A,dx) that will compute the
laplacian of an array A with lattice spacing dx, for this 9-point
stencil:
This is accurate to the same order as the 5-point stencil above, but
better reflects the spherical symmetry of the differential operator.
- Write a function CopyGhost(A) that copies into the edges
of the array A the values of the array from the next row/column in.
This has the effect of establishing a zero normal derivative boundary
condition on all four sides of the array.
- In FitzNag2D.py, fill out the definition of the class
FitzNag2D. This will represent an instance of the simulation
class. You will need to define several methods on this class:
- __init__: This will set up the basic data elements of
the simulation. The fields V and W will be represented
as 2-dimensional scipy arrays. You can set a default grid size to
be 100x100. Initialize V and W
everywhere to be equal to the resting state values V* and W*,
except in a little square of size 10 cells at the center of the grid,
where you will apply a stimulus to V of magnitude delta=3.0.
In addition, you will create an
instance of the DynamicLattice object to enable you to visualize
the dynamics of the field V. (You could construct two of them
if you also want to display W.)
- rhs: This computes the instantaneous right-hand side (rhs)
of the FitzHugh-Nagumo PDE, populating two arrays self.dV_dt
and self.dW_dt describing the time derivatives of the fields
V and W.
- step(dt): This steps the PDE forward by one time step of
size dt, using a forward Euler method. Before computing the right hand
side (rhs) of the PDE, the CopyGhost method should be
applied to enforce the derivative boundary conditions on the field V.
- run(T, dt): This calls step(dt) repeatedly until a
total time increment of size T is reached. In order to steer
and visualize the dynamics, you will want to interact with
the DynamicLattice (self.DL) object previously created. The
DL object can be queried to see if a rectangular box has been selected
(in which case a stimulus pulse to V should be applied in
the box), and the current value of the V field can be
visualized using the display method. You will not want to
interact with self.DL at every time step, but can do so every
DISPLAYINTERVAL steps with a suitable if statement.
- Once you have the basic simulation working, you can explore the
dynamics of the FitzHugh-Nagumo PDE.
- Create an instance of a FitzNag2D object: f = FitzNag2D(). (This will require that you defined default parameter values for the
grid size and system parameters. Alternatively, you can pass them in
explicitly.)
- Advance the simulation forward in time for a specified amount
with a given timestep, e.g., f.run(200, 0.1). If you've seeded
the center of the grid with a pulse in V as instructed, you should
see a disturbance propagate out from the center of the grid.
- Think back to the single-cell FitzHugh-Nagumo model, and the
Stimulate function you wrote. Each cell in the center of
the grid has been stimulated by the pulse applied, and proceeds on
an excursion in the V-W phase plane similar to that studied
previously. The coupling between cells provides a mechanism for cells
to stimulate their neighbors.
- NOTE: If at any time you reach the end of the specified run time,
you can simply issue another command to run again: f.run(200, 0.1),
or however long you want it to go.
- Let the ring spread until it reaches the edges. Can you see evidence
of the zero-derivative boundary conditions as the ring passes out of
the grid?
- Once the ring passes, the simulation is not very exciting, since the
system has more or less returned to its resting state (V*,W*).
Using the mouse, drag a rectangle in the center of the grid of approximately
the same size as the original pulse, which should initiate another stimulus
pulse within that box. You should see another ring spreading
out from the center.
- Now use the mouse again to stimulate a second box that cuts across
a portion of the spreading ring. The intent is to "tear" the circular
ring so as to induce a spiral wave.
- Note that the response of the system
to the stimulus is rather different in the leading edge and the trailing
edge of the ring: why is that?
- Note that the spiral wave rotates at its own frequency, exciting cells
in its neighborhood. Because this excitation is not necessarily
commensurate or in phase with external forcing from the pacemaker cells,
the heart does not beat coherently.
- Introduce more stimuli to produce multiple spiral waves.
- Use the mouse to "shock" the entire system (as one might do with a
defibrillator), and note how spiral waves are erased.
Spatial extensions of the FitzHugh-Nagumo PDE model
For the PDE model described above, the parameters epsilon, gamma and
beta were defined as scalars, i.e., they had the same value for all
cells throughout the system. (Also homogeneous was the implicit coupling
constant between neighboring cells, i.e., the prefactor to the laplacian
of V.) By making these parameters spatially varying, we can
explore other interesting phenomena in the FitzHugh-Nagumo model.
Many of these extensions are described in:
Niels F. Otani,
Further exploration of the FitzHugh-Nagumo model (Project Topics, Section 5,
p. 24).
- Copy your spatially homogeneous code to a new file so that we can
make changes: cp FitzNag2D.py SE_FitzNag2D.py,
where "SE" stands for "spatial extensions". Open SE_FitzNag2D.py
in an editor.
- Modify the __init__() method in your FitzNag2D class
to make eps, gamma and beta two-dimensional arrays of the same shape
as the V and W arrays. They should be initialized to
be uniform, with the values specified by the input parameters eps, gamma,
and beta, respectively. Also introduce a new array (self.D) which
will serve as a diffusion coefficient; initialize it to 1.0 across the
entire array.
- Modify the rhs method to include the new spatially-varying
diffusion coefficient in the equation for dV/dt, which should now be
self.D*del2(self.v, self.dx). Note that no other changes need to
be made to the equations of motion, because the syntax for elementwise
multiplication of two arrays is the same as for multiplication of an
array by a scalar.
- Spontaneous oscillations: Rather than introducing a stimulus
to initiate an excitation, we can reduce the parameter beta locally to
induce spontaneous oscillations (as in
an earlier exercise). In a 10x10 box somewhere in the grid (e.g.,
the center), set the beta field locally to 0.4, run the simulation for
a period of time, and use the mouse to introduce spiral waves: e.g.
f = FitzNag2D()
f.beta[f.N/2-5:f.N/2+5, f.N/2-5:f.N/2+5] = 0.4
f.run(100, 0.1)
- Cross field stimulation, dead tissue, and compartmentalization:
Following the discussion in
Further exploration of the FitzHugh-Nagumo model (Project Topics, Section 5,
p. 24), explore other phenomena by suitably altering the spatial
parameters in SE_FitzNag2D.py.
- Note: the term "gap junction"
refers to the coupling between neighboring cells, i.e., what we have
referred to as the "diffusion coefficient" previously. With the way we
have implemented the gap junction (as an array multiplying the laplacian),
you can disconnect a cell from all its neighbors (D=0) but not selectively
disconnect it from only some neighbors.
Files
References
Links
Notes
James P. Sethna,
Christopher R. Myers.
Niels F. Otani.
Last modified: August 26, 2008