Three-dimensional simulation of calcium waves and contraction in cardiomyocytes using the finite element method

Jun-ichi Okada,1,2 Seiryo Sugiura,2 Satoshi Nishimura,3 and Toshiaki Hisada2

1Core Research for Evolutional Science and Technology of the Japan Science and Technology Agency, Saitama; 2Computational Biomechanics Division, Institute of Environmental Studies, Graduate School of Frontier Sciences, The University of Tokyo, Tokyo; and 3Department of Cardiovascular Medicine, Graduate School of Medicine, The University of Tokyo, Tokyo, Japan

Submitted 1 June 2004 ; accepted in final form 15 October 2004


    ABSTRACT
 TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX: FORMULATION OF THE...
 GRANTS
 REFERENCES
 
To investigate the characteristics and underlying mechanisms of Ca2+ wave propagation, we developed a three-dimensional (3-D) simulator of cardiac myocytes, in which the sarcolemma, myofibril, and Z-line structure with Ca2+ release sites were modeled as separate structures using the finite element method. Similarly to previous studies, we assumed that Ca2+ diffusion from one release site to another and Ca2+-induced Ca2+ release were the basic mechanisms, but use of the finite element method enabled us to simulate not only the wave propagation in 3-D space but also the active shortening of the myocytes. Therefore, in addition to the dependence of the Ca2+ wave propagation velocity on the sarcoplasmic reticulum Ca2+ content and affinity of troponin C for Ca2+, we were able to evaluate the influence of active shortening on the propagation velocity. Furthermore, if the initial Ca2+ release took place in the proximity of the nucleus, spiral Ca2+ waves evolved and spread in a complex manner, suggesting that this phenomenon has the potential for arrhythmogenicity. The present 3-D simulator, with its ability to study the interaction between Ca2+ waves and contraction, will serve as a useful tool for studying the mechanism of this complex phenomenon.

cardiac muscle cell; excitation-contraction coupling; mechanoelectrical feedback; spiral wave; arrhythmia


CARDIAC MUSCLE CONTRACTION is regulated by rhythmic changes in intracellular Ca2+ concentration ([Ca2+]i). Under normal conditions, transsarcolemmal Ca2+ influx activated by membrane depolarization triggers synchronous Ca2+ release from the sarcoplasmic reticulum (SR) to bring about a uniform rise in [Ca2+]i (Ca2+-induced Ca2+ release, or CICR) (3). SR Ca2+ release can occur spontaneously without membrane depolarization to cause a local elevation of [Ca2+]i that propagates throughout the cell in a wavelike pattern under certain conditions (8, 9, 30). Besides their importance in basic physiology, Ca2+ waves also have clinical relevance because a focal increase in [Ca2+]i could activate a transient inward current and membrane depolarization, thus constituting a potentially arrhythmogenic event (2, 6). Accordingly, several studies have attempted to clarify the characteristics and underlying mechanisms of Ca2+ waves using single-cell (8, 15), multicellular (24), and whole heart (17) preparations.

Simulation studies have also been performed to integrate the experimental findings and provide a mechanistic explanation for them (1, 10, 16, 19, 29, 31). Basically, all of these models assumed that Ca2+ diffusion from one release site to another and CICR were the basic mechanisms, although various modifications were made by incorporating novel experimental findings, such as stochastic opening of ryanodine receptors (RyR) (16, 19) and anisotropy in diffusion (31). In addition, to reproduce the evolution of the characteristic wave pattern, including the spiral wave, simulations have been extended to two-dimensional (2-D) space (10, 16, 31). However, there are virtually no 3-D models that can be used for detailed analyses. Indeed, a recent study using a novel confocal microscope clearly demonstrated the 3-D nature of Ca2+ wave propagation (15), thus necessitating the development of a competitive 3-D model.

Another important step in excitation-contraction coupling that is not taken into consideration by current simulation models is myocyte contraction. Initially identified as a focal contraction (7), a Ca2+ wave definitely changes the distance between adjacent Ca2+ release sites as it propagates longitudinally along a myocyte. Furthermore, force development may change the affinity of troponin C (TnC) for Ca2+, a major buffering system in the cytoplasm (13). Both of these aspects could potentially modulate Ca2+ wave propagation but have not yet been investigated fully.

In the present study, we have developed a 3-D simulator of Ca2+ wave propagation and contraction in cardiac myocytes in which the sarcolemma, myofibril, and Z line with Ca2+ release sites were modeled as separate structures using the finite element method. The wave front in 3-D space can be visualized in an arbitrary 2-D plane, facilitating detailed comparisons with the results obtained in earlier experimental and model studies. Furthermore, the effect of contraction on Ca2+ waves can be evaluated by simulation for the first time.


    MATERIALS AND METHODS
 TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX: FORMULATION OF THE...
 GRANTS
 REFERENCES
 
CICR and Ca2+ diffusion model. The basic principle of the present simulation model was similar to those of previous reports (1, 10, 31) and is illustrated schematically in Fig. 1. In each sarcomere, Ca2+ release channels are located at the Z lines. Ca2+ released from the junctional SR (JSR) through a release channel (Irel) is buffered by calmodulin (CaM) in the cytosol and TnC on the thin filament or sequestered by the nonjunctional SR (NSR) (Iup) and translocated to the JSR to replenish it (Itr). The Ca2+ leak current (Ileak) was also modeled. The remaining free Ca2+ diffuses into the sarcomere space, and if the [Ca2+]i at the adjacent Z line exceeds a certain threshold, regenerative Ca2+ release takes place. In this analysis, we adopted a deterministic rule for triggering Ca2+ release (22), instead of a stochastic model (16, 19), for the following reasons: 1) although a stochastic model may be a prerequisite for analyzing the evolution of Ca2+ from a Ca2+ spark, the present analysis exclusively studied the wave propagation; and 2) implementation of a stochastic process heavily increases the computational burden. We admit that this model is not necessarily the best, but it should represent a reasonable starting point. For the same reasons, transsarcolemmal ion exchange processes were eliminated.



