PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of transaThe Royal Society PublishingPhilosophical Transactions AAboutBrowse by SubjectAlertsFree Trial
 
Philos Trans A Math Phys Eng Sci. 2008 September 28; 366(1879): 3381–3409.
Published online 2008 July 4. doi:  10.1098/rsta.2008.0112
PMCID: PMC2778066
NIHMSID: NIHMS111934

From mitochondrial ion channels to arrhythmias in the heart: computational techniques to bridge the spatio-temporal scales

Abstract

Computer simulations of electrical behaviour in the whole ventricles have become commonplace during the last few years. The goals of this article are (i) to review the techniques that are currently employed to model cardiac electrical activity in the heart, discussing the strengths and weaknesses of the various approaches, and (ii) to implement a novel modelling approach, based on physiological reasoning, that lifts some of the restrictions imposed by current state-of-the-art ionic models. To illustrate the latter approach, the present study uses a recently developed ionic model of the ventricular myocyte that incorporates an excitation–contraction coupling and mitochondrial energetics model. A paradigm to bridge the vastly disparate spatial and temporal scales, from subcellular processes to the entire organ, and from sub-microseconds to minutes, is presented. Achieving sufficient computational efficiency is the key to success in the quest to develop multiscale realistic models that are expected to lead to better understanding of the mechanisms of arrhythmia induction following failure at the organelle level, and ultimately to the development of novel therapeutic applications.

Keywords: multiscale modelling, cardiac electrical activity, monodomain, bidomain, ordinary differential equation integration

1. Introduction

Computer simulations of electrical behaviour in the whole ventricles (Xie et al. 2004; Rodriguez et al. 2005; Potse et al. 2006; Ten Tusscher et al. 2007; Ashihara et al. 2008) or atria (Harrild & Henriquez 2000; Vigmond et al. 2001; Virag et al. 2002; Seemann et al. 2006) have become commonplace during the last few years. However, even when powerful state-of-the-art computational resources are used, attempts to integrate behaviour from the protein scale of ion channels to the organ scale of cardiac arrhythmias remain enormously challenging, and typically include significant trade-offs between representations of different types of complexities (subcellular processes versus structural complexities). To arrive at a computational efficiency in modelling the (patho)physiological processes in the heart, such that it permits the exploration of the parameter space of interest, the simplifications typically made are as follows.

  1. The geometry of the organ is represented in a stylized fashion (‘one heart geometry fits all’ approach; Xie et al. 2004; Rodriguez et al. 2005; Ashihara et al. 2008); only parts of the heart are modelled, such as slices across the ventricles (Meunier et al. 2002; Trayanova & Eason 2002; Hillebrenner et al. 2004) or wedges (Burton et al. 2006; Plank et al. 2008); and idealized geometries, such as myocardial slabs (Vigmond & Leon 1999; Cherry et al. 2003; Plank et al. 2005), sheets (Beaumont et al. 1998; Skouibine et al. 1999; Anderson & Trayanova 2001; Samie et al. 2001; Kneller et al. 2002; Weiss et al. 2005) or strands (Thomas et al. 2003; Qu et al. 2006), are used.
  2. In constraining the degrees of freedom, the choice of the computational mesh discretization often leads to (i) under-representation of (to the degree of fully ignoring) the finer details of the cardiac anatomy, such as endocardial trabeculations or papillary muscles, and (ii) the necessity to adjust ad hoc the tissue conductivity tensors in order to avoid the artificial scaling of the wavelength, thus compensating for the dependence of conduction velocities on grid granularity. That is, as the grid is coarsened with all other model parameters remaining unchanged, conduction velocity becomes reduced and thus the wavelength is diminished, with conduction block occurring above a certain spatial discretization limit.
  3. The myocardial mass is treated as a homogeneous continuum, without representing intramyocardial discontinuities such as vascularization, cleavage planes, or patches of fat or collagen.
  4. The specialized cardiac conduction system, i.e. the sinuatrial (SA) and atrioventricular (AV) nodes and the Purkinje network, is typically not represented in the whole-organ simulations, although a few exceptions exist (Berenfeld & Jalife 1998; Vigmond & Clements 2007; Ten Tusscher & Panfilov 2008).
  5. Myocardial membrane ion transport kinetics are modelled in a simplified fashion (Rogers & McCulloch 1994; Ten Tusscher & Panfilov 2006a). Reduced models preserve salient features such as excitability, refractoriness, electrical restitution, etc.; these are usually represented phenomenologically at the scale of the cell (so that they can be easily manipulated). Such approaches have led to important insights into the mechanisms by which action potential characteristics control the stability of electrical propagation (Weiss et al. 2006). Their limitation is, however, that phenomenologically represented parameters do not directly correspond to actual molecular structures or processes, and are thus incapable of accounting for many potentially arrhythmogenic mechanisms, such as propagation instabilities induced by instabilities in calcium cycling or by the altered metabolic state of the cell.

Clearly, the ultimate goal of modelling is to accurately represent the interplay between the subsystems responsible for the primary functions of the heart, including the electrophysiological, Ca2+ handling, contractile and energetic components. This necessitates a significant level of complexity of the cell models, often achieved by implementing highly nonlinear control systems. Moreover, the ability to use simulations to gain a better understanding of cardiac pathophysiology requires a representation of the bidirectional feedback loops connecting the various subsystems. This presents additional computational concerns with respect to the coupling of processes that span very different temporal and spatial scales, described by equations with varying numerical stiffness.

A case in point is a recently described cell model, which combines mitochondrial bioenergetics with previously developed electrophysiological, Ca2+ handling and contractile models, referred to as the excitation–contraction coupling and mitochondrial energetics (ECME) model (Cortassa et al. 2006). This cell model was developed as a framework for studying how the failure of the mitochondrial network of the cardiomyocyte can lead to action potential shortening or complete electrical inexcitability as a result of the depletion of the high-energy phosphate pool and the activation of the ATP-sensitive K+ channels in the sarcolemma. Through a mechanism of reactive oxygen species (ROS)-induced ROS release (Cortassa et al. 2004), it was shown that the uncoupling of oxidative phosphorylation can induce oscillations of background K+ currents (Aon et al. 2003) that were proposed to contribute to electrical heterogeneity and arrhythmias in whole hearts exposed to ischaemia-reperfusion (Akar et al. 2005). Compounds that could inhibit mitochondrial depolarization also prevented the post-ischaemic arrhythmias, demonstrating that a failure at the level of the mitochondrion scales to cause macroscopic desynchronization of cardiac electrical propagation. Thus, there is a strong motivation for developing multiscale models to study excitation–contraction–bioenergetic coupling, spanning from the organelle to the entire organ. Because the ECME model takes into account both fast Markovian ion channel gating kinetics (microsecond relaxation times) and slow enzymatic reactions associated with the tricarboxylic acid cycle (second-long relaxation times), an ordinary differential equation (ODE) solver, which explicitly accounts for the wide range of time-scales involved, would be ideally suited to speed up computational time. The technical challenges of scaling up a model with a temporal range spanning six orders of magnitude are formidable.

The goals of this article are (i) to review the techniques that are currently employed to model cardiac electrical activity in the heart, discussing the strengths and weaknesses of the various approaches, and (ii) to implement a novel modelling approach, based on physiological reasoning rather than on mathematical rigour, that lifts some of the restrictions imposed by ionic models of the most current generation, such as the ECME model. This article thus presents an efficient multiscale model of ventricular arrhythmogenesis that incorporates mitochondrial ionic channels and mitochondrial energetics. Achieving sufficient computational efficiency is the key to success in the quest to develop multiscale realistic models that are expected to lead to better understanding of the mechanisms of arrhythmia induction following failure at the organelle level and, ultimately, to the development of novel therapeutic applications.

2. Review of tissue- and organ-level modelling techniques

(a) The bidomain equations

Besides a few exceptions (Wang et al. 1996), the cardiac bidomain equations, or a simplified version referred to as the monodomain equations, are typically used to model cardiac bioelectric phenomena at the tissue and organ level. The bidomain equations (Plonsey 1988) represent the homogenization of cardiac tissue, replacing discrete components of the intracellular and extracellular tissue matrix, such as cells and gap junctions, with a continuum. Each point in space is considered to exist in two domains, the intra- and extracellular; the domains are superimposed in space and separated from each other by the cellular membrane.

The following equations relate the intracellular potential ϕi to the extracellular potential ϕe, through the transmembrane current density Im:

equation image
(2.1)
equation image
(2.2)
equation image
(2.3)
equation image
(2.4)

where An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ25.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ26.jpg are the intracellular and extracellular conductivity tensors, respectively; β is the surface-to-volume ratio of the cardiac cells; Itrans is the current density of the transmembrane stimulus; Cm is the membrane capacitance per unit area; Vm is the transmembrane voltage, defined as ϕiϕe; and Iion is the current density of the ionic current, which depends on the transmembrane voltage and other state variables represented by υ.

The bidomain conductivity tensors An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ27.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ28.jpg are the result of a homogenization procedure (Henriquez 1993), through which the discrete nature of the tissue matrix is translated into a macroscopic representation of the conductive tissue properties. The conductivity tensors account, in an averaged macroscopic sense, for the observation that conduction velocity is faster along the direction of the fibre orientation, and slower in a direction transverse to it. In the intracellular space, this is a reflection of the discrete coupling network, where cell-to-cell connections (gap junctions) are found most frequently at the ends of a cell (along the long axis of the cell). Furthermore, in both intracellular and extracellular spaces, a preferred direction of conduction is a consequence of the cell geometry, which constrains current flow in both spaces. This property is referred to as anisotropy. The anisotropy ratios between the two domains have been shown to be unequal (Clerc 1976; Roberts et al. 1979; Roberts & Scher 1982; Roth 1997). This property has been demonstrated to be of major importance in the interaction between cardiac tissue and externally applied current (Sepulveda et al. 1989).

