Simulations of Complex and Microscopic Models of Cardiac Electrophysiology Powered by Multi-GPU Platforms

Abstract

Key aspects of cardiac electrophysiology, such as slow conduction, conduction block, and saltatory effects have been the research topic of many studies since they are strongly related to cardiac arrhythmia, reentry, fibrillation, or defibrillation. However, to reproduce these phenomena the numerical models need to use subcellular discretization for the solution of the PDEs and nonuniform, heterogeneous tissue electric conductivity. Due to the high computational costs of simulations that reproduce the fine microstructure of cardiac tissue, previous studies have considered tissue experiments of small or moderate sizes and used simple cardiac cell models. In this paper, we develop a cardiac electrophysiology model that captures the microstructure of cardiac tissue by using a very fine spatial discretization (8 μm) and uses a very modern and complex cell model based on Markov chains for the characterization of ion channel’s structure and dynamics. To cope with the computational challenges, the model was parallelized using a hybrid approach: cluster computing and GPGPUs (general-purpose computing on graphics processing units). Our parallel implementation of this model using a multi-GPU platform was able to reduce the execution times of the simulations from more than 6 days (on a single processor) to 21 minutes (on a small 8-node cluster equipped with 16 GPUs, i.e., 2 GPUs per node).

1. Introduction

Heart diseases are responsible for one third of all deaths worldwide [1]. Cardiac electrophysiology is the trigger to the mechanical deformation of the heart. Therefore, the knowledge of cardiac electrophysiology is essential to understand many aspects of cardiac physiological and physiopathological behavior [2]. Computer models of cardiac electrophysiology [3, 4] have become valuable tools for the study and comprehension of such complex phenomena, as they allow different information acquired from different physical scales and experiments to be combined to generate a better picture of the whole system functionality. Not surprisingly, the high complexity of the biophysical processes translates into complex mathematical and computational models. Modern cardiac models are described by nonlinear system of partial differential equations (PDEs) that may result in a problem with millions of unknowns.

Mathematical models for cell electrophysiology are a key component of cardiac modeling. They serve both as standalone research tools, to investigate the behavior of single cardiac myocytes, and as an essential component of tissue and organ simulation based on the so-called bidomain or monodomain models [4]. The cell models can be written as a general non-linear system of ordinary differential equations (ODEs) and may vary in complexity from simple phenomenological models [5] (based on two variables) to complex models describing a large number of detailed physiological processes [6] (based on 40 to 80 differential variables). Simple models focus on the genesis of action potential (AP), that propagates from cell to cell and generates an electric wave that propagates on the heart. Complex models account not only for the genesis of AP but also describe how this phenomenon is related to cardiac homeostasis and to different sub-cellular components, such as cell membrane’s ion channels. Advances in genetics, molecular biology, and electrophysiology experiments have provided new data and information related to the structure and function of ion channels. The Markov Chain (MC) model formalism has been increasingly used to describe both function and structure of ion channels. MC-based models have enabled simulations of structural abnormalities due to genetic diseases and drug-biding effects on ion channels [7–9]. Unfortunately, these modern cardiac myocyte models pose different challenges to both numerical methods, due to the stiffness of the ODEs introduced by MCs, and to high performance computing, due to the size of the problems, since the number of differential variables rises from a couple to near a hundred [10].

On the tissue level, the bidomain model [4] is considered to be the most complete description of the electrical activity. This nonlinear system of PDEs can be simplified to the so-called monodomain model, which may be less accurate but less computationally demanding than the bidomain model. Unfortunately, large scale simulations, such as those resulting from the discretization of an entire heart, remain a computational challenge. In addition, key aspects of cardiac electrophysiology, such as slow conduction, conduction block, and saltatory or sawtooth effects, demand sub-cellular discretization for the solution of the PDEs and nonuniform, heterogeneous tissue electric conductivity. These aspects of cardiac electrophysiology are strongly related to cardiac arrhythmia, reentry, fibrillation or defibrillation, and have been the research topic of many studies [11–20].

However, the demand of sub-cellular discretization for the solution of the PDEs and nonuniform, heterogeneous tissue electric conductivity have prevented the study of the aforementioned phenomena on large-scale tissue simulations. In addition, due to the high computational costs associated with the simulations of these microscopic models of cardiac tissue, previous works have adopted simple myocyte models, instead of modern MC-based models [6, 10].

In this work, we present a solution for this problem based on multi-GPU platforms (clusters equipped with graphics processing units) that allows fast simulations of microscopic tissue models combined with modern and complex myocyte models. The solution is based on merging two different high-performance techniques. We have previously investigated for cardiac modeling: cluster computing based on message passing communications (MPI) [21–24] and GPGPU (General-purpose computing on graphics processing units) [25–30]. We developed a two-dimensional model that is based on the previous work of Spach and collaborators [11, 17] that accounts for the microstructure of cardiac tissue, gap junction heterogeneous distribution, and discretizations of 8 μm. This microscopic tissue model was combined with the model of Bondarenko et al. [6] which is a modern and complex myocyte model based on MCs. Our parallel implementation of this model using a multi-GPU platform was able to reduce the execution times of the simulations from more than 6 days (on a single processor) to 21 minutes (on a small 8-node cluster equipped with 16 GPUs, that is, 2 GPUs per node). As a result, using this very fast parallel implementation we were able to simulate the formation of spiral waves, a form of self-sustained reentrant activity strongly associated with cardiac arrhythmia. To the best of our knowledge, this is the first time spiral waves are simulated using a cardiac model that accounts for both the microstructure of cardiac tissue and a modern and complex myocyte model.