View larger version (30K):
[in this window]
[in a new window]
 
Fig. 1. Simulation model. A: conceptual framework of the model. The cylindrical myocyte is occupied by myofibrils, and Ca2+ release sites [junctional sarcoplasmic reticulum (JSR) nodes] are located at the Z line (top). In each sarcomere, Ca2+ is released from the JSR at the Z line (Irel) and diffuses into the cytoplasm, where it becomes bound by a buffering system [troponin C (TnC) and calmodulin (CaM)]. Free (unbound) Ca2+ is sequestered by the nonjunctional sarcoplasmic reticulum (NSR) (Iup) and translocates into the JSR (Itr). Ileak: Ca2+ leakage from the NSR. B, left: 4-state model of the regulatory and contractile protein unit (T). Yi (i = 1, 2, 3, 4, or d) and Xi (i = 1, 2, 3, 4, or d) are rate constants governing the transition between the states; right: physical entities of the states. TCa* and T* are force-bearing states.

 
The mathematical formulation for [Ca2+]i dynamics is represented by the reaction-diffusion equation described below for the y-axis coordinate in the longitudinal direction and the x- and z-axis coordinates for the transverse directions (Fig. 1):

(1)
where dt is development over time, Dc is a diagonal matrix describing the diffusivity of Ca2+, and f([Ca2+]i) describes the kinetics of Ca2+ transport into and out of the cytoplasm. In this analysis, the diffusion coefficient for Ca2+ (diagonal elements of Dc) was set to 1.0 µm2/ms for the longitudinal (y-axis) direction and 0.5 µm2/ms for the transverse (x- and z-axis) directions on the basis of previous experimental studies (16, 27). Modifications were made depending on the location to yield the following equations.

At the nodes where NSR is facing cytoplasm (NSR node in Fig. 1A, top):

(2)

At the nodes where JSR is facing cytoplasm (JSR node in Fig. 1A, top):

(3)
where (TnC – Ca2+) is the Ca2+ buffered by TnC and (CaM – Ca2+) is the Ca2+ buffered by CaM. At nodes corresponding to the nucleus, permeability to Ca2+ diffusion was altered to examine its effect. We assumed that the nuclear membrane does not act as a diffusion barrier (4).

Muscle contraction model. To relate the local [Ca2+]i to cross-bridge kinetics and force generation, we adopted the theoretical formalism proposed by Negroni and Lascano (25). Briefly, the TnC state (regardless of whether bound to Ca2+) and the cross-bridge state (attached or detached) are combined into four states as shown in Fig. 1B: 1) TnC is bound to Ca2+, and the cross bridge is attached (TnC*); 2) TnC is not bound to Ca2+, but the cross bridge remains attached (T*); 3) TnC is bound to Ca2+, but the cross bridge is detached (TCa); and 4) TnC is not bound to Ca2+, and the cross bridge is detached (T). The transitions between these four states are governed by the evolutional equations described in the APPENDIX. In the analysis examining the effect of the affinity of TnC for Ca2+, we changed the rate constant for Ca2+ binding to TnC (Y1). Of the four states, TCa* and T* contribute to force generation such that the active force per unit length (Fb) is as follows:

(4)
where A is a constant and h is the cross-bridge elongation. The sarcomere shortening concomitant with the detachment of a certain proportion of cross bridges further decreases the number of attached cross bridges through two additional paths whereby TCa* and T* decrease depending on the velocity of the shortening. This mechanism is known as shortening-induced deactivation (13).

The passive property is characterized by the force developed by the parallel elastic component (Fp):

(5)
where K is a constant, L is the length of a half-sarcomere (L), and L0 is the unstressed length of L.

The total muscle force (F) can be expressed as the sum of the active and passive forces as follows:

(6)

Cell geometry and finite element modeling. We assumed the geometry of the cell to be a cylinder with a diameter of 16 µm and a height of 104 µm. The Z lines, each of which is represented by a truss element network, are spaced at 2-µm intervals. The myofibrils were modeled by 113 vertical truss elements within a sarcomere, each of which had a diameter of 1 µm. They occupied 47% of the cross section. These values were estimated from Figs. 1 and 2 in Lipp and Niggli (21) and Fig. 2 in Ishida et al. (15). The sarcolemma including the cytoskeleton was represented by mixed interpolation of tensorial component shell elements (11), and the cytoplasm was represented by hexahedral bilinear solid elements. The nucleus had a diameter of 5 µm and a height of 10 µm (15). Its center was located 15 µm from the end in the longitudinal direction and 4 µm from the center in the cross section. We also estimated these values from data in the literature (15, 21). There were 15,168 solid elements, 16,884 truss elements, and 5,248 shell elements, and the total number of degrees of freedom was 64,401. As a constitutive law for these finite elements, we assumed an isotropic St. Venant's hyperelastic model, in which the strain energy function is calculated as follows (5):

(7)
where E is the Green-Lagrange strain tensor, E is the Young's modulus, and {nu} is the Poisson's ratio. The colon denotes the scalar products of two second-order tensors (23). Differentiation of the above equation with respect to E gives the second Piola-Kirchhoff stress tensor as follows:

