Computing the stability of steady-state solutions of mathematical models of the electrical activity in the heart

https://doi.org/10.1016/j.compbiomed.2011.05.011Get rights and content

Abstract

Instabilities in the electro-chemical resting state of the heart can generate ectopic waves that in turn can initiate arrhythmias. We derive methods for computing the resting state for mathematical models of the electro-chemical process underpinning a heartbeat, and we estimate the stability of the resting state by invoking the largest real part of the eigenvalues of a linearized model. The implementation of the methods is described and a number of numerical experiments illustrate the feasibility of the methods. In particular, we test the methods for problems where we can compare the solutions with analytical results, and problems where we have solutions computed by independent software. The software is also tested for a fairly realistic 3D model.

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], vt=δΔvI(v,r,g),rt=R(v,r,g),gt=G(v,g).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 (g1,,gn), and r denotes various concentrations (r1,,rm). We assume that the function G has n components, and that it is semi-diagonal in the sense that each component takes the formgi,t=αi(v)(1gi)βi(v)gi,i=1,,n.The 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 vt=I(v,r,g),rt=R(v,r,g),gt=G(v,g).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: δΔvI(v,r,g)=0,R(v,r,g)=0,G(v,g)=0,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 u=F(u).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)

  • CellML Model Repository....
  • T.J. Hund et al.

    Rate dependence and regulation of action potential and calcium transient in a canine cardiac ventricular cell model

    Circulation Research

    (2004)
  • A.T. Winfree

    When Time Breaks Down—The Three Dimensional Dynamics of Electrochemical Waves of Cardiac Arrhythmia

    (1987)
  • D.X. Tran et al.

    Vulnerability to re-entry in simulated two-dimensional cardiac tissue: effects of electrical restitution and stimulation sequence

    Chaos

    (2007)
  • F. Xie et al.

    Effects of simulated ischemia on spiral wave stability

    American Journal of Physiology - Heart and Circulatory Physiology

    (2001)
  • Y. Gong et al.

    Mechanism underlying initiation of paroxysmal atrial flutter/atrial fibrillation by ectopic foci: a simulation study

    Circulation

    (2007)
  • H.W.J. Kirsten et al.

    Organization of ventricular fibrillation in the human heart

    Circulation Research

    (2007)
  • Cited by (1)

    View full text