2. Methods

2.1. Modeling Cardiac Microstructure

We developed a two-dimensional model that is based on the previous work of Spach and collaborators [11, 17] that accounts for the microstructure of cardiac tissue, gap junction heterogeneous distribution, and discretizations of 8 μm × 8 μm. A basic template for myocyte connections was developed and is presented in Figure 1. This basic unit accounts for the connection of a total of 32 cardiac myocytes with different shapes and numbers of neighboring cells. The mean and SD (standard deviation) values for cell length and width are 120.9±27.8 μm and 18.3±3.5 μm, respectively. These values are close to those reported in the literature: [31] (length = 140 μm and width = 19 μm), [32] (length = 134 μm and width = 18 μm), and [20] (length = 100 μm and width = 17.32 μm). On average, each cell connects to other 6 neighboring myocytes. Our two-dimensional model considers a homogeneous depth 𝑑=10 μm [11, 17].

Figure 1 

Basic unit of cardiac myocyte distribution based on a total of 32 cells. Cells are displayed in different alternating colors along the 𝑥-axis. The basic unit spans a total of 648 μm in the longitudinal direction versus 144 μm in the transversal direction.

This basic unit was created in such a way that it allows the generation of larger tissue preparations via the connections of multiple instances of it. Figure 2 presents how this can be achieved.

Figure 2 

Six basic units being combined to form a larger tissue.

Figure 3 presents an example of how the connections between different myocytes can be arranged. The code was developed in a flexible way, so that it allows the user to set up for each discretized volume Vol_𝑖,𝑗 (witharea=ℎ×ℎ) conductivity or conductance values for the north (𝜎_{𝑥𝑖,𝑗+1/2}), south (𝜎_{𝑥𝑖,𝑗−1/2}), west (𝜎_{𝑥𝑖−1/2,𝑗}), and east (𝜎_{𝑥𝑖+1/2,𝑗}) volume faces. These can be any nonnegative values. In this work, we set the discretization ℎ to 8 μm. In addition, based on the work of Spach and collaborators [11, 17], we chose only 5 possible types of connections between neighboring volumes that are membrane (𝜎_𝑚=0.0), cytoplasm (𝜎_𝑐=0.4𝜇S/𝜇m), gap junction plicate (𝐺_𝑝=0.5𝜇S), interplicate (𝐺_𝑖=0.33𝜇S), and combined plicate (𝐺_𝑐=0.062𝜇S), where we use 𝜎 for conductivity and 𝐺 for conductance. For the simulations presented in this work, the distribution of the different gap junctions within the 32 myocytes was not randomly generated. Instead, the gap junction distribution of the basic template unit was manually chosen to reproduce the distribution presented before in [11, 17]. With this setup and conductivity values we found that conduction velocity along the fibers was around 410 μm/ms (LP) and was 130 μm/ms transversal to fiber direction (TP). This results in a ratio LP/TP of 0.32, which is close to the conduction ratio reported in [11].

Figure 3 

In this work, there are only 5 possible types of connections between neighboring volumes that are membrane, which indicated no-flux between neighboring volumes; cytoplasm, which indicates that the neighboring volumes are within the same cell; three possible types of gap junctions, plicate, interplicate, and combined plicate. For the simulations presented in this work, the gap junction distribution of the basic template unit was manually chosen to reproduce the distribution presented before in [11, 17]. This figure presents an example of how different gap junctions are distributed in three neighboring myocytes that belong to the basic unit.

2.2. The Heterogeneous Monodomain Model

Action potentials propagate through the cardiac tissue because the intracellular space of cardiac cells is electrically coupled by gap junctions. In this work, we do not consider the effects of the extracellular matrix. Therefore, the phenomenon can be described mathematically by a reaction-diffusion type partial differential equation (PDE) called monodomain model, given by 𝛽𝐶_𝑚𝜕𝑉𝑥,𝑦,𝑡𝜕𝑡+𝛽𝐼_ion𝑉𝝈𝑥,𝑦,𝑡,𝜼𝑥,𝑦,𝑡=∇⋅𝑥,𝑦∇𝑉𝑥,𝑦,𝑡+𝐼_stim𝑥,𝑦,𝑡,𝜕𝜼𝑥,𝑦,𝑡𝑉,𝜕𝑡=𝐟𝑥,𝑦,𝑡,𝜼𝑥,𝑦,𝑡 (1) where 𝑉 is the variable of interest and represents the transmembrane potential, that is, the difference between intracellular to extracellular potential; 𝜼 is a vector of state variables that also influences the generation and propagation of the electric wave and usually includes the intracellular concentration of different ions (K⁺, Na⁺, Ca²⁺) and the permeability of different membrane ion channels; 𝛽 is the surface-volume ratio of heart cells; 𝐶_𝑚 is the membrane capacitance, 𝐼_ion the total ionic current, which is a function of 𝑉 and a vector of state variables 𝜼; 𝐼_stim is the current due to an external stimulus, 𝝈 is the monodomain conductivity tensor. We assume that the boundary of the tissue is isolated, that is, no-flux boundary conditions (𝐧⋅𝜎∇𝑉=0 on 𝜕Ω).