where I is the second-order unit tensor, and (I {otimes} I) and l are the fourth-order tensors that operate on E as (I {otimes} I) : E = (trE)I and l : E = E, respectively. C is the fourth-order elasticity tensor that results in a constant due to the quadratic form of W, and S is the second Piola-Kirchhoff stress. These equations were defined using the coordinate system shown in Fig. 1. Because of the lack of numeric data on the material properties of the subcellular components, these values were adjusted to reproduce the cardiac muscle properties determined at the tissue level (20). The model-adjusted material properties used in this simulation are summarized in Tables 14.



View larger version (51K):
[in this window]
[in a new window]
 
Fig. 2. A: time-lapse images of a three-dimensional (3-D) simulation showing intracellular Ca2+ concentration ([Ca2+]i) changes and myocyte contraction during normal excitation. B: Ca2+ signal (fluo-3) and contraction of a single rat ventricular myocyte. [Ca2+]i rises and then decays rapidly throughout the cell.

 

View this table:
[in this window]
[in a new window]
 
Table 1. Mechanical properties of myocytes

 

View this table:
[in this window]
[in a new window]
 
Table 4. Model parameters

 
Ca2+ diffusion was analyzed by a Galerkin method-based finite element method in which bilinear hexahedral solid elements are used with the same mesh used for the cytoplasm. It has been shown that a cluster of SR Ca2+ channels is involved in the formation of a Ca2+ spark (i.e., Ca2+ release unit, or CRU) (21, 27). It also has been shown that these CRUs are discretely located at Z lines (27). CRUs were modeled at the nodes on the cross-sectioned planes corresponding to the Z lines (Fig. 1A, top). Ca2+ was sequestered or released on the basis of Eq. 2 and 3, and the diffusion was simulated using the finite element method, in which a time step of {Delta}t = 0.01 ms was used. The Ca2+ concentration thus computed was applied to Negroni and Lascano's model (25) to evaluate the contraction force. Next, the total muscle force was calculated using Eqs. 46 such that the new internal force of the truss element was determined for the finite element deformation analysis of the cell model. The resultant deformation of the myofibril was returned to the finite element model of the Ca2+ diffusion analysis and Negroni and Lascano's model. Excitation-contraction coupling was thus realized. Excitation-contraction coupling analyses were performed for every 200 steps of the Ca2+ diffusion analysis.

Experiments. Although originally recognized as spontaneous myofilament oscillation within cells (7), sarcomere dynamics during Ca2+ waves are not well understood. We monitored the sarcomere length in isolated rat cardiomyocytes during Ca2+ waves. Hearts were removed from adult male Wistar rats (200–300 g) that were under pentobarbital sodium anesthesia (50 mg/kg), and the left ventricular myocytes were isolated using enzymatic dissociation as described previously (35). Myocytes were suspended in 1.8 mmol/l Ca2+ HEPES-Tyrode solution (in mmol/l: 137 NaCl, 5.4 KCl, 1.8 CaCl2, 0.5 MgCl2, 0.33 NaH2PO4, 5 HEPES, and 5 glucose, adjusted to pH 7.4 with NaOH at 22°C), transferred to a silicone-coated (Sigmacote; Sigma) glass chamber, and viewed using an inverted microscope (x40 objective, IX71; Olympus). During spontaneous Ca2+ waves, we recorded the sarcomere length in the middle region of the sarcomere using online Fourier analysis of digitized myocyte images (SarcLen; IonOptix). From simultaneously recorded myocyte images, we measured the cell length using NIH Image software (National Institutes of Health, Bethesda, MD). To visualize Ca2+ waves, myocytes were loaded with fluo-3 by incubation in the 1.8 mmol/l Ca2+ HEPES-Tyrode solution containing 5 mol/l fluo-3 AM (Molecular Probes, Eugene, OR) for 45 min at room temperature. Myocytes were observed using a confocal microscope (CSU22; Yokogawa) equipped with a charge-coupled device camera (EVM285SPD; Texas Instruments).


    RESULTS
 TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX: FORMULATION OF THE...
 GRANTS
 REFERENCES
 
Simulation of normal excitation and Ca2+ waves. Simulated cell shortening and color-coded [Ca2+]i during normal excitation are shown as time-lapse images in Fig. 2. In this case, normal excitation was simulated by homogeneously raising the [Ca2+]i above the threshold (0.2 µmol/l; Refs. 18, 19) at time 0. Ca2+ release was evoked uniformly and instantaneously (0.1 s) along the whole cell length, such that the cell contracted quickly (0.1–0.3 s). When a localized Ca2+ spark occurred (simulated as Ca2+ release from a single release site; 0.05 s), it evolved into a Ca2+ wave and spread in opposite directions at equal velocity (Fig. 3). However, in this case, the Ca2+ wave could induce only weak contraction (0.1–0.3 s), consistent with our experimental observations (Figs. 2B and 3B) and those in another report (24). We also simulated a case in which two independent Ca2+ waves collided (Fig. 4). Waves were initiated at both ends of the cell and propagated toward the center. Because the propagation velocity was the same, the two waves collided in the middle (0.4 s) and then disappeared (0.6 s). This simulation also reproduced our experimental findings (Fig. 4B) and those described in a previous report (8).



View larger version (53K):
[in this window]
[in a new window]
 