Different approaches to recast the bidomain equations have been suggested (Hooke et al. 1994). A particularly popular linear transformation is to add equations (2.1) and (2.2) and to use the definition of Vm according to equation (2.4) to retain ϕe and Vm, the two experimentally measurable quantities of interest, as the independent variables (Pollard et al. 1992):

equation image
(2.5)
equation image
(2.6)

In the most general case, when the heart is immersed in a conductive fluid (e.g. a tissue bath in an experimental context or a surrounding torso to model in vivo scenarios), a Poisson problem has to be additionally solved,

equation image
(2.7)

where σb is the conductivity of the bath medium; and Ie is the volume density of the extracellular stimulus current, which is injected or withdrawn through electrodes located in the bath.

Boundary conditions at the tissue–bath interface enforce continuity of the normal component of the extracellular current and continuity of ϕe, whereas the intracellular domain is isolated there (the normal component of the intracellular current is zero). At the boundaries of the bath, the normal component of the extracellular current is zero. Without enforcing Dirichlet boundary conditions, the elliptic partial differential equation (PDE) components, i.e. equation (2.7), and the elliptic portion of the bidomain equations, i.e. equation (2.5), are singular, which may pose numerical difficulties. Typically, grounding electrodes are used in the extracellular domain, i.e. nodes in the mesh are chosen where ϕe is fixed at zero value, to circumvent this problem. For most scenarios of practical interest this is not a limitation since the validation of simulation results with experimental data requires the use of grounding electrodes to match an experimental set-up. With certain iterative solver techniques, for instance Krylov subspace methods such as conjugate gradients (CG), this is not necessarily required (Potse et al. 2006).

(b) Discretization schemes and issues

Different spatial discretization techniques, such as the finite-difference method (FDM; Skouibine et al. 2000), the finite-volume method (FVM; Harrild & Henriquez 1997; Trew et al. 2005), the interconnected cable model (Leon & Roberge 1990), and the finite-element method (FEM; Rogers & McCulloch 1994; Vigmond et al. 2002; Rodriguez et al. 2005; Ashihara et al. 2008), have been applied to the cardiac bidomain problem. Spatial discretization is mainly governed by the spatial extent of the depolarization wavefront. Assuming that the upstroke of the action potential lasts for approximately 1 ms and that under physiological conditions the wavefront travels at a speed between 0.2 (transverse conduction) and 0.7 m s−1 (longitudinal conduction), the spatial extent of the wavefront is in the range of 0.2–0.7 mm. The grid resolution has to be smaller than the extent of the wavefront, otherwise spatial undersampling effects will influence the simulation results. In the following sections, ‘fine discretization’ will be used to refer to discretizations where the spatial resolution is well below the spatial extent of the wavefront, and ‘coarse’ if the resolution is well above it. A discretization step of 250 μm could be considered as the limit between the two; however, it is not strict. Under pathological conditions, for instance, when a very slow decremental conduction is modelled, even a resolution of 100 μm may be too coarse.

A more formal argument to estimate spatio-temporal discretization constraints is found by considering the Courant–Friedrichs–Lewy (CFL) condition (Courant et al. 1928). Briefly, assuming straight fibres and equal anisotropy ratio, the CFL condition, which governs the maximum possible time-step for stable integration using the forward Euler method, is approximated by

equation image
(2.8)

where An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ33.jpg with An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ34.jpg. Assuming a standard value of 1400 cm−1 for β and the conductivity values as published by Clerc (1976), the CFL condition predicts a maximum possible time-step of 45.5 μs for a 100 μm grid. For a 10 μm grid Δt decreases to 0.45 μs, which renders the forward Euler integration for the parabolic PDE computationally inefficient compared with an implicit scheme such as the Crank–Nicholson scheme discussed below.

Numerically, the bidomain equations can be solved as a coupled system (Vigmond et al. 2002), where transmembrane voltage and extracellular potential are solved for simultaneously. Alternatively, operator-splitting techniques are applied (Keener & Bogar 1998) to decouple the computing scheme into three components, such as an elliptic PDE, a parabolic PDE and a set of nonlinear ODEs (equations (2.9)–(2.12)). It has been shown that the decoupled scheme converges quickly against the coupled scheme by employing a block Gauss–Seidel iteration (Pennacchio & Simoncini 2002). However, in a number of studies, the components were essentially treated as independent (Weber dos Santos et al. 2004a,b; Sundnes et al. 2005; Austin et al. 2006; Plank et al. 2007); solutions were then found by leapfrogging between the decoupled components, where either Vm in equation (2.5) or ϕe in equation (2.6) are considered as constant. In Vigmond et al. (2002), it was found that, with small error tolerance, the difference between coupled and decoupled approaches is negligible.

Discretizing the decoupled bidomain equations leads to a three-step scheme, which involves the solution of the parabolic PDE, the elliptic PDE and the nonlinear system of ODEs at each time-step. In the simplest case, both the parabolic PDE and the nonlinear ODE systems can be solved with an explicit forward Euler scheme (Vigmond et al. 2002),

equation image
(2.9)
equation image
(2.10)
equation image
(2.11)
equation image
(2.12)

where An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ39.jpg is the discretized An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ40.jpg operator, with ξ being either i (intracellular) or e (extracellular); Δt is the time-step; and Vk, An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ41.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ42.jpg are the temporal discretizations of Vm, ϕe and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ43.jpg, respectively, at the time instant of kΔt.

When the computational mesh is of fine discretization and thus the integration time-step of the parabolic PDE has to be reduced substantially to maintain stability, it is advantageous to employ the more expensive semi-implicit Crank–Nicholson scheme instead of the forward Euler scheme to solve the parabolic PDE, with

equation image
(2.13)

replacing equations (2.9) and (2.10).

(i) Linear solver strategies

The PDEs are solved most efficiently with direct methods; however, such methods are limited to small grids only (Vigmond et al. 2002; Plank et al. 2007) because otherwise memory demands increase quickly, which, in turn, significantly increases the required number of operations per solver step. Although direct methods have been implemented to run in parallel environments (Amestoy et al. 2001; Li & Demmel 2003), typically they are harder to parallelize due to the fine-grained parallelism, which is communication intense. For large systems, iterative methods are mandatory because the memory requirements increase as An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ45.jpg and the complexity, in terms of required operations, increases as An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ46.jpg, whereas the requirements for iterative solvers are An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ47.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ48.jpg, respectively, where N is the system size and ϵ is the error tolerance (Meurant 1999). In general, linear systems with more than a few hundreds of thousands of unknowns can be considered as large. For instance, in Plank et al. (2007) a rabbit ventricular geometry immersed in a bath was discretized, resulting in 0.87 million unknowns. The linear system associated with this set-up could not be factorized on a desktop computer with 8 GB of RAM available.

When executing bidomain simulations on sequential computers, the main computational burden can be attributed to the solution of the elliptic problem and the set of ODEs (i.e. the ionic model). The complexity of an ionic model is determined by the number of state variables in the model and the stiffness of the ODE system. With simple ionic models, among which are phenomenological models such as FitzHugh–Nagumo (FitzHugh 1969) and its derivatives as well as electrophysiology-based (Hodgkin–Huxley-type) ionic models with a small number of state variables and without or with limited intracellular compartments and fluxes between them (Drouhard & Roberge 1987), the elliptic problem contributes more than 90 per cent to the overall workload. By contrast, for recent ionic models involving a large number (20–60) of state variables and very stiff ODEs as well as Markov state models (Iyer et al. 2004; Saucerman et al. 2004; Cortassa et al. 2006; Flaim et al. 2007; Mahajan et al. 2008), the ODE solution may even dominate the computations.

The parabolic problem is typically less of a concern. On coarser meshes, simple forward Euler steps can be employed to update Vm. In this case, the contributions of the diffusional component (i.e. the PDE) and the local membrane component to changes in Vm can be updated separately, which renders the PDE linear. On finer grids, semi-implicit Crank–Nicholson schemes perform well. Even when relatively cheap iterative solvers are employed, the parabolic portion can be updated efficiently due to the diagonal dominance of the linear system.

For large systems, parallel computing approaches are necessary to reduce execution times. The parallel computing context alleviates the problem of solving the set of ODEs. State variables in an ionic model do not diffuse, which qualifies the ODEs as an embarrassingly parallel problem. No communication between processors is required to update the state variable and thus the parallel scaling of the ODE portion is linear. The parabolic problem is efficiently solved in parallel as well. Either only a forward Euler step is required (essentially a matrix–vector product, for which good scalability is expected), or the well-posed diagonally dominant linear system is solved efficiently with relatively cheap iterative methods, such as preconditioned CG. Typically, with an incomplete LU (ILU) preconditioner for the iterative CG solver, the parabolic problem can be solved in less than 10 iterations.

The elliptic PDE is the most challenging problem. Standard iterative solvers such as ILU–CG typically require several hundreds of iterations to converge (Plank et al. 2007), which makes this solution significantly more expensive than that of the parabolic system, although both systems share the same sparsity pattern. The parallel scaling of standard iterative solvers is fairly good (Plank et al. 2007); for instance, a parallel ILU–CG solver, where the system is decomposed by a block Jacobi preconditioner with ILU(0), i.e. an incomplete LU factorization with zero fill-in levels that preserves the sparsity pattern of the original matrix, used as a sub-block preconditioner, exhibits good parallel scaling. An example of good parallel scaling is presented in Plank et al. (2007), where parallel efficiency of approximately 92 per cent was achieved over the test range of processor numbers (up to 64 processors). With fewer number of processors, ILU(N) with N levels of fill-in tends to be more efficient; however, with an increasing number of processors, the efficiency of the preconditioning deteriorates since the preconditioner is applied only to the main diagonal block. This deterioration can be avoided by employing overlapping block preconditioners, such as alternating Schwarz methods; however, this increases the burden on the network interconnect owing to the increase in communication caused by the increased number of off-processor elements involved (Toselli & Widlund 2005). If the hardware lacks low-latency interconnects, this is undesirable since parallel efficiency levels off quickly with increasing number of processors.