In this work, the modern and complex Bondarenko et al. model [6] that describes the electrical activity of left ventricular cells of mice was considered to simulate the kinetics of 𝐼_ion in (1). The Bondarenko et al. model (BDK) was the first model presented for mouse ventricular myocytes [6]. The ionic current term 𝐼_ion in this model consists of the sum of 15 transmembrane currents. In short, Bondarenko’s model is based on a ordinary differential equation (ODE) with 41 differential variables that control ionic currents and cellular homeostasis. In this model, most of the ionic channels are represented by Markov chains (MCs).

2.3. Numerical Discretization in Space and Time

The finite volume method (FVM) is a mathematical method used to obtain a discrete version of partial differential equations. This method is suitable for numerical simulations of various types of conservation laws (elliptical, parabolic, or hyperbolic) [33]. Like the finite element method (FEM), the FVM can be used in several types of geometry, using structured or unstructured meshes, and generates robust numerical schemes. The development of the method is intrinsically linked to the concept of flow between regions or adjacent volumes, that is, it is based on the numerical calculation of net fluxes into or out of a control volume. For some isotropic problems discretized with regular spatial meshes, the discretization obtained with the FVM is very similar to the one obtained with the standard finite difference method (FDM).

This section presents a brief description of the FVM application to the time and spatial discretization of the heterogeneous monodomain equations. Detailed information about the FVM applied to the solution of monodomain can be found in [34, 35].

2.3.1. Time Discretization

The reaction and diffusion parts of the monodomain equations were split by employing the Godunov operator splitting [36]. Therefore, each time step involves the solution of two different problems: a nonlinear system of ODEs 𝜕𝑉=1𝜕𝑡𝐶_𝑚−𝐼_ion(𝑉,𝜼)+𝐼_stim,𝜕𝜼(),𝜕𝑡=𝑓𝑉,𝜼 (2) and a parabolic PDE 𝛽𝐶_𝑚𝜕𝑉()𝜕𝑡=∇⋅𝝈∇𝑉. (3)

Since the spatial discretization of our model, ℎ, is extremely small, the CFL [37] condition that assures numerical stability is very restrictive. Therefore, for the PDE we used the unconditionally stable implicit Euler scheme. The time derivative presented in (3), which operates on 𝑉 is approximated by a first-order implicit Euler scheme as follows: 𝜕𝑉=𝑉𝜕𝑡^𝑛+1−𝑉^𝑛Δ𝑡_𝑝, (4) where 𝑉^𝑛 represents the transmembrane potential at time 𝑡_𝑛 and Δ𝑡_𝑝 is the time step used to advance in time the partial differential equation.

For the discretization of the nonlinear system of ODEs, we note that its stiffness demands very small time steps. For simple models based on Hodgkin-Huxley formulation, this problem is normally overcome by using the Rush-Larsen (RL) method [38]. However, for the most modern and complex models that are highly based on MCs, the RL method seems to be ineffective in terms of allowing larger time steps during the numerical integration. For the case of the Bondarenko et al. model, we tested both methods, Euler and RL, and both demanded the same time step, Δ𝑡_𝑜=0.0001ms for stability issues. Since the RL method is more expensive per time step than the Euler method, in this work, we used the simple explicit Euler method for the discretization of the nonlinear ODEs.

However, as already indicated above, we use different time steps for the solution of the two different uncoupled problems, the PDE and the ODEs. Since we use an unconditionally stable method for the PDE, the time step Δ𝑡_𝑝 could be much larger than that used for the solution of the nonlinear system of ODEs, Δ𝑡_𝑜=0.0001ms. In this work, we use Δ𝑡_𝑝=0.01ms, that is, a hundred times larger than Δ𝑡_𝑜. This has not introduced any significant numerical error. We calculated the L2 relative error for the transmembrane potential between a solution that uses the same time step for both the ODE and the PDE, Δ𝑡_𝑜=Δ𝑡_𝑝=0.0001ms, 𝑉_𝑚ref and a solution that uses Δ𝑡_𝑜=0.0001ms and Δ𝑡_𝑝=0.01ms, 𝑉 as follows: error=∑^𝑛𝑡_𝑖=1∑^𝑛𝑣_𝑗=1𝑉𝑖,𝑗−𝑉_𝑚ref𝑖,𝑗²∑^𝑛𝑡_𝑖=1∑^𝑛𝑣_𝑗=1𝑉_𝑚ref𝑖,𝑗², (5) where 𝑛𝑡 is the number of time steps and 𝑛𝑣 is the total number of discretized volumes. For the simulation of a tissue of size 0.5×0.5cm during 20ms (stimulus at the center of the tissue), the error found was 0.01%.

2.3.2. Spatial Discretization

The diffusion term in (3) must be discretized in space. For this we will consider the following: 𝐉=−𝝈∇𝑉, (6) where 𝐉 (μA/cm²) expresses the density of intracellular current flow and ∇⋅𝐉=−𝐼_𝑣. (7) In this equation, 𝐼_𝑣 (μA/cm³) is a volumetric current and corresponds to the left-hand side of (3), serving as the base for this finite volume solution.