Fig. 3. A: time-lapse images of a 3-D simulation during a Ca2+ wave. The cell shortening is moderate in length (shortening fraction 10%) but long lasting compared with normal excitation (Fig. 2; shortening fraction 20%). B: Ca2+ wave observed in the same myocyte shown in Fig. 2B. The cell shortening is moderate compared with Fig. 2B.

 


View larger version (48K):
[in this window]
[in a new window]
 
Fig. 4. A: time-lapse images of a 3-D simulation showing the collision and disappearance of 2 independent Ca2+ waves. B: collision of the Ca2+ waves observed in the experiment.

 
Wave propagation velocity. In this study, the Ca2+ wave propagation velocity was 140 µm/s under the control conditions (resting [Ca2+]i = 0.13 µmol/l, [Ca2+]JSR = 1.86 mmol/l). Because previous experimental studies (24) suggested a dependence of the propagation velocity on Ca2+ loading of the cell, we evaluated the effect of the SR Ca2+ content on the propagation velocity to investigate whether they had a linear relationship with each other (Fig. 5A). The dependence was confirmed when either the myocyte was allowed to shorten (solid line) or its length was fixed (broken line). We changed the Ca2+ content of the JSR because Miura et al. (24), in the relevant experiment, estimated SR Ca2+ loading on the basis of the released Ca2+, which presumably was stored in the JSR. We examined the effect of NSR Ca2+ store to find a similar result with JSR Ca2+ store (data not shown). We think this is because the Ca2+ content of NSR and JSR are closely related. The difference in velocity between the two conditions is discussed below. When we changed the affinity of TnC for Ca2+ by altering the rate constant for Ca2+ binding with TnC (Y1 in Fig. 1B), the velocity decreased as the affinity increased (Fig. 5B). These results are consistent with previous experimental findings (8, 24) and therefore validate the basic assumption of the current simulation model.



View larger version (17K):
[in this window]
[in a new window]
 
Fig. 5. A: relationships between the SR Ca2+ load and the Ca2+ wave propagation velocity under shortening (solid line) and isometric (broken line) conditions. B: relationship between the affinity of TnC for Ca2+ and the Ca2+ wave propagation velocity under shortening conditions.

 
Effect of contraction on Ca2+ waves. In all of the analyses described above, the myocyte was allowed to shorten in the unloaded condition. Under physiological conditions, however, each myocyte is connected to adjacent myocytes and thus under constraint. Accordingly, we simulated a case in which the myocyte length was held constant under load (isometric condition). In Fig. 6, the propagations of both the Ca2+ (left column) and contraction (sarcomere strain pattern; right column) waves are color coded and shown for the unloaded (Fig. 6A) and isometric (Fig. 6B) conditions with time on the x-axis and longitudinal location on the y-axis. The Ca2+ waves propagated smoothly under both conditions, but the unloaded Ca2+ wave reached the end of the cell earlier because of cell shortening. On the other hand, the contraction (sarcomere shortening) wave exhibited different patterns of propagation depending on the mode of contraction. In contrast to the very smooth propagation in the unloaded condition (Fig. 6A, right), isometric contraction produced a heterogeneous progression of sarcomere length distribution (Fig. 6B, right). Such heterogeneity in the strain distribution cannot be explained by the [Ca2+]i distribution or the stress equilibrium on the basis of the static stress-strain relationship. In Negroni and Lascano's contraction model (25), shortening deactivation is implemented by the force deficit terms Qd1 and Qd2, formulated as follows (see also the APPENDIX):

(9)

(10)
where dX/dt is the strain rate. We hypothesized that the complex pattern of contraction was induced by the strain rate, depending on the heterogeneity of the stress-strain relationship. In other words, in the proximity of active shortening (blue region), the sarcomere becomes softer because of the high (dX/dt)2 value and local stretching, thus leading to the complex strain pattern in Fig. 6B. To test this hypothesis, we repeated the simulation under the isometric condition without the force deficit terms Qd1 and Qd2 and compared the result with the original calculation (Fig. 6C). The results clearly revealed that the complex pattern disappeared, thus supporting our hypothesis. A similar tendency was observed experimentally in isolated cardiomyocytes. During a Ca2+ wave, each myocyte exhibited prolonged shortening to a small extent (<20 µm) (Fig. 7, A and B, broken lines). However, because the intensity of the cell attachment to the glass surface varied considerably, thereby applying different loads to the contracting cell, the shortening length of each myocyte differed significantly. In a myocyte that shortened by nearly 20 µm (Fig. 7A, broken line), the sarcomere shortened smoothly by ~0.2 µm (Fig. 7A, solid line). On the other hand, in a myocyte shortened to a lesser extent (<10 µm), the sarcomere stretched first before it shortened (Fig. 7B, solid line). In the simulation, heterogeneity in the sarcomere strain had an influence on the propagation velocity of the Ca2+ wave. Although the color-coded presentations appeared the same (Fig. 6, left), close examination of the isometric Ca2+ wave revealed local heterogeneity in the propagation velocity (Fig. 8A, broken line), in clear contrast to that of unloaded shortening (Fig. 8A, solid line). However, when averaged over its entire course, the propagation velocity of the Ca2+ wave was slightly faster in the isometric condition (Fig. 8B, thin solid line), while the velocity of the contraction wave was faster in the unloaded condition (Fig. 8B, thick broken line).



View larger version (55K):
[in this window]
[in a new window]
 