It has been demonstrated in several recent studies (Weber dos Santos et al. 2004a,b; Austin et al. 2006; Plank et al. 2007) that multilevel preconditioners for CG methods both significantly improve the overall performance and exhibit reasonable parallel efficiency (better than 80%) for up to 128 processors. A generally applicable algebraic multigrid preconditioner (AMG) in conjunction with an iterative Krylov solver reduces the number of iterations per solver step by almost two orders of magnitude compared with ILU–CG (Plank et al. 2007). Although a single iteration with AMG is significantly more expensive than with ILU, the reduction in the number of iterations clearly favours a multilevel approach. For instance, in Plank et al. (2007) a speedup of 6 was reported. Using AMG–CG constitutes to date the most efficient method for solving the elliptic portion of the bidomain equations. Besides its computational efficiency, a major advantage of the AMG approach is that it is applied straightforwardly to unstructured grids, avoiding the complex task of explicitly generating nested representations of the geometry at different levels of resolution, as required by the geometric versions of multilevel preconditioners. Unstructured grids allow smooth representations of the organ boundaries, which is of particular importance in simulations of cardiac defibrillation, where geometric discretizations with jagged boundaries inevitably evoke artificial surface polarizations. A more extensive review of numerical techniques to solve the cardiac bidomain equations can be found in Vigmond et al. (2008).

(c) Adaptivity and domain decomposition techniques

Substantial gains in computational efficiency can be achieved by explicitly taking into account the fact that high spatial and temporal resolution is needed only in the immediate vicinity of a propagating wavefront, where fast transitions in time are translated into steep gradients in space. Typically, steep wavefronts encompass a small portion of the overall domain, although this portion increases substantially when re-entrant activation patterns, particularly fibrillatory activity, are being modelled. Spatially adaptive methods that dynamically adjust grid resolution hold promise in increasing computational efficiency. Cherry et al. (2003) reported significant performance gains with an adaptive mesh refinement technique; however, the method, although very well suited for FDM and explicit stencils, is not easily modifiable to accommodate FEM or FVM discretizations and implicit stencils. A new class of methods loosely referred to as mortar FEM techniques (Pennacchio 2004) will possibly allow the use of FEM discretizations with spatial adaptivity; however, implementations using unstructured grids have not been reported yet.

Alternatively, simpler domain composition methods can be employed, where the region of interest is divided into several subdomains; each subdomain is then integrated at a different rate, the latter being a function of the ongoing activity in the domain itself and in the adjacent domains. However, such techniques introduce a significant bookkeeping overhead. Whether the reduction in floating point operations compensates for these extra costs, and to what degree, remains inconclusive. For instance, in Quan et al. (1998) a speedup of 17 was reported, whereas in Vigmond & Leon (1999) a speedup of only 3 could be achieved. Using different temporal resolutions on a per-domain base inevitably leads to load balancing issues for codes executed in parallel. That is, even though in a particular domain computations are executed quickly, the slowest domain dominates the overall execution time. A possible but non-trivial solution to address this problem is to migrate vertices between parallel partitions, although to our knowledge there are no simulator codes available today that are capable of achieving this in a robust manner.

3. Review of ionic model computation techniques

(a) Integration in standard form

The dynamics of the processes that underlie the action potential in a cardiac cell are typically described by a set of ODEs. These are initial-value problems of the form

equation image
(3.1)

where y is a vector and f is a vector-valued function. The solution to the above equation is typically found as

equation image
(3.2)

Different approaches can be used to approximate the integral over the nonlinear vector function f. Here, all the components of the vector y are updated simultaneously with all the components of An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ51.jpg integrated simultaneously. General-purpose packages incorporating a large set of standard integrators are currently available and allow a convenient integration in a vector form (Cohen & Hindmarsh 1996). Standard integration techniques include explicit (or forward) and implicit (or backward) methods. Explicit methods are popular since they are easy to implement; however, the order of this class of methods is one, which often results in insufficient accuracy. Approaches to overcome this weakness include the use of several previously computed solutions (multistep methods) or additional intermediate solutions in the interval An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ52.jpg (Runge–Kutta methods) to update the current solution. The more sophisticated implicit backward methods such as backward differentiation formula (BDF) or implicit Runge–Kutta methods (e.g. Rosenbrock methods; Rosenbrock 1963) have superior stability properties and allow larger time-steps; however, they are computationally expensive and, in general, robust implementation of these methods is difficult. The use of BDF methods requires the solution of a linear system of equations iteratively, by either fixed-point methods or versions of the Newton–Raphson methods. In this case, Jacobian matrices have to be either analytically determined and repetitively evaluated, or numerically approximated (Hairer & Wanner 2004). Many currently used integrators incorporate advanced features such as variable time-stepping and error control, where a time-step is chosen such that the local error per step is below a prescribed tolerance level.

For the particular purpose of integrating ODEs representing the myocardial action potential, standard methods have not been established and many different techniques are currently employed. Overall, although it is generally believed that implicit methods, owing to the stiffness of membrane kinetic variables, are advantageous and should be preferred over forward methods, the vast majority of implementations rely on forward methods.

(b) Component-wise integration

Although the use of standard integrators is convenient for single-cell action potential simulations, it becomes cumbersome when ionic model integrators are coupled with PDE solvers for calculation of the potentials at the tissue and organ level. Therefore, as described in Maclachlan et al. (2007), in the field of cardiac modelling, it is common to split the vector formulation in equation (3.1) into a set of equations, where each ODE is integrated separately. In this component form, a system of coupled ODEs with N state variables decouples into N ODEs of the form

equation image
(3.3)

where An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ54.jpg and fi refers to the component i of the vector function An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ55.jpg (Maclachlan et al. 2007). Any standard method discussed above is well suited for component-wise integration as well. The application of implicit BDF methods is easier with component-wise integration since only scalar Newton iterations are required.

Possibly the most popular component-wise integration approach was proposed by Rush & Larsen (1978), and it has been widely applied (Vigmond & Leon 1999; Vigmond et al. 2002; Ten Tusscher & Panfilov 2006b) ever since. When the set of ODEs is decoupled and integrated sequentially in a component form, each state variable yi is updated, with all other state variables yj (ji) held constant. Many of the ODEs comprising an ionic model, typically all gating equations in Hodgkin–Huxley-type models, but also ODEs in Markov state formulations, can be written in the form

equation image
(3.4)

where all variables belonging to the subvector An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ57.jpg are considered constant during an integration step. Hence, the nonlinear ODE reduces to a linear ODE with constant coefficients, for which an analytic solution can be found. With A and B being constant, the analytical solution is

equation image
(3.5)

This analytical solution is then used to compute the solution at the next time-step, An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ59.jpg. Numerical analysis has revealed that the Rush–Larsen method offers significant stability and accuracy benefits over the forward Euler method (Maclachlan et al. 2007). Any portions fi of the vector function An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ60.jpg that cannot be written in a linear form need to be integrated in the fully nonlinear form. In the original Rush–Larsen formulation, the forward Euler method was applied to solve the fully nonlinear form. Therefore, the Rush–Larsen method has been classified as a non-standard FDM with forward Euler (NSFD w/FE; Maclachlan et al. 2007). Several suggestions have been made to further improve the method and overcome its weaknesses (Maclachlan et al. 2007). First, stability issues due to the use of the forward Euler method for the nonlinear components could be addressed by the use of implicit methods. Second, the accuracy of the method, originally first order, could be improved to second by a more expensive scheme involving the computation of midpoints. Although it has been shown that such elaborate approaches yield additional benefits for some ionic models, for instance the Courtemanche model of the human atrial cell (Courtemanche et al. 1998), where stiffness properties are less critical, no performance benefits or even worse performance has been observed with models that have pronounced stiffness properties, such as the canine ventricular ionic model by Winslow et al. (1999).

(c) Acceleration techniques

Regardless of the particular method, the expense in ionic model integration is in the evaluation of An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ61.jpg, requiring the calculation of hundreds of computationally expensive expressions (mostly exponential and logarithmic functions). Many of these expressions depend on a single variable or they can be recast into products of terms that depend on a single variable only. The physiological range of these variables is typically well known a priori, thus it is convenient to precompute and tabulate the expensive terms and then simply look up their values when required. For instance, assuming that in equation (3.5) the terms A/B and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ62.jpg depend only on one state variable yj and that no adaptive time-stepping scheme is employed, i.e. Δt is constant throughout the simulation, then the functions An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ63.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ64.jpg can be precomputed and used for table lookup (TLU). As shown in Vigmond et al. (2003), recasting equation (3.5) leads to a more efficient TLU scheme,

equation image
(3.6)

where one yj-dependent lookup operation is required to determine the values of An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ66.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ67.jpg. The solution of the ODE can then be updated with one multiplication and one sum operation only (see algorithm 3.1).

NSFD w/FE integration of the linear part of the vector function f using TLU.

  • for k=0, …, TEndt do
    • for n=0, …, Ncells do
      • 1=TLU(yj);
      • yi[n]=yi[n]*1[0]+1[1];
    • end for
  • end for

The TLU, in its most efficient form, requires one multiplication and one typecast operation only (nearest-neighbour TLU) or two multiplications and two sum operations (TLU with interpolation). Furthermore, if data structures are chosen properly for lookup tables and state variables, excellent cache performance can be achieved with this method.

(d) Tissue-level aspects of ionic model computation

On a current standard desktop computer, single-cell ionic model simulations, no matter how complex the particular ionic model is, can be executed faster than real time and memory requirements are negligible. With regard to organ-level simulations, however, real time is orders of magnitude too slow for a satisfactory performance. For instance, assuming that a human heart is sufficiently finely discretized by 10 million vertices, with real-time performance the simulation of one second of activity would last for approximately 116 days. Using supercomputers with tens or hundreds of CPUs scales the simulation time down to a few hours. However, considering that this is only the ODE part of a tissue-level simulation and that only a time interval of one second is simulated, it becomes apparent that, even on supercomputers, real-time performance is much too slow to allow meaningful computational studies, where a large parameter space is explored and where simulations of several seconds to minutes of electrophysiological activity are desired.