For the space discretization, we will consider a two-dimensional uniform mesh, consisting of regular quadrilaterals (called “volumes”). Located in the center of each volume is a node. The quantity of interest 𝑉 is associated with each node of the mesh.

After defining the mesh geometry and dividing the domain in control volumes, the specific equations of the FVM can be presented. Equation (7) can be integrated spatially over an individual volume 𝑉_𝑖,𝑗 of size ℎ²𝑑, leading to _Ω∇⋅𝐉𝑑𝑣=−_Ω𝐼_𝑣𝑑𝑣. (8) Applying the divergence theorem yields _Ω∇⋅𝐉𝑑𝑣=_𝜕Ω𝐉⋅^→𝜉𝑑𝑠, (9) where ^→𝜉 is the unitary normal vector to the boundary 𝜕Ω. Then, we have _𝜕Ω𝐉⋅^→𝜉𝑑𝑠=−_Ω𝐼_𝑣𝑑𝑣. (10)

Finally, assuming that 𝐼_𝑣 represents an average value in each particular quadrilateral, and substituting (3) in (10), we have 𝛽𝐶_𝑚𝜕𝑉𝜕𝑡_(𝑖,𝑗)=−∫_𝜕Ω𝐉_𝑖,𝑗⋅^→𝜉𝑑𝑠ℎ²𝑑. (11)

For this particular two-dimensional problem, consisting of a uniform grid of quadrilaterals with side ℎ, the calculation of 𝐉_𝑖,𝑗 can be subdivided as a sum of flows on the following faces: _𝜕Ω𝐉_𝑖,𝑗⋅^→𝜉𝑑𝑠=𝐼_{𝑥𝑖+1/2,𝑗}−𝐼_{𝑥𝑖−1/2,𝑗}+𝐼_{𝑦𝑖,𝑗+1/2}−𝐼_{𝑦𝑖,𝑗−1/2}, (12) where 𝐼_{𝑥𝑚,𝑛} and 𝐼_{𝑦𝑚,𝑛} are calculated at faces ((𝑚,𝑛)=(𝑖+1/2,𝑗), (𝑖−1/2,𝑗), (𝑖,𝑗+1/2), or (𝑖,𝑗−1/2)) as follows. For the case in which we have defined a conductivity value at face (𝑚,𝑛), for instance the intracellular, or cytoplasm conductivity, 𝜎_𝑐, as described in Section 2.1, we have 𝐼_{𝑥𝑚,𝑛}=−𝜎_𝑐(𝑚,𝑛)𝜕𝑉𝜕𝑥_(𝑚,𝑛)𝐼ℎ𝑑,_{𝑦𝑚,𝑛}=−𝜎_𝑐(𝑚,𝑛)𝜕𝑉𝜕𝑦_(𝑚,𝑛)ℎ𝑑. (13)

For the case in which we have defined a conductance value at face (𝑚,𝑛), for instance a gap junction conduction 𝐺, as describes in Section 2.1, we have: 𝐼_{𝑥𝑚,𝑛}=−𝐺(𝑚,𝑛)Δ_𝑥𝑉|_(𝑚,𝑛),𝐼_{𝑦𝑚,𝑛}=−𝐺(𝑚,𝑛)Δ_𝑦𝑉|_(𝑚,𝑛). (14)

Using centered finite difference, we have for (13) 𝜕𝑉𝜕𝑥_{(𝑖+1/2,𝑗)}=𝑉_𝑖+1,𝑗−𝑉_𝑖,𝑗ℎ,𝜕𝑉𝜕𝑦_{(𝑖,𝑗+1/2)}=𝑉_𝑖,𝑗+1−𝑉_𝑖,𝑗ℎ. (15) For (14), we have Δ_𝑥𝑉|_{(𝑖+1/2,𝑗)}=𝑉_𝑖+1,𝑗−𝑉_𝑖,𝑗,Δ_𝑦𝑉|_{(𝑖,𝑗+1/2)}=𝑉_𝑖,𝑗+1−𝑉_𝑖,𝑗. (16)

Equations for 𝜕𝑉/𝜕𝑥∣_{(𝑖−1/2,𝑗)}, 𝜕𝑉/𝜕𝑦∣_{(𝑖,𝑗−1/2)}, Δ_𝑥𝑉∣_{(𝑖−1/2,𝑗)} and Δ_𝑦𝑉∣_{(𝑖,𝑗+1/2)} can be obtained analogously.