Fig. 6. A: color-coded presentation of the spatiotemporal distribution of [Ca2+]i (left) and the change in sarcomere length (L/L0) (right) during Ca2+ wave propagation in the unloaded condition. B: color-coded presentation of the spatiotemporal distribution of [Ca2+]i (left) and the change in sarcomere length (L/L0) (right) during Ca2+ wave propagation under the constraint of constant cell length. C: color-coded presentation of the spatiotemporal distribution of [Ca2+]i (left) and the change in sarcomere length (L/L0) (right) during Ca2+ wave propagation under the isometric condition without the force deficit terms Qd1 and Qd2.

 


View larger version (19K):
[in this window]
[in a new window]
 
Fig. 7. Cell length (broken line) and sarcomere length (solid line) under low (A) and high (B) load during Ca2+ wave propagation observed experimentally in isolated rat cardiomyocytes.

 


View larger version (17K):
[in this window]
[in a new window]
 
Fig. 8. A: local variation in the Ca2+ wave propagation velocity under unloaded (solid line) and isometric (broken line) conditions. B: propagation of the contraction (thick lines) and Ca2+ (thin lines) waves under unloaded (broken lines) and constant cell length (solid lines) conditions.

 
Evolution of spiral waves in 3-D space. We examined the effect of the nucleus on Ca2+ wave propagation by initiating Ca2+ sparks in the proximity of the nucleus. Each nucleus was treated not only as a region lacking a releasable Ca2+ pool but also as a buffering pool (21). In the short-axis slice, the Ca2+ wave spread around the nucleus, creating a spiral wave (Fig. 9A, top). In the longitudinal plane (middle and bottom), however, the evolution of the wave differed significantly depending on the depth of the slice. As shown in Fig. 9 (middle), the spiral wave was observed only in the slice including the nucleus. This result is consistent with the experimental observations in Ishida et al. (15), thus substantiating the importance of 3-D simulation. The propagation velocity of a spiral wave was 76 µm/s and thus slower than a planar wave. When the buffering power of the nucleus was excluded, formation of a spiral wave was not evident (data not shown) and the propagation velocity became slightly faster (81 µm/s).



View larger version (35K):
[in this window]
[in a new window]
 
Fig. 9. The numbers below the long arrow indicate the elapsed time after the initial Ca2+ release. Top: short-axis view including the nucleus. A Ca2+ wave originating from the proximity of the nucleus spreads around the nucleus to form a spiral wave. Middle (A): longitudinal section simulating confocal microscopy for a slice including the nucleus (A in short-axis view). The formation of a spiral wave is observed around the nucleus. Bottom (B): longitudinal section for a slice without the nucleus (B in short-axis view). The Ca2+ wave propagates homogeneously.

 

    DISCUSSION
 TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX: FORMULATION OF THE...
 GRANTS
 REFERENCES
 
Assumption of rapid Ca2+ buffering, The role of buffering in Ca2+ waves has been studied in detail by Keizer and colleagues (28, 33). On the basis of a simulation model, they showed that assumption of rapid buffering is valid if the equilibration time of the buffer is much smaller than the time required for Ca2+ diffusion across a region of the size of a typical gradient. Quantitatively, their criteria are expressed as follows:

(11)
where {tau} is the equilibration time constant for the buffer, k and k+ are the rate constants for the binding and dissociation between Ca2+ and the buffer, respectively, [B] is the concentration of the buffer, L is the characteristic length of the [Ca2+]i profile, and Dc is the diffusion constant of Ca2+. Calculating {tau} for CaM with the values used in our study ([Ca2+]i = 1 µmol/l, [CaM – Ca] = 50 µmol/l) and those in the previous report by Smith et al. (28), the {tau} value for CaM is ~0.00004. This value is much smaller than the L2/Dc value calculated in the present study (L2/Dc: 0.001, L = 1 µm, and Dc = 1,000 µm2/s), thus indicating that the assumption of rapid buffering can be applied to the present simulation.

In this simulation, we treated CaM as stationary for the sake of simplicity. According to Smith et al. (28), however, the existence of a mobile buffer reduces the effective diffusion coefficient for Ca2+ (Deff) as follows:

(12)
where {beta} is the differential fraction of free to bound Ca2+ (< 1), and

(13)
where Km is the dissociation constant of the mobile buffer, [Bm]T is the total concentration of the mobile buffer, and Dm is the diffusion constant for the mobile buffer. Again, calculating {gamma}mDm for CaM, the result is ~330 µm2/s for the conditions used in this simulation (Km =2.4 µmol/l, [Bm]T = 50 µmol/l, [Ca2+]i = 1 µmol/l, Dm = 32 µm2/s). This would amount to ~33% of the Dc used in this study, thus leading to an underestimation of the propagation velocity.

Comparisons with previous studies. Studies of Ca2+ waves have been promoted by a number of breakthroughs in experimental technique, such as the application of Ca2+ indicators (26, 32), laser confocal microscopy, and digital image processing (8, 15, 24, 34). Along with these developments, simulation studies also have made a significant contribution to the understanding of this complex phenomenon involving various intracellular dynamics (1, 10, 16, 19, 29). Most of these studies are based on a mechanistic model in which the diffusion of Ca2+ released by Ca2+ sparks induces a regenerative Ca2+ release from the adjacent SR (8). Our simulation model also adopted this conceptual framework and used the mathematical formulation of SR Ca2+ uptake, relocation, release, and intracellular binding proposed by Luo and Rudy (22). Although transsarcolemmal ion flux was ignored, our model could successfully reproduce the bidirectional propagation and disappearance after collision of Ca2+ waves together with the changes in cell length (Figs. 24). Furthermore, dependence of the Ca2+ wave propagation velocity on intracellular Ca2+ loading and affinity binding was also confirmed (24).

