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:

• Niels F. Otani, FitzHugh-Nagumo equation analysis
• Niels F. Otani, Circuit modeling of caridac physiology
• Cardiac Dynamics Exercise

Single-cell analysis of the FitzHugh-Nagumo model

1. Download the file FitzNagHints.py (Hints for FitzHugh-Nagumo ODE model), and rename it as "FitzNag.py". Either load it into a text editor or the ipython dashboard.
2. 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):

   

3. 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.
4. 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.
5. 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.
6. Test you fixed point solution (V*,W*) by integrating using that as an initial condition. Does the trajectory deviate from the initial condition?
7. 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.
8. 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.)
9. 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.
10. 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

1. Download the file FitzNag2DHints.py (Hints for FitzHugh-Nagumo PDE model), and rename it as "FitzNag2D.py".
2. Now consider the extension to the FitzHugh-Nagumo model which now couples the transmembrane potential among neighboring cells :

   

3. 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 Cardiac Dynamics Exercise 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.
4. 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.

5. 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.
   • See Cardiac Dynamics Exercise for discussion of the use of ghost cells in finite-difference methods.
6. 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.
7. 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).

1. 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.
2. 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.
3. 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.
4. 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)
   

5. 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.