Modeling effects of L-type ca(2+) current and na(+)-ca(2+) exchanger on ca(2+) trigger flux in rabbit myocytes with realistic T-tubule geometries.

The transverse tubular system of rabbit ventricular myocytes consists of cell membrane invaginations (t-tubules) that are essential for efficient cardiac excitation-contraction coupling. In this study, we investigate how t-tubule micro-anatomy, L-type Ca(2+) channel (LCC) clustering, and allosteric activation of Na(+)/Ca(2+) exchanger by L-type Ca(2+) current affects intracellular Ca(2+) dynamics. Our model includes a realistic 3D geometry of a single t-tubule and its surrounding half-sarcomeres for rabbit ventricular myocytes. The effects of spatially distributed membrane ion-transporters (LCC, Na(+)/Ca(2+) exchanger, sarcolemmal Ca(2+) pump, and sarcolemmal Ca(2+) leak), and stationary and mobile Ca(2+) buffers (troponin C, ATP, calmodulin, and Fluo-3) are also considered. We used a coupled reaction-diffusion system to describe the spatio-temporal concentration profiles of free and buffered intracellular Ca(2+). We obtained parameters from voltage-clamp protocols of L-type Ca(2+) current and line-scan recordings of Ca(2+) concentration profiles in rabbit cells, in which the sarcoplasmic reticulum is disabled. Our model results agree with experimental measurements of global Ca(2+) transient in myocytes loaded with 50 μM Fluo-3. We found that local Ca(2+) concentrations within the cytosol and sub-sarcolemma, as well as the local trigger fluxes of Ca(2+) crossing the cell membrane, are sensitive to details of t-tubule micro-structure and membrane Ca(2+) flux distribution. The model additionally predicts that local Ca(2+) trigger fluxes are at least threefold to eightfold higher than the whole-cell Ca(2+) trigger flux. We found also that the activation of allosteric Ca(2+)-binding sites on the Na(+)/Ca(2+) exchanger could provide a mechanism for regulating global and local Ca(2+) trigger fluxes in vivo. Our studies indicate that improved structural and functional models could improve our understanding of the contributions of L-type and Na(+)/Ca(2+) exchanger fluxes to intracellular Ca(2+) dynamics.

Introduction

In cardiac ventricular myocytes, invaginations of the cell membrane, known as t-tubules, promote rapid propagation of the action potential (AP) in the cell interior (Savio-Galimberti et al., 2008; Orchard et al., 2009; Smyrnias et al., 2010). The AP activates and modulates sarcolemmal Ca2+ fluxes, including fluxes through the L-type Ca2+ channels (LCC), Na+/Ca2+ exchangers (NCX), Ca2+ ATPase pumps, and background sarcolemmal Ca2+ leak (Bers, 2001). The entry of Ca2+ via LCC and NCX triggers the sarcoplasmic reticulum (SR) Ca2+ release via ryanodine receptors (RyRs). The SR provides Ca2+ for the troponin C (TnC) myofilament protein, thereby activating and regulating myocyte contraction (Bridge et al., 1990; Bers, 2001). Controversy, however, surrounds our understanding of whether the openings of LCCs can activate allosteric Ca2+-binding sites on NCX and how the clustering of LCCs, the allosteric catalysis of NCX, and cell surface shape modulate Ca2+ trigger flux controlling SR Ca2+ release (Sipido et al., 1997; Litwin et al., 1998; Egger and Niggli, 1999; Inoue and Bridge, 2003; Bers and Ginsburg, 2007; Cheng et al., 2010).

To investigate relationships between ion fluxes via LCCs and NCXs at voltages corresponding to the early AP plateau, Sobie et al. (2008) recently measured the time-dependent Ca2+ concentration profiles ([Ca2+]) in isolated rabbit ventricular myocytes. Pharmaceutical disruption of SR activity enabled them to examine the contributions of LCC and NCX to Ca2+ trigger flux, which otherwise would have been masked by SR-bound fluxes from RyRs and the SR Ca2+ ATPase (SERCA). The study suggests that at positive voltages (Vm = +30 mV), the trigger flux of Ca2+ is greater than estimates from a simple summation of the fluxes through LCCs and Ca2+ entry via reverse NCX mode. The authors hypothesized that openings of LCCs increase local ([Ca2+]) near the NCX protein complex, which activates Ca2+-binding sites on NCX and results in an increase of Ca2+ influx via the exchanger.