Rearranging and substituting the discretizations of (4) and (12) in (11) and decomposing the operators as described by (2), and (3) yields 𝜎_{𝑖+1/2,𝑗}+𝜎_{𝑖−1/2,𝑗}+𝜎_{𝑖,𝑗+1/2}+𝜎_{𝑖,𝑗−1/2}+𝛼𝑉^∗_𝑖,𝑗−𝜎_{𝑖,𝑗−1/2}𝑉^∗_{𝑖,𝑗−1}−𝜎_{𝑖+1/2,𝑗}𝑉^∗_𝑖+1,𝑗−𝜎_{𝑖,𝑗+1/2}𝑉^∗_𝑖,𝑗+1−𝜎_{𝑖−1/2,𝑗}𝑉^∗_{𝑖−1,𝑗}=𝛼𝑉^𝑛_𝑖,𝑗,𝐶_𝑚𝑉^𝑛+1_𝑖,𝑗−𝑉^∗_𝑖,𝑗Δ𝑡_𝑜=−𝐼_ion𝑉^∗_𝑖,𝑗,𝜼^𝑛,𝜼^𝑛+1−𝜼^𝑛Δ𝑡_𝑜=𝑓𝜼^𝑛,𝑉^∗,𝑡, (17)(18) where 𝛼=(𝛽𝐶_𝑚ℎ²)/Δ𝑡_𝑝, 𝑛 is the current step, ∗ is an intermediate step, and 𝑛+1 is the next time step. In addition 𝜎 can stand for any of the gap junction conductance (𝐺_𝑝, 𝐺_𝑖, 𝐺_𝑐) divided by the depth 𝑑 or for any conductivity value (𝜎_𝑐, 𝜎_𝑚) defined for each volume face as described in Section 2.1. This defines the equations for each finite volume Vol_𝑖,𝑗. First we solve the linear system associated with (17) to advance time by Δ𝑡_𝑝 and then we solve the nonlinear system of ODEs associated with (18) 𝑁_𝑜 times until we have 𝑁_𝑜Δ𝑡_𝑜=Δ𝑡_𝑝.

2.4. Parallel Numerical Implementations

Large scale simulations, such as those resulting from fine spatial discretization of a tissue, are computationally expensive. For example, when an 8 μm discretization is used in a 1 cm × 1 cm tissue and the Bondarenko et al. model (BDK), which has 41 differential variables, is used as cardiac cell model, a total of 1250 × 1250 × 41 = 64,062,500 unknowns must be computed at each time step. In addition, to simulate 100 ms of cardiac electrical activity, 64 millions of unknowns of the nonlinear systems of ODEs must be computed one million times (with Δ𝑡_𝑜=0.0001ms) and the PDE with 1.5 million of unknowns must be computed ten thousand times.

To deal with this high computational cost, two distinct tools for parallel computing were used together: MPI and GPGPU.

2.4.1. Cluster Implementation

The cluster implementation is a parallel implementation tailored to cluster of CPUs. The cluster implementation uses the PETSc [39] and MPI [40] libraries. It uses a parallel conjugate gradient preconditioned with ILU(0) (with block Jacobi in parallel) to solve the linear system associated with the discretization of the PDE of the monodomain model. More details about this parallel implementation can be found in our previous works on this topic [21–24].

To solve the non-linear systems of ODEs, the explicit Euler method was used. This is an embarrassingly parallel problem. No dependency exists between the solutions of the different systems of ODEs of each finite volume Vol_𝑖,𝑗. Therefore, it is quite simple to implement a parallel version of the code: each MPI process is responsible for computing a fraction 𝑁𝑝 of the total number of volumes of the simulation, where 𝑁𝑝 is the number of processes involved in the computation.

2.4.2. Multi-GPU Implementation

In our multi-GPU implementation, we have decided to keep the cluster approach for the solution of the linear system associated with the discretization of the PDE of the monodomain model. Therefore, multi-GPU also solves the discretized PDE with the parallel conjugate gradient preconditioned with ILU(0) (with block Jacobi in parallel) available in the PETSc library.

However, we have accelerated the solution of the systems of ODEs by using multiple GPUs. This is a different strategy from those we have used before when the full Bidomain equations (elliptic PDE, parabolic PDE, and systems of ODEs) were completely implemented in a single GPU [30], or the full Monodomain equations (parabolic PDE and system of ODEs) were completely implemented in a single GPU [25, 26, 29].

The motivation for choosing a different strategy is based on several reasons. As presented in [29], the monodomain model can be accelerated using a single GPU by 35-fold when compared to a parallel OpenMP [41] implementation running on a quad-core computer. However, this final speedup obtained by the GPU comes from a near 10-fold speedup for the solution of the PDE and a near 450-fold speedup for the solution of the nonlinear systems of ODEs. Nowadays, as manycore architecture evolves, one may easily find in the market a single computer equipped with 64 processing cores. Therefore, we believe that solving the PDE on these new machines with traditional MPI or OpenMP-based parallel implementations may outperform a single GPU implementation. On the other side, for the parallel solution of the nonlinear systems of ODEs a single GPU still easily outperforms these new manycore-based computers. This bring, us to focus GPU implementations to the parallel solutions of the millions of nonlinear systems of ODEs. A second motivation is related to the preconditioners that can be easily and efficiently implemented for the conjugate gradient method in GPUs. For the bidomain equations, efficient geometric multigrid preconditioners [30] were implemented in a single GPU, and sophisticated algebraic multigrid preconditioners [42] were implemented in a multi-GPU platform. However, both implementations are only viable for the solution of the linear system associated with the elliptic PDE of the bidomain equations. Multigrid preconditioners are too expensive and turns out to be an inefficient option for the solution of the parabolic PDE, which is the PDE type of the monodomain model. Until now, the cheap but inefficient 𝑤-Jacobi preconditioner has been the best choice for GPU implementations when it concerns the solution of the parabolic PDE [29, 42]. However, it is well known that incomplete LU (ILU) preconditioners combined with block Jacobi or additive Schwarz domain decomposition methods [23] greatly outperform Jacobi-like preconditioners on cluster computing for the solution of the PDE of the monodomain model. This argument favors cluster-like implementations as the best choice for the parallel solution of the parabolic PDE of the monodomain model (see [43] and the references cited therein). Finally, a third and last motivation is related to the particular problem we propose to investigate in this work: models that reveal the microstructure of cardiac tissue. Another recent work presented an implementation for the bidomain model for multi-GPU platforms [44]. Both PDEs and systems of ODEs were implemented on GPUs using explicit methods, Jacobi relaxation, and explicit Euler, respectively. We note that for our particular microscopic tissue model with spatial discretization of 8 μm, the approach of using an explicit and cheap solver for the PDE would be very inefficient due to the severe stability restrictions imposed by the CFL conditions [37]. Therefore, once more, this argument also favors cluster- like implementations based on implicit methods for the parallel solution of the parabolic PDE of the monodomain model.

