## Abstract

Intracellular Ca^{2+} dynamics provides excitation–contraction coupling in cardiac myocytes. Under pathological conditions, spontaneous Ca^{2+} release events can lead to intracellular Ca^{2+} travelling waves, which can break, giving transitory or persistent intracellular re-entrant Ca^{2+} scroll waves. Intracellular Ca^{2+} waves can trigger cellular delayed after-depolarizations of membrane potential, which if they occur in a cluster of a few hundred neighbouring myocytes may lead to cardiac arrhythmia. Quantitative prediction of the initiation and propagation of intracellular Ca^{2+} waves requires the dynamics of Ca^{2+}-induced Ca^{2+} release, and the intracellular spatial distribution of Ca^{2+} release units (CRUs). The spatial distribution of ryanodine receptor clusters within a few sarcomeres was reconstructed directly from confocal imaging measurements. It was then embedded into a three-dimensional ventricular cell model, with a resting membrane potential and simple stochastic Ca^{2+}-induced Ca^{2+} release dynamics. Isotropic global Ca^{2+} wave propagation can be produced within the anisotropic intracellular architecture, by isotropic local Ca^{2+} diffusion, and the branching Z-disc structure providing inter Z-disc pathways for Ca^{2+} propagation. The branching Z-disc provides a broader spatial distribution of ryanodine receptor clusters across Z-discs, which reduces the likelihood of wave initiation by spontaneous Ca^{2+} releases. Intracellular Ca^{2+} dynamics during catecholaminergic polymorphic ventricular tachycardia (CPVT) was simulated phenomenologically by increasing the Ca^{2+} sensitivity factor of the CRU, which results in an increased rate of Ca^{2+} release events. Flecainide has been shown to prevent arrhythmias in a murine model of CPVT and in patients. The modelled actions of flecainide on the time course of Ca^{2+} release events prevented the initiation of Ca^{2+} waves.

## 1. Introduction

Biological processes generally involve nonlinear interactions between a large number of variables and operate over a wide range of time scales and so their local dynamics are modelled by high-order, nonlinear, stiff systems of equations. They occur in an organized hierarchy of spatially extended systems of molecules, organelles, cells and tissues, and where coarse graining allows a gradient representation, may be modelled by partial differential equations, such as excitable media (Holden *et al.* 1991) models of cardiac tissue (Holden & Panfilov 1997). Spatio-temporal behaviours can be obtained by numerical solutions. As more and more biological details of the mechanisms and anisotropic and even orthotropic structures are included, the computational load increases. Quantitative prediction in medical applications, say to cardiac arrhythmias (Bernus *et al.* 2009), leads to a demand for real-time (or faster) solutions. In spite of technological advances in computing power and speed, this remains an unresolved grand challenge. The basic science dissection of mechanisms does not require real-time solutions, and can be explored in reduced dimensionality (Benson *et al.* 2008).

Intracellular Ca^{2+} regulates a diverse range of cellular processes concurrently, and this is achieved by its spatio-temporal characteristics (Berridge *et al.* 2000). In cardiac myocytes, the spatial properties include the localization and colocalization of Ca^{2+}-regulating proteins (Orchard *et al.* 2009), such as receptors, pumps, exchangers and binders, and the temporal properties including the time course of channel activation and inactivation, binding and off-binding rates of Ca^{2+} buffers. Both the spatial and temporal properties of intracellular Ca^{2+} can span across a large scale, from nanometres to a few hundred micrometres, from milliseconds to hours (Hinch *et al.* 2006). Solving such a multi-scale problem requires high-end computing.

