Numerical Simulation of Calcium Waves in Heart Cells
Matthias K. Gobbert and Bradford E. Peercy
Department of Mathematics and Statistics
University of Maryland, Baltimore County
This page can be reached via the homepage
http://www.math.umbc.edu/~gobbert.
Overview
Calcium ions play a crucial role in driving the heart beat. Calcium is
released into the heart cell at a lattice of discrete positions throughout
the cell, known as calcium release units (CRUs). The probability that
calcium will be released from a CRU depends on the calcium concentration
at that CRU. When a CRU releases calcium (it 'fires'), the local
concentration of calcium ions increases sharply and the calcium that
diffuses raises the probability for release at neighboring sites. As a
sequence of CRUs begins to release calcium throughout the cell, the
release can selforganize into a wave of increasing concentration. If the
wave is triggered amongst the normal physiological signaling (i.e., the
cardiac action potential) then it can lead to irregular heart beats and
possibly a life threatening ventricular fibrillation.
A mathematical model for this effect has been developed by Leighton T. Izu
and coworkers
[Biophysical Journal, pp. 88102 and pp. 103120, 2001]. It
consists of a system of three timedependent partial differential
equations, socalled reactiondiffusion equations, which are coupled
through nonlinear reaction terms. The most important component in the
driving terms of the equation for the calcium concentration is the
superposition of Dirac delta distributions that models the CRUs as point
sources.
Realistic simulations can only be obtained using mathematical models
that discretize the the threedimensional cell domain with
a highresolution numerical mesh.
Therefore, a specialpurpose computer code is needed that
pays great attention to the efficient use of computer memory
and that is able to pool the memory from several nodes
in a parallel computing cluster.
Moreover, large final times are required to reach laboratory time scales
and to simulate for long enough to allow several waves to propagate
through the cell.
The numerical solution of this problem is made additionally
challenging by the discontinuities caused by the Dirac delta
distributions, for which the classical finite element theory does not
apply. This means that standard numerical analysis cannot
guarantee the convergence of the numerical method,
which is necessary to trust the computed solutions.
Initial work by graduate student Alexander Hanhart [M.S. thesis, 2002],
funded by a UMBC DRIF grant, demonstrated that a parallel implementation
can achieve the required high resolution of the numerical mesh.
Moreover, it confirmed computationally that a finite element method with
trilinear nodal basis functions converges for this problem
[Journal of Computational and Applied Mathematics,
pp. 431458, 2004].
Work involving several students, including
the undergraduate student Kevin Allen and the high school senior
Greg McGlynn, created a full specialpurpose code for the
application that uses stateoftheart timestepping and
matrixfree linear solvers. Its confirmed numerical convergence,
scalable parallel performance,
ability to achieve spontaneous waveinitiation,
and correct representation of vital physical effects such as mass conservation
make it suitable for longtime simulations on the scales
needed for the application
[Gobbert, SIAM Journal on Scientific Computing,
pp. 29222947, 2008].
Such simulations are made possible by the computing cluster tara
with stateoftheart Intel Nehalem processors
and extended memory of 24 GB per node that is part of the
UMBC High Performance Computing Facility (www.umbc.edu/hpcf),
funded in part by the National Science Foundation
through MRI and SCREMS grants with
additional significant support from UMBC.
Very recent work in the department
has allowed us to solve two key issues:
The convergence of the finite element method is now rigorously established
together with departmental colleague Thomas I. Seidman
[Seidman, Gobbert, Trott, and Kruzik, Finite Element Approximation for
TimeDependent Diffusion with MeasureValued Source, submitted (2011)],
in addition to previous computational evidence.
However, the original values of the experimental parameters for the code
resulted in unrealistic physiological behavior.
This issue was now resolved by graduate student Zana Coulibaly,
funded by a National Science Foundation and a UMBC DRIF grant,
and his advisor Bradford E. Peercy.
A directed search of parameter space indicated by a lower order model has
established that the 3D model and its implementation can yield waves
with recovery [Coulibaly, Peercy, and Gobbert, 3D Cardiac Cell Based on
Analysis of a 1D Deterministic Model, in preparation (2011)].
Moreover, further studies produced spontaneous
spiral waves [Coulibaly, Peercy, and Gobbert, Spontaneous Spiral Wave
Initiation in a 3D Cardiac Cell, in preparation (2011)],
which are an experimentally observed phenomenon.
The first plot shows a simulated confocallike image of the calcium
concentration with a spiral shape that is rotating clockwise.
This image of a simulation result is designed to mirror the
appearance of an experimental result.
Simulated confocallike image of the calcium concentration
with a spiral wave.
The second plot shows an iso surface that
brings out the threedimensional shape of the region of
highly elevated calcium concentration and which also shows
the bending spiral shape in the right half of the cell.
Threedimensional iso surface plot of
elevated calcium concentration.
This result demonstrates the significance of the accomplishments
in that our simulator produces physiologically correct behavior
in full threedimensional simulations
and can provide a better understanding of physiological processes
than available from other simulators to date.
With this ability, we hope to provide input
to help guide in what ranges of parameters to focus experiments,
thus making research significantly more effective.
Work on improving the numerical efficiency of the code continues,
which will help make runs faster and thus allow for more
rapid parameter studies as would be necessary for the support
of experiments. In collaboration with
graduate student Yu Wang, supported by UMBC as HPCF RA,
and Marc Olano from
the Department of Computer Science and Engineering at UMBC,
we are working on leveraging cuttingedge GPGPUs on each node
in the parallel computing cluster
to improve the computational kernel of the simulator.
People Involved
 Xuan Huang, graduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County
 Tylynn Pettrey, undergraduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County
 Zana Coulibaly, graduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County
 Michael Muscedere, graduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County
 Dr. Thomas I. Seidman,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County
 Jonas Schäfer, graduate student,