Our multi-GPU implementation uses CUDA [45] to implement the numerical solution of the BDK cardiac cell model. The CUDA model extends the C programming language with a set of abstractions to express parallelism, that is, CUDA includes C software development tools and libraries to hide the GPGPU hardware details from programmers that can focus on important issues of the parallelism of their code rather than dealing with unfamiliar and complicated concepts from computer graphics in order to explore the computational power of GPUs for general purpose computation.

In order to run an application, the programmer must create a parallel function called kernel. A kernel is a special C function callable from the CPU but executed on the GPU simultaneously by many threads. Each thread is run by a GPU stream processor. They are grouped into blocks of threads or just blocks. The blocks can be one-, two-, or three-dimensional. A set of blocks of threads form a grid, that can be one- or two-dimensional. When the CPU calls the kernel, it must specify how many blocks and threads will be created at the GPU to execute the kernel. The syntax that specifies the number of threads that will be created to execute a kernel is formally known as the execution configuration and is flexible to support CUDA’s hierarchy of threads, blocks of threads, and grids of blocks. Since all threads in a grid execute the same code, a unique set of identification numbers is used to distinguish threads and to define the appropriate portion of the data they must process. These threads are organized into a two-level hierarchy composed by blocks and grids and two unique coordinates, called blockId and threadId, are assigned to them by the CUDA runtime system. These two built-in variables can be accessed within the kernel functions and they return the appropriate values that identify a block and thread, respectively. All the threads within a single block are allowed to synchronize with each other via a special barrier operator, called syncthread, and have access to a high-speed, per-block shared memory which allows interthread communication. Threads from different blocks in the same grid can coordinate their execution only through the use of atomic global memory operations. No assumptions are made about the execution order of thread blocks, which means that a kernel must execute correctly no matter the order in which blocks are scheduled by the hardware to run.

Some additional steps must be followed to use the GPU: (a) the device must be initialized; (b) memory must be allocated in the GPU and data transferred to it; (c) the kernel is then called. After the kernel have finished its execution, results are transferred back to the CPU.

Two kernels have been developed to solve each of the systems of ODEs related to BDK model. The first kernel is responsible for setting the initial conditions of the systems of ODEs, whereas the second one integrates the systems of ODEs at each time step.

Both kernel implementations were optimized in many different ways. The state variables of 𝑀 cardiac cells were stored in an array called SV, whose size is equal to 𝑀𝑁_eq, where 𝑁_eq is the number of differential equations of the ionic model (in this work, 𝑁_eq is equal to 41). The SV array was organized in such a way that the first 𝑀 entries correspond to the first state variable, followed by 𝑀 entries of the next state variable, and so on. Moreover, for all ionic models, the first 𝑀 entries of the SV array correspond to the transmembrane potential 𝑉. During the solution of the systems of PDEs, after the integration of the ODEs systems, the transmembrane potential of each node should be passed to the PETSC solver. Due to the memory organization chosen for the SV array, this is a straightforward task since, as stated before, the 𝑀 first entries of the array correspond to the transmembrane potential 𝑉 of each node. This organization allows us to avoid extra memory transactions between CPU and GPU, improving performance. Another implementation choice that impact performance positively was the way the SV array has been allocated. The SV array was allocated in global GPU memory using the cudaMallocPitch routine from the CUDA API. This routine may pad the allocation in order to ensure that corresponding memory addresses of any given row will continue to meet the alignment requirements for the coalescing operations performed by the hardware. In short, a strict coalescing requires that thread 𝑗 out of 𝑛 threads has to access data 𝑢[𝑗] if 𝑢[0] is accessed by thread 0, that is, each thread should perform data access by stride 𝑛. Therefore, in the first kernel, to set the initial conditions, each thread sets the values of all its state variables. The kernel that solves the system of ODEs operates similarly, that is, each thread computes and updates its state variables writing to the right position in memory that corresponds to their variables. In addition, the second kernel was optimized to use as much as local memory operations as possible.

