# Numerical Simulation of Calcium Waves in a Heart Cell

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

1. 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.
2. 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.
3. 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.
4. 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.
5. 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.
6. 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.
7. 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.
8. 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.
9. 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.
10. 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.
11. 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.
12. 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.

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

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

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

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

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

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

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

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

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

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

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

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

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

26. Kevin P. AllenU and Matthias K. Gobbert. A Matrix-Free Conjugate Gradient Method for Cluster Computing. Technical Report, University of Maryland, Baltimore County, 2003.

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

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

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

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