Mathematical modeling complements experimental studies by enabling the examination of Ca2+ signaling and excitation-contraction coupling (ECC) in cellular and sub-cellular environments. To this end, whole-cell models have suggested an intimate relationship between Ca2+, Na+, K+ ionic fluxes, and Ca2+ transient in rabbit ventricular myocytes (Shannon et al., 2004; Mahajan et al., 2008; Aslanidi et al., 2010). These models permitted examination of the contribution of various cellular components to the evolution of the Ca2+ transient under normal and certain pathological conditions. Furthermore, recent approaches have included approximate representations of the sub-cellular geometry to introduce spatial control of the predicted local and global Ca2+ transients (Langer and Peskoff, 1996; Michailova et al., 2002; Izu et al., 2006; Koh et al., 2006; Lines et al., 2006; Means et al., 2006; Lu et al., 2009; Soeller et al., 2009; Cheng et al., 2010, 2011, 2012a,b; Louch et al., 2010; Hatano et al., 2011, 2012; Sato and Bers, 2011; Yu et al., 2011; Hake et al., 2012). Lu et al. (2009) introduced a cylindrical representation of a t-tubule of rat ventricular myocytes. Cheng et al. (2010) extended this approach to use recently published t-tubule structures imaged in mice; these t-tubules displayed a large degree of branching unlike the simple cylindrical representation used by Lu et al.

In contrast to mice, the transverse tubular system (t-system) in rabbit ventricular myocytes exhibits a simple topology that more closely resembles the t-system of canine and human ventricular myocytes (Hayashi et al., 2009; Crossman et al., 2011; Sachse et al., 2012). Scanning confocal microscopy has yielded sub-micrometer resolution details of 3D structure of the rabbit t-system (Sachse et al., 2008, 2009) including wider t-tubule cross-sections compared to rodents, and structural variations such as constrictions and flattening (Savio-Galimberti et al., 2008). These large structural differences of the t-system, as well as variations in the density and expression of proteins, suggest that the local and global Ca2+ dynamics during myocyte contraction and relaxation may be species-dependent (Bers, 2001).

In this study we extended the model of Cheng et al. (2010), originally developed for rodent myocytes to incorporate structural data and model parameters specific to rabbit ventricular myocytes. We examined local and global Ca2+ dynamics under the influence of two applied transmembrane voltages (0 and +50 mV). We included equations describing the voltage-dependent background Ca2+ leak and the sarcolemmal Ca2+ ATPase pump (Lu et al., 2009). Furthermore, we evaluated a model of NCX allostery proposed by Weber et al. (2001) against experimental data (Sobie et al., 2008) and tested the hypothesis that the allosteric activation of NCX augments the global trigger flux for SR Ca2+ release in rabbits. Preliminary results of this work have been presented to the Biophysical Society in abstract form (Kekenes-Huskey et al., 2011).

Materials and Methods

A summary of important model properties is presented here. Detailed methods, descriptions, definitions of variables and abbreviations, parameter values, and model equations are provided in the Supplementary Material. The model code is available to download from http://mccammon.ucsd.edu/smol and http://www.fetk.org, respectively.

Geometric model

The model geometry was derived from the published structural data in rabbit ventricular myocytes (Sachse et al., 2008, 2009; Savio-Galimberti et al., 2008). The model contains one repeating unit inside the cell that includes: realistic surface sarcolemma; realistic t-tubule and its surrounding half-sarcomeres (Figures 1A,B). The surrounding half-sarcomeres were modeled as a rectangular-shaped box of 2.34 μm × 2.58 μm in the plane of external sarcolemma and 5.76 μm in depth.

Figure 1