Pure domain decomposition was used for parallelism. The tissue domain was linearly decomposed on 𝑁𝑝 nonoverlapping subdomains (or 𝑁𝑝 tasks, 𝑇₁ to 𝑇_𝑁𝑝, see Figure 4), where 𝑁𝑝 is the number of MPI processes or processing cores. The parallel solution of the PDE is implemented via PETSc (see [21]), with each processing core 𝑝 responsible for updating the variables associated to subdomain 𝑇_𝑝. In our computational environment each machine or node has more CPU cores (8) than GPUs (2). Therefore, for the solution of the ODEs each GPU device will be responsible for processing more than one task. The tasks assigned with one node are distributed to the GPUs in a round-robin fashion. For example, if 𝑁𝑝=16 and we have two machines (each with 8 cores and 2 GPU devices), Figure 4 presents how the tissue domain will be partitioned. Four tasks would be assigned to each GPU device. For instance, at node 0, GPU 0 would process tasks 𝑇₁, 𝑇₃, 𝑇₅, and 𝑇₇, GPU 1 the tasks 𝑇₂, 𝑇₄, 𝑇₆, and 𝑇₈.

Figure 4 

Linear parallel decomposition of tissue. Example for the case of two nodes (each with 8 cores and 2 GPU devices, that is, a total of 16 CPU cores and 4 GPU devices). Each CPU core processes one task. Four tasks are assigned to each GPU device. For instance, at node 0, GPU 0 processes tasks 𝑇1, 𝑇3, 𝑇5, and 𝑇7, and GPU 1 the tasks 𝑇2, 𝑇4, 𝑇6, and 𝑇8.

For the solution of the ODEs, both sequential and parallel (CUDA) codes used single precision. For the solution of the PDE we have used double precision. For the case of monodomain simulations, we have shown in [25] that the use of single precision in CUDA does not affect the numerical precision of the solver.

3. In Silico Experiments, Computational Environment, and Metrics

The simulations were performed using the microscopic model with spatial discretization of 8 μm and heterogeneous conductivity values as described in Section 2.1. The values used for 𝛽 and 𝐶_𝑚 were set to 0.14 cm⁻¹ and 1.0 μF/cm², respectively. The time step used to solve the linear system associated with (17) was set to Δ𝑡_𝑝=0.01 ms and to solve the nonlinear system of ODEs associated to (18) was set with Δ𝑡_𝑜=0.0001 ms.

Three different tissue setups were used to test our model and parallel implementations: a cardiac tissue of 0.5 cm × 0.5 cm size that was stimulated in the center and was executed for 10 ms, a cardiac tissue of 1.0 cm × 1.0 cm size that was stimulated in the center and was executed for 10 ms, and a cardiac tissue of 1.0 cm × 1.0 cm that was stimulated using the S1-S2 protocol to generate a spiral wave, a form of self-sustained reentrant activity strongly associated with cardiac arrhythmia.

Our experiments were performed on a cluster of 8 SMP computers. Each computer contains two Intel E5620 Xeon quad-core processors and 12 GB of RAM. All nodes run Linux version 2.6.18 – 194.17.4.el5. The codes were compiled with gcc 4.1.2 and CUDA 3.2. Each node contains two Tesla C1060. The Tesla C1060 card has 240 CUDA cores and 4 GB of global memory.

All tests were performed three times. The average execution time, in seconds, is then used to calculate the speedup, defined as the sequential execution time divided by the parallel execution time.

4. Results

Figure 5 presents the propagation of a central stimulus on the tissue of size 1 cm × 1 cm for different time instants. As expected, macroscopically, the propagation looks very smooth and continuous. However, when highlighting a smaller region of size 1 mm × 1 mm, see Figure 6, we can already observe the discrete nature of propagation, that is, the influence of the cardiac microstructure on the propagation of action potentials.

Figure 5 

Action potential propagation (transmembrane potential) after a central stimulus on a tissue of size 1 cm × 1 cm for different time instants. (a) 𝑡=1 ms, (b) 𝑡=7 ms, (c) 𝑡=13 ms, and (d) 𝑡=20 ms.

Figure 6 

(a) Transmembrane potential at 𝑡=7 ms after a central stimulus on a tissue of size 1 cm × 1 cm. (b) Microscopic details revealing the discrete nature of AP propagation, that is, the influence of cardiac microstructure on a large tissue size simulation.

Table 1 presents the results obtained by the parallel implementations for the experiment with a square tissue of 0.5 cm × 0.5 cm. As one can observe, the time spent solving the ODEs is responsible for near 90% of the execution time. It can also be observed that although the obtained speedups with the cluster are respectable and almost linear (near 61 with 64 cores), the total execution time remains high. With respect to the multi-GPU implementation, the results are much better. It must be stressed that although 64 cores where used in this simulation, only 16 GPGPU devices where available for executing the simulation, so 8 processes share 2 GPGPU devices per machine. As one can observe, the obtained speedup was huge, about 343 times faster than a single core processor. The execution time drops from 1.6 days (using one processing core) to only 6.7 minutes (using the 8-node multi-GPU platform).

Table 1 

Average execution time and speedup of parallel implementations for a tissue of 0.5 cm × 0.5 cm. The execution times are presented in seconds.

Table 2 presents the results obtained by the parallel implementations for the experiment with a square tissue of 1.0 cm × 1.0 cm. Once again, the speedups obtained with the cluster implementation, were almost linear (61 with 64 cores). With respect to the multi-GPU implementation the results are much better. The speedup was huge, about 420 times faster than a single core processor. The execution time drops from more than 6 days (using one processing core) to only 21 minutes (using the 8-node multi-GPU platform). We can also observe that the multi-GPU implementation was near 7 times faster than the cluster implementation when running on the 8 computers.

