Numerical Simulation of Calcium Waves in a Heart Cell
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.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 self-organize 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 co-workers
[Biophysical Journal, vol. 80, pp. 88-102 and pp. 103-120, 2001].
It consists of a system of three time-dependent partial differential
equations, so-called reaction-diffusion equations, which are coupled
through non-linear 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 three-dimensional cell domain with
a high-resolution numerical mesh.
Therefore, a special-purpose 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
tri-linear nodal basis functions converges for this problem
[Journal of Computational and Applied Mathematics,
vol. 169, no. 2, pp. 431-458, 2004].
Work involving several students, including
the undergraduate student Kevin Allen and the high school senior
Greg McGlynn, created a full special-purpose code for the
application that uses state-of-the-art time-stepping and
matrix-free linear solvers. Its confirmed numerical convergence,
scalable parallel performance,
ability to achieve spontaneous wave-initiation,
and correct representation of vital physical effects such as mass conservation
make it suitable for long-time simulations on the scales
needed for the application
[Gobbert, SIAM Journal on Scientific Computing,
vol. 30, no. 6, pp. 2922-2947, 2008].
Such simulations are made possible by the computing cluster tara
with state-of-the-art 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
[Numerische Mathematik, vol. 122, no. 4, pp. 709-723, 2012],
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 3-D model and its implementation can yield waves
with recovery [International Journal of Computer Mathematics,
published online 2014].
Moreover, further studies produced spontaneous
spiral waves [Numerical Methods Partial Differential Equations,
vol. 31, no. 1, pp. 143-167, 2015],
which are an experimentally observed phenomenon.
The first plot shows a simulated confocal-like 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 confocal-like image of the calcium concentration
with a spiral wave.
The second plot shows an iso surface that
brings out the three-dimensional shape of the region of
highly elevated calcium concentration and which also shows
the bending spiral shape in the right half of the cell.
Three-dimensional 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 three-dimensional 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 [Xuan Huang et al.,
Order Investigation of Scalable Memory-Efficient Finite Volume Methods for
Parabolic Advection-Diffusion-Reaction Equations with Point Sources,
In preparation].
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 cutting-edge GPGPUs on each node
in the parallel computing cluster
to improve the computational kernel of the simulator
[23rd High Performance Computing Symposium (HPC 2015),
accepted (2015)].
The material of this project was used as a
CSE Success Story
for the SIAM Report on the State of Computational Science and Engineering.
People Involved
- Jonathan Graf, graduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County
- Samuel Khuvis, graduate 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
- 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
- Xuan Huang, graduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County,
Ph.D. May 2015;
- Matthew W. Brewster, undergraduate student,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County,
B.A. May 2014;
went to private industry
- Jonas Schäfer, graduate student,
Universität Kassel, Germany,
Dipl.-Math. 2013;
went to private industry
- 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 self-initiated 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 mpeg-file);
Movie
Movie for iPad
of isosurfaces of the calcium concentration with isolevel 65 uM
(93 MB mpeg-file);
Movie
Movie for iPad
of simulated confocal-like images of the calcium concentration
(31 MB mpeg-file).
- Spontaneous spiral wave:
The confocal-like 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 confocal-like images of the calcium concentration
(48 MB mpeg-file);
Movie
of the open CRUs in the cell during the spiral waves
(81 MB mpeg-file);
Movie
of isosurfaces of the calcium concentration with isolevel 65 uM
(86 MB mpeg-file).
- Waves with recovery:
Movie
of simulated confocal-like images of the calcium concentration
(4 MB mpeg-file);
Movie
of the open (blue) and closed (red) CRUs in the cell
(129 MB mpeg-file).
Long-Time Simulations on High Resolution Meshes to Model Calcium Waves
Long-time simulation of calcium flow in a cell
of size 12.8-by-12.8-by-64.0 um with a 15-by-15-by-31 lattice of CRUs
using a 64-by-64-by-256 element mesh
from time 0 ms to 1000 ms without any initial stimulation.
The simulations use the parameters and boundary conditions of the
three-dimensional model as specified in Gobbert (2008).
-
Current through CRU = 10 pA:
No wave self-organizes.
Movie
of the open CRUs in the cell (24 MB mpeg-file);
notice output time step size is 1 ms.
Movie
of the calcium concentration with isolevel 65 uM (2 MB mpeg-file);
notice output time step size is 10 ms.
-
Current through CRU = 20 pA:
Two waves self-organize and travel through the cell repeatedly.
Movie
of the open CRUs in the cell (27 MB mpeg-file);
notice output time step size is 1 ms.
Movie
of the calcium concentration with isolevel 65 uM (3 MB mpeg-file);
notice output time step size is 10 ms.
Parallel Performance Studies for Scalar Test Problem
Scalability studies for a time-dependent 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 distributed-memory cluster
kali
with progressively finer meshes up to final time 1000 ms.
The test problem and more details are contained in
Gobbert (CiSE, 2005).
Three-Dimensional 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 self-organizes 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:
non-linear three-equation model, realistic numerical resolution,
but unphysical domain size and CRU placement.
Isosurface of the calcium concentration at one particular time step:
Grid 64-by-32-by-32,
Grid 128-by-64-by-64.
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 this project; please go to the
publication list
on my homepage for additional papers.
Student co-authors are indicated by superscripts,
with U for an undergraduate and
G for a graduate student.
-
Xuan HuangG, Matthias K. Gobbert, Bradford E. Peercy,
Stefan Kopecz, Philipp Birken, and Andreas Meister.
Order Investigation of Scalable Memory-Efficient Finite Volume Methods for
Parabolic Advection-Diffusion-Reaction Equations with Point Sources.
In preparation.
-
Matthew W. BrewsterU, Jonathan S. GrafG,
Xuan HuangG, Zana CoulibalyG,
Matthias K. Gobbert, and Bradford E. Peercy.
Calcium Induced Calcium Release with Stochastic Uniform Flux Density
in a Heart Cell.
47th Summer Computer Simulation Conference 2015 (SCSC 2015),
6 pages, 2015.
Preprint in PDF-format.
-
Xuan HuangG and Matthias K. Gobbert.
Long-time Simulation of Calcium Induced Calcium Release
In A Heart Cell using Finite Element Method on a Hybrid CPU/GPU Node.
23rd High Performance Computing Symposium (HPC 2015),
8 pages, 2015.
Preprint in PDF-format.
-
Zana A. CoulibalyG, Bradford E. Peercy, and
Matthias K. Gobbert.
Insight into Spontaneous Recurrent Calcium Waves in a
3-D Cardiac Cell Based on Analysis of a 1-D Deterministic Model.
International Journal of Computer Mathematics,
vol. 92, no. 3, pp. 591-607, 2015.
Link to this article at Taylor and Francis online
and preprint in PDF-format.
-
Jonas SchäferG, Xuan HuangG,
Stefan Kopecz, Philipp Birken, Matthias K. Gobbert, and Andreas Meister.
A Memory-Efficient Finite Volume Method for Advection-Diffusion-Reaction
Systems with Non-Smooth Sources.
Numerical Methods Partial Differential Equations,
vol. 31, no. 1, pp. 143-167, 2015.
Link to this article at the Wiley Online Library
and preprint in PDF-format.
-
Matthew W. BrewsterU.
The Influence of Stochastic Parameters
on Calcium Waves in a Heart Cell.
Senior Thesis,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County, May 2014.
Reprint in PDF-format.
-
Xuan HuangG and Matthias K. Gobbert.
Parallel Performance Studies for a Three-Species Application Problem
on the Cluster maya.
Technical Report number HPCF-2014-8,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2014.
Reprint in PDF-format.
-
Jonathan GrafG and Matthias K. Gobbert.
Parallel Performance Studies for a Parabolic Test Problem
on the Cluster maya.
Technical Report number HPCF-2014-7,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2014.
Reprint in PDF-format.
-
Samuel KhuvisG and Matthias K. Gobbert.
Parallel Performance Studies for an Elliptic Test Problem
on the Cluster maya.
Technical Report number HPCF-2014-6,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2014.
Reprint in PDF-format.
-
Thomas I. Seidman, Matthias K. Gobbert, David W. TrottG, and
Martin Kružík.
Finite Element Approximation for Time-Dependent
Diffusion with Measure-Valued Source.
Numerische Mathematik,
vol. 122, no. 4, pp. 709-723, 2012.
Link to this article at SpringerLink
and preprint in PDF-format.
-
Yu WangG, Marc Olano, Matthias Gobbert, and
Wesley GriffinG.
Parallel Computing for Long-Time Simulations of Calcium Waves in a Heart Cell.
Proceedings in Applied Mathematics and Mechanics (PAMM),
vol. 12, pp. 637-638, 2012.
Link to this article at the Wiley Online Library and
Reprint in PDF-format.
-
David W. TrottG and Matthias K. Gobbert.
Finite Element Convergence for Time-Dependent 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 PDF-format.
David W. TrottG and Matthias K. Gobbert.
Parallel Performance Studies for a Three-Species Application Problem
on the Cluster tara,
Technical Report number HPCF-2010-11,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2010.
Reprint in PDF-format.
-
David W. TrottG and Matthias K. Gobbert.
Finite Element Convergence Studies using COMSOL 4.0a and LiveLink for MATLAB.
Technical Report number HPCF-2010-8,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2010.
Reprint in PDF-format.
-
Michael MuscedereG, Andrew M. RaimG,
and Matthias K. Gobbert.
Parallel Performance Studies for a Parabolic Test Problem on the Cluster tara.
Technical Report number HPCF-2010-4,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2010.
Reprint in PDF-format.
-
Andrew M. RaimG and Matthias K. Gobbert.
Parallel Performance Studies for an Elliptic Test Problem on the Cluster tara.
Technical Report number HPCF-2010-2,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2010.
Reprint in PDF-format.
-
Michael MuscedereG and Matthias K. Gobbert.
Parameter Study of a Reaction-Diffusion System Near the Reactant Coefficient
Asymptotic Limit.
Dynamics of Continuous, Discrete and Impulsive Systems
Series A Supplement,
pp. 29-36, 2009.
Preprint in PDF-format.
-
Zana CoulibalyU, Michael MuscedereG,
Matthias K. Gobbert, and Bradford E. Peercy.
Long-Time 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 Self-Organizing Wave Behavior.
Technical Report number HPCF-2009-6,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2009.
Reprint in PDF-format.
-
Matthias K. Gobbert and Shiming YangG.
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 PDF-format;
see also the
COMSOL area
of my homepage.
-
Shiming YangG and Matthias K. Gobbert.
Convergence Order Studies for Elliptic Test Problems with COMSOL Multiphysics.
Technical Report number HPCF-2008-4,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2008.
Reprint in PDF-format.
-
Michael MuscedereG and Matthias K. Gobbert.
Parallel Performance Studies for a Parabolic Test Problem.
Technical Report number HPCF-2008-2,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2008.
Reprint in PDF-format.
-
Matthias K. Gobbert.
Parallel Performance Studies for an Elliptic Test Problem.
Technical Report number HPCF-2008-1,
UMBC High Performance Computing Facility,
University of Maryland, Baltimore County, 2008.
Reprint in PDF-format.
-
Matthias K. Gobbert.
Long-Time Simulations on High Resolution Meshes
to Model Calcium Waves in a Heart Cell.
SIAM Journal on Scientific Computing,
vol. 30, no. 6, pp. 2922-2947, 2008.
[Special issue on Computational Science & Engineering.]
Direct link to the abstract
and
reprint in PDF-format.
-
Matthias K. Gobbert.
Configuration and Performance of a Beowulf Cluster for
Large-Scale Scientific Simulations.
Computing in Science and Engineering (CiSE),
vol. 7, no. 2, pp. 14-26, March/April 2005.
Link to IEEE's Xplore site for this article
and reprint in PDF-format.
-
Kevin P. AllenU.
Efficient Parallel Computing for Solving Linear Systems of Equations.
UMBC Review: Journal of Undergraduate Research and Creative Works,
vol. 5, pp. 8-17, 2004.
-
Kevin P. AllenU.
A Parallel Matrix-Free 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. AllenU and Matthias K. Gobbert.
A Matrix-Free Conjugate Gradient Method for Cluster Computing.
Technical Report, University of Maryland, Baltimore County, 2003.
-
Alexander L. HanhartG, Matthias K. Gobbert,
and Leighton T. Izu.
A Memory-Efficient Finite Element Method for Systems
of Reaction-Diffusion Equations with Non-Smooth Forcing.
Journal of Computational and Applied Mathematics,
vol. 169, no. 2, pp. 431-458, 2004.
Link to Elsevier's ScienceDirect Service for this article
and reprint in PDF-format.
-
Kevin P. AllenU and Matthias K. Gobbert.
Coarse-Grained Parallel Matrix-Free Solution of
a Three-Dimensional 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. 290-299,
Springer-Verlag, 2003.
Preprint in PDF-format.
-
Alexander L. HanhartG.
Coarse-Grained Parallel Solution of a Three-Dimensional Model
for Calcium Concentration in Human Heart Cells.
M.S. thesis,
Department of Mathematics and Statistics,
University of Maryland, Baltimore County, 2002.
-
Alexander L. HanhartG, Matthias K. Gobbert,
and Leighton T. Izu.
Coarse-Grained 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 © 2000-2015 by Matthias K. Gobbert. All Rights Reserved.
This page version 9.1, May 2015.