Model geometry and diagram illustrating Ca. (A) The cardiac sarcolemma, including external and t-tubule membranes were visualized using scanning confocal microscopy and labeled in Blender. Localized aggregates of L-type Ca2+ channels (red spots) were placed randomly within t-tubule membrane. (B) T-tubule mesh and its surrounding half-sarcomeres (upper panel); external membrane and t-tubule mouth (lower panel). (C) Schematic drawing of Ca2+ entrance and extrusion via the sarcolemma and Ca2+ buffering and diffusion inside the myocyte: LCC, L-type Ca2+ current; NCX, Na+/Ca2+ exchanger; Pump, membrane Ca2+ ATPase pump; Leak, background sarcolemmal Ca2+ leak; SR, sarcoplasmic reticulum; TnC, troponin C; CaM, calmodulin; ATP, adenosine triphosphate; Fluo-3, fluorescent dye. In all numerical experiments: LCC and NCX current densities were ninefold and threefold higher, respectively, in the t-tubule membrane; Ca2+ leak and pump were uniformly distributed along the sarcolemma; LCC clusters (diameter of (∼200 nm) had the same current density in the outer and t-tubular sarcolemma. The line-scan was positioned at 200 nm away from the t-tubule mouth [yellow line and yellow spot in (B)]. Local Ca2+ transients were extracted at two featured spots along the scanning line (black spots) and along the t-tubule membrane (green spots).

The t-tubule diameter varied from between 0.39 and 0.62 μm and the length was ∼4.6 μm. The constrictions occurred every 1.87 ± 1.09 μm along the principal axis of the t-tubule and cross-sectional area of these constrictions was reduced to an average of 57.7 ± 27.5% (Savio-Galimberti et al., 2008). The volume of the model compartment was estimated to be 0.0282 pL. The compartment membrane area was 15.9 μm2 where the cell membrane within t-tubule was 7.8 μm2 and within the external membrane 8.1 μm2.

Reaction-diffusion model

The effects of four exogenous and endogenous Ca2+ buffers (Fluo-3, ATP, calmodulin, and TnC) were considered (Figure 1C). The endogenous stationary buffer TnC was distributed uniformly throughout the cytosol, but not on the sarcolemma or the sub-sarcolemmal space (40–50 nm in depth; Shannon et al., 2004; Cheng et al., 2010). The free Ca2+ and mobile buffers (Fluo-3, ATP, and calmodulin) diffuse and react throughout the cytoplasm. The outer sarcolemma and sarcomere box faces were subject to reflective boundary conditions.

Model currents

We examined the contributions of four ionic currents; ILCC, INCX, IPump, and ILeak. Formulation of INCX was presented as the product of an electrochemical (ΔE) and an allosteric factor (Allo; Weber et al., 2001). The maximum NCX rate value (Vmax_NCX = 0.207 μM/ms) used here was from Shannon et al. (2004). The NCX allosteric constant (KmCaAct) was fitted to 0.29 vs. 0.256 μM in Shannon et al. (2004). The KmCaAct value was adjusted to approximate the curvature of +50 mV whole-cell Ca2+ transient data from Sobie et al. (2008).

Immunohistochemical studies have demonstrated that most of the LCCs are concentrated in the t-tubules (from 3 to 9 times more than on the external sarcolemma; Scriven et al., 2000; Soeller et al., 2009). In this study, in agreement with the reported data, the LCC current density was assumed to be ninefold higher in the t-tubular membrane than in the outer cell surface. Measurements of Ca2+ sparks in rabbits and other species imply that a cluster of LCCs is likely involved in gating a cluster of RyRs (Lipp and Niggli, 1998; Bridge et al., 1999; Scriven et al., 2000, 2010; Takagishi et al., 2000; Inoue and Bridge, 2003; Altamirano and Bers, 2007; Dan et al., 2007; Poláková et al., 2008; Sobie and Ramay, 2009; Louch et al., 2010). Furthermore, 3D visualizations of fragments from rabbit ventricular myocytes have revealed that the majority of RyR clusters are adjacent to the t-system and that these clusters are irregularly distributed along t-tubules (Dan et al., 2007; Sachse et al., 2009). Thus here three patches (∼200 nm in diameter with the same LCC current density) were placed randomly along the t-tubule (Figure 1A). Data also suggest a smaller number of RyR clusters to co-localize with LCCs on the external membrane (Franzini-Armstrong et al., 1999; Chen-Izu et al., 2006; Dan et al., 2007; Baddeley et al., 2011; Sachse et al., 2012). Dan et al. (2007) measured ∼2 μm longitudinal periodicity of RyR clusters on the cell surface and Sachse et al. (2009, 2012) observed irregular cluster distributions in transverse sheets. As surface RyR cluster positions with respect to t-tubule mouth remain controversial, in our model we assumed that the LCC current density was continuous along the outer membrane. Because in adult rabbits a lesser degree of clustering for NCX has been demonstrated (Lin et al., 2009; Gershome et al., 2011), we also assumed a continuous NCX distribution with a threefold higher density along the t-tubule (Neco et al., 2010). Here Ca2+ pump and leak currents were uniformly distributed along the model surface (Shannon et al., 2004; Lu et al., 2009; Brini and Carafoli, 2011).

In this study, each simulation started with basal cytosolic ([Ca2+]) of 0.1 μM, and buffers in equilibrium. The extracellular Ca2+ and Na+ concentrations were 2 and 140 mM, respectively, and remained constant. The voltage-clamp protocols (holding potential −50 mV, voltage step to 0 or +50 mV for 200 ms) were derived from data in rabbits (Sobie et al., 2008). Each current density (ILCC, INCX, IPump, and ILeak) was converted to Ca2+ flux (see Eq. A7, Supplementary Material) by using the experimentally suggested surface to volume ratio (Cm/Vcell ∼ 4.54 pF/pL) in adult rabbit ventricular myocytes (Bers, 2001). To ensure the total Ca2+ flux through L-type channels to be the same as measured at given voltage, the model protocols for whole-cell LCC current were fitted vs. data reported in rabbits (Sobie et al., 2008). To justify the model predictions a solution convergence analysis has been performed. Multiple tests, including refining the original mesh or changing the original time-step size for integration (0.5 ms), demonstrated that the used mesh and time step were correctly chosen (see Figure A1 in Appendix).

Figure A1

Solution convergence analysis. Upper panel: refining the original mesh. Lower panel: changing the original time step of 0.5 ms. Global Ca2+ transients experimentally measured (symbols). Membrane voltage at 0 mV (red lines and symbols) and at +50 mV (blue lines and symbols). [Na+] 10 mM.

Results

Global Ca2+ signals

We validated the model against whole-cell measurements of the L-type Ca2+ current and Ca2+ transient at two voltages (0 and +50 mV) in the presence of 50 μM Fluo-3 and [Na+] of 10 mM (Sobie et al., 2008). The voltage-clamp protocols and LCC currents are shown in Figures 2A,B. Global NCX time-courses are shown in Figure 2C. Calcium pump and leak fluxes are not shown here due to their small contribution. The model predicts at 0 mV a steep increase in global [Ca2+] that tapers off at ∼40 ms and converges to 0.19 μM as shown experimentally (Figure 2E). At +50 mV, a more gradual accumulation of Ca2+ is predicted until converging to ∼0.2 vs. 0.19 μM experimentally. These results indicate that our model is a reliable representation of whole-cell Ca2+ dynamics as measured in rabbit ventricular myocytes (Sobie et al., 2008).

Figure 2

Membrane currents, calcium signals, and trigger fluxes in the presence of 50 μM Fluo-3. (A) Voltage-clamp protocols. (B) Whole-cell LCC current fitted and plotted (solid lines) vs. data reported in rabbits with SR blocked (triangles and dots). (C,D) Predicted global and sub-sarcolemmal Na+/Ca2+ fluxes. (E) Global Ca2+ transients (solid and dash lines) vs. experimentally measured (triangles and dots). (F) Predicted trigger fluxes from global Ca2+ transient. Inset shows normalized whole-cell LCC currents (black lines) and normalized trigger fluxes at 0 and +50 mV. (G,H) Local Ca2+ trigger fluxes at 1.5 and 5.3 μm (green spots in Figure 1B). Membrane voltage at 0 mV (red lines and symbols) and at +50 mV (blue lines and symbols). [Na+] 10 mM (solid lines) and [Na+] 0 mM (dash lines).

Global and sub-sarcolemmal NCX fluxes

In agreement with experiment (Sobie et al., 2008), the model predicts that global NCX flux (INCXglobal; computed by averaging local [Ca2+] levels for the entire compartment) is in reverse mode at both applied voltages with 10 mM [Na+] (Figure 2C). For [Na+] of 0 mM, the NCX reverse mode is inactivated and an outward Ca2+ flux is predicted. The contribution to the global Ca2+ transient due to NCX was largest at +50 mV and [Ca2+] monotonically increases with time.

Figure 2D shows NCX flux calculated by averaging local [Ca2+] levels in the sub-sarcolemmal space (INCXsarc). In the presence of 10 mM [Na+], our results demonstrate that: (1) the overall scale of INCXsarc was ∼5.6-fold greater than INCXglobal; (2) the increase in INCXsarc during ILCC upstroke relative to INCXglobal was much faster; (3) while INCXglobal monotonically increased over the entire simulation, local peaks in INCXsarc are predicted at 15 ms. Moreover, in absence of [Na+], the model predicts a sharp reversal in INCXsarc at 0 mV, whereas INCXsarc had a more gradual and monotonically decreasing flux.

Global and local Ca2+ trigger fluxes

The whole-cell trigger flux reported by Sobie et al. (defined as dF/dtmax) quantifies the maximum inward Ca2+ flux. By normalizing dF/dtmax and ILCC to 1.0, it was suggested that the relative contribution of global NCX flux can be estimated (Sobie et al., 2008). Here we first calculated our global trigger flux (d[Ca2+]global/dt, Figure 2F) from the predicted global Ca2+ transient. We then normalized ILCC and d[Ca2+]global/dt based on their maximum values, which both occur at 15 ms for 0 mV. In the presence of 10 mM [Na+], our results show that at +50 mV, d[Ca2+]global/dtmax is ∼30% of the value at 0 mV while experimental value is reported as 45%. Furthermore, ILCC constitutes 70% of the trigger flux at +50 mV in comparison to 50% measured (Figure 2F Inset vs. Figure 2 in Sobie et al.). For zero [Na+], in agreement with experiment (Figure 2F Inset vs. Figure 3 in Sobie et al.), there is negligible impact on d[Ca2+]global/dtmax at 0 mV (∼10% reduction) and at +50 mV, a ∼30% reduction in d[Ca2+]global/dtmax is predicted relative to a measured 50% decline.

The model also provides estimates of local Ca2+ trigger fluxes within the sub-sarcolemmal space. At 0 mV, the model predicts an approximately threefold increase in d[Ca2+]/dtmax near the t-tubule mouth relative to the global trigger flux (Figures 2F,G). Moreover, the local d[Ca2+]/dtmax at the cell exterior and distal end of t-tubule (adjacent to LCC cluster) were very different: 3 vs. 9 μM/s at 0 mV; 1.3 vs. 2.5 μM/s at +50 mV (Figures 2G,H). At +50 mV a roughly 30% drop in d[Ca2+]mouth/dtmax was predicted in the absence of [Na+].

Local Ca2+ signals

In Sobie’s et al. (2008) experiment in rabbits the sarcolemma were not labeled and the fluorescence signal was recorded along the single scan lane at unknown orientations (Sobie personal communication). Due to these experimental limitations, we assumed the scanned line oriented in transverse cell direction (Figure 1 in Sobie et al., Figure (1B yellow line) to probe contributions of the realistic surface and t-tubule shape to local Ca2+ profiles. The calculated line-scan images and local Ca2+ time-courses are shown in Figures 3A–H. In agreement with experiment with 10 mM [Na+] non-uniform Ca2+ distributions for both voltages are predicted (Figures 3E,G vs. Figure 1 in Sobie et al.). The contribution of a LCC cluster located ∼2 μm away from upper surface is evident as a spike in the line-scan images. At +50 mV a slight Ca2+ gradient is predicted in transverse cell direction (Figure 3G). In addition, results in Figures 3A,D show that local Ca2+ transients at both locations along the scanning line (1.5 μm red and blue circles; 5.3 μm red and blue triangles) follow the same trends as the global Ca2+ transients.

Figure 3

Local calcium signals in the presence of 50 μM Fluo-3. (A–D) Local Ca2+ transients taken at two featured spots along the scanning line (red and blue lines and symbols) and the t-tubule membrane (black lines and symbols): circle −1.5 μm; triangles −5.3 μm. Membrane voltage at 0 mV [red and black lines in (A,B)] and at +50 mV [blue and black lines in (C,D)]. [Na+] 10 mM (solid lines) and [Na+] 0 mM (dash lines). (E–H) Calcium concentrations visualized as line-scan images in transverse cell direction for the cases described in (A–D).

The model also predicts that local Ca2+ transients in the proximity of a t-tubule (Figure 1B green spots) differ considerably with respect to local cytosolic Ca2+ transients. Adjacent to the LCC cluster located at the t-tubule distal end (Figure 1A), the local Ca2+ transients closely resembles the ILCC profiles with times to peak ∼15 ms followed by a decay to the cytosolic [Ca2+] by 200 ms (Figures 3A–D black triangles). Toward the cell surface, the trends depend on the applied voltage (Figures 3A–D black circles). At 0 mV, the [Ca2+] peak at the cell surface is nearly half of the transient at the LCC cluster and occurs ∼10 ms later. At +50 mV, the sub-membrane [Ca2+] initially increases faster than the cytosolic analog (for t < 20 ms) and thereafter [Ca2+] increases at the same rate as local cytosolic. We should mention here, that such local Ca2+ signals are difficult to resolve experimentally due to optical blurring and noise. Thus, our modeling study is one more example that computational models may serve as powerful tools for prediction and analysis on how local Ca2+ dynamics is regulated.

The most compelling evidence of the voltage and [Na+]-influenced [Ca2+] transient is shown in Movie S1 in Supplementary Material. Consistent in all movies is the predicted large and steep Ca2+ gradient in the narrow sub-sarcolemmal region. In these movies, the LCC clusters are clearly evident as localized regions of enhanced local [Ca2+]. Comparison of the 0 vs. +50 mV cases (top and bottom rows) demonstrates that at both voltages [Ca2+] increases heterogeneously in transverse cell direction as suggested by experiment (Sobie et al., 2008). A spontaneous increase in sub-membrane [Ca2+] at +50 mV propagating within the cell late in the simulation is predicted also, while initiation and propagation of Ca2+ gradients at both voltages has not been observed during experiment. Additional interesting results in the presence 10 mM [Na+] are that: (1) Ca2+ gradient traveling from the external membrane to the cell interior is predicted at both voltages when LCCs were continuously distributed along the t-tubule (Movie S2 in Supplementary Material, right column); and (2) lesser inwardly propagating Ca2+ gradient was observed at +50 mV when using the non-allosteric model of NCX (e.g., Allo = 1; Movie S3 in Supplementary Material, lower left panel).

The interaction between adjacent Ca2+ release units (CRU, local functional unit where LCC and RyR clusters reside) has been suggested to be critically dependent on the distance between one unit and its immediate neighbor (Franzini-Armstrong et al., 1999; Scriven et al., 2000; Izu et al., 2006; Dan et al., 2007; Hayashi et al., 2009). To test this, we held fixed the LCC cluster placed ∼2 μm away from upper surface while placing the other two clusters at various inter-cluster spacings. These spacings include ∼0.57, ∼0.8, 1.07, and ∼1.8 μm, whereas 0.78 ± 0.21 μm for the nearest end-to-end CRU distance was measured in rabbits (Dan et al., 2007; Savio-Galimberti et al., 2008). We found that predicted spatial [Ca2+] distributions are highly sensitive to the spacing between LCC clusters (Movie S4 in Supplementary Material).

Finally, to gain further insights on the role of NCX flux in regulating local Ca2+ dynamics, we tested how the changes in NCX allosteric constant (KmCaAct) and maximum exchanger rate (Vmax_NCX) affect the results. Our data indicate that local Ca2+ profiles are highly sensitive to KmCaAct and Vmax_NCX alterations (Movie S5 in Supplementary Material vs. Movie S1 in Supplementary Material left lower panel).

Discussion

Model geometry and L-type Ca2+ channels clustering

Numerical studies of Ca2+ signaling in rodent cardiomyocytes that consider non-trivial geometries have shown that the spatial Ca2+ dynamics depends on myocyte micro-structure, including the branching of t-tubules (Lu et al., 2009; Soeller et al., 2009; Cheng et al., 2010; Hake et al., 2012; Hatano et al., 2012). Therefore, we hypothesized that details specific to rabbits, e.g., the linear t-tubule structure with varying diameter and eccentricity, provide unique control of the Ca2+ transient (Sachse et al., 2008; Savio-Galimberti et al., 2008). Our studies indicate that these factors may contribute to a spatially non-uniform [Ca2+] distribution (Sobie et al., 2008; Bridge and Sachse, personal communication) while in rats, the measured [Ca2+] profiles were more evenly distributed with SR activity disabled (Cheng et al., 1994).

In this study, we also assumed a clustered distribution for LCCs along the t-tubule based on a random distribution. The model predicts that clustering of LCCs resulted in more uniform [Ca2+] profiles along the transverse cell direction relative to a continuous LCCs distribution. This greatly reduced the amplitude of the outer sarcolemmal compared to the continuous LCCs distribution (Movie S2 in Supplementary Material). New findings are also that local Ca2+ levels are highly sensitive to LCC cluster positions along the t-tubule (Movie S4 in Supplementary Material). In the model LCCs were also uniformly distributed on the cell surface, yet RyR distribution data suggest a small number of LCC clusters may also co-localize on the external membrane (Chen-Izu et al., 2006; Dan et al., 2007; Sachse et al., 2009, 2012). Thus, we placed two LCC clusters within the surface membrane near the t-tubule mouth but found that the LCC cluster placement intensified the under-membrane [Ca2+] non-uniformity (data not shown). We speculate here that the localization of two LCC clusters along the cell surface within the half-sarcomere micro-domain probably overestimates the outer sarcolemma contribution (two outer LCC clusters vs. three t-tubule LCC clusters). In addition, our assumption of identical current profiles for each LCC cluster may be inappropriate, since the cluster shape and size, and the number of L-type channels involved in single spark triggering will certainly modulate the local current and local [Ca2+] profiles (Franzini-Armstrong et al., 1999; Chen-Izu et al., 2006; Hayashi et al., 2009; Louch et al., 2010; Scriven et al., 2010). Although our concept of LCC clustering is quite rudimentary, this is a first attempt to examine the effects of LCC clusters, incorporated in more realistic membrane shapes, on local Ca2+ dynamics in light of evidence that LCC clusters exist in junctional clefts and form part of the couplon. A more appropriate description may require modeling LCC patches comprised of LCC “sub-clusters,” given their suggested involvement in the triggering of RyRs in the dyadic junction (Louch et al., 2010). Labeling techniques capable of resolving the localized positions of LCC clusters would permit a more detailed analysis of LCC and NCX contributions to Ca2+ trigger flux (Jayasinghe et al., 2009).

Allosteric catalysis of Na+-Ca2+ exchanger and Ca2+ trigger fluxes

A unique feature of the presented model was its ability to directly examine the role of catalytic Ca2+-binding sites on NCX in regulating local and global Ca2+ trigger fluxes. Here our studies indicate that upon depolarization, the non-allosteric NCX rapidly entered reverse mode and contributed a substantial constant inward flux during the entire simulation with 10 mM [Na2+] [INCXglobal(0mV) ∼ 3 μM/s, INCXglobal(+50mV) ∼ 8 μM/s, INCXsarc(0mV) ∼ 9 μM/s, INCXsarc(+50mV) ∼ 22 μM/s]. At 0 mV, when the LCC flux was dominant, a modest decrease in d[Ca2+]mouth/dtmax (2.5 vs. 3 μM/s) was predicted. At +50 mV, when NCX contributes a significantly larger fraction of the total Ca2+ transient, the effects of NCX allostery were more evident. Assumption of a non-allosteric NCX model resulted in a approximately twofold drop in d[Ca2+]mouth/dtmax and lesser inwardly propagating Ca2+ gradient at +50 mV (Movie S3 in Supplementary Material lower left panel vs. Movie S1 in Supplementary Material lower left panel). These data reveal that the allosteric catalysis of NCX can augment the local trigger fluxes as well (Litwin et al., 1996; Ramirez et al., 2011). No visible changes, however, were detected in d[Ca2+]end/dtmax and d[Ca2+]global/dtmax at both applied voltages (data not shown). A possible reason for the negligible predicted effects of NCX allostery on d[Ca2+]end/dtmax and d[Ca2+]global/dtmax, which is at odds with experimental estimates of the global trigger fluxes (Figure 2F), might be the assumed continuous NCX distribution along the cell membrane. We speculate that the NCX allosteric effect would be much more pronounced, if reverse-mode Ca2+ entry were localized to a smaller region typical of a cluster, which could then activate Ca2+ binding to NCX to a greater degree (Jayasinghe et al., 2009). Moreover, the relative spacing between an NCX cluster and LCC and alterations in normal NCX maximum rate and allosteric constant (Movie S5 in Supplementary Material) could modulate the amplitude of the NCX reverse mode.

Limitations of the model

Although we demonstrate good correlation with whole-cell experimental data using a relatively small sub-cellular geometric domain, this assumption implies that all domains are identical in shape and flux distribution. T-tubule data in rabbits from Savio-Galimberti et al. (2008), however, suggest considerable variety both in terms of the diameter, shape, and arrangement of tubules, as well as the shape of the cell exterior. Including such details could permit investigation into coupling between adjacent tubules in rabbits.

Furthermore, we restrict our model to Ca2+ and Ca2+-buffer dynamics, whereas additional ions, namely Na+ and K+, play a significant role in modulating ECC coupling (Despa et al., 2003; Bers and Despa, 2009; Torres et al., 2010). In particular, our results suggest that changes in [Na+] modulate NCX activity independent from any Ca2+-dependent allosteric effects (Movie S3 in Supplementary Material). By explicitly including the primary fluxes for Na+ and K+ (INa, INaK), intracellular Na+ diffusion and buffering, and distributions of Na+ channels and Na+/K+ pumps along the sarcolemma we may further examine the role of NCX allostery due to LCC opening, especially early in the AP. We further anticipate that the allosteric interaction between NCX and LCCs will be more nuanced during an AP, given that the fast Na+ promotes activation of the exchanger before many LCCs have opened (Torres et al., 2010). In this study also, the mitochondrial Ca2+ fluxes were omitted (Dash et al., 2009; Pradhan et al., 2010). These are likely to have important impact on predicted Ca2+ profiles as well (Nguyen et al., 2007; Cortassa and Aon, 2012; Dedkova and Blatter, 2012).

Conclusion

We developed and validated a 3D reaction-diffusion model of Ca2+ signaling in rabbit ventricular myocytes. The model incorporates realistic t-tubule and cell surface topologies, as well as clusters of LCCs along the t-tubule membrane. The key findings are: (1) the linear t-tubule topology and the punctuate spatial distribution of LCC along the sarcolemma function to inhibit inwardly propagating Ca2+ gradients; (2) local trigger fluxes of Ca2+ are at least threefold to eightfold higher than whole-cell Ca2+ trigger flux; and (3) the activation of allosteric Ca2+-binding sites on Na+/Ca2+ exchanger could provide a mechanism that regulates local and global Ca2+ trigger fluxes in vivo. We concluded that improved models representing the localized positions and the functional diversity of LCC clusters could help to improve our understanding of ILCC and INCX contributions to global and local Ca2+ trigger fluxes.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at http://www.frontiersin.org/Computational_Physiology_and_Medicine/10.3389/fphys.2012.00351/abstract

Effects of voltage and intracellular Na.

Click here for additional data file.

Effects of continuous LCC current distribution along the t-tubule on spatial Ca.

Click here for additional data file.

Effects of voltage and intracellular Na.

Click here for additional data file.

Effects of LCC cluster placement along t-tubule on spatial Ca.

Click here for additional data file.

Effects of changes in NCX allosteric constant and maximum NCX rate on spatial Ca.

Click here for additional data file.

Table A1

Physical constants and cell geometry parameters.

Symbol	Definition	Value	Reference
PHYSICAL CONSTANTS
F	Faraday constant	96.5°C mmol−1	Physical constant
T	Temperature	295 K	Physical constant
R	Universal gas constant	8.314 J mol−1 K−1	Physical constant
WHOLE-CELL GEOMETRY
Vcell	Cell volume	30.4 pL	Satoh et al. (1996)
Cm	Cell capacitance	138 pF	Satoh et al. (1996)
COMPARTMENT AND t-TUBULE GEOMETRIES
Vms	Compartment volume	2.82e−2 pL	Estimated
Smc	Compartment surface	15.9 μm2	Estimated
Smc-tubule	T-tubule surface area	7.8 μm2	Estimated
dlng	Longitudinal cell direction	2.58 μm	Estimated
dax	Axial cell direction	2.34 μm	Estimated
dtr	Transverse cell direction	5.76 μm	Estimated
rtubule	T-tubule radius	0.2–0.31 μm	Savio-Galimberti et al. (2008)
ltubule	T-tubule length	4.6 μm	Savio-Galimberti et al. (2008)

Table A2

Calcium and buffer reaction-diffusion parameters.

Symbol	Definition	Value	Reference
Ca2+, Na+, AND BUFFER CONCENTRATIONS
[Ca2+]o	Extracellular Ca2+ concentration	2000 μM	Sobie et al. (2008)
[Ca2+]i0	Resting Ca2+ concentration	0.1 μM	Sobie et al. (2008)
[Na+]o	Extracellular Na+ concentration	140 mM	Sobie et al. (2008)
[Na+]i0	Resting Na+ concentration	10 mM	Sobie et al. (2008)
[TN]	Total troponin concentration	70 μM	Shannon et al. (2004)
[ATP]	Total free ATP concentration	260 μM	Cheng et al. (2010)
[CaM]	Total calmodulin concentration	24 μM	Shannon et al. (2004)
[Fluo]	Total Fluo-3 concentration	50 μM	Sobie et al. (2008)
DIFFUSION CONSTANTS (AT 22°C)
Dca	Diffusion coefficient for Ca2+	0.39 μm2 ms−1	Cheng et al. (2010)
DCaFluo	Diffusion coefficient for CaFluo	0.1 μm2 ms−1	Cheng et al. (2010)
DCaATP	Diffusion coefficient for CaATP	0.168 μm2 ms−1	Cheng et al. (2010)
DCaCam	Diffusion coefficient for CaCaM	0.025 μm2 ms−1	Cheng et al. (2010)
RATE COEFFICIENTS AND DISSOCIATION CONSTANTS (AT 22°C)
k + CaTN	Ca2+ on-rate constant for TN	0.04 μM−1 ms−1	Shannon et al. (2004)
k + CaTN	Ca2+ off-rate constant for TN	0.04 ms−1	Shannon et al. (2004)
kDTN	Ca2+ dissociation constant for TN	1 μM	Shannon et al. (2004)
k + CaATP	Ca2+ on-rate constant for CaATP	0.225 μM−1 ms−1	Cheng et al. (2010)
k-CaATP	Ca2+ off-rate constant for CaATP	45 ms−1	Cheng et al. (2010)
kDATP	Ca2+ dissociation constant for ATP	200 μM	Cheng et al. (2010)
k + CaCaM	Ca2+ on-rate constant for Cal	0.125 μM−1 ms−1	Shannon et al. (2004)
k-CaCaM	Ca2+ off-rate constant for Cal	0.2975 ms−1	Shannon et al. (2004)
kDCaM	Ca2+ dissociation constant for Cal	2.38 μM	Shannon et al. (2004)
k + CaFluo	Ca2+ on-rate constant for CaFluo	0.23 μM−1 ms−1	Cheng et al. (2010)
k-CaFluo	Ca2+ off-rate constant for CaFluo	0.17 ms−1	Cheng et al. (2010)
kDFluo	Ca2+ dissociation constant for Fluo	0.739 μM	Cheng et al. (2010)

Table A3

Membrane calcium fluxes parameters.

Symbol	Definition	Value	Reference
L-TYPE Ca2+ CURRENT
a(0mV)	Constant	3.28e−08	Estimated
b(0mV)	Constant	−1.388e−05	Estimated
c(0mV)	Constant	0.002	Estimated
d(0mV)	Constant	−0.136	Estimated
e(0mV)	Constant	3.711	Estimated
a(+50mV)	Constant	3.366e−09	Estimated
b(+50mV)	Constant	−1.426e−06	Estimated
c(+50mV)	Constant	0.000	Estimated
d(+50mV)	Constant	−0.019	Estimated
e(+50mV)	Constant	0.981	Estimated
βCa	Scaling constant	28.5	Estimated
Na+/Ca2+ EXCHANGE CURRENT
Vmax_NCX	Maximum NCX rate	0.207 μM ms−1	Shannon et al. (2004)
KmCaAct	Allosteric constant	0.29 μM	Estimated
K_mCa_o	Extracellular Ca2+ dissociation constant	1.3e−3 μM	Shannon et al. (2004)
K_mCa_i	Intracellular Ca2+ dissociation constant	3.59 μM	Shannon et al. (2004)
K_mNa_o	Extracellular Na+ dissociation constant	87.5e−3 μM	Shannon et al. (2004)
K_mNa_i	Intracellular Ca2+ dissociation constant	12.29e−3 μM	Shannon et al. (2004)
η	Voltage-dependent factor	0.35	Shannon et al. (2004)
ksat	Low potential saturation factor	0.27	Shannon et al. (2004)
βNCX	Scaling constant	28.5	Estimated
SARCOLEMMAL Ca2+ ATPase PUMP
Vmax_Pump	Maximum pump rate	2.2e−3 μM ms−1	Shannon et al. (2004)
KmPump	Half-saturation constant	0.5 μM	Shannon et al. (2004)
nHill	Hill coefficient	1.6	Shannon et al. (2004)
βPump	Scaling constant	28.5	Estimated
SARCOLEMMAL Ca2+ LEAK
Vmax⁡_Leak10mM[Na^+]_i	Conductance	0.001984 μM ms−1	Estimated
Vmax⁡_Leak0mM[Na^+]_i	Conductance	0.004382 μM ms−1	Estimated
βLeak	Scaling constant	28.5	Estimated