Computing the stability of steady-state solutions of mathematical models of the electrical activity in the heart
Introduction
In each heartbeat an electrical wave traverses the entire cardiac muscle in order to initiate a well-organized contraction of the muscle. The synchronized contraction is crucial for efficient pumping of blood and thus is a necessary condition for the well being of the body. Cardiac arrhythmias are perturbations to the normal wave which may represent a major risk; about 300,000 US-citizens suffer from sudden cardiac death every year (see e.g. Weiss et al. [1]). Mathematical models of the electrical activity in the heart have been an active area of research for almost 50 years; see Keener and Sneyd [2] for an introduction. It is commonly believed that the electrical wave traversing the cardiac muscle is adequately modeled by the Monodomain model [2], Here v denotes the transmembrane potential, denotes the spatial diffusion coefficient, I denotes the ionic current density scaled by the cell membrane capacitance, g denotes a vector of gating variables , and r denotes various concentrations . We assume that the function G has n components, and that it is semi-diagonal in the sense that each component takes the formThe function R is assumed to be a general non-linear function with m components. The model is defined on a domain and v satisfies a no-flux condition at the boundary of .
Suppose for a moment that we study a single cell. Then the system above reduces to a system of ordinary differential equations of the form Beginning with the seminal paper by Hodgkin and Huxley from 1952 (see [3]), a series of mathematical models of this form have been developed for excitable cells. The web page [4] maintains a comprehensive survey of cell models.
It is well known that computational arrhythmias can be generated in models of the form (1). For illustration we have graphed the solution of the Hund-Rudy [5] model at time t=100, 500 and 1000 ms after the onset of the second stimulus; see Fig. 1, and see [6] for further details on this specific computation. First a planar wave is initiated, then after a critical time a second planar, and orthogonal wave is started. Without the extra stimulus, the solution would simply be a traveling wave solution traversing the computational domain. But the introduction of an extra stimulus at a specific time leads to a very complex solution. It has been known for a long time that an extra stimulus can generate arrhythmias in computational models of cardiac muscles; see e.g. Winfree [7].
Over the past 20 years, numerous studies have verified the existence of fibrillatory computational solutions based on variants of the system (1). For examples of computational arrhythmias consult e.g. Aliev and Panfilov [8], Tran et al. [9], Xie et al. [10], Gong et al. [11], Reumann et al. [12], or ten Tusscher et al. [13].
As mentioned above, the common trigger for computational arrhythmias is an extra stimulus. But as pointed out by Keener [14] it is far from obvious how such an extra stimulus can be generated in the heart. There are many possible explanations of this phenomenon; one is the existence of ectopic waves triggering activity from other places than the sinoatrial node. The cells of the sinoatrial node are special in the sense that if they were allowed to act as a single cell, they would create a periodic solution (see Keener and Sneyd [2]). For that reason they are commonly referred to as pacemaker cells. The common cardiac cell is stable in the sense that nothing happens if the cell is not stimulated. Therefore, the resting state of a normal cardiac cell is inherently stable in contrast to a pacemaker cell. It has, however, been observed that pacemaker cells also exist in other parts of the cardiac muscle; see e.g. Khan [15]. In particular, it is known that pacemaker cells can be found around the pulmonary veins (see e.g. [16]), and that the amount and electrical coupling of such cells are essential for the stability of collections of pacemaker cells (see e.g. [17]). It is therefore important to be able to analyze a state where pacemaker cells are surrounded by stable (normal) cells and see how the stability of such a blend depends on the amount of pacemaker cells and the electrical coupling (diffusion). Mathematically, this problem is closely related to the problem of estimating the critical size of the sinoatrial node (see Keener and Sneyd [2]). The resting state for a single cell was carefully analyzed by Jacquemet [18] and an efficient method for computing the resting state in tissue was derived by Jacquemet and Henriquez [19]. Based on mathematical and computational arguments, we have derived a condition for setting off ectopic waves in models of excitable cells where the dynamics is governed by a version of the Luo-Rudy I model; see [20]. In that paper we studied the problem consisting of a set of pacemaker cells surrounded by stable cells in 1D and on the unit square. The simplified geometry and mathematical models allowed us to do all the computations in Matlab [21]. It is observed in [20] that for simplified models, the stability of pacemaker cells surrounded by stable cells can be assessed by considering only a collection of pacemaker cells equipped with Dirichlet boundary conditions. This approximation is analyzed in the paper [22] in the case of fairly realistic dynamics defined on simple geometries. Our aim is now to pursue this route of reasoning to the more interesting case of ectopic activity in the vicinity of the pulmonary veins. In order to do so, we need a much more powerful computational tool.
It is the purpose of the present paper to present a computational framework suitable for analyzing the stability of models of pacemaker cells surrounded by stable cells. For mathematical models of the form (1), we will discuss how to compute an equilibrium solution of the system. For a given equilibrium solution, we will derive a linearized model and study the stability of the linear model by computing the largest real part of the eigenvalues of an associated discrete problem. We want to have a versatile tool allowing us to easily change mathematical models and geometry, and therefore we extensively apply automatic tools for handling the complicated models involved. Another possible method of investigating the stability would be to compute a stationary solution of the problem and then assesses the stability computationally by simply perturbing the steady-state solution and see if a wave is initiated. Albeit simple, this approach has a serious limitation since it can only determine whether the steady-state solution is stable with respect to a specific perturbation. Another perturbation may yield another conclusion. By using the eigenmode approach we know precisely in what manner a steady-state solution is stable or unstable.
The paper is organized as follows: In the next section we identify a number of software components needed to compute the eigenvalue with the largest real part of linearized versions of (1). In Section 3 we briefly describe the implementation of these software components, and in Section 4 we present experiments displaying the properties and applicability of the software system.
Section snippets
Methods
The purpose of the present paper is to present a method for analyzing the stability of equilibrium solutions of systems written on the form (1). An equilibrium solution of (1) is a solution to the following system: defined on the domain , and where, again, the transmembrane potential v satisfies a homogeneous Neumann type boundary condition.
Implementation
We have implemented a software module that enables analysis of systems of ordinary differential equations on the form The software, based on Swiginac [25], utilizes a symbolic representation of the system, and it allows us to automatically identify equations of the form (2). More specifically, the main motivation for utilizing symbolic representation of the system is the possibility to detect linear and semi-diagonal terms. In the single cell case the knowledge of such terms is used to
Results
The method described above is implemented in terms of a series of software modules that individually are quite complex and thus the combination of them may carry many errors. Testing is thus very important in order to apply the software for problems of interest from a physiological point of view. We have tested the software for two systems of reaction-diffusion equations. The first test case is based on a special version of the FitzHugh–Nagumo model (see [22]). For that model we derive an
Conclusion
We have developed a method for computing and analyzing a steady-state solution of a system of partial differential equations modeling a blend of stable and unstable cardiac cells. The accuracy of the computed results has been tested for a modified FizhHugh–Nagumo model where eigenvalues can be computed analytically. Furthermore, the system has been applied to a simplified model of the outlet of a pulmonary vein where the inner part consists of unstable cells and the outer part consists of
Conflict of interest statement
None declared.
Acknowledgement
This research was supported by a Center of Excellence grant from the Research Council of Norway to the Center for Biomedical Computing at Simula Research Laboratory.
References (34)
- et al.
A note on a method for determining advantageous properties of an anti-arrhythmic drug based on a mathematical model of cardiac cells
Mathematical Biosciences
(2009) - et al.
A simple two-variable model of cardiac excitation
Chaos, Solitons & Fractals
(1996) - et al.
Multiple wavelets, rotors, and snakes in atrial fibrillation—a computer simulation study
Journal of Electrocardiology
(2007) - et al.
Genesis of ectopic waves: role of coupling, automaticity and heterogeneity
Biophysical Journal
(2005) Steady-state solutions in mathematical models of atrial cell electrophysiology and their stability
Mathematical Biosciences
(2007)- et al.
A condition for setting off ectopic waves in computational models of excitable cells
Mathematical Biosciences
(2008) Impulses and physiological states in theoretical models of nerve membrane
Biophysical Journal
(1961)- et al.
From pulsus to pulseless: the saga of cardiac alternans
Circulation Research
(2006) - et al.
Mathematical Physiology
(1998) - et al.
A quantitative description of membrane current and its application to conduction and excitation in nerve
Journal of Physiology
(1952)
Rate dependence and regulation of action potential and calcium transient in a canine cardiac ventricular cell model
Circulation Research
When Time Breaks Down—The Three Dimensional Dynamics of Electrochemical Waves of Cardiac Arrhythmia
Vulnerability to re-entry in simulated two-dimensional cardiac tissue: effects of electrical restitution and stimulation sequence
Chaos
Effects of simulated ischemia on spiral wave stability
American Journal of Physiology - Heart and Circulatory Physiology
Mechanism underlying initiation of paroxysmal atrial flutter/atrial fibrillation by ectopic foci: a simulation study
Circulation
Organization of ventricular fibrillation in the human heart
Circulation Research
Cited by (1)
Instabilities of the resting state in a mathematical model of calcium handling in cardiac myocytes
2012, Mathematical Biosciences