Clearly, performance plays a key role in tissue/organ simulations whereas it is of minor importance for single-cell simulations. Therefore, ionic model implementations used in tissue-/organ-level simulations are typically not straightforward extensions of single-cell codes. Ionic state variables in tissue models become vectors, memory demands increase, and the flexibility needed to allow implementation of spatial heterogeneities in model parameters or the use of different ionic models with different parameter settings in one organ model requires a careful software engineering design.

At the tissue/organ level, there are specific issues that prevent concepts successfully applied at the single-cell level from being carried over to tissue simulations without major adjustments. These issues can be summarized as follows.

  1. Organ-level simulations are almost exclusively carried out in parallel computing environments in order to achieve reasonable performance, whereas parallelization is pointless for single-cell codes. Since state variables do not diffuse, the solution of the system of ODEs qualifies as an embarrassingly parallel problem. That is, parallelization is easily achieved by simply partitioning the problem into evenly sized chunks, with no communication required. In tissue-level codes, the ODE solver is coupled to the solution of the PDEs, adding complexity. In the most complex case, where a bidomain formulation including a bath is discretized on an unstructured grid and where globally defined quantities other than Vm (for instance, strain or interstitial concentrations) are linked to sets of different ODE systems representing different types of tissues (SA and AV node, Purkinje network, ventricular tissue, with spatial heterogeneities in transmural and apicobasal directions), the partitioning of the ODE system and its interaction with the global quantities becomes a truly challenging problem. Addressing parallelization and domain decomposition issues will become even more important in the near future because any standard desktop computer will be essentially a parallel computer with several computational cores. Approaches that do not address these issues will not benefit from future generations of computers and thus will be of very limited use.
  2. In single-cell simulations, computational efficiency can be improved by allowing adaptivity in time-stepping, with dt being a function of yi or dyi. In tissue-level simulations this becomes problematic because ODE time-stepping has to be synchronized with the update intervals of the PDEs. Similar to the approach of solving bidomain equations in the decoupled form, where the elliptic portion is solved and hence the extracellular potential ϕe is updated at a slower rate than the parabolic portion updating Vm, the ODEs can be updated at a different rate than the PDEs. This is likely to work for small time-steps, with the ODE time-step chosen to be smaller than that of the PDE. Such an approach is advantageous if the stiffness associated with the fast-changing internal state variables enforces a time-step well below the time-step limit imposed by the mesh ratio of the parabolic portion of the bidomain equations. Furthermore, for time-adaptive ODE methods, the solver time depends largely on the current state of the ionic model. Therefore, solution time is larger in portions of the tissue grid where activation is ongoing, and shorter in portions where the tissue is quiescent or at plateau/repolarization. In tissue simulations, particularly when re-entrant phenomena are simulated, all phases of the action potential coexist at any time. Therefore, adaptive ODE time-stepping introduces a load imbalance in parallel codes, where partitions with the fastest processes and thus the smallest time-steps become the limiting factor preventing linear parallel scaling of the ODE system.
  3. For current computer architectures, cache performance is the key factor. Efforts to reduce the overall number of flops are not nearly as efficient as the careful layout of data structures, which increases the cache performance of the code. For single-cell simulators, code efficiency is of minor importance. No matter how complex an ionic model is, it will always fit into cache.
  4. Single-cell codes aim to increase the maximum time-step while preserving stability and accuracy. Requirements are, however, very different in tissue-/organ-level simulations. The maximum time-step that allows stable integration could be limited by the mesh ratio and not the ODEs. For instance, using the NSFD w/FE technique, the Courtemanche model of the human atrial action potential (Courtemanche et al. 1998) can be integrated with a time-step of 400 μs. However, when incorporated in a tissue model of a fine discretization, the time-step has to be reduced well below this number; the particular value depends on the method used to solve the parabolic PDE (Weber dos Santos et al. 2004a,b). That is, the potential of expensive ODE integration schemes, which allow large time-steps, cannot be fully exploited since PDE time-step constraints limit Δt to values for which cheaper ODE integrators perform well without any stability problems. On the other hand, it has been demonstrated that methods that increase the stability of an integration scheme allow for larger time-steps; however, they lead to increasingly inaccurate solutions with increasing Δt. For instance, in Whiteley (2006, 2007) tissue simulations were conducted with time-steps from 100 to 500 μs. Although the numerical scheme was stable, oscillations were observed immediately after the action potential upstroke (Whiteley (2006); figure 1c,d). Such oscillations may lead to spurious transients in internal state variables, which, in turn, will be reflected in Vm. For instance, owing to the tight coupling between Vm and intracellular calcium handling, oscillations in Vm at the end of an action potential upstroke would give rise to artificial calcium transients; these are purely numerical artefacts and are not related to any intrinsic properties of the model. We demonstrate this in figure 1. Hence, such large time-steps have to be used with great care. In Whiteley (2006), the largest time-step used in a fully implicit ODE solver that resulted in reasonably accurate solutions (a small artefactual ripple in Vm during early repolarization was visible, though) turned out to be 100 μs. The same simulation can be carried out with the NSFD w/FE method with the same time-step without any stability problems. In addition, a single NSFD w/FE step implemented efficiently could take only a fraction of the execution time required for an implicit BDF solver step, thus allowing one to execute tens to hundreds of NSFD w/FE solver steps within the time-frame required for a single BDF step. Thus, it can be concluded that, in most cases of practical interest, BDF methods are not likely to help in reducing the computational burden.
    Figure 1
    Stable integration schemes for the parabolic PDE allow larger time-steps, but, as illustrated in Whiteley (2006, 2007), may lead to oscillations in (a) transmembrane potential Vm. Such oscillations are numerical artefacts. Although they do not affect ...

4. A new modelling paradigm for complex ionic models based on temporal multiscale decoupling

(a) The need for a new paradigm

Although computationally efficient implementations of versions of the NSFD w/FE method continue to be an excellent choice for many recent ionic models (Mahajan et al. 2008), for other recent models such as the ECME ionic model (Cortassa et al. 2006), their application becomes impracticable. For instance, with the fastest extended NSFD method implementation to date (Maclachlan et al. 2007), it took 25.77 s to execute a simulation of over 500 ms of activity in a single cell using the Winslow et al. (1999) canine ventricular model (no specific detail regarding the implementation was provided). The same study showed that fully implicit methods indeed allow much larger time-steps, but the overall performance was not better and the accuracy was strikingly inferior to that of NSFD w/FE.

The reason why some recent models have imposed a much higher burden on the ODE integrators than others stems from the current trend to develop more realistic (and hence, predictive) mechanistic formulations. This has led to the incorporation of processes whose dynamics occur on a faster time-scale and require a large number of state variables to faithfully reproduce experimentally observed phenomena. This is the case for formulations where Hodgkin–Huxley-type ionic channel representations have been replaced by Markov process descriptions (Clancy & Rudy 1999). Taking this approach has led to the development of a series of integrative cardiac myocyte models that can explain the mechanisms underlying observations, such as action potential prolongation in heart failure (Winslow et al. 1999), the coexistence of high gain and graded sarcoplasmic reticulum Ca release (Greenstein & Winslow 2002), the role of beta-adrenergic stimulation on excitation–contraction coupling (Greenstein et al. 2004), and how properties of the cardiac dyad influence excitation–contraction coupling gain (Greenstein et al. 2006).

The wider range of time-scales captured in such ionic models increases the overall model stiffness primarily owing to the very fast time constants that govern intermediate states of these Markov state models. In a number of these ionic models (Jafri et al. 1998; Winslow et al. 1999; Cortassa et al. 2006; Greenstein et al. 2006), the fastest processes are linked to formulations describing the intrinsically fast dynamics of the L-type calcium channel and the ryanodine receptor, RyR, as well as the calcium dynamics in the dyadic space, i.e. the restricted space between the membrane of the T-tubules, where the L-type calcium channels reside predominantly, and the co-localized RyR in the membrane of the junctional sarcoplasmic reticulum. The positive feedback loop between the L-type calcium currents and the calcium release controlled by RyR is referred to as calcium-induced calcium release (CICR). Descriptions of CICR processes stiffen the ODE system considerably, as illustrated in figure 2. The time-scales governing CICR are often in the sub-microsecond range. Using the NFSD w/FE method enforces a time-step as small as 0.2 μs to successfully integrate these equations.

Figure 2
(a) (i) CICR-related currents and fluxes in the ECME model (blue, ICa; green, INa; red, Jrel; light green, Jxfer; violet, Jtr; dark yellow, Jup). Fluxes Jrel and Jxfer act at a substantially faster time-scale during CICR. The sodium current, INa, which ...

To avoid issues with ODE integration that are related to the stiffness of the system of equations resulting from detailed descriptions of CICR, a simplified phenomenological description could be used. Such descriptions employ a simple set of equations that reproduce the profile of a typical calcium-release signal with the time course of the release ‘hard coded’ into the equations. However, because CICR is the process that is at the very core of myocyte contractile function, and because understanding how alterations in CICR impact integrative function in both normal and diseased cells is a fundamental goal of building mathematical models, the use of phenomenological descriptions of CICR or any subcellular process hinders progress towards elucidating the integrative mechanisms at play in the cardiac myocyte. In order to be considered predictive, ionic models must be developed that accurately capture the underlying biophysical mechanisms of experimentally observed phenomena.