Table 2 

Average execution time and speedup of parallel implementations for a tissue of 1.0 cm × 1.0 cm. The execution times are presented in seconds.

As a result, using this very fast parallel implementation, we were able to simulate the formation of spiral waves, a form of self-sustained reentrant activity strongly associated with cardiac arrhythmia, see Figure 7. To the best of our knowledge, this is the first time spiral waves are simulated using a cardiac model that accounts for both the microstructure of cardiac tissue and a modern and complex myocyte model. After a couple of tries using the S1-S2 protocol to find the correct vulnerable window, we managed to generate a sustained spiral wave using this cardiac model that accounts for both the microstructure of cardiac tissue and a modern and complex myocyte model. The whole process took less than one day (around 13 hours with each simulation taking between 3 and 7 hours). Without our multi-GPU parallel implementation, this process would have taken 227 days using a single core computer or near 4 days using our cluster implementation running with 64 cores but without the GPUs.

Figure 7 

Spiral wave formation after an S1-S2 stimulus protocol at different time instants: (a) 𝑡=80 ms, (b) 𝑡=100 ms, (c) 𝑡=112 ms, and (d) 𝑡=120 ms.

5. Discussion and Future Works

Our results show that our multi-GPU parallel implementation described in Section 2.4.2 was able to significantly accelerate the numerical solution of a cardiac electrophysiology model that captures the microstructure of cardiac tissue (using a very fine spatial discretization) and is based on a very modern and complex cell model (with Markov chain formulation that has been extensively used for the characterization of ion channels). Speedups around 420 times were obtained, reducing execution times from more than 6 days (using one processing core) to only 21 minutes (using the 8-node multi-GPU platform). The hybrid Multi-GPU parallel implementation presented in this work is even more attractive if one considers that the architectures of GPUs and multicore processors continue to evolve on a fast pace.

Nevertheless, we believe our parallel implementation can be further improved. For instance, in the current implementation, the CPU cores are idle while waiting for the results of the nonlinear ODEs that are being computed by the GPU devices. For future work, we intend to evaluate different load balancing techniques to better distribute the parallel tasks between GPU devices and CPU cores and make a more efficient use of all the computational resources. Another possible improvement is related to the multilevel parallelism introduced for the solution of the bidomain equations [24] that combines task parallelism (via pipeline) and data parallelism (via data decomposition). We believe a similar combination of data and task parallelism could be also exploited for the solution of the monodomain equations to further enhance the parallel efficiency of our algorithms.

Recent studies that focus on the discrete or discontinuous nature of AP propagation have avoided the computational challenges that arise from microscopic models via the development and use of discrete models, where each cardiac myocyte is represented by a single point connected with the neighboring myocytes by different conductivities [46, 47]. This description has allowed the study of the effects of randomly distributed conductivities in the conduction velocity and on the formation of reentry patterns on cardiac tissue. Discrete models were introduced by Keener in [48] to describe the electrical propagation in a 1𝐷 cable of 𝑛𝑐 connected cells for the case of low gap-junctional coupling. In this model, the cells are assumed to be isopotential. Therefore, only gap junction conductances are considered for the connection of neighboring myocytes, that is, cytoplasmic resistance is considered to be insignificant. Recently, we have compared discrete and microscopic models for a 1𝐷 cable of connected cells [49]. We have shown that the numerical results obtained by the discrete model are similar to those obtained by the heterogeneous microscopic model for the case of low gap-junctional coupling (1%–10% of normal coupling). However, the discrete model failed for the case of normal gap-junctional coupling or moderate reduced gap-junctional coupling (50%–100% of normal coupling). The two-dimensional microscopic model developed in this work will allow us to further compare these two approaches (detailed microscopic models versus discrete models) and to better understand the benefits and limitations of each one of them. In addition, we hope that our microscopic model may also suggest ways to better develop discrete models, which are computationally less expensive than the detailed microscopic ones.

6. Conclusion

In this paper, we developed a cardiac electrophysiology model that captures the microstructure of cardiac tissue by using a very fine spatial discretization and uses a very modern and complex cell model based on Markov chains for the characterization of ion channel’s structure and dynamics. To cope with the computational challenges, the model was parallelized using a hybrid approach: cluster computing and GPGPUs. Different in silico tissue preparations were used in this work for the performance tests. We have shown that in all cases, our parallel multi-GPU implementation was able to significantly reduce the execution times of the simulations, for instance, from more than 6 days (on a single processor) to 21 minutes (on a small 8-node cluster equipped with 16 GPUs, that is, 2 GPUs per node). We believe that this new parallel implementation paves the way for the investigation of many open questions associated with the complex and discrete propagation nature of action potentials on cardiac tissue.

Authors' Contribution

B. G. de Barros and R. S. Oliveira contributed equivalently in this paper. “Therefore they can be both considered as first authors appearing in alphabetical order”.

Acknowledgments

The authors would like to thank CAPES, CNPq, FAPEMIG, FINEP, UFMG, UFSJ, and UFJF for supporting this work.