Universität Kassel, Germany
 Dr. Andreas Meister,
Universität Kassel, Germany
 Yu Wang, graduate student,
Department of Computer Science and Electrical Engineering,
University of Maryland, Baltimore County
 Dr. Marc Olano,
Department of Computer Science and Electrical Engineering,
University of Maryland, Baltimore County
Students Formerly Involved
 David W. Trott, graduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County,
M.S. May 2011;
continued to Ph.D. program
 Kyle Stern, graduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County,
M.S. May 2010;
went to private industry
 Greg E. McGlynn, Catonsville High School,
graduated Spring 2007;
went to Carnegie Mellon University
 Kevin P. Allen, undergraduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County,
B.S. May 2003;
continued to M.S. program
 Alexander L. Hanhart, graduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County,
M.S. May 2002;
continued to Ph.D. studies at the
School of Mathematics, University of Minnesota
Spontaneous Calcium Waves
For accurate representation of the experiments,
it is vital that computationally created waves are
both selfinitiated and sustained,
while also returning the calicum concentrations to basal levels
in between waves.
This is borne out by the waves with recovery below.
Spontaneous spiral waves are phenomenon that can be seen experimentally.
It is extremely exciting that our simulation tool
has produced this type of wave below!
Both results by Zana Coulibaly and Dr. Peercy
(publications in preparation (2011)).
 Simulations using the finite volume method:
Movie
Movie for iPad
of the open CRUs in the cell during the spiral waves
(79 MB mpegfile);
Movie
Movie for iPad
of isosurfaces of the calcium concentration with isolevel 65 uM
(93 MB mpegfile);
Movie
Movie for iPad
of simulated confocallike images of the calcium concentration
(31 MB mpegfile).
 Spontaneous spiral wave:
The confocallike image shows the spiral shape of the wave most clearly.
The movie of the open CRUs shows the repeated waves moving through the cells.
The isosurface plots show the recovery to basal levels in between waves
(contrary to older simulations farther below on this page).
The output time step size is 1 ms.
Movie
of simulated confocallike images of the calcium concentration
(48 MB mpegfile);
Movie
of the open CRUs in the cell during the spiral waves
(81 MB mpegfile);
Movie
of isosurfaces of the calcium concentration with isolevel 65 uM
(86 MB mpegfile).
 Waves with recovery:
Movie
of simulated confocallike images of the calcium concentration
(4 MB mpegfile);
Movie
of the open (blue) and closed (red) CRUs in the cell
(129 MB mpegfile).
LongTime Simulations on High Resolution Meshes to Model Calcium Waves
Longtime simulation of calcium flow in a cell
of size 12.8by12.8by64.0 um with a 15by15by31 lattice of CRUs
using a 64by64by256 element mesh
from time 0 ms to 1000 ms without any initial stimulation.
The simulations use the parameters and boundary conditions of the
threedimensional model as specified in Gobbert (2008).