In the following sections, we present a novel ODE integration scheme that builds on the NSFD w/FE method but overcomes its limitations to improve computational efficiency and thus enable the use of ionic models such as the ECME model in tissue- and organ-level simulations. The needed gains in efficiency are obtained by taking advantage of two key observations.

  1. There is a wide range of time constants involved in the system of ODEs, ranging from the sub-microsecond range to hundreds of milliseconds (figure 2). Decoupling the ODE system into a suitable number of subsystems that are governed by similar time constants, and using different time-steps and/or different integrators for each subsystem could lead to substantial performance gains. We refer to this approach as temporal multiscale decoupling (TMSD).
  2. A salient feature of ionic models such as the ECME model is the extreme stiffness of the system of ODEs, which is, however, only present during very brief phases of the action potential. By making simplifications to reduce this stiffness, the integrator can be relieved from having to track this extremely fast dynamics. We refer to this technique as dynamic reformulation of the ODEs.

(b) Formalization of the TMSD approach

In this section, we present the general TMSD approach, which we later apply to the ECME model to provide a specific illustration of its implementation. The essence of the approach is that fast state variables are considered constant as long as the difference in time-scales with the rest of the system is sufficiently large. In applying the temporal multiscale decoupling approach, state variables that act at a similar time-scale are grouped together, and blockwise updates of the different groups are performed using different update frequencies without any impact on the overall solution. However, special care is taken of equations where fi depends on fluxes. In this case, fluxes are integrated and adjusted to a particular update frequency when used outside of the group in which they were computed. Often such groups can be mapped to subsystems of the cell, for instance the mitochondria in the ECME model (Cortassa et al. 2006). Temporal changes of state variables in this subsystem occur at time-scales that are at least three orders of magnitude slower than those of the RyR subsystem and thus infrequent updates of the state variables in this subsystem can be executed without any loss of information. Since updates of slow variables require the same expensive evaluations of a group of functions fn of the vector function f, not updating all state variables at the same rate (figure 2) leads to significant savings in terms of execution time.

The wide range of time-scales involved in these ionic models can be exploited in a systematic manner to arrive at more efficient integration schemes, where different components of the state vector are decoupled according to a temporal multiscale modelling approach. In the simplest case, the state vector y is split into M portions and each partial state vector ym is integrated either with the same time-step Δt, when fixed time-step integration schemes are employed, or if an adaptive integration method with a varying internal time-step is used, the integrator output is synchronized at a desired Δt. With respect to the ODE integration frameworks presented in the review above, each subvector can then be integrated either in the standard form or component-wise in the non-standard form. Such decoupling has several benefits.

  1. Each group An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ68.jpg is integrated with the method that is best suited. Stiff partial state vector systems with very fast components can be integrated with an expensive ODE solver, whereas partial state vectors with slow components can be integrated with much cheaper methods.
  2. When implicit backward integrators in the standard form are employed, the nonlinear systems that have to be solved are much smaller and thus easier to solve. For instance, solving the ODEs of the ECME model in the standard form involves the solution of a system with a sparse Jacobian of dimension 50×50. By decoupling and employing the method described here to the RyR subsystem, only a 4×4 system has to be solved. The evaluation of the Jacobian is much simpler in this case because evaluation of the portions of the vector function f that are not pertinent to the subsystem fm can be skipped.

A formal description of the TMSD method for three groups is provided in algorithm 4.1. Each group is integrated using a different time-step with every outer time-step being a multiple of the time-step of the innermost loop.

TMSD integration of a system of ODEs using time-steps Δt0, Δt1 and Δt2. State variables An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ69.jpg are considered constant. Primed and double-primed fluxes J′ and J″ are used to update state variables with a different time-step in an outer loop.

  • for k=0, …, TEndt2 do
    • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ70.jpg;
    • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ71.jpg;
    • t2=kΔt2;
    • for j=0 to Δt2t1 do
      • t1=t2+jΔt1;
      • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ72.jpg;
      • for i=0, …, Δt1t0 do
        • t0=t1+iΔt0;
        • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ73.jpg;
        • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ74.jpg;
      • end for
      • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ75.jpg;
      • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ76.jpg;
      • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ77.jpg;
    • end for
    • An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ78.jpg;
  • end for

(c) Dynamic reformulation of the ODEs

The main reason for the increased computational cost of the current state-of-the-art ionic models that use Markov chains is their stiffness and, to a much lesser extent, the large number of state variables. With Markov state models, the stiffness of the governing equations may vary substantially depending on the state vector y, since the Jacobian is not a constant but a function of y. That is, during a particular interval, a problem may be stiff, but it may be non-stiff during other intervals. For instance, with the ECME model, stiffness seems to be a problem only for a few microseconds during the onset and upstroke of the calcium transient in the dyad; other than that, the ECME model can be treated as other models are, with the minor difference that the number of state variables is two to three times higher than that of most standard ionic models. During the phases of pronounced stiffness, the integration time-step has to be reduced to values as small as 0.2 μs, even when expensive implicit ODE integrators are employed (Maclachlan et al. 2007).

In other words, this problem can be circumvented neither with more elaborate ODE integrators nor with the TMSD approach, since the inner loop of the TMSD integrator is tied to the same small time-step. Depending on the number of Markov state variables relative to the overall size of the set of ODEs, the inner loop may impose a computational burden that renders organ-level applications computationally intractable or at least extremely expensive. In such cases, we suggest that the equations themselves be modified to reduce the stiffness of the system without changing the solution. This is again achieved by taking advantage of the separation of time-scales using the rationale as follows.

  1. Processes that are sufficiently faster than the minimum time-scale of interest can be assumed to occur instantaneously. In this case, an ODE An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ79.jpg can be replaced by an algebraic equation by assuming that An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ80.jpg, which yields An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ81.jpg.
  2. In a Markov state model, rapid equilibrium assumptions can be made if the transition between two connected states is substantially faster than other transition rates in the Markov network. For instance, in the Markov state model of RyR in the ECME model (Jafri et al. 1998; Cortassa et al. 2006), for a transition between the closed state An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ82.jpg and the open state An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ83.jpg (figure 3), one can assume that An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ84.jpg, where An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ85.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ86.jpg are rate constants and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ87.jpg is the calcium concentration in the dyad, and use the combined state An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ88.jpg instead of the individual states An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ89.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ90.jpg. This approach has been used in a stochastic implementation of a Markov state RyR model (Greenstein & Winslow 2002), and in the development of a simplified (yet mechanistically detailed) model of CICR in the cardiac myocyte (Greenstein et al. 2006); however, unlike here, an integrator in standard form was used.
    Figure 3
    An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ9.jpg-dependent toggling of RyR subsystem models. Depending on a threshold calcium level An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ10.jpg in the dyad, different formulations to describe the RyR subsystem were used. For low An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ11.jpg concentrations with An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ12.jpg, when rate coefficients are slow, the full-blown formulation ...

(d) Implementation of a TMSD ODE integrator with dynamic reformulation of ODEs for the ECME model

As part of this study, the suggested TMSD ODE numerical integration technique with dynamic reformulation of ODEs was implemented for the ECME model and validated against the well-established integrator used in the original paper (Cortassa et al. 2006); the latter relies on an implicit BDF method provided through the CVODE library (Cohen & Hindmarsh 1996). Our validated TMSD code used three different groups, which we integrated at time-steps of 2.5, 20 and 100 μs. The smallest time-step was applied to integrate the equations related to An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ91.jpg and the Markov state models representing RyR and the L-type calcium current. Standard Hodgkin–Huxley-type gating variables and other ionic fluxes were integrated at a time-step of 20 μs using the NSFD w/FE integrator. The very slow processes linked to mitochondrial metabolism and force generation were integrated with the largest time-step of 100 μs.

Besides using the TMSD approach, dynamic reformulation was applied to the equations for RyR and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ92.jpg, but not to the Markov state model representing the L-type calcium current. Details of the procedure are illustrated in figure 3. Formulations for the RyR subsystem were switched based on the transition rates of the An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ93.jpg level. For An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ94.jpg below a threshold level, An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ95.jpg, the full-blown Markov state model of the RyR (formulation A, figure 3a) was employed and the ODE that governs the calcium flux balance in this compartment was integrated to update An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ96.jpg. Since there are no An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ97.jpg transients, large time-steps could be implemented to speed up execution time. Fast transients in the RyR states are triggered by the quick rise in An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ98.jpg. As soon as An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ99.jpg rose above the threshold level An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ100.jpg, rapid equilibrium assumptions were made, and a simplified model was employed to update the RyR states (formulation B, figure 3b). The combined state An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ101.jpg substituted the individual states An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ102.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ103.jpg, where the probability of An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ104.jpg relative to the new combined state was determined by An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ105.jpg with α being determined from An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ106.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ107.jpg, i.e. An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ108.jpg. The threshold value for switching between the formulations was chosen in such a way that the rate coefficients governing the transitions between An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ109.jpg and An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ110.jpg were at least 10 times faster than all other processes. Specifically, we used a threshold An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ111.jpg=100 m s−1, which translated into a threshold concentration of An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ112.jpg=0.095 mM.

When switching to formulation B, further numerical stability benefits were gained by replacing, with an algebraic equation, the ODE that governs the An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ113.jpg transient as a function of the fluxes into and out of the dyad. This was accomplished by assuming that An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ114.jpg equilibrates rapidly, i.e. the Ca2+ influx into the dyadic space and the Ca2+ efflux out of it are equal. This is equivalent to assuming that An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ115.jpg holds true, resulting in the reduction of the ODE governing the transients into an algebraic equation. This assumption is well justified because Ca2+ diffusion out of the dyad occurs far more rapidly than the L-type Ca channel and RyR gating processes (Hinch et al. 2004; Greenstein et al. 2006). As evidenced by figure 3c, this modification does not lead to any differences in results (compare traces A and B in figure 3c), but allowed a much larger time-step, 2.5 μs. Otherwise, for stability reasons, the integrator was forced into a small time-step of 0.2 μs. However, for An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ116.jpg close to the base level, employing the algebraic equation to update An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ117.jpg results in small deviations. In this case, the onset of the An external file that holds a picture, illustration, etc.
Object name is rsta20080112equ118.jpg transient occurs approximately 5 μs earlier than that in the original model (compare traces A and B with B* in figure 3c).