On the other hand, we did not incorporate the stochastic nature of Ca2+ sparks (16, 19), owing to limitations in computer power. In this sense, our simulator is ineffective for studying the evolution process of Ca2+ waves from partial Ca2+ sparks but can be used to investigate the coupling of a Ca2+ wave with local sarcomere contraction and complete simulation in 3-D space, which were not achieved by previous simulators.

Effect of contraction on Ca2+ waves. Using the cardiac contraction model of Negroni and Lascano (25), we coupled [Ca2+]i with contraction to evaluate its influence on Ca2+ waves. To our knowledge, this is the first attempt to couple these aspects of cardiomyocytes using full 3-D simulation. Currently, it requires ~37 h to complete the computation for a single contraction (2 s) using a CPU running at 3.2 GHz with 2-GB memory. Although the limitations in computational power forced us to adopt this relatively simple model, the analysis revealed interesting findings. As expected on the basis of experimental observations, there was significant internal shortening and stretching along the course of Ca2+ wave propagation while the cell length was kept constant (Fig. 6B, right). These changes in sarcomere length can influence a Ca2+ wave in at least two ways. First, strain applied to the thin filament increases the affinity of TnC for Ca2+ (14) and thus potentially slows the propagation. This property was incorporated into our model as a decrease in Ca2+ from TnC as the sarcomere shortened. Second, sarcomere shortening can decrease the diffusion distance to accelerate the propagation. Through these mechanisms, isometric contraction is expected to slow the propagation. In our analysis, however, isometric contraction exhibited a complex velocity pattern (Fig. 7A). To clarify these apparently contradictory results, further studies are required in both the experimental and simulation fields.

Nucleus and Ca2+ waves. Although 3-D simulations of Ca2+ waves have been reported (16, 31), most have been simple extrapolations of 2-D simulations. Considering the axisymmetric and repetitive structure of cardiomyocytes, this type of simplification may be validated. However, a recent report by Ishida et al. (15) showing the evolution of spiral waves in 3-D space highlights the necessity for complete 3-D simulation, especially in the presence of the nucleus as an obstacle to propagation. The role of the nucleus in the generation of spiral waves was investigated by Dupont et al. (10), but their 2-D simulation could not provide an answer to the vertical heterogeneity in wave propagation reported by Ishida et al. (15). Although our simulation model successfully reproduced the previously reported generation, velocity, and vertical heterogeneity of spiral waves (15, 21), we could not model the oscillatory waves emanating from the spiral waves (21). Generation of repetitive waves may require a different model for CICR. Because recent studies have revealed spreading of Ca2+ waves over trabeculae (24) and the whole heart (17), oscillatory waves may be an origin of fatal arrhythmia. Further studies are required to clarify the conditions for generating oscillatory waves.

Study limitations. As stated repeatedly, we excluded the transsarcolemmal ion flux from the analysis because of limitations in computer power. It has been shown that spontaneous Ca2+ release can cause membrane depolarization that sometimes induces an action potential (6). To analyze this potentially arrhythmogenic event, more comprehensive models that include ion channels and exchangers are needed.

Because of a lack of experimental data, the nucleus was treated as a region lacking a Ca2+ release store and a site for Ca2+ buffering. However, judging by the complex wave pattern around the nucleus reported in the literature (15, 21), the nucleus may have another role in cytosolic Ca2+ regulation. Further experimental evaluation of not only this aspect but also anisotropy in diffusion (12, 31) and Ca2+ buffering will surely contribute to the development of a complete model of cardiomyocytes.

Summary. To simulate Ca2+ wave propagation in cardiomyocytes, a finite element simulation program incorporating Ca2+ diffusion and excitation-contraction coupling mechanisms was developed. The results clearly indicate the 3-D nature of this phenomenon, the modulatory effect of contraction, and the role of the nucleus in the evolution of complex wave patterns.


    APPENDIX: FORMULATION OF THE MODEL
 TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX: FORMULATION OF THE...
 GRANTS
 REFERENCES
 
1) Ca2+ dynamics

At the NSR node:





At the JSR node:






For both the NSR and JSR nodes:

2) Ca2+ buffering in the cytoplasm

CMDN:

TnC:












At nodes including the nucleus:


3) Mechanical properties of myofibrils