Modelling intracellular Ca^{2+} waves provides an example of how localized, stochastic, transient elementary events can interact with each other, and build up global and sustained signals (Ishida *et al.* 1999; Izu *et al.* 2001*b*; Falcke & Malchow 2003; Okada *et al.* 2005). In cardiac myocytes, during membrane depolarization, Ca^{2+} influx via L-type Ca^{2+} channel (LCC) located at T-tubule network evokes major Ca^{2+} release from the nearby clusters of ryanodin receptors (RyR) expressed on the surface of sarcoplasmic reticulum (SR). This process is referred to as Ca^{2+}-induced Ca^{2+} release, and these discrete Ca^{2+} release events from SR as calcium sparks (Cheng *et al.* 1993). The local machinery linking LCC to RyR clusters in a confined subspace forms the calcium-release unit (CRU). Under some pathological conditions, such as increased sensitivity of CRUs or calcium overload in SR, spontaneous calcium release occurs more frequently. These more frequent calcium events can lead to initiation and propagation of Ca^{2+} waves, which can trigger delayed after-depolarizations (DADs) and cardiac arrhythmia (Eisner *et al.* 2006). However, the detailed mechanism of initiation and propagation of Ca^{2+} waves remains unclear. Computational models can be particularly useful in identifying and assessing the quantitative roles of each of the Ca^{2+} regulating components involved, such as SR content, Ca^{2+} diffusion, activation and inactivation time courses of release current, and the spatial distribution of CRUs (Soeller *et al.* 2009).

An integrative approach was adopted to combine both high-resolution experimental measurement of spatial distribution of RyR clusters and high-performance computing, to reproduce spatial and temporal properties of intracellular Ca^{2+} waves and provide insights into their initiation and propagation. Detailed spatial distribution of RyR clusters is directly reconstructed from confocal imaging to define the three-dimensional computational domain for computer simulations. Parallel computing with combined MPI and OpenMP programming is used for computer simulation. The aims are to identify the effects of branching in the Z-disc structure on the symmetry of Ca^{2+} sparks and isotropy of Ca^{2+} wave propagation, and to identify the cellular mechanisms of flecainide suppression of catecholaminergic polymorphic ventricular tachycardia.

## 2. Methods

### (a) Reconstruction of spatial distribution of RyR clusters

Rat cardiac myocytes were prepared using protocols described in Evans & Cannell (1997), and immediately fixed. Briefly, Wistar rats weighing approximately 250 g were given a lethal injection of pentobarbitone and their hearts quickly mounted on a Langendorff perfusion system. After 10–12 min of enzyme treatment, the ventricles were cut free, chopped into small pieces and triturated. The resulting cell suspension was filtered and the Ca^{2+} concentration raised to 1 mM over a period of 5 min. Cells were fixed for 10 min at room temperature in phosphate-buffered saline (PBS) containing 2 per cent paraformaldehyde and washed for a further 10 min in fresh PBS (Soeller *et al.* 2009).

Mouse monoclonal anti-RyR_{2} antibodies were used to label RyR_{2} clusters in isolated rat cardiac myocytes. The Z-discs were labelled using a mouse monoclonal antibody against alpha-actinin. A Zeiss LSM410 laser scanning confocal microscope using a Zeiss 63X NA 1.25 oil-immersion objective was used to obtain fluorescent images.

The distribution of RyR clusters was determined using a cluster detection algorithm, and provided for computational simulations in the format of Cartesian coordinates. The number of RyRs per cluster was estimated according to the methods described in Soeller *et al.* (2007).

The reconstructed experimental dataset has the dimensions of 24 μm×24 μm×7.2 μm, and spans across about three and a half Z-discs (figure 1*a*). It was discretized into a three-dimensional computational domain with 120×120×36 voxels, with a spatial resolution of 0.2 μm. This study focused on calcium interactions among spatially distributed RyR receptor clusters rather than local calcium signalling within an individual cluster. The spatial resolution is well below the nearest spacing between CRUs, e.g. 1–2 μm, and thus sufficient for this study. For global Ca^{2+} wave simulation, this experimentally reconstructed geometric model was then extended to the full length of the cell, about 114 μm, by stacking itself repetitively along the long axis.

### (b) Mathematical model for Ca^{2+} excitation

A spatially extended, stochastic model of intracellular calcium dynamics was used (Izu *et al.* 2001*a*; Li & Holden 2009; Li *et al.* 2009), with calcium excitation based on Keizer *et al.* (1998). The model consists of a set of reaction–diffusion equations, with details of stochastic triggering of CRUs, calcium pumps, leak current and cytoplasmic buffers, and can be described as2.1where *D*_{c} denotes the isotropic diffusion coefficient for myoplasm Ca^{2+}.