Finally, longer single-cell simulations, over 1 min, were executed to ensure the long-term stability of the scheme. It was established, for all 50 state variables of the model, that the results were indiscernible from those obtained with the reference implementation.

5. Implementation of the new paradigm in organ-level simulations with the ECME model

The gain in execution speed obtained by using TMSD with the dynamic reformulation of several ODEs versus the classical NSFD w/FE scheme turned out to be at least two orders of magnitude. Although the exact speedup depended on several factors, such as platform, compiler, optimization and implementation details, differences in execution time between a straightforward NSFD w/FE implementation of the ECME model, using an integration time-step of 0.2 μs, and the TMSD integrator with dynamic reformulation reached a factor of up to 285. This is mainly attributed to gains in the maximum stable time-step size that can be used to integrate the RyR state equations, while less frequent updates of state variables related to the metabolic state and force generation helped to reduce execution times further by reducing the number of expensive evaluations of the function f.

To assess to what degree the gain in efficiency in the ECME model was sufficient to allow simulations at the tissue and organ levels, we carried out benchmark simulations in tissue consisting of one million cells. The benchmark was executed without coupling the cells, i.e. the parabolic PDE was not solved and no spatial domain was discretized. This problem size was chosen (i) by considering that one million degrees of freedom suffice to discretize any small mammalian heart up to the size of a rabbit at a sufficiently fine spatial resolution (approx. 200 μm) and (ii) to ensure that realistic organ-level simulation is tested where cache performance will play an important role.

Execution times to simulate one second of activity were measured on a standard desktop computer with four CPUs for the ECME model and several other widely used models, from a simple ventricular model with seven state variables, the MBRDR model (Drouhard & Roberge 1987), through the Courtemanche model of the human atrial action potential (Courtemanche et al. 1998) and the Luo–Rudy guinea-pig ventricular model in its latest revision (Faber & Rudy 2000; Luo & Rudy 2006), to the most recent Markov state model of the rabbit ventricular cell by Mahajan et al. (2008). For all models other than the ECME, we used the NSFD w/FE method with the acceleration techniques described in the review portion of this article. Although some of these models could be integrated with significantly larger time-steps in a single-cell mode, all models were integrated with a time-step of 20 μs, assuming that this value is a constraint imposed by the CFL conditions of the parabolic PDE.

Benchmark results are summarized in table 1. Even with all the optimizations implemented, the ECME model is, by a wide margin, computationally significantly more expensive than all the other models in this benchmark. Nonetheless, the performance is sufficient to investigate phenomena occurring on the time-scale of several seconds to minutes when simulations are executed on a supercomputer. Since the ODE computations scale linearly on a parallel computer, the time spent on solving the ODEs can be scaled down to execute the computations within a reasonable time. For instance, using 512 processors, the execution time will be reduced by a factor of 128 relative to this benchmark. That is, in 2.34 min, 1 s of activity can be simulated; 1 min of activity would be executed in 140 min.

Table 1
Summary of execution times Texec, in minutes, to simulate 1 s of activity in a tissue model with one million cells running on a standard desktop computer for various ionic models. (A global time-step of 20 μs was used for all models.) ...

Finally, we coupled the ECME model with the implemented TMSD ODE integration scheme with dynamic reformulation to tissue-/organ-level bidomain codes to conduct true organ-level simulations. A guinea-pig-size ventricular geometry model, which represents the rabbit ventricular geometry used in our previous studies (Trayanova et al. 2002; Rodriguez et al. 2005; Ashihara et al. 2008), scaled half in size, was chosen to facilitate the execution of the test on a standard desktop computer with reasonable performance. The discretization of the myocardial volume at an average element size of approximately 200 μm led to approximately 75e3 degrees of freedom. The guinea-pig ventricles were paced at the apex to initiate a simple activation sequence (figure 4). Using a current desktop computer (Dual Opteron X2), the simulation of 1 ms of activity was executed in approximately 1.2 s for a monodomain formulation and in approximately 5 s for a bidomain formulation. That is, on a standard desktop computer in 1 hour execution time, approximately 3 s of activity could be simulated with a monodomain formulation and 0.72 s with a bidomain formulation. These numbers clearly demonstrate the potential of the new methods presented in this study. Again, when parallel computers are employed to execute simulations, much longer observation periods, of the order of tens of seconds to minutes, can be studied with ionic models as complex as the ECME model within a reasonable time.

Figure 4
Application of the TMSD ODE integrator with dynamic reformulation of ODEs to the ECME model in an organ-level simulation. An apical pacing sequence was simulated in a monodomain guinea-pig-size ventricular model (approx. 75e3 degrees of freedom) using ...

6. Concluding remarks on the new paradigm

This article has presented a review of various numerical techniques to model electrical activity in the heart, spanning the spatial scales from subcellular mechanisms to the entire organ. Furthermore, a new version of the NSFD w/FE method was developed, which takes advantage of the very different time-scales of the processes represented in recent ionic models of very high complexity, such as the ECME model, to obtain substantial increases in computational efficiency. In Maclachlan et al. (2007), it was reported that the simulation of 500 ms activity in a single cell modelled by the Winslow et al. (1999) canine ventricular model lasted for 25 s. The ECME model used in this study employed the same set of equations to represent the L-type calcium current and RyR subsystem, which are by far the computationally most expensive components of the model. Although direct comparison is difficult, the ECME model can be judged as computationally more expensive than the canine model by Winslow et al. The performance we obtained in this study for the ECME model translates into 36 ms to simulate 500 ms activity in a single cell and, thus, is significantly more efficient than the standard NSFD w/FE method.

This is the first attempt at an implementation based on the TMSD ODE integration method with dynamic reformulation of ODEs. Although the performance is satisfactory and opens new avenues to employ highly detailed cell models in organ-level simulations, this is a first step only and further significant performance improvements can be expected. The choice of decomposing the set of ODEs into groups and the maximum possible time-step for each group cannot be determined automatically and needs manual tuning. In general, this is a significant disadvantage of the TMSD technique since intimate knowledge of the model is required to tune up the method for best performance. Decoupling an ionic model into different subsystems requires a non-trivial effort and careful examination of the maximum rates of change of different state variables. The effort to develop ionic models based on the suggested methods is further complicated by the fact that at least two implementations of an ionic model are required, a very accurate reference implementation such as an implicit BDF integrator in the standard form and the TMSD implementation itself. At the current stage of development, the method is sufficiently efficient to allow parameter sweeps in organ modelling studies of small mammalian hearts, even when high-performance computing (HPC) environments are not available. To simulate larger mammalian hearts or for questions that require the use of the bidomain formulation, the use of HPC hardware is mandatory. Today, with the availability of HPC environments, there are virtually no limitations in using highly detailed Markov chain ionic models, and the testing of various simulation scenarios as well as explorations of parameter spaces indeed appear feasible.

Acknowledgments

This work was supported by Marie Curie Fellowship MC-OIF 040190 funded by the European Commission and SFB grant F3210-N18 from the Austrian Science Fund (FWF) to G.P.; NIH grants R37-HL54598, R33-HL87345 and P01-HL081427 to B.O'R.; NIH grants R33 HL87345, N01-HV28180, PO1 HL081427 and PO1-HL077180 to R.L.W.; and NIH grants R01-HL063195, R01-HL082729 and R01-HL067322 and NSF grant CBET-0601935 to N.A.T.

Footnotes

One contribution of 11 to a Theme Issue ‘The virtual physiological human: building a framework for computational biomedicine II’.