Current through CRU = 10 pA:
No wave selforganizes.
Movie
of the open CRUs in the cell (24 MB mpegfile);
notice output time step size is 1 ms.
Movie
of the calcium concentration with isolevel 65 uM (2 MB mpegfile);
notice output time step size is 10 ms.

Current through CRU = 20 pA:
Two waves selforganize and travel through the cell repeatedly.
Movie
of the open CRUs in the cell (27 MB mpegfile);
notice output time step size is 1 ms.
Movie
of the calcium concentration with isolevel 65 uM (3 MB mpegfile);
notice output time step size is 10 ms.
Parallel Performance Studies for Scalar Test Problem
Scalability studies for a timedependent linear scalar test problem
with realistic diffusion coefficients on the domain and on
time scales representative of the application problem,
studies up to 32 processors using one CPU per node
on the distributedmemory cluster
kali
with progressively finer meshes up to final time 1000 ms.
The test problem and more details are contained in
Gobbert (CiSE, 2005).
ThreeDimensional Results Using the Finite Element Method

Isosurface plots of the concentration of calcium ions at times
0 ms,
1 ms,
2 ms,
3 ms,
4 ms,
5 ms,
6 ms,
7 ms,
8 ms,
9 ms,
10 ms.
The calcium release selforganizes into a wave.
Click
here
for a complete movie of the effect.
Results by Alexander L. Hanhart.

Comparison of two numerical solutions for a test problem:
nonlinear threeequation model, realistic numerical resolution,
but unphysical domain size and CRU placement.
Isosurface of the calcium concentration at one particular time step:
Grid 64by32by32,
Grid 128by64by64.
This shows that the numerical method converges.
Results by Alexander L. Hanhart.
Publications Resulting From This Research
In reverse chronological order
The following list is a partial list with a focus on
student involvement in the work; please go to the
publication list
on my webpage
for additional papers.

Zana A. Coulibaly^{G}, Bradford E. Peercy, and
Matthias K. Gobbert.
Spontaneous Spiral Wave Initiation in a 3D Cardiac Cell.
In preparation.

Zana A. Coulibaly^{G}, Bradford E. Peercy, and
Matthias K. Gobbert.
Insight Into Spontaneous Recurrent Calcium Waves in a
3D Cardiac Cell Based on Analysis of a 1D Deterministic Model.
In preparation.

Jonas Schäfer^{G}, Xuan Huang^{G},
Stefan Kopecz, Philipp Birken, Matthias K. Gobbert, and Andreas Meister.
A MemoryEfficient Finite Volume Method for AdvectionDiffusionReaction
Systems with NonSmooth Sources.
Submitted.
Preprint in PDFformat.

Thomas I. Seidman, Matthias K. Gobbert, David W. Trott^{G}, and
Martin Kružík.
Finite Element Approximation for TimeDependent
Diffusion with MeasureValued Source.
Numerische Mathematik,
vol. 122, no. 4, pp. 709723, 2012.
Link to this article at SpringerLink
and preprint in PDFformat.

David W. Trott^{G} and Matthias K. Gobbert.
Finite Element Convergence for TimeDependent PDEs
with a Point Source in COMSOL 4.
In: Yeswanth Rao, editor,
Proceedings of the COMSOL Conference 2011, Boston, MA,
6 pages, 2011.
Reprint in PDFformat.
David W. Trott^{G} and Matthias K. Gobbert.
Parallel Performance Studies for a ThreeSpecies Application Problem
on the Cluster tara,
Technical Report number HPCF201011,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2010.
Reprint in PDFformat.

David W. Trott^{G} and Matthias K. Gobbert.
Finite Element Convergence Studies using COMSOL 4.0a and LiveLink for MATLAB.
Technical Report number HPCF20108,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2010.
Reprint in PDFformat.

Michael Muscedere^{G}, Andrew M. Raim^{G},
and Matthias K. Gobbert.
Parallel Performance Studies for a Parabolic Test Problem on the Cluster tara.
Technical Report number HPCF20104,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2010.
Reprint in PDFformat.

Andrew M. Raim^{G} and Matthias K. Gobbert.
Parallel Performance Studies for an Elliptic Test Problem on the Cluster tara.
Technical Report number HPCF20102,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2010.
Reprint in PDFformat.

Michael Muscedere^{G} and Matthias K. Gobbert.
Parameter Study of a ReactionDiffusion System Near the Reactant Coefficient
Asymptotic Limit.
Dynamics of Continuous, Discrete and Impulsive Systems
Series A Supplement,
pp. 2936, 2009.
Preprint in PDFformat.