At a CRU, the reaction term *J*([Ca^{2+}]) is represented as2.2where is the scaling factor to Ca^{2+} release magnitude. *N*_{RyR} denotes the number of RyR receptors per cluster, and ranges from 8 to 500 with a mean value of 174 and s.d. of 144. *s* is the scale factor to release current according to the number of open channels in a cluster. Here, *s*=2 is used as suggested by Thul & Falcke (2004). *i*_{Ca} is the release current and *F* is the Faraday constant. *τ* is the release duration. *δ*(*r*) represents the spatial discretion of RyR clusters. denotes the *N*th sparks at site *M*. *Θ* is the Heaviside step function.

The probability to release Ca^{2+} from a CRU is governed by a function *P*, which is dependent on both local Ca^{2+} concentration and Ca^{2+} sensitivity factor *K*_{poss}. A CRU releases Ca^{2+} when2.3where *P*_{max} is the maximum open probability of the CRU, *n* is the Hill coefficient and *u*_{rand} is a uniformly distributed random number between 0 and 1. *τ*_{R} is the CRU refractory duration, during which no more release events can be activated.

In the myoplasm surrounding RyR clusters, *J*([Ca^{2+}]) is described as2.4where *J*_{buffers} represents immobile cytoplasmic buffers:

In the above equation, are the given association and dissociation rate constants, and [*B*_{1}]_{T} and [*B*_{2}]_{T} are the total concentrations of the two buffers. The values of Ca*B*_{1} and Ca*B*_{2} are described by two ordinary differential equations:

Moreover, the *J*_{pump} term in equation (2.4) is given bywhere *v*_{pump}, *K*_{pump} and *m* are constants. Finally, we have a constant SR leak current, *J*_{leak}=*J*_{pump}(*c*_{0}) in equation (2.4). Values for all parameters are listed in table 1.

### (c) Numerical method and its parallelization

The three-dimensional computational domain for the reaction–diffusion model is directly reconstructed from confocal imaging data, with no-flux boundary conditions. The time derivative ∂[Ca^{2+}]/∂*t* is approximated by the forward Euler method, making the entire numerical scheme fully explicit.

The computational domain has an irregular shape bounded by the membrane. In order to allow finite difference approximation of the diffusion term, by a standard seven-point stencil, we ‘immerse’ the three-dimensional irregular domain into a box-shaped uniform mesh. Each mesh point that lies outside the solution domain is marked as non-computational, and the remaining computational points are marked either as RyR clusters or myoplasm.

Finite differences are used in both time and space to discretize the reaction–diffusion equations. This is for the sake of simplicity and efficiency. An alternative is to use finite elements in the spatial discretization, which can better model the irregularly shaped domain. The drawback, however, is that this approach will make the overall numerical strategy implicit, even if the forward Euler time integrator is used (unless the mass matrix is lumped in the finite-element discretization). The finite-element approach will thus require intensive computations in the form of solving linear systems at every time step. On the other hand, the finite difference approach is fully explicit when the forward Euler time integrator is used, making it very computationally effective. Our numerical experiments also confirm that the efficient finite difference approach produces satisfactory numerical solutions.

Owing to high computational requirement, e.g. hundreds of thousands of time steps and millions of computational points, parallel computing is necessary for carrying out simulations of interest.

A general three-dimensional curved partitioning of the computational points can optimize work load balance, but irregularly shaped subdomains are not convenient for finite difference computations. We thus chose (Li *et al.* 2009) a one-dimensional partitioning of the three-dimensional box-shaped computational mesh by using cutting planes orthogonal to the longitudinal axis.

The resulting subdomains were all box-shaped, each consisting of a number of two-dimensional rectangular planes from the three-dimensional global box-shaped mesh. The advantage of such a straightforward one-dimensional partitioning scheme is that the communication pattern between the subdomains is simple. A major disadvantage arises when the number of subdomains is large, making it difficult to evenly divide the longitudinal stack of two-dimensional rectangular planes.

