Modeling Effects of L-Type Ca2+ Current and Na+-Ca2+ Exchanger on Ca2+ 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 Ca2+ channel (LCC) clustering, and allosteric activation of Na+/Ca2+ exchanger by L-type Ca2+ current affects intracellular Ca2+ 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+/Ca2+ exchanger, sarcolemmal Ca2+ pump, and sarcolemmal Ca2+ leak), and stationary and mobile Ca2+ 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 Ca2+. We obtained parameters from voltage-clamp protocols of L-type Ca2+ current and line-scan recordings of Ca2+ concentration profiles in rabbit cells, in which the sarcoplasmic reticulum is disabled. Our model results agree with experimental measurements of global Ca2+ transient in myocytes loaded with 50 μM Fluo-3. We found that local Ca2+ concentrations within the cytosol and sub-sarcolemma, as well as the local trigger fluxes of Ca2+ crossing the cell membrane, are sensitive to details of t-tubule micro-structure and membrane Ca2+ flux distribution. The model additionally predicts that local Ca2+ trigger fluxes are at least threefold to eightfold higher than the whole-cell Ca2+ trigger flux. We found also that the activation of allosteric Ca2+-binding sites on the Na+/Ca2+ exchanger could provide a mechanism for regulating global and local Ca2+ 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+/Ca2+ exchanger fluxes to intracellular Ca2+ 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 Orchard et al., 2009;Smyrnias et al., 2010). The AP activates and modulates sarcolemmal Ca 2+ fluxes, including fluxes through the L-type Ca 2+ channels (LCC), Na + /Ca 2+ exchangers (NCX), Ca 2+ ATPase pumps, and background sarcolemmal Ca 2+ leak . The entry of Ca 2+ via LCC and NCX triggers the sarcoplasmic reticulum (SR) Ca 2+ release via ryanodine receptors (RyRs). The SR provides Ca 2+ for the troponin C (TnC) myofilament protein, thereby activating and regulating myocyte contraction (Bridge et al., 1990;. Controversy, however, surrounds our understanding of whether the openings of LCCs can activate allosteric Ca 2+ -binding sites on NCX and how the clustering of LCCs, the allosteric catalysis of NCX, and cell surface shape modulate Ca 2+ trigger flux controlling SR Ca 2+ release (Sipido et al., 1997;Litwin et al., 1998;Egger and Niggli, 1999;Inoue and Bridge, 2003;Bers and Ginsburg, 2007;. To investigate relationships between ion fluxes via LCCs and NCXs at voltages corresponding to the early AP plateau,  recently measured the time-dependent Ca 2+ concentration profiles ([Ca 2+ ] i ) in isolated rabbit ventricular myocytes. Pharmaceutical disruption of SR activity enabled them to examine the contributions of LCC and NCX to Ca 2+ trigger flux, which otherwise would have been masked by SR-bound fluxes from RyRs and the SR Ca 2+ ATPase (SERCA). The study suggests that at positive voltages (V m = +30 mV), the trigger flux of Ca 2+ is greater than estimates from a simple summation of the fluxes through LCCs and Ca 2+ entry via reverse NCX mode. The authors hypothesized that openings of LCCs increase local ([Ca 2+ ] i ) near the NCX protein complex, which activates Ca 2+ -binding sites on NCX and results in an increase of Ca 2+ influx via the exchanger.
Mathematical modeling complements experimental studies by enabling the examination of Ca 2+ signaling and excitationcontraction coupling (ECC) in cellular and sub-cellular environments. To this end, whole-cell models have suggested an intimate relationship between Ca 2+ , Na + , K + ionic fluxes, and Ca 2+ transient in rabbit ventricular myocytes Mahajan et al., 2008;Aslanidi et al., 2010). These models permitted examination of the contribution of various cellular components to the evolution of the Ca 2+ 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 Ca 2+ 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;Cheng et al., , 2011Cheng et al., , 2012aLouch et al., 2010;Hatano et al., 2011Hatano et al., , 2012Sato 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.  extended this approach to use recently published ttubule 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 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., , 2009 including wider t-tubule crosssections compared to rodents, and structural variations such as constrictions and flattening . 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 Ca 2+ dynamics during myocyte contraction and relaxation may be species-dependent . In this study we extended the model of , originally developed for rodent myocytes to incorporate structural data and model parameters specific to rabbit ventricular myocytes. We examined local and global Ca 2+ dynamics under the influence of two applied transmembrane voltages (0 and +50 mV). We included equations describing the voltage-dependent background Ca 2+ leak and the sarcolemmal Ca 2+ ATPase pump (Lu et al., 2009). Furthermore, we evaluated a model of NCX allostery proposed by  against experimental data  and tested the hypothesis that the allosteric activation of NCX augments the global trigger flux for SR Ca 2+ 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., , 2009. 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.
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 crosssectional area of these constrictions was reduced to an average of 57.7 ± 27.5% . The volume of the model compartment was estimated to be 0.0282 pL. The compartment membrane area was 15.9 µm 2 where the cell membrane within t-tubule was 7.8 µm 2 and within the external membrane 8.1 µm 2 .

REACTION-DIFFUSION MODEL
The effects of four exogenous and endogenous Ca 2+ 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 subsarcolemmal space (40-50 nm in depth; . The free Ca 2+ 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; I LCC , I NCX , I Pump , and I Leak . Formulation of I NCX was presented as the product of an electrochemical (∆E) and an allosteric factor (Allo; . The maximum NCX rate value (V max_NCX = 0.207 µM/ms) used here was from . The NCX allosteric constant (K mCaAct ) was fitted to 0.29 vs. 0.256 µM in . The K mCaAct value was adjusted to approximate the curvature of +50 mV whole-cell Ca 2+ transient data from .
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;. 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 Ca 2+ 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., 2000Scriven et al., , 2010Takagishi Frontiers in Physiology | Computational Physiology and Medicine The cardiac sarcolemma, including external and t-tubule membranes were visualized using scanning confocal microscopy and labeled in Blender. Localized aggregates of L-type Ca 2+ 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 Ca 2+ entrance and extrusion via the sarcolemma and Ca2 + buffering and diffusion inside the myocyte: LCC, L-type Ca2 + current; NCX, Na + /Ca 2+ exchanger; Pump, membrane Ca2 + ATPase pump; Leak, background sarcolemmal Ca 2+ 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; Ca 2+ 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 Ca 2+ transients were extracted at two featured spots along the scanning line (black spots) and along the t-tubule membrane (green spots). 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. (2009Sachse et al. ( , 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 Ca 2+ pump and leak currents were uniformly distributed along the model surface Lu et al., 2009;Brini and Carafoli, 2011).
In this study, each simulation started with basal cytosolic ([Ca 2+ ] i ) of 0.1 µM, and buffers in equilibrium. The extracellular Ca 2+ 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 . Each current density (I LCC , I NCX , I Pump , and I Leak ) was converted to Ca 2+ flux (see Eq. A7, Supplementary Material) by using the experimentally suggested surface to volume ratio (C m /V cell ∼ 4.54 pF/pL) in adult rabbit ventricular myocytes . To ensure the total Ca 2+ 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 . 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).

GLOBAL Ca 2+ SIGNALS
We validated the model against whole-cell measurements of the L-type Ca 2+ current and Ca 2+ transient at two voltages (0 and +50 mV) in the presence of 50 µM Fluo-3 and [Na + ] i of 10 mM www.frontiersin.org . 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 [Ca 2+ ] i 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 Ca 2+ 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 Ca 2+ dynamics as measured in rabbit ventricular myocytes .

GLOBAL AND SUB-SARCOLEMMAL NCX FLUXES
In agreement with experiment , the model predicts that global NCX flux (I NCXglobal ; computed by averaging local [Ca 2+ ] i levels for the entire compartment) is in reverse mode at both applied voltages with 10 mM [Na + ] i ( Figure 2C). For [Na + ] i of 0 mM, the NCX reverse mode is inactivated and an outward Ca 2+ flux is predicted. The contribution to the global Ca 2+ transient due to NCX was largest at +50 mV and [Ca 2+ ] i monotonically increases with time. Figure 2D shows NCX flux calculated by averaging local [Ca 2+ ] i levels in the sub-sarcolemmal space (I NCXsarc ). In the presence of 10 mM [Na + ] i , our results demonstrate that: (1) the overall scale of I NCXsarc was ∼5.6-fold greater than I NCXglobal ; (2) the increase in I NCXsarc during I LCC upstroke relative to I NCXglobal was much faster; (3) while I NCXglobal monotonically increased over the entire simulation, local peaks in I NCXsarc are predicted at 15 ms. Moreover, in absence of [Na + ] i , the model predicts a sharp reversal in I NCXsarc at 0 mV, whereas I NCXsarc had a more gradual and monotonically decreasing flux.

GLOBAL AND LOCAL Ca 2+ TRIGGER FLUXES
The whole-cell trigger flux reported by Sobie et al. (defined as dF /dt max ) quantifies the maximum inward Ca 2+ flux. By normalizing dF /dt max and I LCC to 1.0, it was suggested that the relative contribution of global NCX flux can be estimated . Here we first calculated our global trigger flux (d[Ca 2+ ] global /dt, Figure 2F) from the predicted global Ca 2+ transient. We then normalized I LCC and d[Ca 2+ ] global /dt based on their maximum values, which both occur at 15 ms for 0 mV. In the presence of 10 mM [Na + ] i , our results show that at +50 mV, d[Ca 2+ ] global /dt max is ∼30% of the value at 0 mV while experimental value is reported as 45%. Furthermore, I LCC 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 + ] i , in agreement with experiment ( Figure 2F Inset vs. Figure 3 in Sobie et al.), there is negligible impact on d[Ca 2+ ] global /dt max at 0 mV (∼10% reduction) and at +50 mV, a ∼30% reduction in d[Ca 2+ ] global /dt max is predicted relative to a measured 50% decline.
The model also provides estimates of local Ca 2+ trigger fluxes within the sub-sarcolemmal space. At 0 mV, the model predicts an approximately threefold increase in d[Ca 2+ ]/dt max near the t-tubule mouth relative to the global trigger flux (Figures 2F,G). Moreover, the local d[Ca 2+ ]/dt max 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[Ca 2+ ] mouth /dt max was predicted in the absence of [Na + ] i .

LOCAL Ca 2+ SIGNALS
In  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 Ca 2+ profiles. The calculated line-scan images and local Ca 2+ time-courses are shown in Figures 3A-H 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 Ca 2+ gradient is predicted in transverse cell direction ( Figure 3G). In addition, results in Figures 3A,D show that local Ca 2+ 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 Ca 2+ transients.
The model also predicts that local Ca 2+ transients in the proximity of a t-tubule ( Figure 1B green spots) differ considerably with respect to local cytosolic Ca 2+ transients. Adjacent to the LCC cluster located at the t-tubule distal end (Figure 1A), the local Ca 2+ transients closely resembles the I LCC profiles with times to peak ∼15 ms followed by a decay to the cytosolic [Ca 2+ ] i 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 [Ca 2+ ] i 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 [Ca 2+ ] i initially increases faster than the cytosolic analog (for t < 20 ms) and thereafter [Ca 2+ ] i increases at the same rate as local cytosolic. We should mention here, that such local Ca 2+ 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 Ca 2+ dynamics is regulated.
The most compelling evidence of the voltage and [Na + ] iinfluenced [Ca 2+ ] i transient is shown in Movie S1 in Supplementary Material. Consistent in all movies is the predicted large and steep Ca 2+ gradient in the narrow sub-sarcolemmal region. In these movies, the LCC clusters are clearly evident as localized regions of enhanced local [Ca 2+ ] i . Comparison of the 0 vs. +50 mV cases (top and bottom rows) demonstrates that at both voltages [Ca 2+ ] i increases heterogeneously in transverse cell direction as suggested by experiment . A spontaneous increase in sub-membrane [Ca 2+ ] i at +50 mV propagating within the cell late in the simulation is predicted also, while initiation and propagation of Ca 2+ gradients at both voltages has not been observed during experiment. Additional interesting results in the presence 10 mM [Na + ] i are that: (1) Ca 2+ 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 Ca 2+ 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 Ca 2+ 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 endto-end CRU distance was measured in rabbits (Dan et al., 2007;. We found that predicted spatial [Ca 2+ ] i 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 Ca 2+ dynamics, we tested how the changes in NCX allosteric constant (K mCaAct ) and maximum exchanger rate (V max_NCX ) affect the results. Our data indicate that local Ca 2+ profiles are highly sensitive to K mCaAct and V max_NCX alterations (Movie S5 in Supplementary Material vs. Movie S1 in Supplementary Material left lower panel).

MODEL GEOMETRY AND L-TYPE Ca 2+ CHANNELS CLUSTERING
Numerical studies of Ca 2+ signaling in rodent cardiomyocytes that consider non-trivial geometries have shown that the spatial Ca 2+ dynamics depends on myocyte micro-structure, including the branching of t-tubules (Lu et al., 2009;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 Ca 2+ transient . Our studies indicate that these factors may contribute to a spatially non-uniform [Ca 2+ ] i distribution Bridge and Sachse, personal communication) while in rats, the measured [Ca 2+ ] i 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 [Ca 2+ ] i 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 Ca 2+ 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., 2009Sachse et al., , 2012. Thus, we placed two LCC clusters within the surface membrane near the ttubule mouth but found that the LCC cluster placement intensified the under-membrane [Ca 2+ ] i 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 [Ca 2+ ] i 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 Ca 2+ 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 Ca 2+ trigger flux ).

ALLOSTERIC CATALYSIS OF Na + -Ca 2+ EXCHANGER AND CA 2+ TRIGGER FLUXES
A unique feature of the presented model was its ability to directly examine the role of catalytic Ca 2+ -binding sites on NCX in regulating local and global Ca 2+ 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 [Na 2+ ] i [I NCXglobal(0mV) ∼ 3 µM/s, I NCXglobal(+50mV) ∼ 8 µM/s, I NCXsarc(0mV) ∼ 9 µM/s, I NCXsarc(+50mV) ∼ 22 µM/s]. At 0 mV, when the LCC flux was dominant, a modest decrease in d[Ca 2+ ] mouth /dt max (2.5 vs. 3 µM/s) was predicted. At +50 mV, when NCX contributes a significantly larger fraction of the total Ca 2+ transient, the effects of NCX allostery were more evident. Assumption of a non-allosteric NCX model resulted in a approximately twofold drop in d[Ca 2+ ] mouth /dt max and lesser inwardly propagating Ca 2+ 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) (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 Ca 2+ entry were localized to a smaller region typical of a cluster, which could then activate Ca 2+ binding to NCX to a greater degree . 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 , 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 Ca 2+ and Ca 2+ -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 + ] i modulate NCX activity independent from any Ca 2+ -dependent allosteric effects (Movie S3 in Supplementary Material). By explicitly including the primary fluxes for Na + and K + (I Na , I NaK ), 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 Ca 2+ fluxes were omitted (Dash et al., 2009;Pradhan et al., 2010). These are Frontiers in Physiology | Computational Physiology and Medicine likely to have important impact on predicted Ca 2+ 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 Ca 2+ 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 Ca 2+ gradients; (2) local trigger fluxes of Ca 2+ are at least threefold to eightfold higher than whole-cell Ca 2+ trigger flux; and (3) the activation of allosteric Ca 2+ -binding sites on Na + /Ca 2+ exchanger could provide a mechanism that regulates local and global Ca 2+ 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 I LCC and I NCX contributions to global and local Ca 2+ trigger fluxes.

ACKNOWLEDGMENTS
The authors thank Eric Sobie, Ryan Szypowski, and Zeyun Yu for helpful suggestions and discussion on the experimental and numerical tools. This work is supported by the National Biomedical Computational Resource (NIH grant 8P41 GM103426-19), an NHLBI Systems Biology Collaboration 1R01HL105242-01 (McCulloch). Dr. McCammon acknowledges additional support from NIH, NSF, Howard Hughes Medical Institute (HHMI), and the Center for Theoretical Biological Physics (CTBP). We acknowledge also funding by NIH HLBI R01 HL094464, the Richard A. and Nora Eccles Fund for Cardiovascular Research, and the Nora Eccles Treadwell Foundation. The work was also partly supported by Center of Excellence grant from the Research Council of Norway to the Center for Biomedical Computing at Simula Research Laboratory.

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 Movie S1 | Effects of voltage and intracellular Na + on spatial Ca 2+ dynamics in the presence of allosteric catalysis and 50 µM Fluo-3.

APPENDIX EXPERIMENTAL AND SIMULATION PROTOCOLS CONFOCAL IMAGING, POST-PROCESSING, AND GEOMETRIC MODELING
Myocytes were superfused with membrane-impermeant dextran conjugated to fluorescein and imaged with a BioRad MRC-1024 laser-scanning confocal microscope . Image data were devonvolved, corrected for background signals, and depth dependent attenuation. Single or groups of ttubules were segmented using region-growing and stored together with an associated segment of outer sarcolemma. A modified Marching cube algorithm was used to create an initial triangulation from the segmented data (Heiden et al., 1993). A GAMer plugin to Blender was used to assign markers for the cell surface, t-tubule, and interior boundaries (Yu et al., 2008a). Tet-Gen was then used to convert the surface mesh into a tetrahedralized mesh suitable for finite-element solution (http://wiasberlin.de/software/tetgen/). The created geometric model contained realistic surface sarcolemma, realistic t-tubule, and its surrounding half-sarcomeres.

CURRENTS AND CALCIUM CONCENTRATION SIMULATION PROTOCOLS
The voltage-clamp protocols (holding potential −50 mV, voltage step to 0 or +50 mV for 200 ms) and whole-cell LCC currents were derived from data in rabbit ventricular myocytes with RyR Ca 2+ release and SERCA Ca 2+ uptake disabled pharmacologically .
Model simulation protocols for [Ca 2+ ] i resemble the protocols applied previously in experimental studies where the cells, incubated with ryanodine and thapsigargin, were loaded with 50 µM Fluo-3, and confocal line-scan recordings of fluorescence ratio (F /F o ) reflecting the changes of Ca 2+ concentration in time were made concurrently . This allowed simultaneous measurements of whole-cell LCC current and the total transmembrane flux of Ca 2+ under conditions promoting Ca 2+ entry via NCX; positive potentials and [Na + ] i of 10 mM Sobie's et al. data indicate: (1) that at V m ≥ 20 mV, dF /dt max was larger with 10 mM [Na + ] i than that measured with 0 mM [Na + ] i (zero Na + prohibits Ca 2+ entry through NCX), suggesting that reverse-mode NCX augments Ca 2+ trigger flux at positive voltages; and (2) that under control conditions (10 mM [Na + ] i and V m ≥ +30 mV) dF /dt max is much greater than a simple sum of fluxes due to LCC current and reverse-mode NCX, suggesting possible activation of an allosteric NCX Ca 2+ -binding site.

Reaction-diffusion equations
where [B m ] represents the concentration of mobile buffer Fluo-3, ATP, or calmodulin; [B s ] is the concentration stationary buffer troponin C.
In the model we assume: (1) isotropic diffusion for Ca 2+ and all mobile buffers ); (2) Ca 2+ binds to Fluo-3, calmodulin, ATP, and TN without cooperativity; (3) the initial total concentrations of the mobile buffers are spatially uniform; (4) the diffusion coefficients of Fluo-3, ATP, or calmodulin with bound Ca 2+ are equal to the diffusion coefficients of free Fluo-3, ATP, or calmodulin.

MODEL CURRENTS
The total Ca 2+ flux across the t-tubule and external membrane is described as: Each current density was converted to Ca 2+ flux by using the experimentally suggested surface to volume ratio C m /V cell ∼4.54 pF/pL in adult rabbit ventricular myocytes (Satoh et al., 1996;: where V ms and S ms are total surface and volume of the real or ideal model compartment and β i scaling constant.

Na + /Ca 2+ exchanger current (I NCX )
The steady-state relation between [Ca 2+ ] i and I NCX is described as in : The three factors are: Allo-allosteric regulation by [Ca 2+ ] i , which is not transported; Bind-competitive binding of Na + and Ca 2+ as substrates; Elect-effects of membrane potential. The latter two combined give an electrochemical term, i.e., ∆E = Bind × Elect.

Sarcolemmal Ca 2+ ATPase pump (I Pump )
The equation and describing I Pump dependence on [Ca 2+ ] i is from :

Sarcolemmal Ca 2+ leak (I Leak )
The equation for I Leak dependence on [Ca 2+ ] i is as in : At rest the Ca 2+ influx via background Ca 2+ leak is adjusted to match the Ca 2+ efflux via NCX and Ca 2+ ATPase so that no net movement across the cell membrane to occurred.

NUMERICAL ALGORITHMS AND SOFTWARE
In finite-element methods, a complex domain needs to be discretized into a number of small elements (such as triangles or tetrahedra). This process is usually referred to as mesh generation (Yu et al., 2008a,b). Although different types of meshes may be generated depending on the numerical solvers to be employed, we restrict ourselves to triangular (surface) and tetrahedral (volumetric) mesh generation as commonly used in biomedical simulation. In the present simulation, the number of finite-element nodes and tetrahedral elements are 4274 and 13,673, respectively.
The non-linear reaction-diffusion system was solved by a finite difference method in time and finite-element method in space using our SMOL software tool (Smoluchowski Solver, http://mccammon.ucsd.edu/smol/) with the time step of 0.5 ms. It takes around 60 min to run 200 ms snapshots with a single Intel Xeon X5355 processor. The SMOL program utilizes libraries from the finite-element tool kit (FETK), which previously has been used in several molecular level studies (Cheng et al., 2007a(Cheng et al., ,b, 2008. One bottleneck for dynamic 3D simulation of non-linear reaction-diffusion system is the computing complexity involved in solving the problem. Here we successfully extended SMOL to solve multiple coupled partial differential equations with non-linear ordinary equations. Multiple tests demonstrate that our SMOL program is quite robust and flexible for various boundary and initial conditions. The simulation results were visualized using GMV mesh viewers (Ortega, 1996). Post-processing and data analyses were implemented by customized Python, MATLAB 2008b (The MathWorks, Natick, MA, USA) scripts and Xmgrace software (Vaught, 1996). A version control system, subversion, was used to monitor the development of software (Collins-Sussman et al., 2002). www.frontiersin.org