Definitions of variables: [Ca2+]i, free Ca2+ concentration in the cytoplasm (mmol/l); [Ca2+]total, total cytoplasmic Ca2+ concentration (mmol/l); [CaM – Ca2+], CaM-buffered Ca2+ concentration (mmol/l); [TnC – Ca2+], TnC-buffered Ca2+ concentration (mmol/l); [Ca2+]NSR, NSR Ca2+ concentration (mmol/l); [Ca2+]JSR, JSR Ca2+ concentration (mmol/l); [Ca2+]JSR,total, total JSR Ca2+ concentration including buffered Ca2+ (mmol/l); [CSQN – Ca2+], calsequestrin-buffered Ca2+ concentration (mmol/l); Iup, Ca2+ uptake from the cytoplasm to the NSR (mmol/l/ms); Ileak, Ca2+ leakage from the NSR to the cytoplasm (mmol/l/ms); Itr, translocation current of Ca2+ from the NSR to the JSR (mmol/l/ms); IrelCICR, Ca2+ release from the JSR to the cytoplasm due to CICR (mmol/l/ms); Inucleus, Ca2+ uptake from the cytoplasm to the nucleus (mmol/l/ms); [TCa2+], thin filament site with TnC bound to Ca2+ (mmol/l); [TCa2+]eff, effective [TCa2+] (mmol/l); [TCa*], thin filament site with TnC bound to Ca2+ and an attached cross bridge (mmol/l); [T], thin filament site with free TnC (mmol/l); [T*], thin filament site with an attached cross bridge and without Ca2+ bound to TnC (force generator) (mmol/l); Qb, net rate of Ca2+ binding to T (mmol/l/ms); Qa, net rate of cross-bridge attachment; Qr, net rate of Ca2+ release from TCa* (mmol/l/ms); Qd, rate of T* detachment to give Ca2+ and T (mmol/l/ms); Qd1, additional rate of T* detachment during filament sliding (mmol/l/ms); Qd2, additional rate of TCa* detachment during filament sliding (mmol/l/ms); Fb, equivalent cross-bridge force normalized with respect to the muscle cross-sectional area (kPa); Fp, elastic force in parallel with Fb normalized with respect to the muscle cross-sectional area (kPa); F, force developed by a muscle unit normalized with respect to the muscle cross-sectional area (kPa); Km,up, half-concentration saturation of Iup; L, half-sarcomere length (µm); h, cross-bridge elongation (µm); X, length composed of half of the thick filament and the free portion of the thin filament (µm); VNSRnode and VJSRnode, NSR and VSR volume for each FE node, respectively; Vmyonode, cytoplasm volume for each FE node.


    GRANTS
 TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX: FORMULATION OF THE...
 GRANTS
 REFERENCES
 
This work was supported in part by grants from Core Research for Evolutional Science and Technology of the Japan Science and Technology Agency and National Cardiovascular Research Center, The Program for Promotion of Fundamental Studies in Health Sciences of the Organization for Pharmaceutical Safety and Research, The Research Grant for Cardiovascular Disease from the Ministry of Health, Labor and Welfare Support, the Suzuken Memorial Foundation, and the Vehicle Racing Commemorative Foundation.


View this table:
[in this window]
[in a new window]
 
Table 2. Cell geometry

 

View this table:
[in this window]
[in a new window]
 
Table 3. Initial conditions of stimulation

 

    FOOTNOTES
 

Address for reprint requests and other correspondence: J. Okada, Computational Biomechanics Division, Institute of Environmental Studies, Graduate School of Frontier Sciences, The Univ. of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo 113-0033, Japan (E-mail: okada{at}sml.k.u-tokyo.ac.jp)

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.


    REFERENCES
 TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX: FORMULATION OF THE...
 GRANTS
 REFERENCES
 
1. Backx PH, de Tombe PP, Van Deen JHK, Mulder BJM, and ter Keurs HEDJ. A model of propagating calcium-induced calcium release mediated by calcium diffusion. J Gen Physiol 93: 963–977, 1989.[Abstract]

2. Berlin JR, Cannell MB, and Lederer WJ. Cellular origins of the transient inward current in cardiac myocytes: role of fluctuations and waves of elevated intracellular calcium. Circ Res 65: 115–126, 1989.[Abstract]

3. Bers DM. Control of cardiac contraction by SR and sarcolemmal Ca fluxes. Excitation-Contraction Coupling and Cardiac Contractile Force (2nd ed.). Dordrecht, The Netherlands: Kluwer Academic, 2001, p. 245–272.

4. Bers DM. Sources and sinks of activator calcium. Excitation-Contraction Coupling and Cardiac Contractile Force (2nd ed.). Dordrecht, The Netherlands: Kluwer Academic, 2001, p. 39–62.

5. Bonet J and Burton AJ. A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations. Comput Methods Appl Mech Eng 162: 151–164, 1998.[CrossRef][ISI]

6. Capogrossi MC, Houser SR, Bakinski A, and Lakatta EG. Synchronous occurrence of spontaneous localized calcium release from the sarcoplasmic reticulum generates action potentials in rat cardiac ventricular myocytes at normal resting membrane potential. Circ Res 61: 498–503, 1987.[Abstract]

7. Capogrossi MC and Lakatta EG. Frequency modulation and synchronization of spontaneous oscillations in cardiac cells. Am J Physiol Heart Circ Physiol 248: H412–H418, 1985.[Abstract/Free Full Text]

8. Cheng H, Lederer MR, Lederer WJ, and Cannell MB. Calcium sparks and [Ca2+]i waves in cardiac myocytes. Am J Physiol Cell Physiol 270: C148–C159, 1996.[Abstract/Free Full Text]

9. Cheng H, Lederer WJ, and Cannell MB. Calcium sparks: elementary events underlying excitation-contraction coupling in heart muscle. Science 262: 740–744, 1993.[ISI][Medline]

10. Dupont G, Pontes J, and Goldbeter A. Modeling spiral Ca2+ waves in single cardiac cells: role of the spiral heterogeneity created by the nucleus. Am J Physiol Cell Physiol 271: C1390–C1399, 1996.[Abstract/Free Full Text]

11. Dvorkin EN and Batha KJ. A continuum mechanics based four-node shell element for general nonlinear analysis. Eng Comput 1: 77–88, 1984.

12. Engel J, Fechner M, Sowerby AJ, Finch AAE, and Stier A. Anisotropic propagation of Ca2+ waves in isolated cardiomyocytes. Biophys J 66: 1756–1762, 1994.[Abstract]

13. Gordon AM, Homsher E, and Regnier M. Regulation of contraction in striated muscle. Physiol Rev 80: 853–924, 2000.[Abstract/Free Full Text]