Zana Coulibaly^{UG}, Michael Muscedere^{G},
Matthias K. Gobbert, and Bradford E. Peercy.
LongTime Simulation of Calcium Waves in a Heart Cell to Study the
Effects of Calcium Release Flux Density and of Coefficients in
the Pump and Leak Mechanisms on SelfOrganizing Wave Behavior.
Technical Report number HPCF20096,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2009.
Reprint in PDFformat.

Matthias K. Gobbert and Shiming Yang^{G}.
Numerical Demonstration of Finite Element Convergence for Lagrange Elements
in COMSOL Multiphysics.
In: Vineet Dravid, editor,
Proceedings of the COMSOL Conference 2008, Boston, MA,
6 pages, 2008.
Reprint in PDFformat;
see also the
COMSOL area
of my homepage.

Shiming Yang^{G} and Matthias K. Gobbert.
Convergence Order Studies for Elliptic Test Problems with COMSOL Multiphysics.
Technical Report number HPCF20084,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2008.
Reprint in PDFformat.

Michael Muscedere^{G} and Matthias K. Gobbert.
Parallel Performance Studies for a Parabolic Test Problem.
Technical Report number HPCF20082,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2008.
Reprint in PDFformat.

Matthias K. Gobbert.
Parallel Performance Studies for an Elliptic Test Problem.
Technical Report number HPCF20081,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2008.
Reprint in PDFformat.

Matthias K. Gobbert.
LongTime Simulations on High Resolution Meshes
to Model Calcium Waves in a Heart Cell.
SIAM Journal on Scientific Computing,
vol. 30, no. 6, pp. 29222947, 2008.
Direct link to the abstract
and
reprint in PDFformat.

Matthias K. Gobbert.
Configuration and Performance of a Beowulf Cluster for
LargeScale Scientific Simulations.
Computing in Science and Engineering (CiSE),
vol. 7, no. 2, pp. 1426, March/April 2005.
Link to IEEE's Xplore site for this article
and reprint in PDFformat.

Kevin P. Allen^{UG}.
Efficient Parallel Computing for Solving Linear Systems of Equations.
UMBC Review: Journal of Undergraduate Research and Creative Works,
vol. 5, pp. 817, 2004.

Kevin P. Allen^{UG}.
A Parallel MatrixFree Implementation of the Conjugate Gradient Method
for the Poisson Equation.
Senior Thesis,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County, May 2003.

Kevin P. Allen^{UG} and Matthias K. Gobbert.
A MatrixFree Conjugate Gradient Method for Cluster Computing.
Technical Report, University of Maryland, Baltimore County, 2003.

Alexander L. Hanhart^{G}, Matthias K. Gobbert,
and Leighton T. Izu.
A MemoryEfficient Finite Element Method for Systems
of ReactionDiffusion Equations with NonSmooth Forcing.
Journal of Computational and Applied Mathematics,
vol. 169, no. 2, pp. 431458, 2004.
Link to Elsevier's ScienceDirect Service for this article
and reprint in PDFformat.

Kevin P. Allen^{UG} and Matthias K. Gobbert.
CoarseGrained Parallel MatrixFree Solution of
a ThreeDimensional Elliptic Prototype Problem.
In: Vipin Kumar, Marina L. Gavrilova, Chih Jeng Kenneth Tan,
and Pierre L'Ecuyer, editors,
Computational Science and Its Applications  ICCSA 2003, Part II,
Lecture Notes in Computer Science, vol. 2668, pp. 290299,
SpringerVerlag, 2003.
Preprint in PDFformat.

Alexander L. Hanhart^{G}.
CoarseGrained Parallel Solution of a ThreeDimensional Model
for Calcium Concentration in Human Heart Cells.
M.S. thesis,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County, 2002.

Alexander L. Hanhart^{G}, Matthias K. Gobbert,
and Leighton T. Izu.
CoarseGrained Parallel Solution for a System of Reaction Diffusion Equations.
Technical Report, University of Maryland, Baltimore County, 2002.
Disclaimer
The information provided on this page is related to the numerical
solution of a system of partial differential equation that is
designed to model physiological phenoma in human heart cells.
This information does not constitute medical advice. Please see
your physician or other medical professional, if you believe to
have health problems.
Copyright © 20002013 by Matthias K. Gobbert. All Rights Reserved.
This page version 8.0, July 2013.