To reduce the load imbalance, here we used a mixed parallelization approach, keeping the number of subdomains to be relatively small, while introducing thread-based parallelism within each subdomain. This mixed approach is particularly suitable for the hardware architecture of clusters that are made up of multicore-based compute nodes. Each compute node is assigned with one subdomain, where MPI calls are used for the inter-subdomain communication. OpenMP directives are used for the intranode parallelism. Such a mixed parallelization approach has the advantage of having a relatively large size of the subdomains. The number of MPI messages needed is dramatically reduced. Parallelism between the cores on the same compute node is automatically enabled by an OpenMP capable compiler. For further performance enhancement, we have used non-blocking MPI send and receive calls (instead of blocking MPI calls) to overlap communication with computation, thus hiding the communication overhead.

Figure 1 compares the speed of three parallelization approaches for carrying out a 2.5 s simulation of calcium wave propagation. The time step is Δ*t*=10^{−5} s, the box-shaped global computational mesh has 570 points in the longitudinal direction, and 120 points in each of the two other spatial directions. Time usage measurements are obtained on a cluster using Gigabit Ethernet as the interconnect. Each compute node has two Xeon(R) E5420 2.5 GHz quad-core processors, i.e. eight cores per node. We can see from figure 1 that the mixed MPI–OpenMP approach without communication hiding is only advantageous over the flat MPI approach from Li *et al.* (2009) when the number of subdomains is large. The second approach of mixed MPI–OpenMP parallelization, which overlaps communication with computation, consistently outperforms the other two parallelization approaches.

## 3. Results

### (a) Isotropic Ca^{2+} sparks generate isotropic Ca^{2+} waves

Previous simulation studies assumed both anisotropic Ca^{2+} diffusion, with Ca^{2+} diffusion coefficient in the longitudinal direction twice that in the transversal direction, and a 2 : 1 ratio of longitudinal and transversal spacing between neighbouring CRUs (Izu *et al.* 2001*b*). However, although high-resolution imaging of RyRs and Z-disc in rat ventricular myocytes shows the 2 : 1 ratio in peripheral region, structural bifurcation of Z-discs is observed in the centre of the cell (figure 2*a*). RyR cluster distribution is colocalized with the Z-discs and followed this branching structure (Soeller *et al.* 2007, 2009). Further, Ca^{2+} sparks were recently reported to be circular (i.e. isotropic diffusion) rather than elongated (Banyasz *et al.* 2001). These lead to a different explanation for isotropic Ca^{2+} waves (figure 2*a*). The bifurcating structure of Z-disc provides an inter Z-disc pathway for calcium propagation.

For comparison with previous studies, a pseudo-geometry (figure 2*b*), in which all RyR clusters were placed on idealized flat and non-branching Z-discs, was also constructed, and to mimic the general geometry of the experimental data, these flat planes representing Z-discs were oriented at an angle of 7.6° to the transverse axis. Because of the repetitive nature of RyR distribution along the longitudinal axis, an extended geometric model (realistic) with dimensions of 120×120×570 containing 10 830 RyR clusters was constructed by cascading the original experimentally measured dataset throughout the full length of the cell for more detailed investigation of global behaviours.

In figure 2*c*, RyR clusters in a small region in the centre of the cell were activated to release Ca^{2+} at *t*=0 ms, using both realistic and pseudo-geometries. The blue semi-transparent background indicated the cellular membrane, and the Ca^{2+} waves were colour-coded according to the local Ca^{2+} concentration. Under same conditions, the realistic geometry with branching structures produced isotropic Ca^{2+} waves, with estimated conduction velocity of about 50 μm s^{−1}, which is within the range of experimental measurements (e.g. 50–90 μm s^{−1}; Cheng *et al.* 1996). The pseudo-geometry tends to slow down the longitudinal propagation due to a larger spacing between the neighbouring Z-discs, while the transversal Ca^{2+} propagation was less affected, and leading to an anisotropic wave profile.