14. Hannon JD, Martyn DA, and Gordon AM. Effects of cycling and rigor crossbridges on the conformation of cardiac troponin C. Circ Res 71: 984–991, 1992.[Abstract]

15. Ishida H, Genka C, Hirota Y, Nakazawa H, and Barry WH. Formation of planar and Ca2+ waves in isolated cardiac myocytes. Biophys J 77: 2114–2122, 1999.[Abstract/Free Full Text]

16. Izu LT, Wier WG, and Balke CW. Evolution of cardiac calcium waves from stochastic calcium sparks. Biophys J 80: 103–120, 2001.[Abstract/Free Full Text]

17. Kaneko T, Tanaka H, Oyamada M, Kawata S, and Takamastu T. Three distinct types of Ca2+ waves in Langendorff-perfused rat heart revealed by real-time confocal microscopy. Circ Res 86: 1093–1099, 2000.[Abstract/Free Full Text]

18. Keizer J and Levine L. Ryanodine receptor adaptation and Ca2+-induced Ca2+ release-dependent Ca2+ oscillations. Biophys J 71: 3477–3487, 1996.[Abstract]

19. Keizer J and Smith GD. Spark-to-wave transition: saltatory transmission of calcium waves in cardiac myocytes. Biophys Chem 72: 87–100, 1998.[CrossRef][ISI][Medline]

20. Lin DHS and Yin FCP. A multiaxial constitutive law for mammalian left ventricular myocardium in steady-state barium contracture or tetanus. J Biomech Eng 120: 504–517, 1998.[ISI][Medline]

21. Lipp P and Niggli E. Microscopic spiral waves reveal positive feedback in subcellular calcium signaling. Biophys J 65: 2272–2276, 1993.[Abstract]

22. Luo CH and Rudy Y. A dynamic model of the cardiac ventricular action potential: I. simulations of ionic currents and concentration changes. Circ Res 74: 1071–1096, 1994.[Abstract]

23. Malvern LE. Introduction to the Mechanics of a Continuous Medium. Englewood Cliffs, NJ: Prentice-Hall, 1969, p. 34–35.

24. Miura M, Boyden PA, and ter Keurs HEDJ. Ca2+ waves during triggered propagated contractions in intact trabeculae: determinants of the velocity of propagation. Circ Res 84: 1459–1466, 1999.[Abstract/Free Full Text]

25. Negroni JA and Lascano EC. A cardiac muscle model relating sarcomere dynamics to calcium kinetics. J Mol Cell Cardiol 28: 915–929, 1996.[CrossRef][ISI][Medline]

26. Orchard CH, Eisner DA, and Allen DG. Oscillations of intracellular Ca2+ in mammalian cardiac muscle. Nature 304: 735–738, 1983.[CrossRef][ISI][Medline]

27. Parker I, Zang WJ, and Wier G. Ca2+ sparks involving multiple Ca2+ release sites along Z-lines in rat heart cells. J Physiol 497: 31–38, 1996.[Abstract]

28. Smith GD, Wagner J, and Keizer J. Validity of rapid buffering approximation near a point source of calcium ions. Biophys J 70: 2527–2539, 1996.[Abstract]

29. Sneyd J, Girard S, and Clapham D. Calcium wave propagation by calcium-induced calcium release: an unusual excitable system. Bull Math Biol 55: 315–344, 1993.[ISI][Medline]

30. Stuyvers BD, Boyden PA, and ter Keurs HEDJ. Calcium waves: physiological relevance in cardiac function. Circ Res 86: 1016–1018, 2000.[Free Full Text]

31. Subramanian S, Viatchenko-Karpinski S, Lukyanenko V, Györke S, and Wiesner TF. Underlying mechanisms of symmetric calcium wave propagation in rat ventricular myocytes. Biophys J 80: 1–11, 2001.[Abstract/Free Full Text]

32. Takamastu T and Wier WG. Calcium waves in mammalian heart: quantification of origin, magnitude, and velocity. FASEB J 4: 1519–1525, 1990.[Abstract/Free Full Text]

33. Wagner J and Keizer J. Effects of rapid buffers on Ca2+ diffusion and Ca2+ oscillations. Biophys J 67: 447–456, 1994.[Abstract]

34. Wier WG, ter Keurs HEDJ, Marban E, Gao WD, and Balke CW. Ca2+ "spark" and waves in intact ventricular muscle resolved by confocal imaging. Circ Res 81: 462–469, 1997.[Abstract/Free Full Text]

35. Yasuda SI, Sugiura S, Kobayakawa N, Fujita H, Yamashita H, Katoh K, Saeki Y, Kaneko H, Suda Y, Nagai R, and Sugi H. A novel method to study contraction characteristics of a single cardiac myocyte using carbon fibers. Am J Physiol Heart Circ Physiol 281: H1442–H1446, 2001.[Abstract/Free Full Text]





This Article
Abstract
Full Text (PDF)
Supplemental Movies
All Versions of this Article:
288/3/C510    most recent
00261.2004v1
Alert me when this article is cited
Alert me if a correction is posted
Citation Map
Services
Email this article to a friend
Similar articles in this journal
Similar articles in ISI Web of Science
Similar articles in PubMed
Alert me to new issues of the journal
Download to citation manager
Google Scholar
Articles by Okada, J.-i.
Articles by Hisada, T.
Articles citing this Article
PubMed
PubMed Citation
Articles by Okada, J.-i.
Articles by Hisada, T.


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Visit Other APS Journals Online
Copyright © 2005 by the American Physiological Society.