References

  • Akar F.G, Aon M.A, Tomaselli G.F, O'Rourke B. 2005. The mitochondrial origin of postischemic arrhythmias. J. Clin. Invest. 115, 3527–3535doi:10.1172/JCI25371 [PMC free article] [PubMed]
  • Amestoy, P., Duff, I. S., L'Excellent, J. Y. & Koster, J. 2001 Mumps: a general purpose distributed memory sparse solver. In Para '00: Proc. 5th Int. Workshop on Applied Parallel Computing, New Paradigms for HPC in Industry and Academia, pp. 121–130. London, UK: Springer.
  • Anderson C, Trayanova N.A. 2001. Success and failure of biphasic shocks: results of bidomain simulations. Math. Biosci. 174, 91–109doi:10.1016/S0025-5564(01)00076-1 [PubMed]
  • Aon M.A, Cortassa S, Marban E, O'Rourke B. 2003. Synchronized whole cell oscillations in mitochondrial metabolism triggered by a local release of reactive oxygen species in cardiac myocytes. J. Biol. Chem. 278, 44 735–44 744doi:10.1074/jbc.M302673200 [PubMed]
  • Ashihara T, Constantino J, Trayanova N.A. 2008. Tunnel propagation of postshock activations as a hypothesis for fibrillation induction and isoelectric window. Circ. Res. 102, 737–745doi:10.1161/CIRCRESAHA.107.168112 [PMC free article] [PubMed]
  • Austin T.M, Trew M.L, Pullan A.J. 2006. Solving the cardiac bidomain equations for discontinuous conductivities. IEEE Trans. Biomed. Eng. 53, 1265–1272doi:10.1109/TBME.2006.873750 [PubMed]
  • Beaumont J, Davidenko N, Davidenko J.M, Jalife J. 1998. Spiral waves in two-dimensional models of ventricular muscle: formation of a stationary core. Biophys. J. 75, 1–14 [PubMed]
  • Berenfeld O, Jalife J. 1998. Purkinje-muscle reentry as a mechanism of polymorphic ventricular arrhythmias in a 3-dimensional model of the ventricles. Circ. Res. 82, 1063–1077 [PubMed]
  • Burton R.A, et al. 2006. Three-dimensional models of individual cardiac histoanatomy: tools and challenges. Ann. NY Acad. Sci. 1080, 301–319doi:10.1196/annals.1380.023 [PMC free article] [PubMed]
  • Cherry E.M, Greenside H.S, Henriquez C.S. 2003. Efficient simulation of three-dimensional anisotropic cardiac tissue using an adaptive mesh refinement method. Chaos. 13, 853–865doi:10.1063/1.1594685 [PubMed]
  • Clancy C.E, Rudy Y. 1999. Linking a genetic defect to its cellular phenotype in a cardiac arrhythmia. Nature. 400, 566–569doi:10.1038/23034 [PubMed]
  • Clerc L. 1976. Directional differences of impulse spread in trabecular muscle from mammalian heart. J. Physiol. 255, 335–346 [PubMed]
  • Cohen S, Hindmarsh A. 1996. CVODE, a stiff/nonstiff ODE solver in C. Comput. Phys. 10, 138–143
  • Cortassa S, Aon M.A, Winslow R.L, O'Rourke B. 2004. A mitochondrial oscillator dependent on reactive oxygen species. Biophys. J. 87, 2060–2073doi:10.1529/biophysj.104.041749 [PubMed]
  • Cortassa S, Aon M.A, O'Rourke B, Jacques R, Tseng H.J, Marban E, Winslow R.L. 2006. A computational model integrating electrophysiology, contraction, and mitochondrial bioenergetics in the ventricular myocyte. Biophys. J. 91, 1564–1589doi:10.1529/biophysj.105.076174 [PubMed]
  • Courant R, Friedrichs K, Lewy H. 1928. Über die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann. 100, 32–74doi:10.1007/BF01448839
  • Courtemanche M, Ramirez R.J, Nattel S. 1998. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am. J. Physiol. 275, H301–H321 [PubMed]
  • Drouhard J.P, Roberge F.A. 1987. Revised formulation of the Hodgkin–Huxley representation of the sodium current in cardiac cells. Comput. Biomed. Res. 20, 333–350doi:10.1016/0010-4809(87)90048-6 [PubMed]
  • Faber G.M, Rudy Y. 2000. Action potential and contractility changes in [Na+]i overloaded cardiac myocytes: a simulation study. Biophys. J. 78, 2392–2404 [PubMed]
  • FitzHugh R.A. 1969. Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466 [PubMed]
  • Flaim S, Giles W, McCulloch A. 2007. Arrhythmogenic consequences of Na+ channel mutations in the transmurally heterogeneous mammalian left ventricle: an analysis of the I1768V SCN5A mutation. Heart Rhythm. 4, 768–778doi:10.1016/j.hrthm.2007.02.009 [PubMed]
  • Greenstein J.L, Winslow R.L. 2002. An integrative model of the cardiac ventricular myocyte incorporating local control of Ca2+ release. Biophys. J. 83, 2918–2945 [PubMed]
  • Greenstein J.L, Tanskanen A.J, Winslow R.L. 2004. Modeling the actions of beta-adrenergic signaling on excitation–contraction coupling processes. Ann. NY Acad. Sci. 1015, 16–27doi:10.1196/annals.1302.002 [PMC free article] [PubMed]
  • Greenstein J.L, Hinch R, Winslow R.L. 2006. Mechanisms of excitation–contraction coupling in an integrative model of the cardiac ventricular myocyte. Biophys. J. 90, 77–91doi:10.1529/biophysj.105.065169 [PubMed]
  • Hairer, E. & Wanner, G. 2004 Solving ordinary differential equations II: stiff and differential-algebraic problems, 2nd edn. Springer Series in Computational Mathematics, vol. 14. Berlin, Germany: Springer.
  • Harrild D.M, Henriquez C.S. 1997. A finite volume model of cardiac propagation. Ann. Biomed. Eng. 25, 315–334doi:10.1007/BF02648046 [PubMed]
  • Harrild D.M, Henriquez C.S. 2000. A computer model of normal conduction in the human atria. Circ. Res. 87, E25–E36 [PubMed]
  • Henriquez C.S. 1993. Simulating the electrical behavior of cardiac tissue using the bidomain model. Crit. Rev. Biomed. Eng. 21, 1–77 [PubMed]
  • Hillebrenner M.G, Eason J.C, Trayanova N.A. 2004. Mechanistic inquiry into decrease in probability of defibrillation success with increase in complexity of preshock reentrant activity. Am. J. Physiol. Heart. Circ. Physiol. 286, H909–H917doi:10.1152/ajpheart.00492.2003 [PubMed]
  • Hinch R, Greenstein J.L, Tanskanen A.J, Xu L, Winslow R.L. 2004. A simplified local control model of calcium-induced calcium release in cardiac ventricular myocytes. Biophys. J. 87, 3723–3736doi:10.1529/biophysj.104.049973 [PubMed]
  • Hooke N, Henriquez C.S, Lanzkron P, Rose D. 1994. Linear algebraic transformations of the bidomain equations: implications for numerical methods. Math. Biosci. 120, 127–145doi:10.1016/0025-5564(94)90049-3 [PubMed]
  • Iyer V, Mazhari R, Winslow R.L. 2004. A computational model of the human left-ventricular epicardial myocyte. Biophys. J. 87, 1507–1525doi:10.1529/biophysj.104.043299 [PubMed]
  • Jafri M.S, Rice J.J, Winslow R.L. 1998. Cardiac Ca2+ dynamics: the roles of ryanodine receptor adaptation and sarcoplasmic reticulum load. Biophys. J. 74, 1149–1168 [PubMed]
  • Keener J, Bogar K. 1998. A numerical method for the solution of the bidomain equations in cardiac tissue. Chaos. 8, 234–241doi:10.1063/1.166300 [PubMed]
  • Kneller J, Zou R, Vigmond E.J, Wang Z, Leon L.J, Nattel S. 2002. Cholinergic atrial fibrillation in a computer model of a two-dimensional sheet of canine atrial cells with realistic ionic properties. Circ. Res. 90, E73–E87doi:10.1161/01.RES.0000019783.88094.BA [PubMed]
  • Leon L.J, Roberge F.A. 1990. A new cable model formulation based on Green's theorem. Ann. Biomed. Eng. 18, 1–17doi:10.1007/BF02368414 [PubMed]
  • Li X, Demmel J. 2003. SuperLU_DIST: a scalable distributed-memory sparse direct solver for unsymmetric linear systems. ACM Trans. Math Software (TOMS). 29, 110–140doi:10.1145/779359.779361
  • Luo C.H, Rudy Y. 2006. A model of the ventricular cardiac action potential depolarization, repolarization, and their interaction. Circ. Res. 68, 1501–1526 [PubMed]
  • Maclachlan M.C, Sundnes J, Spiteri R.J. 2007. A comparison of non-standard solvers for ODEs describing cellular reactions in the heart. Comput. Methods Biomech. Biomed. Eng. 10, 317–326doi:10.1080/10255840701259301 [PubMed]
  • Mahajan A, et al. 2008. A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophys. J. 94, 392–410doi:10.1529/biophysj.106.98160 [PubMed]
  • Meunier J.M, Eason J.C, Trayanova N.A. 2002. Termination of reentry by a long-lasting AC shock in a slice of canine heart: a computational study. J. Cardiovasc. Electrophysiol. 13, 1253–1261doi:10.1046/j.1540-8167.2002.01253.x [PubMed]
  • Meurant, G. 1999 Computer solution of large linear systems Studies in Mathematics and its Applications, vol. 28. Amsterdam, The Netherlands: North Holland.
  • Pennacchio M. 2004. The mortar finite element method for the cardiac ‘bidomain’ model of extracellular potential. J. Sci. Comp. 20, 191–210doi:10.1023/B:JOMP.0000008720.85771.d0
  • Pennacchio M, Simoncini V. 2002. Efficient algebraic solution of reaction–diffusion systems for the cardiac excitation process. J. Comp. Appl. Math. 145, 49–70doi:10.1016/S0377-0427(01)00535-0
  • Plank G, Leon L.J, Kimber S, Vigmond E.J. 2005. Defibrillation depends on conductivity fluctuations and the degree of disorganization in reentry patterns. J. Cardiovasc. Electrophysiol. 16, 205–216doi:10.1046/j.1540-8167.2005.40140.x [PubMed]
  • Plank G, Liebmann M, Weber dos Santos R, Vigmond E.J, Haase G. 2007. Algebraic multigrid preconditioner for the cardiac bidomain model. IEEE Trans. Biomed. Eng. 54, 585–596doi:10.1109/TBME.2006.889181 [PubMed]
  • Plank G, Prassl A, Hofer E, Trayanova N.A. 2008. Evaluating intramural virtual electrodes in the myocardial wedge preparation: simulations of experimental conditions. Biophys. J. 94, 1904–1915doi:10.1529/biophysj.107.121343 [PubMed]
  • Plonsey R. 1988. Bioelectric sources arising in excitable fibers (ALZA lecture). Ann. Biomed. Eng. 16, 519–546doi:10.1007/BF02368014 [PubMed]
  • Pollard A.E, Hooke N, Henriquez C.S. 1992. Cardiac propagation simulation. Crit. Rev. Biomed. Eng. 20, 171–210 [PubMed]
  • Potse M, Dube B, Richer J, Vinet A, Gulrajani R.M. 2006. A comparison of monodomain and bidomain reaction–diffusion models for action potential propagation in the human heart. IEEE Trans. Biomed. Eng. 53, 2425–2435doi:10.1109/TBME.2006.880875 [PubMed]
  • Qu Z, Garfinkel A, Weiss J.N. 2006. Vulnerable window for conduction block in a one-dimensional cable of cardiac cells, 1: single extrasystoles. Biophys. J. 91, 793–804doi:10.1529/biophysj.106.080945 [PubMed]
  • Quan W, Evans S.J, Hastings H.M. 1998. Efficient integration of a realistic two-dimensional cardiac tissue model by domain decomposition. IEEE Trans. Biomed. Eng. 45, 372–385doi:10.1109/10.661162 [PubMed]
  • Roberts D.E, Scher A.M. 1982. Effect of tissue anisotropy on extracellular potential fields in canine myocardium in situ. Circ. Res. 50, 342–351 [PubMed]
  • Roberts D.E, Hersh L.T, Scher A.M. 1979. Influence of cardiac fiber orientation on wavefront voltage, conduction velocity, and tissue resistivity in the dog. Circ. Res. 44, 701–712 [PubMed]
  • Rodriguez B, Li L, Eason J.C, Efimov I.R, Trayanova N.A. 2005. Differences between left and right ventricular chamber geometry affect cardiac vulnerability to electric shocks. Circ. Res. 97, 168–175doi:10.1161/01.RES.0000174429.00987.17 [PMC free article] [PubMed]
  • Rogers J.M, McCulloch A.D. 1994. A collocation–Galerkin finite element model of cardiac action potential propagation. IEEE Trans. Biomed. Eng. 41, 743–757doi:10.1109/10.310090 [PubMed]
  • Rosenbrock H.H. 1963. Some general implicit processes for the numerical solution of differential equations. Comput. J. 5, 329–331doi:10.1093/comjnl/5.4.329
  • Roth B.J. 1997. Electrical conductivity values used with the bidomain model of cardiac tissue. IEEE Trans. Biomed. Eng. 44, 326–328doi:10.1109/10.563303 [PubMed]
  • Rush S, Larsen H. 1978. A practical algorithm for solving dynamic membrane equations. IEEE Trans. Biomed. Eng. 25, 389–392doi:10.1109/TBME.1978.326270 [PubMed]
  • Samie F.H, Berenfeld O, Anumonwo J, Mironov S.F, Udassi S, Beaumont J, Taffet S, Pertsov A.M, Jalife J. 2001. Rectification of the background potassium current: a determinant of rotor dynamics in ventricular fibrillation. Circ. Res. 89, 1216–1223doi:10.1161/hh2401.100818 [PubMed]
  • Saucerman J.J, Healy S.N, Belik M.E, Puglisi J.L, McCulloch A.D. 2004. Proarrhythmic consequences of a KCNQ1 AKAP-binding domain mutation: computational models of whole cells and heterogeneous tissue. Circ. Res. 95, 1216–1224doi:10.1161/01.RES.0000150055.06226.4e [PubMed]
  • Seemann G, Höper C, Sachse F.B, Dössel O, Holden A.V, Zhang H. 2006. Heterogeneous three-dimensional anatomical and electrophysiological model of human atria. Phil. Trans. R. Soc. A. 364, 1465–1481doi:10.1098/rsta.2006.1781 [PubMed]
  • Sepulveda N.G, Roth B.J, Wikswo J.P., Jr 1989. Current injection into a two-dimensional anisotropic bidomain. Biophys. J. 55, 987–999 [PubMed]
  • Skouibine K.B, Trayanova N.A, Moore P.K. 1999. Anode/cathode make and break phenomena in a model of defibrillation. IEEE Trans. Biomed. Eng. 46, 769–777doi:10.1109/10.771186 [PubMed]
  • Skouibine K.B, Trayanova N.A, Moore P.K. 2000. A numerically efficient model for simulation of defibrillation in an active bidomain sheet of myocardium. Math. Biosci. 166, 85–100doi:10.1016/S0025-5564(00)00019-5 [PubMed]
  • Sundnes J, Lines G.T, Tveito A. 2005. An operator splitting method for solving the bidomain equations coupled to a volume conductor model for the torso. Math. Biosci. 194, 233–248doi:10.1016/j.mbs.2005.01.001 [PubMed]
  • Ten Tusscher K.H, Panfilov A.V. Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions. Phys. Med. Biol. 51, 2006a. 6141–6156doi:10.1088/0031-9155/51/23/014 [PubMed]
  • Ten Tusscher K.H, Panfilov A.V. Alternans and spiral breakup in a human ventricular tissue model. Am. J. Physiol. Heart Circ. Physiol. 291, 2006b. H1088–H1100doi:10.1152/ajpheart.00109.2006 [PubMed]
  • Ten Tusscher K.H, Panfilov A.V. 2008. Modelling of the ventricular conduction system. Prog. Biophys. Mol. Biol. 96, 152–170doi:10.1016/j.pbiomolbio.2007.07.026 [PubMed]
  • Ten Tusscher K.H, Hren R, Panfilov A.V. 2007. Organization of ventricular fibrillation in the human heart. Circ. Res. 100, e87–e101doi:10.1161/CIRCRESAHA.107.150730 [PubMed]
  • Thomas S.P, Kucera J.P, Bircher-Lehmann L, Rudy Y, Saffitz J.E, Kleber A.G. 2003. Impulse propagation in synthetic strands of neonatal cardiac myocytes with genetically reduced levels of connexin43. Circ. Res. 92, 1209–1216doi:10.1161/01.RES.0000074916.41221.EA [PMC free article] [PubMed]
  • Toselli, A. & Widlund, O. 2005 Domain decomposition methods—algorithms and theory Springer Series in Computational Mathematics, vol. 34. Berlin, Germany: Springer.
  • Trayanova N, Eason J. 2002. Shock-induced arrhythmogenesis in the myocardium. Chaos. 12, 962–972doi:10.1063/1.1483955 [PubMed]
  • Trayanova N, Eason J, Aguel F. 2002. Computer simulations of cardiac defibrillation: a look inside the heart. Comput. Visual. Sci. 4, 259–270doi:10.1007/s00791-002-0082-8
  • Trew M, Le Grice I, Smaill B, Pullan A. 2005. A finite volume method for modeling discontinuous electrical activation in cardiac tissue. Ann. Biomed. Eng. 33, 590–602doi:10.1007/s10439-005-1434-6 [PubMed]
  • Vigmond E.J, Clements C. 2007. Construction of a computer model to investigate sawtooth effects in the Purkinje system. IEEE Trans. Biomed. Eng. 54, 389–399doi:10.1109/TBME.2006.888817 [PubMed]
  • Vigmond E.J, Leon L.J. 1999. Computationally efficient model for simulating electrical activity in cardiac tissue with fiber rotation. Ann. Biomed. Eng. 27, 160–170doi:10.1114/1.160 [PubMed]
  • Vigmond E.J, Ruckdeschel R, Trayanova N.A. 2001. Reentry in a morphologically realistic atrial model. J. Cardiovasc. Electrophysiol. 12, 1046–1054doi:10.1046/j.1540-8167.2001.01046.x [PubMed]
  • Vigmond E, Aguel F, Trayanova N. 2002. Computational techniques for solving the bidomain equations in three dimensions. IEEE Trans. Biomed. Eng. 49, 1260–1269doi:10.1109/TBME.2002.804597 [PubMed]
  • Vigmond E.J, Hughes M, Plank G, Leon L.J. 2003. Computational tools for modeling electrical activity in cardiac tissue. J. Electrocardiol. 36, Suppl.69–74doi:10.1016/j.jelectrocard.2003.09.017 [PubMed]
  • Vigmond E.J, Weber Dos Santos R, Prassl A.J, Deo M, Plank G. 2008. Solvers for the cardiac bidomain equations. Prog. Biophys. Mol. Biol. 96, 3–18doi:10.1016/j.pbiomolbio.2007.07.012 [PMC free article] [PubMed]
  • Virag N, Jacquemet V, Henriquez C.S, Zozor S, Blanc O, Vesin J.M, Pruvot E, Kappenberger L. 2002. Study of atrial arrhythmias in a computer model based on magnetic resonance images of human atria. Chaos. 12, 754–763doi:10.1063/1.1483935 [PubMed]
  • Wang S, Leon L.J, Roberge F.A. 1996. Interactions between adjacent fibers in a cardiac muscle bundle. Ann. Biomed. Eng. 24, 662–674doi:10.1007/BF02684179 [PubMed]
  • Weber dos Santos R, Plank G, Bauer S, Vigmond E. Parallel multigrid preconditioner for the cardiac bidomain model. IEEE Trans. Biomed. Eng. 51, 2004a. 1960–1968doi:10.1109/TBME.2004.834275 [PubMed]
  • Weber dos Santos, R., Plank, G., Bauer, S. & Vigmond, E. 2004b Preconditioning techniques for the bidomain equations Lecture Notes in Computational Science and Engineering, no. 40, pp. 571–580. Berlin, Germany: Springer.
  • Weiss J.N, Qu Z, Chen P.S, Lin S.F, Karagueuzian H.S, Hayashi H, Garfinkel A, Karma A. 2005. The dynamics of cardiac fibrillation. Circulation. 112, 1232–1240doi:10.1161/CIRCULATIONAHA.104.529545 [PubMed]
  • Weiss J.N, Karma A, Shiferaw Y, Chen P.S, Garfinkel A, Qu Z. 2006. From pulsus to pulseless: the saga of cardiac alternans. Circ. Res. 98, 1244–1253doi:10.1161/01.RES.0000224540.97431.f0 [PubMed]
  • Whiteley J.P. 2006. An efficient numerical technique for the solution of the monodomain and bidomain equations. IEEE Trans. Biomed. Eng. 53, 2139–2147doi:10.1109/TBME.2006.879425 [PubMed]
  • Whiteley J.P. 2007. Physiology driven adaptivity for the numerical solution of the bidomain equations. Ann. Biomed. Eng. 35, 1510–1520doi:10.1007/s10439-007-9337-3 [PubMed]
  • Winslow R.L, Rice J, Jafri S, Marban E, O'Rourke B. 1999. Mechanisms of altered excitation–contraction coupling in canine tachycardia-induced heart failure, II: model studies. Circ. Res. 84, 571–586 [PubMed]
  • Xie F, Qu Z, Yang J, Baher A, Weiss J.N, Garfinkel A. 2004. A simulation study of the effects of cardiac anatomy in ventricular fibrillation. J. Clin. Invest. 113, 686–693doi:10.1172/JCI17341 [PMC free article] [PubMed]

Articles from Philosophical transactions. Series A, Mathematical, physical, and engineering sciences are provided here courtesy of The Royal Society