However, these structural effects are sensitive to the value of diffusion coefficient, since a larger diffusion coefficient can eliminate the structural effects, and a smaller one may amplify the effects. The relationship between structural effects and diffusion coefficient was nicely studied in Dawson *et al.* (1999). When *D*_{c}*τ*/*d*^{2} (*d* refers to the spacing between neighbouring CRUs) is larger than 1, wave propagation becomes continuous, and more dependent on *D*_{c} instead of *d*. In Izu *et al.* (2001*b*), 150 μm^{2} s^{−1} (transversal) and 300 μm^{2} s^{−1} (longitudinal) were used for anisotropic calcium diffusion; here, the isotropic calcium diffusion coefficient of 150 μm^{2} s^{−1} was used to test the hypothesis that branching structure of Z-disc can compensate for the reduced longitudinal calcium diffusion.

### (b) Branching structure of Z-discs may suppress sustained Ca^{2+} wave events

The rate of spontaneous Ca^{2+} release events is determined by the probability function, through both local Ca^{2+} concentration and Ca^{2+} sensitivity factor. With a given Ca^{2+} sensitivity factor, local Ca^{2+} elevation seems to be the only factor to increase the firing probability of a RyR cluster. However, there are other less explicit factors, such as the differences of Ca^{2+} release intensity owing to the size of RyR clusters, and distance to the nearest neighbouring RyR cluster.

In figure 3*a*, only localized Ca^{2+} sparks were observed during 500 ms simulation within the realistic geometry. However, with the same condition, sustained global wave propagation was initiated by 200 ms within the pseudo-geometry. The Ca^{2+} spark rate was similar, since the Ca^{2+} sensitivity factors were identical; however, local Ca^{2+} releases have a greater chance to evolve into a travelling wave within the pseudo-geometry. This may be owing to a higher in-plane density of RyR clusters, which decreases the nearest neighbouring distance among CRUs.

### (c) Flecainide can suppress spontaneous calcium activities effectively

Catecholaminergic polymorphic ventricular tachycardia (CPVT) is an inherited syndrome in which the electrocardiogram appears normal, but potentially fatal polymorphic ventricular tarchycardia is induced by adrenergic stress (Watanabe *et al.* 2009). CPVT-related mutations in genes *RYR2* (encoding RyR receptor) and *CASQ2* (encoding cardiac calsequestrin) cause increased rates of spontaneous Ca^{2+} release, which may trigger DADs, and ectopic action potentials. The adrenergic intervention of the ventricles releases transmitter at varicosities, and so clusters of neighbouring cells would be exposed to elevated noradrenaline during adrenergic drive, enabling ectopic action potentials in a cluster of cells to trigger propagation.

Recently, the Na^{+} channel blocker flecainide has been shown to alter Ca^{2+} release kinetics, and to be effective in treating CPVT in a mouse model, and in patients. In figure 4, model parameters were changed to mimic CPVT and flecainide effects. In the CPVT model, *K*_{poss} was simply decreased to produce much higher frequency of Ca^{2+} release events; when applying flecainide at 50 μM, parameters were modified according to Watanabe *et al.* (2009) (see table 2 for details).

The simulation results show clearly that the effect of flecainide is to reduce the likelihood of global spontaneous Ca^{2+} waves, even though the rate of Ca^{2+} release events (determined by *K*_{poss}) is still elevated.

## 4. Discussion

A simplified stochastic model of calcium release and diffusion was embedded into the geometry of a ventricular myocyte and spatio-temporal intracellular Ca^{2+} behaviours (sparks and propagating waves) computed for normal, pathological and pharmacologically modified pathological case studies. This approach could be generalized and modified for calcium dynamics in other cells, such as smooth muscle cells and skeletal muscle cells, since various cell types share similar intracellular calcium regulation pathways.

The intracellular three-dimensional spatial distribution of RyR clusters of rat ventricular myocytes was reconstructed from high-resolution confocal imaging of a few sarcomeres, and extended to cover the full length of a typical myocyte, about 100 μm, for whole-cell simulations of intracellular Ca^{2+} phenomena. This cell model has simple dynamics, for example, the voltage-dependent membrane ionic currents are ignored, but the details of the spatial distribution of CRUs are included. One limitation is that CRUs were considered as point sources with no local control mechanisms, so a detailed stochastic calcium dynamics model including local CRU gating complexity and SR geometry needs to be developed. The simple dynamics and complicated geometry contrast with most models of ventricular cell electrochemistry that ignore spatial effects, and are high-order ordinary differential equations for intracellular compartmental ionic concentrations and ion transfer, binding and sequestration dynamics (e.g. Rudy & Silva 2006). It still remains to bring these two approaches together, to combine details of the kinetics of ion channels, exchangers and binding proteins, and their spatial distributions, within the intracellular architecture of sarcomeres and mitochondria. Such a comprehensive ventricular cell model is necessary for understanding the emerging data and ideas on the phosphorylation-based intracellular regulation of excitation–contraction coupling.

The large computational domain to model both local and global Ca^{2+} phenomena within the whole cell requires parallel computing; of three approaches tested, a mixed MPI–OpenMP parallelization that overlaps communication with computation was consistently the most efficient.

The results illustrated in figure 2 showing isotropic calcium sparks can also produce isotropic calcium wave propagation, suggest that the branching Z-discs provide inter Z-disc pathways for calcium propagation. Anisotropic propagation was obtained in the absence of branching in the the pseudo 2 : 1 geometry. This anisotropy can be explained by a denser in-plane distribution, and decreased density of cross-plane distribution of CRUs. However, the extent of speed-up of in-plane propagation and slowdown of cross-plane propagation both depend on the magnitude of the diffusion coefficient chosen. For example, an arbitrarily high diffusion coefficient will eliminate any structural effects on wave propagation, while a low diffusion coefficient that is just enough to propagate along the pseudo-plane might even lead to propagation failure across the neighbouring planes. Further, variation in the size of RyR clusters and spatial distribution of large-sized RyR clusters can also play a role. Since large-sized RyR clusters can release Ca^{2+} at higher amplitude, they can be less sensitive to the structural differences between the realistic and pseudo-geometries. If most of the large-sized RyR clusters were expressed in sub-sarcolemmal region, one may expect that the structural difference may have higher impact in the centre of the cell, where the branching structure of Z-discs is exactly located. Then, instead of a longitudinal slowdown as shown in figure 2*c*, a transversal slowdown may occur as in Li *et al.* (2009).

In our model, spontaneous release event is dependent on both local calcium concentration and sensitivity factor. When both of these two terms are fixed, structural difference between realistic and pseudo-geometries can lead to distinctive results, as shown in figure 3. Although the calcium release events are the same in both cases, local sparks develop into global waves more often in the pseudo-geometry.

The parameters can be phenomenologically modified to simulate pathologies and drug actions. The increased spontaneous Ca^{2+} release rate seen in CPVT can be reproduced by a decrease in *K*_{poss}. The effects of flecainide on Ca^{2+} release dynamics, RyR open probability, mean open time and mean closed time are readily incorporated, as in table 2. If the cellular mechanism of CPVT is via DADs triggering ectopic action potentials, the reduced incidence of Ca^{2+} triggering events illustrated in figure 4 would interact multiplicatively with a reduced likelihood of a DAD triggering an action potential, as flecainide is primarily a Na^{+} channel blocker. Thus, the clinical efficacy of flecainide for CPVT may result from its actions on two independent stages in a causal chain.

To conclude, the combination of high-resolution experimental imaging techniques, to generate three-dimensional structural models, with high-performance computing allows modelling of intracellular dynamics that links protein level data to whole cell physiology and clinical phenomenology.

## Acknowledgements

P.L. was supported by the National Science Foundation (award no. CBET-0929633), at the Center of Bioelectricity and Arrhythmia Center (CBAC), Department of Bioengineering, Washington University in St Louis. W.W. and X.C. were supported by a Center of Excellence grant from the Norwegian Research Council to the Center for Biomedical Computing at Simula Research Laboratory. This work was also supported by grants from the Auckland Medical Research Foundation, the Health Research Council and the Wellcome Trust (C.S. and M.B.C.) and the European Union through the Network of Excellence BioSim, contract no. LSHB-CT-2004-005137 (A.V.H.).

- © 2010 The Royal Society