Next Article in Journal
Novel Surrogates for Membrane Fouling and the Application of Support Vector Machine in Analyzing Fouling Mechanism
Next Article in Special Issue
Acid-Sensing Ion Channel 2: Function and Modulation
Previous Article in Journal
Membrane Distillation of Saline Water Contaminated with Oil and Surfactants
Previous Article in Special Issue
Cardiac Alternans Occurs through the Synergy of Voltage- and Calcium-Dependent Mechanisms
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Stochastic Spatiotemporal Model of Rat Ventricular Myocyte Calcium Dynamics Demonstrated Necessary Features for Calcium Wave Propagation

by
Tuan Minh Hoang-Trong
1,
Aman Ullah
1,
William Jonathan Lederer
2 and
Mohsin Saleet Jafri
1,2,*
1
School of Systems Biology, Krasnow Institute for Advanced Study, George Mason University, Fairfax, VA 22030, USA
2
Center for Biomedical Engineering and Technology, University of Maryland School of Medicine, Baltimore, MD 21201, USA
*
Author to whom correspondence should be addressed.
Membranes 2021, 11(12), 989; https://doi.org/10.3390/membranes11120989
Submission received: 16 November 2021 / Revised: 12 December 2021 / Accepted: 15 December 2021 / Published: 18 December 2021
(This article belongs to the Special Issue Membrane Channel of Cells)

Abstract

:
Calcium (Ca2+) plays a central role in the excitation and contraction of cardiac myocytes. Experiments have indicated that calcium release is stochastic and regulated locally suggesting the possibility of spatially heterogeneous calcium levels in the cells. This spatial heterogeneity might be important in mediating different signaling pathways. During more than 50 years of computational cell biology, the computational models have been advanced to incorporate more ionic currents, going from deterministic models to stochastic models. While periodic increases in cytoplasmic Ca2+ concentration drive cardiac contraction, aberrant Ca2+ release can underly cardiac arrhythmia. However, the study of the spatial role of calcium ions has been limited due to the computational expense of using a three-dimensional stochastic computational model. In this paper, we introduce a three-dimensional stochastic computational model for rat ventricular myocytes at the whole-cell level that incorporate detailed calcium dynamics, with (1) non-uniform release site placement, (2) non-uniform membrane ionic currents and membrane buffers, (3) stochastic calcium-leak dynamics and (4) non-junctional or rogue ryanodine receptors. The model simulates spark-induced spark activation and spark-induced Ca2+ wave initiation and propagation that occur under conditions of calcium overload at the closed-cell condition, but not when Ca2+ levels are normal. This is considered important since the presence of Ca2+ waves contribute to the activation of arrhythmogenic currents.

1. Introduction

Ca2+ is released from the sarcoplasmic reticulum (SR) intracellular Ca2+ store in the form of Ca2+ sparks at discrete microdomains known as calcium release sites of size ~300 nm in width and ~12 nm in height [1]. The large size of the ventricular myocyte (~120 µm in length, ~20 µm in width and ~15 µm in depth) allows calcium signals to propagate as Ca2+ waves if the conditions are right. The links between calcium sparks and calcium waves have been documented before and has been implicated in the origin of cardiac arrhythmias [2,3,4]. During Ca2+ alternans, a wave-like Ca2+ release was observed during systole in ventricular myocytes by triggering the cells with a small stimulus [5,6]. However, the mechanism and conditions under which calcium sparks can form a propagating calcium wave is still under investigation. The relation between Ca2+ sparks and Ca2+ waves is complicated by the cell structure, asymmetric spatial distribution of RyR2 clusters, anisotropic diffusion of Ca2+, Ca2+ sensitivity of CRUs, etc. Although fluctuations in RyR2 gating have been linked to the development of calcium waves in pathological situations, the links between RyR2 activity and propagation of calcium waves are not well understood. Ca2+ wave generation in myocytes is linked to RyR2 gating and sarcoplasmic reticulum Ca2+ overload [7]. Ca2+ waves in rat ventricular myocytes have been observed during diastole under Ca2+ overload conditions [8,9,10].
How the nature and properties of calcium sparks contribute to the generation of Ca2+ waves and arrhythmias has raised intriguing issues in the field since Ca2+ sparks were discovered [1,11,12,13,14]. Among them, there is a major unresolved problem between the full-width at half-maximal concentration (FWHM) of experimentally recorded Ca2+ sparks (which is about 1.8–2.2 µm) and the pseudo-line scan generated Ca2+ spark from computational modeling (which is about 1.1–1.6 µm) [15,16,17,18,19,20,21]. Izu and co-workers developed a 3D model with spherical geometry using FASCIMILE (AEA Technologies, Harwell, UK) for studying Ca2+ sparks. In their model, in order to reproduce FWHM ~2 µm, they created a supercluster with 4 CRUs at 0.4 µm apart and a very large release of Ca2+ was assumed, i.e., from 2 pA to 5–10 pA for each CRU [17]. A recent 3D model developed using FASCIMILE by Kong and co-workers also produced a FWHM of 1.2 μm [22]. Such models are not supported by physiological measurements and anatomic characterizations of the CRUs.
Here, we have developed a more realistic model. The study of Ca2+-spark-induced Ca2+ wave started with simple deterministic models [4,23] or stochastic models using a simple representation of calcium release units. Keizer and co-workers, using a simple single-channel model of CRUs, with no buffers and very high Ca2+ sensitivity (KCa < 1.5 μM), suggested that the threshold for wave propagation is governed by a single, dimensionless parameter Dτ/d2, where τ is the mean time that a site is open, D is the diffusion constant of Ca2+ and the distance d between 2 sites. The wave is continuous when D τ / d 2 1 , or saltatory if D τ / d 2 1 [24,25]. Izu and co-workers used their large CRU model in a 2D setting, with the fixed open time for Ca2+ release and a more realistic Ca2+ sensitivity (KCa = 15 μM), estimated Po,trigger depending on the distance d and the rate and total mass of the source Ca2+ release [26]. Ramay and co-workers created a 3D model to study the Ca2+ wave, using a T-tubule of 200 nm as a diffusion barrier [27]. The model didn’t incorporate the sarcolemma and its ionic current. In 2010, a 2D model to study spark-induced spark was created by Rovetti and co-workers to form a network of 100 × 100 CRUs with each CRU occupying a domain of 5 × 5 grids [28]. However, their model wasn’t able to produce a proper spark model, i.e., the peak [Ca2+]i < 10 μM, while the estimated level should be around 100 μM [1]. A recently developed 3D model for the rabbit ventricular myocyte by Nivala and co-workers [29] extended the previous model by Rovetti et al. [28] to 3D with 100 × 20 × 10 CRUs. In their model, they assumed (1) fast buffering, (2) the transmembrane voltage senses the calcium everywhere, (3) a simple Hill equation to describe the SERCA pump which does not capture true pump dynamics, (4) no flux into the subspace calcium and the calcium in the dyadic subspace is at equilibrium, and (5) Ca2+-dependent inactivation, however, the model cannot reproduce realistic Ca2+ sparks. The model described in this paper overcomes all these issues.
For atrial cells, due to the lack of the transverse T-tubule system or reduced T-tubule presence and altered properties, the physiological activation of the contraction process requires the propagation of Ca2+ in the form of Ca2+ waves from the sarcolemma to the interior of the cell [30,31]. Thul and co-workers developed a 3D model in which the stochastic nature of Ca2+ release was simply modeled using a threshold model in which a random variable is generated to determine whether Ca2+ is released or not [32].
At the whole-cell level, Shiferaw and colleagues [33] split the cells into subregions that each represented an individual sarcomere and combined them into 1D chains of these identical subregions to approximate the myocyte. Li and co-workers [34] developed the first detailed model by describing the cell as a cylinder with space step 200 nm, with cell volume 20.1 pL and modeling the release of calcium as a point source. Thus, they were unable to investigate the stochastic nature of calcium release which is important to the triggering of calcium waves [35].
We first introduce the 3D spatiotemporal model of the rat ventricular myocyte and then focus on the computational methodology used in model development. The major contribution of this study is two-fold: it investigates the probabilities and conditions for Ca2+ spark-induced sparks to occur, and the probabilities and conditions for spark-induced waves to occur.

2. Materials and Methods

2.1. Model Development

Based on a previously published compartmental model for a rat ventricular myocyte, we created a spatial model [36]. The rat ventricular myocyte is represented as a rectangular solid with dimensions of 120   μ m × 20.8   μ m × 10   μ m , yielding a total cell volume of 24.96 pL, which is consistent with the cell volume with 20,000 CRUs [1,37,38]. To capture the local effect of calcium release, the rectangular cell is divided into grid points of size 0.2   μ m × 0.2   μ m × 0.2   μ m , yielding a total of 3,120,000 grid points, as shown in Figure 1.
The traditional electrophysiological approach of using total membrane capacitance to estimate the surface area and as an indirect index of cell volume is considered not accurate due to the variation in degree of membrane folding. Despite differences in cell volumes and capacitance-to-volume ratios between species, data from optical cell sections obtained with laser scanning confocal microscopy (LSCM) revealed a surprisingly constant membrane capacitance-to-cell volume ratio for each species over a wide range of cell volume using Sprague-Dawley rats (with volumes ranging from 23 pL to 63 pL) [38].The capacitance-volume ratio in rats is quite high (8.42 pF/pL), despite the fact that it varies with age (about 6.76 pF/pL in 3-month-old rats and 8.88 pF/pL in 6-month-old rats). In the Wistar rat (with cell volume 23.6 pL and capacitance 200 pF), a comparable number (8.43 pF/pL) was measured [39]. Other rat models, namely 6.25 pF/pL, substantially underestimate this [40,41].

2.2. Sarcolemma(SL) and T-Tubule Membrane

Given that the specific capacitance of cell membranes is consistent across cell types and species, C sc 1 μ F / cm 2 the total membrane surface area was given as A m = 2.00 × 10 10 cm 2   and the equivalent whole-cell membrane capacitance in the model is 200 pF. The surface membrane and the membrane in the T-tubular system are both part of the sarcolemmal membrane. Depending on the cell type, the thickness of a single cell’s plasma membrane is typically in the range of a few to tens of nanometers [42]; consequently, it can be easily fit into a single grid point of the size employed in the model. If we assume that the surface membrane at each grid point is 1.12 times the area of one side of the grid point, and using the total sarcolemma area above, it requires 445,000 grid points. The total number of grid points for the external surface membrane takes only 195,200 grid points, which means more than half of the membrane (i.e., 56%) belongs to the T-tubular system.
The fraction of total surface membrane area that forms T-tubules is known to be species-dependent. Early measurements revealed that it is approximately 33% in rats [43,44], with approximately 20–50% of the T-tubular membrane being part of the junctional complex with the junctional-SR [45]. Electron microscope measurements of cell capacitance after formamide-induced detubulation indicate a comparable value (32%) [46]. Recent optical measurements, however, revealed that the fraction of folded membrane forming the T-tubular system in rats is substantially higher (about 65%) [47,48]. The T-tubule fractions would be in the range of 50 to 65% based on the development stage-dependent volume-ratio data [37,49]. Thus, the value of 56% in the model is in the range. The discrepancy between the capacitance reduction and detubulation measurement can be explained by incomplete detubulation, structural distortions during the preparation of tissue for electron microscopy and even the smaller specific capacitance (0.56 µF/cm2) in the T-tubule due to high cholesterol level as suggested by [49].
Furthermore, new research showed that, despite the fact that T-tubules leave the surface membrane at Z-lines, only around 60% of the tubular volume occurs near the Z-line [47]. The number of grid points along the Z-line that include the T-tubule in our model is 151,000, which means the number of longitudinal grid points that constitute the branching T-tubules is 101,000 (40.6%). This is close to the value (40%) estimated experimentally in rats [47,50]. Due to the fact that the membrane folds not only in the transverse direction but also in the axial direction, we call it a transverse-axial tubular system (TATS or T-Ax), as suggested by [51] or sarcolemmal Z rete (ZRE) [47]. The TATS is also referred to as the sarcolemmal tubule network [52]. Several studies suggest that disruption or loss of the TATS is associated to a variety of cardiac diseases, and that understanding the link between them is critical for treatment development [53,54,55,56,57,58]. This emphasizes the importance of developing a spatiotemporal model for the ventricular myocyte that incorporates the detailed structure of the TATS to improve our understanding of its function.
In summary, the external sarcolemmal membrane spans the 6 surfaces of the rectangular solid. The TATS is composed of two parts. The T-tubular part is assumed to span along the y-axis of the cell, co-localizing at the calcium release site placements. The longitudinal part of the TATS is added to the network forming 40.6% of the membrane folding.

2.3. Calcium Release Site

The dynamics of calcium and calcium-calmodulin (CaCalm) complex at each release site is given by the equations:
d [ Ca ] ds ( i ) dt = ( J ryr ( i ) J efflux ( i ) + J dhpr ( i ) ) λ ds 2 d [ CaCalm ] ( i ) dt d [ CaSL ] ( i ) dt d [ CaSR ] ( i ) dt
d [ CaCalm ] ( i ) dt = k Calm + ( [ Ca ] ds ( i ) ) 2 ( [ Calm ] T [ CaCalm ] ds ( i ) ) k Calm [ CaCalm ] ds ( i )
d [ CaSL ] ( i ) dt = k SL + ( [ Ca ] ds ( i ) ) ( [ SL ] T [ CaSL ] ds ( i ) ) k SL [ CaSL ] ds ( i )
d [ CaSR ] ( i ) dt = k SR + ( [ Ca ] ds ( i ) ) ( [ SR ] T [ CaSR ] ds ( i ) ) k SR [ CaSR ] ds ( i )
where i represents the index of the calcium release sites (i = 1…20,000). Jryr is the flux of calcium release via RyR2 channels. Jefflux is the flux of calcium from the dyadic subspace (ds) into the cytosol. The volume ratio λ ds = V ds / V myo was incorporated into the fluxes to address the volume differences between model compartments. The model used for this study is based on our published work where the reader can find all the details of the ionic currents and other details [36].

2.4. Spatial Placement of CRUs

Detubulation reduces the calcium current (ICa) by 87% in L-type calcium channels, suggesting a smaller fraction on the surface membrane (Kawai, Hussain et al., 1999) [46]. This was incorporated in the model by using 13% of CRU on the surface membrane. There are two strategies for placing CRUs: (1) uniformly along each dimension at a given distance, or (2) non-uniformly based on the given distribution of nearest distance along the T-tubule. In the first case, a typical assumption is 1.8 µm along × and 0.8 µm along the transversal direction. In the second case, an algorithm was developed to place CRUs based on the distribution between adjacent CRUs along each T-tubule measured previously [59]. An example of the generated location of CRUs on a given depth value is shown in Figure 2.

2.5. Spatial Placement of Na+/Ca2+ Exchangers, SERCA Pump, SR and SL Buffers

The main pathway for Ca2+ extrusion is via the Na+/Ca2+ exchanger (NCX). With immunofluorescence and/or immunoelectron microscopy, the distribution of NCX was unclear, with both even distribution between external membrane surface and the T-tubules in rats and guinea pigs [60] and higher distribution in the T-tubular system in guinea pigs [61]. In another study for rat ventricular myocytes, NCX distributed largely in the T-tubule, yet NCX has not been observed in the dyad [62]. In a more recent study, by measuring the rise in extracellular calcium, ref. [63] showed that the detubulated cell has a smaller rise which suggested a higher distribution of NCX in T-tubules in rats, yet they could not quantify the change. Experiments using formamide-induced detubulation, ref. [64] estimated NCX to be 3–3.5-fold more concentrated in the T-tubules. However, in these studies, the decrease in cell surface area in detubulated cells is in the range 25–32%. Therefore, we focus on modeling NCX with 3 times higher in the T-tubules than in the external sarcolemma.
The transport of Ca2+ from the cytosol of the cardiomyocyte to the lumen of the sarcoplasmic reticulum (SR) is the major mechanism of removing Ca2+, therefore, it plays a major role in the contraction-relaxation cycle of the myocardium. There are different isoforms of SERCA pump that have been found in cardiac myocytes; with some being species-specific. In rat cardiac myocytes, SERCA2a is the major cardiac isoform while the level of SERCA2b is small. In mice, unlike SERCA2b which has a preferential localization around the T-tubules, SERCA2a is distributed transversely and longitudinally in the SR membrane [65,66]. Using immunostaining of rat ventricular myocytes with anti-SERCA2a primary antibody, it shows a uniform striated pattern where the brightest regions on the image are the Z-lines of the sarcomeres [16]. Smith and Keizer modeled this using a hypothetical function that is bell-shaped with the peak at the Z-line. The authors of [67] also found the presence of SERCA2a in the perinuclear region of the cardiac cell, yet the signal is weaker than that in the Z-lines. It is also noted that phospholamban (PLB) interact directly with SERCA2a to inhibit the uptake of Ca2+ back to the SR. Phosphorylation of PLB, which relives the inhibition of PLB on SERCA2a, can occur at two sites: Ser16 (PLB-16P) and Thr17 (PLB-17P), via cAMP-dependent kinases (PKA) and Ca2+/Calm-dependent kinase (CaMKII). Thus, the subcellular distribution of PLB-16P and PLB-17P can be of functional significance. The authors of [67] found that PLB-16B has a higher density at the Z-lines and perinuclear region; while PLB-17 is stronger at the Z-lines, in the intercalated disk region, and the external surface membrane. In essence, PLB-17 is found mainly in the region of importance to the EC-coupling. However, here we only focus on SERCA2a. The role of PLB will be a part of future studies.
In addition to calmodulin and troponin, SR, and SL also play an important role as Ca2+ buffers. The membranes near the Z-line have a greater Ca2+ buffering capacity than those close to the M-line [19]. Thus, in our model, 90% of the SR and SL buffers are assumed near the Z-line.

2.6. Diffusion of Ions

The diffusion of Ca2+ in aqueous solution of physiological ionic strength was estimated at about 700–800 μm2/s [68]. Kushmerick and Polodosky estimated that the diffusion coefficient in muscle is 50 times slower than in aqueous solution, i.e., DCafree = 14 μm2/s [69]. However, recent data suggested the diffusion in cytosol of free (unbuffered) Ca2+ is only 2–2.5× [16], with the value 400 μm2/s in smooth muscle cells and 223 μm2/s from Xenopus laevis oocytes [70,71]. Other studies used DCafree = 300 μm2/s [22,27,29]. The diffusion constant DCafree = 270 μm2/s being used in the model is in the range.
Using FRAP (fluorescence recovery after bleach) of intra-SR Ca2+ indicator Fluo-5N, Wu and Bers estimated the diffusion of Ca2+ in the SR of rabbit cardiac myocytes to be about 60 μm2/s [72]. However, another study estimated a much smaller value, 8–9 μm2/s in rat and guinea pig myocytes [73]. Swietach and co-workers suggested that the slow diffusion for Ca2+ in the SR helps to explain the long recovery time for [Ca2+]jSR within 100–200 ms. In addition, the authors suggested that the high value measured by Wu and Bers did not consider the SR Ca2+ leak during measurement. A recent study by Picht and co-workers supported the previous result in Bers’ group and rejected that claim [74]. A study by Sobie and Lederer supported the results of Bers’ group in that their simulation data can reproduce simulation results with DCasr = 60 μm2/s, rather than with DCasr = 20 μm2/s [75]. In order to explain the slow recovery during Ca2+ sparks, e.g., 161 ms by [76], given the fast diffusion of Ca2+ in the SR, one possible hypothesis is that the duration of SR Ca2+ release, time-to-blink nadir, is longer than time-to-spark peak [75,77]. The authors suggested some irregularities in the structure of the RyR2 cluster, e.g., involving rogue RyR2s.
Another factor for the long recovery time is not only the sharp depletion of [Ca2+] in the jSR compared to the network SR, but also the partial depletion of [Ca2+] in the neighboring nSR grid points. The released Ca2+ includes not only the free calcium but also the calcium that binds to CSQ. The amount of Ca2+ that binds to CSQ in the cardiac muscle was estimated to be about 50% [52]; although this value is likely to vary with species and experimental conditions. Terentyev and co-workers applied an acute, 3-fold reduction in CSQ2 in quiescent rats and found about a 2-fold decrease in SR Ca2+; without change in free intra-SR [Ca2+], which suggested 70% of Ca2+ was bound to CSQ2 [78,79]. Given the depletion of [Ca2+]jSR from 1000 μM to ~100 μM during a Ca2+ spark, as suggested by a previous version of our current model [20] and a recent study [19], and the binding constant between CSQ and Ca2+ is about 400 μM, a significant amount of Ca2+ unbound from CSQ are released. In addition, we also observe a significant Ca2+ reduction in the neighboring nSR (~500 μM). Thus, the recovery time for [Ca2+]jSR should be long enough to provide an adequate amount of Ca2+ to refill the free [Ca2+]jSR, to neighboring nSR, including Ca2+ to bind to [CSQ]. This, not including the possible rate-limiting effect between the jSR and nSR due to geometrical structure, can be used to explain the long recovery time, given the diffusion constant of [Ca2+]SR at 60 μm2/s. Thus, in the model, we chose DCasr = 60 μm2/s.
The diffusion of Fluo-3 constant 90 μm2/s was used in the model. This is in agreement with the measurement [80] that is almost 5-fold larger than what has been measured earlier in skeletal muscle [81]. The role of ATP mobility and its serving as a Ca2+ buffer has not been considered in the model [82]. Based on molecular weight, calmodulin (Calm) is expected to diffuse with about one order lower than that of calcium. Even though other studies have modeled it as a mobile buffer with a small diffusion constant: 20–40 μm2/s, calmodulin is modeled as a stationary buffer in our model.

2.7. Model Formulation in the Spatial Cell

At each grid point, the model contains the concentrations of the chemical species in the myoplasm, and in the network SR. Each grid point contains a myoplasmic fraction and SR fraction (Figure 3). In the myoplasmic volume part of a grid point, the chemical species are calcium, calcium bound to fluorescent indicator dye and other calcium-bound buffers. Depending upon the location of the grid point, i.e., if it resides at the external SL or the T-tubule, there can be other species, e.g., SL buffer. Due to the fast diffusion rate of the membrane potential Vm, it is assumed that the transmembrane potential Vm is spatially uniform distributed across the cell. The dynamics of transmembrane potential is derived from the ionic currents in the form:
dV m dt = 1 C sc ( i = ion I i ) + I app
where I app is the stimulus current. At this scope of the study, it is assumed that there is no ionic exchange between the extracellular media and the inside of the cell. However, for completeness of the model description, all the mathematical equations are described.
i = ion I i = I Na + I dhpr T A m + I K 1 + I Kss + I Ktof + I Ktos + I NCX ¯ + I Na / K + I pmca ¯ + I bNa + I bCa ¯ + I bK
where
I dhpr T = i = indexCRU I dhpr ( i )
I NCX ¯ = i = 1 # grids as membrane I NCX ( i ) # grids as membrane
I bCa ¯ = i = 1 # grids as membrane I bCa ( i ) # grids as membrane
I pmca ¯ = i = 1 # grids as membrane I pmca ( i ) # grids as membrane
In the current model, only the spatial distribution of calcium and calcium-bound species are considered. In the myoplasm, calcium is buffered by the endogenous buffers (calmodulin (Calm), troponin-C (Trpn), SL, SR membrane buffers) or exogeneous buffers (Fluo-3 (F)) with the kinetics based on [83]:
J buffer = J CaTrpn + J CaCM + J CaSR + J CaF + J CaSL
with CaTrpn = calcium-bound to troponin, CaCalm = calcium-bound calmodulin, CaSR = calcium-bound to SR buffer, CaF = calcium-bound to Fluo-3. Of these, calmodulin, the SL, SR membrane buffers and troponin are assumed to be stationary and are described by the following ordinary differential equations:
d [ CaCalm ] dt = J Calm = k Calm + ( [ Ca ] myo ) ( [ Calm ] T [ CaCalm ] myo ) k Calm [ CaCalm ] myo
d [ CaSR ] dt = J CaSR = k SR + ( [ Ca ] myo ) ( [ SR ] T [ CaSR ] myo ) k SR [ CaSR ] myo
d [ CaSL ] dt = J CaSL = k SL + ( [ Ca ] myo ) ( [ SL ] T [ CaSL ] myo ) k SL [ CaSL ] myo
d [ CaTrpn ] dt = J Trpn = k Trpn + ( [ Ca ] myo ) ( [ Trpn ] T [ CaTrpn ] myo ) k Trpn [ CaTrpn ] myo
The exogenous buffer Fluo-3 are mobile buffers and are describe by the following differential equations:
[ CaF ] myo t = J CaF + 2 [ CaF ] myo
[ F ] myo t = J CaF + 2 [ F ] myo
where the fluxes are:
J CaF = k F + ( [ Ca ] myo ) [ F ] k F [ CaF ]
Free calcium and mobile buffers passively diffuse from one grid point to the neighboring ones following Fickian diffusion laws. Combined with the above fluxes, the following partial differential equations describe these model variables:
[ Ca 2 + ] myo t = J ryr J buffer J serca + 2 [ Ca 2 + ] myo
[ Ca 2 + ] nsr t = J serca + J refill + 2 [ Ca 2 + ] nsr  
The calcium in the jSR is replenished from the nSR by a diffusive flux ( J refill ) which has been described by [84]. Even though Ca2+ release via RyR2 from the jSR can spread by more than one grid point, the jSR is treated as a single volume. Thus, the changes in [Ca2+]jSR at the ith release site is described by the following ordinary differential equation:
[ Ca 2 + ] jsr i t = β jsr ( J refill J RyR 2 ) / λ jsr
where fast buffering is assumed in the jSR, and λjsr = Vjsr/Vmyo is the volume fraction.
β jsr = 1 / ( 1 + ( B jsr T × Km jsr ) / ( ( Km jsr + [ Ca 2 + ] jsr i ) 2 ) )

2.8. Computational Methods

The model is fully stochastic in terms of the channel gating of RyR2s and LCCs. The program was mainly written in Fortran and CUDA Fortran, using the CUDA programming toolkit to run on Nvidia Fermi GPU, with some parts written in C++. The Euler method was used to solve the partial-differential equations (PDEs) of diffusion-reaction, and other ordinary differential equations (ODEs). The adaptive time-step ranging from 10 ns to 1 µs was used. When there is some activity, due to channel gating, the time-step will be reduced for numerical stability. The units in the systems are as follows: transmembrane potential—mV, membrane currents—µA/cm2, the ionic concentration—µM (defined based on the corresponding volume), time-second (s). In the spatial model, the forward difference in time and central difference in space derivatives yields:
u i,j,k t + Δ t u i,j,k t Δ t = f ( u i,j,k t , t ) + D x ( u i + 1 , j , k t 2 u i,j,k t + u i 1 , j,k t ( Δ x ) 2 ) + D y ( u i , j + 1 , k t 2 u i,j,k t + u i , j 1 , k t ( Δ y ) 2 )   + D z ( u i , j , k + 1 t 2 u i,j,k t + u i , j , k 1 t ( Δ z ) 2 ) u i,j,k t + Δ t u i,j,k t Δ t   = f ( u i,j,k t , t ) + D x ( u i + 1 , j,k t 2 u i,j,k t + u i 1 , j,k t ( Δ x ) 2 ) + D y ( u i , j + 1 , k t 2 u i,j,k t + u i , j 1 , k t ( Δ y ) 2 )   + D z ( u i , j , k + 1 t 2 u i,j,k t + u i , j , k 1 t ( Δ z ) 2 )
u ijk p + 1 u ijk p Δ t = f ( u p , t )   + D ( u i + 1 , j , k p 2 u ijk p + u i 1 , j , k p ( Δ x ) 2 + u i , j + 1 , k p 2 u ijk p + u i , j 1 , k p ( Δ y ) 2 + u i , j , k + 1 p 2 u ijk p + u i , j , k 1 t ( Δ z ) 2 ) u ijk p + 1 u ijk p Δ t   = f ( u p , t ) + D ( u i + 1 , j , k p 2 u ijk p + u i 1 , j , k p ( Δ x ) 2 + u i , j + 1 , k p 2 u ijk p + u i , j 1 , k p ( Δ y ) 2 + u i , j , k + 1 p 2 u ijk p + u i , j , k 1 t ( Δ z ) 2 )
where u represents the species concentration in a single grid point.
u ijk p + 1 = u ijk p + Δ t { f ( u , t ) + D ( u i + 1 , j , k p 2 u ijk p + u i 1 , j , k p ( Δ x ) 2 + u i , j + 1 , k p 2 u ijk p + u i , j 1 , k p ( Δ y ) 2 + u i , j , k + 1 p 2 u ijk p + u i , j , k 1 t ( Δ z ) 2 ) } u ijk p + 1 = u ijk p + Δ t { f ( u , t ) + D ( u i + 1 , j , k p 2 u ijk p + u i 1 , j , k p ( Δ x ) 2 + u i , j + 1 , k p 2 u ijk p + u i , j 1 , k p ( Δ y ) 2 + u i , j , k + 1 p 2 u ijk p + u i , j , k 1 t ( Δ z ) 2 ) }
For a single cell, the Neumann boundary condition is used, i.e., no flux:
u 0 , j , k   p = u 1 , j , k p   and   u X + 1 , j , k   p = u X , j , k p
u i , 0 , k   p = u i , 1 , k p   and   u i , Y + 1 , k   p = u i , Y , k p
u i , j , 0   p = u 1 , j , 1 p   and   u i , j , Z + 1   p = u X , j , Z p
To handle the extremely large data generated by the 3D model, the HDF5 (hierarchical data format) library was used [85], with the data exported in a compressed format that can be extracted later for further data analysis. The IDL language (data visualization software, 2014) was used to perform all data visualization and data analysis. In this study of calcium waves, the two ions (K+, Na+) were kept constant for the duration of the simulation. Similarly, the corresponding currents were also modeled as uniform [63].
In the 3D model, several 3D data arrays need to be created, e.g., data for cytosolic calcium, myoplasmic calcium, calcium-bound to fluorescence, free fluorescence. However, there are certain data that are not ‘everywhere’ but are only present at a certain number of grid points. Many of them are species that reside on the sarcolemma membranes, e.g., SL buffer, ionic channels (K+, Na+, LCC) and pumps/exchangers (NCX, PMCA, Na+/K+). To save memory and to enhance the performance, a single 3D data array called grid point data is used. Each element of this array is a 32-bit value, where we split them into groups of 4 bits. Each group of bits tells us some information about that grid point.
  • Bit 0–3: index to the array that tells NCX distribution
  • Bit 4–7: index to the array that tells SERCA distribution
  • Bit 8–11: index to the array that tells SR buffer distribution
  • Bit 12–15: index to the array that tells Troponin-C distribution
  • Bit 16–19: index to the array that tells the grid-type (MEMBRANE, INNER-GRIDPOINT, OUTER-GRIDPOINT, NUCLEUS, MITO). Currently, we only use MEMBRANE (both SL and T-tubule) and INNER-GRIDPOINT and OUTER-GRIDPOINT (stencil grid point).
  • Bit 20–31: reserved (for future use)
This significantly reduces the memory usage. From six 3D arrays, we now only need the same number of 1D arrays of a few data elements, and a single 3D array. In the current study, only channels, buffers related to Ca2+ are modeled spatially. The arrays that describe the SL buffer distribution, SR buffer distribution, and NCX distribution contain just one (in the case of uniform distribution), or two (in the case we have high and low distribution) data elements. The value in each array represents the concentration of the corresponding species in each grid point. This is important due to the limited memory on CUDA-capable GPUs, and to reduce the overhead of memory access, where reading data from global memory is very expensive compared to reading data on constant memory of much smaller sizes as bit-operations are very fast.

3. Results

3.1. Ca2+ Transient

A line-scan image along the longitudinal direction during the Ca2+ transient is shown in Figure 4. The current model allows the investigation of the level of Ca2+ directly without using Ca2+-bound fluorescent intensities. Even though the value F/F0 is in agreement with experimental data, compared to the result using the back-calculation formula given by [1], the peak of free Ca2+ is about 3 times higher than the value estimated. This means that the calculated calcium from the experiment underestimated the amount of level free calcium in the cell. This emphasizes that the functional role of local Ca2+ near the Ca2+ release site (at the z-line) in the cell due to the restricted space can regulate different cellular signals. Aside from a snapshot of a pseudo-line scan calcium transient during a Ca2+ transient along the longitudinal direction (A), we also demonstrate the dynamics of free calcium content in (B), the dynamics of calcium-bound fluorescence in (C), and the dynamics of total calcium content in (D).

3.2. Ca2+ Spark Induced Ca2+ Sparks

The current model has the following key features: (1) The gating of ion channels is modeled stochastically with the RyR2 kinetics fitted to give the proper fidelity of spark generation that matches the experimental data, i.e., ~100 sparks/cell/s at rest. (2) The non-uniform distribution of SL, SR buffers and NCX were incorporated in the model to reflect experimental reality. (3) The dynamics of calcium buffering in the subspace is modeled explicitly, without using fast binding assumption. (4) The ultra-fast MCMC method allows us to obtain a better statistic of spark-induced spark compared to other studies which only enable the examination one, two, or a few pairs of CRUs at a time. The first simulation series explores the role of different factors, e.g., [Ca2+]myo, [Ca2+]SR, and CRU distances in spark triggering. The result was collected based on 1416 simulation cases.
Under normal conditions, a sharp transition occurs when the CRU distance goes beyond 0.6 µm (Figure 5) with the probability for triggering activation of a CRU (Po,trigger) at 0.8 µm being only 0.14% and the delay in time being 21 ms on average, as shown in Figure 5A,B. This is in agreement with experimental data where Ca2+ sparks are the local release of calcium that span a restricted space, without affecting the neighboring sites [1]. The increase of cytosolic calcium, to 0.1564 µM, brings the Po,trigger for the CRU at distance 0.6 µm from 45.40% to 58.41%. However, at [Ca2+]myo = 0.4 µM, it significantly increases the Po,trigger which is now 89.12% for CRU at 0.6 µm apart and is 3.56% for the CRU at 0.8 µm apart. Along with this, we see a reduction in the delay time, which reduces from 6.6 ms for the CRU at 0.6 µm apart in Figure 5A,B, to 4.8 ms for the CRU at 0.6 µm apart in Figure 5C,D. A similar result can be achieved when increasing [Ca2+]sr (Figure 6). In particular, we tested under two overload conditions with different levels of [Ca2+]nsr. In Figure 6A,B, the high overload conditions are shown ([Ca2+]myo = 0.156 µM, and [Ca2+]nsr = 1.70 mM), and in Figure 6C,D, the low overload condition is shown ([Ca2+]myo = 0.156 µM, and [Ca2+]nsr = 1.30 mM). In each condition, we measure the probabilities of one CRU triggering the neighboring ones at different distances and the delay upon which condition it occurs under. The main cluster can trigger the neighboring cluster at 0.6 µm apart with high fidelity (almost 99%) and very small delay (<8 ms). This strongly suggests the role of small clusters in bridging the gap between CRUs in generating calcium waves across all tested conditions.
To induce higher SR calcium, ref. [86] applied anti-phospholamban (APL) antibodies. The results showed that calcium sparks increase from ~4 events/100 µm to ~7 events/100 µm and remains at this level during the experiment. Using [Ca2+]o = 10 mM, a 4-fold increase in Ca2+ spark was observed [2]. In the cardiac cell, Ca2+ binds to calsequestrin (CSQ) to directly regulate RyR2 gating from the lumenal side. Even though CSQ can exist in multiple forms (monomers, dimers and multimers), the monomer undertakes the regulatory function in RyR2 [87]; while the role of multimers in RyR2 regulation has not been confirmed other than serving as Ca2+ buffers [88].
During Ca2+ overload, a significant amount of SR Ca2+ is bound to CSQ in oligomer forms. In the current model, lumenal dependency is modeled as a function of [Ca2+]jSR, not [Ca2+/CSQ]. The multiple forms of CSQ, i.e., dimers and multimers, which can function as Ca2+ buffers to hold more SR Ca2+ during Ca2+ overload have not been modeled. With a fixed amount of CSQ at the jSR, we suggest that when [Ca2+]jSR exceeds a certain amount, the dependency on [Ca2+]jSR saturates. Thus, the lumenal function is revised to match the spark frequency from the experimental results at the Ca2+ overload condition, while also allowing a higher [Ca2+]SR that can trigger a Ca2+ wave. The lumenal function is changed from:
Θ = k 0 × [ Ca 2 + ] jSR + k 1
to:
Θ = k 0 × min ( Ca max , [ Ca 2 + ] jSR ) + k 1
with Camax = 1.13 mM being the saturation value of [Ca2+]jSR that was selectedwas to be in agreement with the previous experiment.
The difference in Po,trigger between the two cases is small, as shown in Table 1 and Table A1 in Appendix A, suggesting that the Po,trigger is strongly influenced by the Ca2+ diffusing from the neighboring release site. Based on the measurement in Figure 4C of [59], a significant number of transversal CRUs have nearest neighbors that fall into the 0.7 µm range (21%), which may allow multiple CRUs to be activated on the same z-disc.
Experimental studies using high resolution imaging showed that RyR2s are organized in subclusters that can function as single functional CRUs [77,89,90,91]. In this study, by using a big cluster with 36 RyR2 and two smaller satellite clusters of 15 RyR2 at 0.150 µm apart, the realistic FWHM = 1.85 μm was simulated. In Figure 7A, there is a pseudo-line scan image showing a single calcium spark, this is the result of applying Gaussian noise and incorporating the effect of optical blurring to the calculated Δ F / F 0   which is a dynamic variable in the model. The FWHM of the calcium spark is calculated after applying smoothing (3 × 3 window) filtering. Using the computational model, it enables us to record the underlying calcium dynamics which shows the contributions from multiple subclusters. This is, however, invisible if back-calculating the value from the recorded fluorescent signal occurs, as is done in the experimental protocol. This model, for the first time, shed light onto the unresolved question of the FWHM that has not been replicated in previous models. The role of these satellite clusters in spark-induced spark triggering is tested by modeling each CRU with 3 satellite clusters of size 10 RyR2s at 0.2 µm apart. Figure 8 shows that the satellite clusters help boost the chance for spark-induced sparks. Compared to the result in Figure 6B, the CRU can now trigger a neighboring CRU at a further distance, with a higher probability for the ones at the same distances as used in the earlier case. This is in agreement with the study by the authors of [92], in which they suggested that blocking small, non-spark producing clusters of RyR2 using ruthenium red (RuR) are important to the process of Ca2+ wave propagation.

3.3. Ca2+ Spark Induced Ca2+ Waves

In the study by Cheng and co-workers [2], increasing extracellular from 1 mM to 10 mM elicited a 4-fold increase in Ca2+ sparks, with spark amplitude and spark size increasing 4.1-fold and 1.7-fold accordingly. These macro Ca2+ sparks also served as the sites of wave initiation (65%), indicating that the macrosparks may represent aborted waves involving propagation between release sites. Our study, which is in agreement with these studies, also suggests that the Ca2+ spark observed is the result of calcium release, not from a single CRU but multiple CRUs, and thus the mass of Ca2+ release from multiple CRUs is considered strong enough to trigger the wave in most cases.
The general mechanism for Ca2+ wave propagation is the fire-diffuse-fire model [24,25] with an estimated wave velocity of 67 µm/s with CRUs at d = 2.0 µm apart. For sustained Ca2+ wave propagation, there are two important factors: first the diffusing Ca2+ from one Z-line should be able to activate the CRU in the next Z-line, and then calcium release from this CRU should be able to activate the neighboring ones in the same Z-line, and the process repeats. Given the shorter nearest neighbor distances for CRU on the same Z-line, in Figure 4, the CRUs on the same Z-line should be activated first during a Ca2+ wave initiation. This result is in agreement with the prediction that a single spark is unlikely to trigger a spark at the next Z-line [26,93].
The previous section already provided the Po,trigger at which a Ca2+-spark from one CRU can induce the activation of the neighboring CRU in the same Z-line under the different conditions of [Ca2+]myo, [Ca2+]nsr, different distances and of the presence of satellite RyR2. The second factor contributing wave formation along the longitudinal direction at which the distance between the CRUs is much longer at 1.4–2.2 µm. The question is under what condition can this happen? This is determined by two factors: the number of CRUs activated on one Z-line and the mass of Ca2+ release, which depends critically on the level of Ca2+ overload.
A simulation study placing a few CRUs on the same Z-line that are activated, and another CRU at the other end of the sarcomere, was performed. Instead of modeling the activation of a single CRU, in this case, a model at which a few ‘black’ CRUs are activated, allowed the analysis of the Po,trigger of the ‘blue’ CRU (Figure 9).
Simulations with [Ca2+]SR = 1.3 mM with 9 activated CRUs failed to produce enough Ca2+ release to activated nearby release sites. Given the spatial resolution of the confocal microscope is 10 times larger than the size of a junctional SR (0.06 fL vs. 0.008 fL), it is practically impossible to accurately detect the level of Ca2+ in the jSR. In addition, fluo-5N, the dye being used to estimate [Ca2+]jSR, has the Ca2+ affinity ~400 μM, which is 2–3 times smaller than the diastolic level of [Ca2+]SR. Thus, estimating the level of [Ca2+]SR overload is currently impossible [27]. To estimate the total [Ca2+]SR during a Ca2+ overload, it requires detailed modeling of CSQ in multiple forms which is not available in the current model. Therefore, we focus on the number of activated sites. Nevertheless, we were able to see a full and repetitive Ca2+ wave with [Ca2+]SR = 1.7 mM and [Ca2+]myo = 0.156 μM. This seems to provide a “safety factor” where normal [Ca2+]SR cannot produce Ca2+ activation of adjacent sites to form waves.
Izu and co-workers used a 5–7 supercluster [26]. Each super-CRU releases 10–20 pA, compared to the typical 3 pA current [1,94]. This maps to the number of 15–27 CRUs being fired at the same time to achieve the wave velocity of 126 μm/s. In our model, 7–9 CRUs were enough to trigger a Ca2+ wave. Thus, the spark model in our system, the triggering of a Ca2+ wave during Ca2+ overload is more likely to happen. In addition, Izu et al. (2001) used a two-dimensional space which clearly enhances the Ca2+ wave propagation compared to 3D. The general condition to test is [Ca2+]myo = 0.156 µM, [Ca2+]sr = 1.7 mM, and the statistics were collected based on 144 trial cases.
We calculated the propagation velocity as shown in the case of spark-induced spark, when the source of Ca2+ is from a single release site, the velocity is 30–45 µm/s. Using 9 activated CRUs, the wave velocity in the simulation setting, as shown in Figure 10, is 63.6 µm/s at 1.4 µm apart and then reduced to 48.4 µm/s at 1.8 µm apart. Another simulation series was created where one intermediate cluster of RyR2s was added, Figure 11. When an intermediate RyR2 cluster was added at 0.6 µm apart, not only the probability for trigger Po,trigger increased but also the wave speed was 74 µm/s (for Z-line 2.0 µm apart) and 100 µm/s (for Z-line 1.8 µm apart). This suggested that the wave velocity is strongly dependent on the mass of calcium release, i.e., the number of activated CRUs on one Z-line and the distance to the CRU in the next Z-line and suggests an important role of intermediate RyR2 clusters.
The results presented strengthen the positive amplitude-velocity relationship that has been suggested by previous experiments [95,96]. The estimated wave velocities vary significantly between different studies. The early estimation of wave velocity in cardiac myocytes was 100 µm/s [97]. Recently, by testing different diffusion constants of calcium, [98] estimated the value to be around 78 µm/s. Experimental data by [92] in rabbit ventricular myocytes suggested the range 102.2 ± 4.19 µm/s. In other experiments for rats, the authors have showed that the wave speed strongly depends on the amount of calcium overload and temperature [99,100]. Using [Ca2+]o from 2 mM to 15 mM, the wave speed increases from 33 µm/s to 87 µm/s; while at [Ca2+] = 3 mM, the wave speed was 33 µm/s at 23 °C and 74 ± 19 µm/s at 30 °C. Another experiment in rats confirmed the diversity in spatiotemporal properties of Ca2+ waves, modulated by the Ca2+-loading state [101] with different wave velocities: 84 ± 16 µm/s and 116 ± 29.4 µm/s. This is in agreement with our results, that the wave speed varies with the mass calcium release that elevates basal levels of [Ca2+]i. In addition, we showed the important role of off-Z-lines CRUs in wave propagation. A whole-cell simulation of Ca2+ wave is shown in Figure 12A. The model produced repetitive waves similar to Ca2+ waves in experiments under Ca2+ overload conditions, as shown in Figure 12B (Supplemental Videos S1 and S2). Figure 13 shows the simulated fluorescence and network SR Ca2+ changes accompanying the Ca2+ waves at three example times during the simulation. Note that the fall in the network SR calcium occurs slightly behind the Ca2+ wave front as diffusion in the myoplasm is the key driver of wave propagation due to CICR. It is important to note that the irregular spacing of the release sites was essential for wave propagation because in the model with uniform site placement, waves did not propagate beyond a few sarcomeres.
The role of SERCA in the spark-induced Ca2+ wave was tested, as shown in Figure 14. Compared to Figure 10A, with uniform distribution of SERCA along the longitudinal direction (Figure 14A), Ca2+ waves can propagate farther because less Ca2+ will be sequestered at the release site. When SERCA is reduced, the waves can propagate farther because there is more Ca2+ available for diffusion (Figure 14B). This is in agreement with a previous study that found SERCA inhibition should enhance the wave velocity [102]. Given the fact the cytosolic calcium diffuses faster than calcium in the SR, the inhibition of SERCA should increase the diffusion of calcium to the next release site to trigger the opening of RyR2 there. After applying UV-flash photolysis of ‘caged’ Nmoc-DBHQ or the blocking of SERCA affecting the activation of CRUs (a precursor of the SERCA-inhibitor DBHQ [103]), Kelly and co-workers found a decrease in Ca2+ wave velocity. The authors suggested that calcium overloaded at one CRU, due to the activity of the SERCA pump, leads to RyR2 sensitization ahead of the cytosolic Ca2+ wave, which explains the reason why blocking SERCA decreases the Ca2+ wave [104]. However, this does not follow the CICR mechanism. Based on a computational study, Ramay and co-workers concluded that this ‘sensitization’ only occurs at a slow diffusion of SR Ca2+ [27]. However, this slow diffusion was not found in other studies [74]. The activation of one CRU, during a Ca2+ wave, depends on the Ca2+ released from CRUs, not only from the adjacent Z-line, which in turn depends on the distance, but also from activated CRUs on the same Z-line. Thus, one possible explanation for this when SERCA is blocked is when they measure the velocity, the CRU distance is longer for the Ca2+ to reach and the blocking of SERCA reduces the activated CRUs on the same Z-line, which in turn reduce diffusing Ca2+ to activate the release site.

4. Discussion

A 3D spatiotemporal model has been used in this investigation to study spark-induced sparks and spark-induced Ca2+ waves. The model was developed using available data for rat ventricular myocytes from the literature. This model is superior to other models in many ways. The first is that this realistic spark model has been tested to reproduce proper pump/leak balance in rat ventricular myocytes [84]. The second is the detailed 3D spatiotemporal representation of the cell at the resolution of 200 nm. The third is the realistic representation of different ionic currents on the SL, SR membranes, as well as the non-uniform distribution of the endogenous buffers in the cells. The model’s results are in agreement with other studies that the Ca2+ wave velocity varies upon the condition of the setting, in particular the loading of Ca2+. The important role of intermediate clusters has been shown in the model as well, which has not been investigated before. This is in agreement with experimental data showing that T-tubules are a far more complicated structure that branches not only transversally but also longitudinally [47]. Moreover, the role of satellite clusters was also tested, which implied a heterogeneity in RyR2 cluster arrangements. The widely used assumption is that a single connected RyR2 cluster forms a planar array for each release site. However, the experimental data have showed RyR2 arrangement is more heterogeneous and that functional release sites can be the result of one main RyR2 cluster and a few smaller sized RyR2 clusters at a distance of 100–200 nm apart [77,89,90,91]. By using a model with one central cluster with 2 satellite clusters, the model was able to reproduce the proper FWHM of an experimentally observed Ca2+ spark.
Using the local control model, the results showed how SR Ca2+ can induce Ca2+ waves based on stochastic gating of RyR2 channels. In the current model, we suggest that an activation of 7–9 CRUs is required to trigger a wave from one Z-line to another. However, the probability of activation depends on the distance. There are two factors that may contribute to this increase in spark-induced spark activity leading to wave propagation. The first is the amplification of the Ca2+ signal by the opening of RyR2 channels from these satellite clusters. In addition, it has been suggested that clusters of smaller sizes tend to open more often due to the weak allosteric coupling between the channels [84,105]. Thus, we suggested that the high opening activity of these satellites is the second factor that contribute to this rising in spark-induced-spark phenomenon. The overal effect is quite important as we can now see that a CRU can activate another one at a farther distance, forming the so-called ‘macrosparks’ that have been observed in other experiments [2,14,93,106]. The diffusion constant of ATP as well as CaATP is significant, 168 μm2/s [82]. The space consumed by the T-tubule has not been considered, i.e., it was assumed to be zero. Adding these to the model may help to boost Ca2+ diffusion as a fraction of space is now consumed by the T-tubule and the diffusion of calcium is facilitated by Ca-ATP. Thus, it may help to spread Ca2+, which makes wave initiation and propagation even easier. Other factors that may affect Ca2+ wave or arrhythmias that have not been included in the model are the presence of mitochondria and the nucleus [107,108,109].
The T-tubule structure is complex with recent studies showing it varies for different myocytes in the heart [110]. Atrial myocytes have a lower T-tubule density than ventricular myocytes [111]. Remodeling of the T-tubular network has been observed during disease [112]. Furthermore, during disease there is also remodeling of the T-tubular network post-myocardial infarction and is present during cardiac hypertrophy and heart failure [56,113]. In fact, remodeling of the T-tubular network has been shown to disrupt β-adrenergic signaling during heart failure [114]. Computational modeling has been used to help understand the implications of T-tubule remodeling on cardiac calcium dynamics. Wagner et al. demonstrated, with computational modeling, that the T-tubule remodeling causing reduction SR and T-tubule proximity post myocardial infarction led to delayed subcellular Ca2+ release and action potential prolongation [113]. Other work by Song et al., using a spatial computational model, showed that Ca2+ alternans and Ca2+ wave-mediated arrhythmias increased when the percentage of orphaned Ca2+ release sites was in an intermediate range, but reduced in myocytes exhibiting either a well-organized T-tubule network or low T-tubule density [115]. The model presented here demonstrated the importance of the T-tubule location and non-uniformity of placement for the generation of Ca2+ waves because close proximity of some Ca2+ release sites is needed for the Ca2+ wave nucleation and propagation, consistent with these previous findings.
Previous studies have explored the distance between release sites as being an important factor for Ca2+ wave propagation. Izu and co-workers were able to achieve propagation between z-disks in a cardiac myocyte, but required large Ca2+ sparks to achieve this [26]. Nivala and co-workers demonstrated the importance of distance in a simplified model of Ca2+ release [116]. In a simulation site, Coleman and co-workers observed that the distance between release sites was an important factor in Ca2+ propagation [117].
The importance of rogue RyR2s in SR Ca2+ leaks has been previously described [84,105]. They have also been suggested to be critical for Ca2+ wave propagation [118]. However, in this work Ca2+ waves could propagate between sites at normal SR [Ca2+]. A simulation study by Chen and co-workers observed that rogue ryanodine receptors were needed to get propagation between release sites [119]. However, a single spark could not initiate a Ca2+ wave. Four adjacent release sites had to be triggered to initiate a wave. This condition is not required in our model. It could be because the modeling work presented here paid careful attention to achieving accurate Ca2+ spark dynamics and FWHM.
To understand Ca2+ wave propagation it is essential to capture correct Ca2+ spark dynamics and FWHM. Experimentally observed FWHM of recorded Ca2+ sparks is about 1.8–2.2 µm. Previous computational models have not been able to produce realistic FWHM, only reaching up to ~1.2 µm [15,16,17,18,19,20,21,120]. Izu and co-workers developed a 3D model with spherical geometry using FASCIMILE (AEA Technologies, Harwell, UK) for studying Ca2+ sparks. In their model, in order to reproduce FWHM ~2 µm, they created a supercluster with 4 CRUs at 0.4 µm apart and a very large release of Ca2+ was assumed, i.e., from 2 pA to 5–10 pA for each CRU [17]. A recent 3D model developed using FASCIMILE by Kong and co-workers also produced FWHM of 1.2 μm [22]. The computational model presented here was able to simulate realistic FWHM of 1.85 µm which is sufficient to cover the sarcomere with elevated Ca2+ and ensure contraction.
The current model hasn’t incorporated the different oligomers of CSQ, thus the total free SR [Ca2+] during a calcium wave has not been investigated. However, we were able to produce, for the first time, a repetitive Ca2+ wave at [Ca2+]SR = 1.7 mM and [Ca2+]myo = 0.0156 μM. The simulation showed that the local high elevation of cytosolic calcium, which was significantly underestimated using the back-calculation method, may play an important role in regulating the cellular signals. Lastly, the non-uniform placement of CRU is important for calcium wave initiation and propagation under calcium overload.
Keller and co-workers [104] proposed the concept of RyR2 ‘sensitization’ as a possibility for decreasing Ca2+ wave velocity. We did the test by (1) reducing SERCA during Ca2+ waves, (2) by using a different SR Ca2+ diffusion constant, however the result follows CICR. As we pointed out, the Ca2+ velocity depends not only on the CRU distance but also the role of intermediate CRU’s. Thus, it is possible that the increase/decrease in Ca2+ wave velocity is the result of the presence of an intermediate RyR2 cluster or a longer CRU distance. Another possible explanation is that with SERCA blocking, the number of activated CRUs on one Z-line can be smaller, thus reducing the mass Ca2+ release, which in turn reduces the Ca2+ wave velocity. We thus reject the hypothesis that the Ca2+ diffusion constant is higher on the SR side.
Using the Ca2+-sensitive RyR2’s compatible with experimental measurements, the model has successfully explained the difference in wave speed, the mechanism and conditions at which Ca2+ waves are initiated. In an earlier 2D model, Ca2+-bound fluorescence had not declined to the baseline level after the wave [26]. The model, for the first time, is able to reproduce repetitive Ca2+ waves, without any change to RyR2 channel Ca2+-sensitivity and the SERCA pump. The model is the first of its kind at this level of detail and promises to be able to provide a further understanding of different pathophysiological conditions, such as T-tubule remodeling or diverse channel mutations.

5. Conclusions

Our spatiotemporal model of Ca2+ dynamics simulates realistic spark frequency, dynamics, and FWHM. This is the first model to produce realistic FWHM with physiologically realistic parameters. The model suggests that the robustness of Ca2+ wave propagation depends on Ca2+ release site placement, the density of RyR2 clusters, the presence of non-junctional or rogue RyR2 channels, and the level of the SR Ca2+ load. The model shows that repetitive Ca2+ waves can be produced under Ca2+ overload conditions, similar to those observed in recorded experiments.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/membranes11120989/s1, Video S1: Experimental Calcium Waves) and Video S2: Simulated Calcium Waves.

Author Contributions

Conceptualization, M.S.J. and T.M.H.-T.; methodology, M.S.J. and T.M.H.-T.; software, T.M.H.-T.; validation, M.S.J., A.U., W.J.L. and T.M.H.-T.; formal analysis, T.M.H.-T.; investigation, T.M.H.-T.; resources, M.S.J.; data curation, T.M.H.-T.; writing—original draft preparation, T.M.H.-T.; writing—review and editing, M.S.J., A.U., W.J.L. and T.M.H.-T.; visualization, T.M.H.-T. and A.U.; supervision, M.S.J.; project administration, M.S.J.; funding acquisition, M.S.J. and A.U. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the National Institutes of Health (NIH) grant numbers 5R01HL105239 (W.J.L. & M.S.J.), 5R01HL106059 (W.J.L.), 5U01HL116321 (M.S.J. & W.J.L.), and R01 HL142290 (W.J.L.).

Institutional Review Board Statement

No human subjects were involved.

Informed Consent Statement

No human subjects were involved.

Data Availability Statement

Model codes are publicly available at the Mason Archival Repository Service (MARS) at the following link: Available online: http://mars.gmu.edu/handle/1920/12166 (accessed on 14 December 2021).

Acknowledgments

We would like to thank Brian Hagan for the videos of the experimental Ca2+ waves in cardiac ventricular myocytes under Ca2+ overload conditions.

Conflicts of Interest

The authors declare that there is no conflict of interest.

Appendix A

Appendix A.1. Parameters

Table A1. Parameters in the model.
Table A1. Parameters in the model.
ParameterDefinitionValue
VcellCell volume24.96 (pL)
VmyoMyoplasmic volume50% × Vcell
VnsrNetwork SR volume3.2% × Vcell
VjsrJunctional SR volume0.5% × Vcell
Dyad_heightThe distance between jSR and T-tubule12 nm
Couplon_sizeThe size of the subspace300 nm
X_lenCell dimension120 µm
Y_lenCell dimension20.8 µm
Z_lenCell dimension10.0 µm
Dmyo(x,y,z)Free cytosolic calcium diffusion270 µm2/s
Dnsr(x,y,z)Free SR calcium diffusion60 µm2/s
Ddye(x,y,z)Fluorescence and calcium-bound fluorescence diffusion90 µm2/s
FTTotal fluorescence (Fluo-3)50 µM
konBinding constant for fluorescence80 µM−1.s−1
koffUnbinding constant for fluorescence72 s−1
EccChange in free energy between closed RyR2 pairs−0.872 kBT
EooChange in free energy between open RyR2 pairs−1.15 kBT
kBBoltzmann constant1.381 × 10−23 J/K
TTemperature295.15 K
[K+]oExtracellular potassium concentration5.4 mM
[Na+]oExtracellular sodium concentration140 mM
[Ca2+]oExtracellular calcium concentration1.8 mM
[K+]iCytosolic potassium concentration143.72 mM
[Na+]iCytosolic sodium concentration10.2 mM
[Ca2+]iCytosolic calcium concentration0.096 µM
[Ca2+]srSR calcium concentration1.02 mM
BmyoTTotal myoplasmic buffer concentration3.703026 × 10−2 µM
Km,myoDisassociation constant of myoplasm buffer1.1900 µM
[Trpn]Troponin Ca2+ binding sites concentration
(k+ = 2.37 µM−1.s−1, k = 0.032 s−1)
140 µM
[Calm]The dynamics buffers with Kd = 2.38 µM
(k+ = 30 µM−1.s−1, k = 71.4 s−1)
24 µM
[B]SLThe SL buffer
(k+ = 115 µM−1.s−1, k = 1000 s−1)
750 µM
[B]SRThe SR buffer with Kd=0.86µM
(k+ = 115 µM−1.s−1, k = 100 s−1)
47 µM
ApumpConcentration of SERCA pump300 µM
Kp,myoThe binding affinity of cytosolic calcium to SERCA900 µM
Kp,nsrThe binding affinity of SR calcium to SERCA2150 µM
vTrefillTotal refill rate for 20,000 CRUs3 (1/s/(L-cyt))
vTeffluxTotal efflux rate for 20,000 CRUs120 (1/s/(L-cyt))
iryrSingle channel RyR2 current0.2 (pA)
gbCaBackground Ca2+ conductance2.9 × 10−4 (mS/µF)
gbNaBackground Na+ conductance1.066 × 10−4 (mS/µF)
gbKBackground K+ conductance0 (mS/µF)
gNaNa+ conductance13 (mS/µF)
gK1K1 conductance0.2 (mS/µF)
gKssKss conductance4.21 × 10−2 (mS/µF)
gKtofKtof conductance0.0798 (mS/µF)
gKtosKtos conductance6.29 × 10−2 (mS/µF)
AmCell surface area
IpmcaMaximum PMCA current density0.10 (µA/µF)
IncxMaximum NCX current density750 (µA/µF)
INaKMaximum NaK current density0.88 (µA/µF)
eta_RyRHill coefficient of Ca2+-dependent in RyR2.2
eta_LCCHill coefficient of Ca2+-dependent in LCC2.0
EjCoupling strength energy (for 49 RyRs)0.0714
EccCoupling energy between 2 closed channels−0.78 (kBT)
EooCoupling energy between 2 open channels−1.26 (kBT)

References

  1. Cheng, H.; Lederer, W.J.; Cannell, M.B. Calcium Sparks: Elementary Events Underlying Excitation-Contraction Coupling in Heart Muscle. Science 1993, 262, 740–744. [Google Scholar] [CrossRef]
  2. Cheng, H.; Lederer, M.R.; Lederer, W.J.; Cannell, M.B. Calcium sparks and [Ca2+]i waves in cardiac myocytes. Am. J. Physiol. 1996, 270, C148–C159. [Google Scholar] [CrossRef]
  3. Wier, W.G.; ter Keurs, H.E.; Marban, E.; Gao, W.D.; Balke, C.W. Ca2+ ‘sparks’ and waves in intact ventricular muscle resolved by confocal imaging. Circ. Res. 1997, 81, 462–469. [Google Scholar] [CrossRef]
  4. Lukyanenko, V.; Gyorke, S. Ca2+ sparks and Ca2+ waves in saponin-permeabilized rat ventricular myocytes. J. Physiol. 1999, 521, 575–585. [Google Scholar] [CrossRef] [PubMed]
  5. Diaz, M.E.; Eisner, D.A.; O’Neill, S.C. Depressed ryanodine receptor activity increases variability and duration of the systolic Ca2+ transient in rat ventricular myocytes. Circ. Res. 2002, 91, 585–593. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Diaz, M.E.; O’Neill, S.C.; Eisner, D.A. Sarcoplasmic reticulum calcium content fluctuation is the key to cardiac alternans. Circ. Res. 2004, 94, 650–656. [Google Scholar] [CrossRef] [Green Version]
  7. Petrovič, P.; Valent, I.; Cocherová, E.; Pavelková, J.; Zahradníková, A. Ryanodine receptor gating controls generation of diastolic calcium waves in cardiac myocytes. J. Gen. Physiol. 2015, 145, 489–511. [Google Scholar] [CrossRef] [Green Version]
  8. Orchard, C.H.; Eisner, D.A.; Allen, D.G. Oscillations of intracellular Ca2+ in mammalian cardiac muscle. Nature 1983, 304, 735–738. [Google Scholar] [CrossRef] [PubMed]
  9. Wier, W.G.; Blatter, L.A. Ca2+-oscillations and Ca2+-waves in mammalian cardiac and vascular smooth muscle cells. Cell Calcium 1991, 12, 241–254. [Google Scholar] [CrossRef]
  10. Wier, W.G.; Cannell, M.B.; Berlin, J.R.; Marban, E.; Lederer, W.J. Cellular and subcellular heterogeneity of [Ca2+]i in single heart cells revealed by fura-2. Science 1987, 235, 325–328. [Google Scholar] [CrossRef] [PubMed]
  11. Wang, S.Q.; Stern, M.D.; Ríos, E.; Cheng, H. The quantal nature of Ca2+ sparks and in situ operation of the ryanodine receptor array in cardiac cells. Proc. Natl. Acad. Sci. USA 2004, 101, 3979–3984. [Google Scholar] [CrossRef] [Green Version]
  12. Shen, J.X.; Wang, S.; Song, L.S.; Han, T.; Cheng, H. Polymorphism of Ca2+ sparks evoked from in-focus Ca2+ release units in cardiac myocytes. Biophys. J. 2004, 86, 182–190. [Google Scholar] [CrossRef] [Green Version]
  13. Niggli, E.; Shirokova, N. A guide to sparkology: The taxonomy of elementary cellular Ca2+ signaling events. Cell Calcium 2007, 42, 379–387. [Google Scholar] [CrossRef]
  14. Cheng, H.; Lederer, W.J. Calcium sparks. Physiol. Rev. 2008, 88, 1491–1545. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Pratusevich, V.R.; Balke, C.W. Factors shaping the confocal image of the calcium spark in cardiac muscle cells. Biophys. J. 1996, 71, 2942–2957. [Google Scholar] [CrossRef] [Green Version]
  16. Smith, G.D.; Keizer, J.E.; Stern, M.D.; Lederer, W.J.; Cheng, H. A simple numerical model of calcium spark formation and detection in cardiac myocytes. Biophys. J. 1998, 75, 15–32. [Google Scholar] [CrossRef] [Green Version]
  17. Izu, L.T.; Mauban, J.R.H.; Balke, C.W.; Wier, W.G. Large currents generate cardiac Ca2+ sparks. Biophys. J. 2001, 80, 88–102. [Google Scholar] [CrossRef] [Green Version]
  18. Koh, X.; Srinivasan, B.; Ching, H.S.; Levchenko, A. A 3D Monte Carlo analysis of the role of dyadic space geometry in spark generation. Biophys. J. 2006, 90, 1999–2014. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Hake, J.; Edwards, A.G.; Yu, Z.; Kekenes-Huskey, P.M.; Michailova, A.P.; McCammon, J.A.; Holst, M.J.; Hoshijima, M.; McCulloch, A.D. Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit. J. Physiol. 2012, 590, 4403–4422. [Google Scholar] [CrossRef] [Green Version]
  20. Sobie, E.A.; Dilly, K.W.; dos Santos, C.J.; Lederer, W.J.; Jafri, M.S. Termination of Cardiac Ca2+ Sparks: An Investigative Mathematical Model of Calcium-Induced Calcium Release. Biophys. J. 2002, 83, 59–78. [Google Scholar] [CrossRef] [Green Version]
  21. Laver, D.R.; Kong, C.H.T.; Imtiaz, M.S.; Cannell, M.B. Termination of calcium-induced calcium release by induction decay: An emergent property of stochastic channel gating and molecular scale architecture. J. Mol. Cell. Cardiol. 2013, 54, 98–100. [Google Scholar] [CrossRef]
  22. Kong, C.H.T.; Laver, D.R.; Cannell, M.B. Extraction of Sub-microscopic Ca Fluxes from Blurred and Noisy Fluorescent Indicator Images with a Detailed Model Fitting Approach. PLoS Comput. Biol. 2013, 9, e1002931. [Google Scholar] [CrossRef] [Green Version]
  23. Backx, P.H.; de Tombe, P.P.; Van Deen, J.H.; Mulder, B.J.; ter Keurs, H.E. A model of propagating calcium-induced calcium release mediated by calcium diffusion. J. Gen. Physiol. 1989, 93, 963–977. [Google Scholar] [CrossRef] [Green Version]
  24. Keizer, J.; Smith, G.D. Spark-to-wave transition: Saltatory transmission of calcium waves in cardiac myocytes. Biophys. Chem. 1998, 72, 87–100. [Google Scholar] [CrossRef]
  25. Keizer, J.; Smith, G.D.; Ponce-Dawson, S.; Pearson, J.E. Saltatory propagation of Ca2+ waves by Ca2+ sparks. Biophys. J. 1998, 75, 595–600. [Google Scholar] [CrossRef] [Green Version]
  26. Izu, L.T.; Wier, W.G.; Balke, C.W. Evolution of Cardiac Calcium Waves from Stochastic Calcium Sparks. Biophys. J. 2001, 80, 103–120. [Google Scholar] [CrossRef] [Green Version]
  27. Ramay, H.R.; Jafri, M.S.; Lederer, W.J.; Sobie, E.A. Predicting local SR Ca2+ dynamics during Ca2+ wave propagation in ventricular myocytes. Biophys. J. 2010, 98, 2515–2523. [Google Scholar] [CrossRef] [Green Version]
  28. Rovetti, R.; Cui, X.; Garfinkel, A.; Weiss, J.N.; Qu, Z. Spark-induced sparks as a mechanism of intracellular calcium alternans in cardiac myocytes. Circ. Res. 2010, 106, 1582–1591. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Nivala, M.; de Lange, E.; Rovetti, R.; Qu, Z. Computational modeling and numerical methods for spatiotemporal calcium cycling in ventricular myocytes. Front. Physiol. 2012, 3, 114. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Berlin, J.R. Spatiotemporal changes of Ca2+ during electrically evoked contractions in atrial and ventricular cells. Am. J. Physiol. 1995, 269, H1165–H1170. [Google Scholar] [CrossRef] [PubMed]
  31. Kockskamper, J.; Sheehan, K.A.; Bare, D.J.; Lipsius, S.L.; Mignery, G.A.; Blatter, L.A. Activation and propagation of Ca2+ release during excitation-contraction coupling in atrial myocytes. Biophys. J. 2001, 81, 2590–2605. [Google Scholar] [CrossRef] [Green Version]
  32. Thul, R.; Coombes, S.; Roderick, H.L.; Bootman, M.D. Subcellular calcium dynamics in a whole-cell model of an atrial myocyte. Proc. Natl. Acad. Sci. USA 2012, 109, 2150–2155. [Google Scholar] [CrossRef] [Green Version]
  33. Shiferaw, Y.; Karma, A. Turing instability mediated by voltage and calcium diffusion in paced cardiac cells. Proc. Natl. Acad. Sci. USA 2006, 103, 5670–5675. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Li, P.; Lancaster, M.; Holden, A.V. A Three Dimensional Ventricular E-Cell (3Dv E-Cell) with Stochastic Intracellular Ca2+ Handling. In International Conference on Functional Imaging and Modeling of the Heart; Springer: Berlin/Heidelberg, Germany, 2007; pp. 180–189. [Google Scholar]
  35. Soeller, C.; Jayasinghe, I.D.; Li, P.; Holden, A.V.; Cannell, M.B. Three-dimensional high-resolution imaging of cardiac proteins to construct models of intracellular Ca2+ signalling in rat ventricular myocytes. Exp. Physiol. 2009, 94, 496–508. [Google Scholar] [CrossRef]
  36. Hoang-Trong, M.T.; Ullah, A.; Lederer, W.J.; Jafri, M.S. Cardiac Alternans Occurs through the Synergy of Voltage- and Calcium-Dependent Mechanisms. Membranes 2021, 11, 794. [Google Scholar] [CrossRef] [PubMed]
  37. Satoh, H.; Delbridge, L.M.; Blatter, L.A.; Bers, D.M. Surface:volume relationship in cardiac myocytes studied with confocal microscopy and membrane capacitance measurements: Species-dependence and developmental effects. Biophys. J. 1996, 70, 1494–1504. [Google Scholar] [CrossRef] [Green Version]
  38. Delbridge, L.M.; Satoh, H.; Yuan, W.; Bassani, J.W.; Qi, M.; Ginsburg, K.S.; Samarel, A.M.; Bers, D.M. Cardiac myocyte volume, Ca2+ fluxes, and sarcoplasmic reticulum loading in pressure-overload hypertrophy. Am. J. Physiol. 1997, 272, H2425–H2435. [Google Scholar] [CrossRef] [PubMed]
  39. Swift, F.; Stromme, T.A.; Amundsen, B.; Sejersted, O.M.; Sjaastad, I. Slow diffusion of K+ in the T tubules of rat cardiomyocytes. J. Appl. Physiol. 2006, 101, 1170–1176. [Google Scholar] [CrossRef] [Green Version]
  40. Padmala, S.; Demir, S.S. Computational model of the ventricular action potential in adult spontaneously hypertensive rats. J. Cardiovasc. Electrophysiol. 2003, 14, 990–995. [Google Scholar] [CrossRef] [PubMed]
  41. Pandit, S.V.; Clark, R.B.; Giles, W.R.; Demir, S.S. A mathematical model of action potential heterogeneity in adult rat left ventricular myocytes. Biophys. J. 2001, 81, 3029–3051. [Google Scholar] [CrossRef] [Green Version]
  42. Andersen, O.S.; Koeppe, R.E., II. Bilayer thickness and membrane protein function: An energetic perspective. Annu. Rev. Biophys. Biomol. Struct. 2007, 36, 107–130. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Page, E. Quantitative ultrastructural analysis in cardiac membrane physiology. Am. J. Physiol. 1978, 235, C147–C158. [Google Scholar] [CrossRef] [Green Version]
  44. Stewart, J.M.; Page, E. Improved stereological techniques for studying myocardial cell growth: Application to external sarcolemma, T system, and intercalated disks of rabbit and rat hearts. J. Ultrastruct. Res. 1978, 65, 119–134. [Google Scholar] [CrossRef]
  45. Page, E.; Surdyk-Droske, M. Distribution, surface density, and membrane area of diadic junctional contacts between plasma membrane and terminal cisterns in mammalian ventricle. Circ. Res. 1979, 45, 260–267. [Google Scholar] [CrossRef] [Green Version]
  46. Kawai, M.; Hussain, M.; Orchard, C.H. Excitation-contraction coupling in rat ventricular myocytes after formamide-induced detubulation. Am. J. Physiol. 1999, 277, H603–H609. [Google Scholar] [CrossRef] [PubMed]
  47. Soeller, C.; Cannell, M.B. Examination of the Transverse Tubular System in Living Cardiac Rat Myocytes by 2-Photon Microscopy and Digital Image Processing Techniques. Circ. Res. 1999, 84, 266–275. [Google Scholar] [CrossRef]
  48. Brette, F.; Orchard, C.H. No apparent requirement for neuronal sodium channels in excitation-contraction coupling in rat ventricular myocytes. Circ. Res. 2006, 98, 667–674. [Google Scholar] [CrossRef] [Green Version]
  49. Pasek, M.; Brette, F.; Nelson, A.; Pearce, C.; Qaiser, A.; Christe, G.; Orchard, C.H. Quantification of t-tubule area and protein distribution in rat cardiac ventricular myocytes. Prog. Biophys. Mol. Biol. 2008, 96, 244–257. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  50. Brette, F.; Orchard, C. T-tubule function in mammalian cardiac myocytes. Circ. Res. 2003, 92, 1182–1192. [Google Scholar] [CrossRef]
  51. Forbes, M.S.; Hawkey, L.A.; Sperelakis, N. The transverse-axial tubular system (TATS) of mouse myocardium: Its morphology in the developing and adult animal. Am. J. Anat. 1984, 170, 143–162. [Google Scholar] [CrossRef]
  52. Bers, D.M. Excitation-Contraction Coupling and Cardiac Contractile Force, 2nd ed.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2001. [Google Scholar]
  53. Wei, S.; Guo, A.; Chen, B.; Kutschke, W.; Xie, Y.P.; Zimmerman, K.; Weiss, R.M.; Anderson, M.E.; Cheng, H.; Song, L.S. T-tubule remodeling during transition from hypertrophy to heart failure. Circ. Res. 2010, 107, 520–531. [Google Scholar] [CrossRef] [Green Version]
  54. Polakova, E.; Sobie, E.A. Alterations in T-tubule and dyad structure in heart disease: Challenges and opportunities for computational analyses. Cardiovasc. Res. 2013, 98, 233–239. [Google Scholar] [CrossRef] [Green Version]
  55. Heinzel, F.R.; Bito, V.; Biesmans, L.; Wu, M.; Detre, E.; von Wegner, F.; Claus, P.; Dymarkowski, S.; Maes, F.; Bogaert, J.; et al. Remodeling of T-tubules and reduced synchrony of Ca2+ release in myocytes from chronically ischemic myocardium. Circ. Res. 2008, 102, 338–346. [Google Scholar] [CrossRef] [Green Version]
  56. Wagner, E.; Lauterbach, M.A.; Kohl, T.; Westphal, V.; Williams, G.S.B.; Steinbrecher, J.H.; Streich, J.-H.; Korff, B.; Tuan, H.-T.M.; Hagen, B.; et al. Stimulated emission depletion live-cell super-resolution imaging shows proliferative remodeling of T-tubule membrane structures after myocardial infarction. Circ. Res. 2012, 111, 402–414. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Lyon, A.R.; MacLeod, K.T.; Zhang, Y.; Garcia, E.; Kanda, G.K.; Lab, M.J.; Korchev, Y.E.; Harding, S.E.; Gorelik, J. Loss of T-tubules and other changes to surface topography in ventricular myocytes from failing. Proc. Natl. Acad. Sci. USA 2009, 106, 6854–6859. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Guo, A.; Zhang, C.; Wei, S.; Chen, B.; Song, L.S. Emerging mechanisms of T-tubule remodelling in heart failure. Cardiovasc. Res. 2013, 98, 204–215. [Google Scholar] [CrossRef] [Green Version]
  59. Chen-Izu, Y.; Mcculle, S.L.; Ward, C.W.; Soeller, C.; Allen, B.M.; Rabang, C.; Cannell, M.B.; Balke, C.W.; Izu, L.T. Three-dimensional distribution of ryanodine receptor clusters in cardiac myocytes. Biophys. J. 2006, 91, 1–13. [Google Scholar] [CrossRef] [Green Version]
  60. Kieval, R.S.; Bloch, R.J.; Lindenmayer, G.E.; Ambesi, A.; Lederer, W.J. Immunofluorescence localization of the Na-Ca exchanger in heart cells. Am. J. Physiol. 1992, 263, C545–C550. [Google Scholar] [CrossRef] [PubMed]
  61. Frank, J.S.; Mottino, G.; Reid, D.; Molday, R.S.; Philipson, K.D. Distribution of the Na+-Ca2+ exchange protein in mammalian cardiac myocytes: An immunofluorescence and immunocolloidal gold-labeling study. J. Cell Biol. 1992, 117, 337–345. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  62. Scriven, D.R.; Dan, P.; Moore, E.D. Distribution of proteins implicated in excitation-contraction coupling in rat ventricular myocytes. Biophys. J. 2000, 79, 2682–2691. [Google Scholar] [CrossRef] [Green Version]
  63. Yang, Z.; Pascarel, C.; Steele, D.S.; Komukai, K.; Brette, F.; Orchard, C.H. Na+-Ca2+ exchange activity is localized in the T-tubules of rat ventricular myocytes. Circ. Res. 2002, 91, 315–322. [Google Scholar] [CrossRef] [Green Version]
  64. Despa, S.; Brette, F.; Orchard, C.H.; Bers, D.M. Na/Ca exchange and Na/K-ATPase function are equally concentrated in transverse tubules of rat ventricular myocytes. Biophys. J. 2003, 85, 3388–3396. [Google Scholar] [CrossRef] [Green Version]
  65. Periasamy, M.; Huke, S. SERCA pump level is a critical determinant of Ca2+ homeostasis and cardiac contractility. J. Mol. Cell. Cardiol. 2001, 33, 1053–1063. [Google Scholar] [CrossRef]
  66. Greene, A.L.; Lalli, M.J.; Ji, Y.; Babu, G.J.; Grupp, I.; Sussman, M.; Periasamy, M. Overexpression of SERCA2b in the heart leads to an increase in sarcoplasmic reticulum calcium transport function and increased cardiac contractility. J. Biol. Chem. 2000, 275, 24722–24727. [Google Scholar] [CrossRef] [Green Version]
  67. Drago, G.A.; Colyer, J.; Lederer, W.J. Immunofluorescence Localization of SERCA2a and the Phosphorylated Forms of Phospholamban in Intact Rat Cardiac Ventricular Myocytes. Ann. N. Y. Acad. Sci. 1998, 853, 273–279. [Google Scholar] [CrossRef]
  68. Wang, J.H. Tracer-diffusion in Liquids. IV. Self-diffusion of Calcium Ion and Chloride Ion in Aqueous Calcium Chloride Solutions. J. Am. Chem. Soc. 1953, 75, 1769–1770. [Google Scholar] [CrossRef]
  69. Kushmerick, M.J.; Podolsky, R.J. Ionic mobility in muscle cells. Science 1969, 166, 1297–1298. [Google Scholar] [CrossRef] [PubMed]
  70. Kargacin, G.; Fay, F.S. Ca2+ movement in smooth muscle cells studied with one- and two-dimensional diffusion models. Biophys. J. 1991, 60, 1088–1100. [Google Scholar] [CrossRef] [Green Version]
  71. Allbritton, N.L.; Meyer, T.; Stryer, L. Range of messenger action of calcium ion and inositol 1,4,5-trisphosphate. Science 1992, 258, 1812–1815. [Google Scholar] [CrossRef]
  72. Wu, X.; Bers, D.M. Sarcoplasmic reticulum and nuclear envelope are one highly interconnected Ca2+ store throughout cardiac myocyte. Circ. Res. 2006, 99, 283–291. [Google Scholar] [CrossRef] [Green Version]
  73. Swietach, P.; Spitzer, K.W.; Vaughan-Jones, R.D. Ca2+-mobility in the sarcoplasmic reticulum of ventricular myocytes is low. Biophys. J. 2008, 95, 1412–1427. [Google Scholar] [CrossRef] [Green Version]
  74. Picht, E.; Zima, A.V.; Shannon, T.R.; Duncan, A.M.; Blatter, L.A.; Bers, D.M. Dynamic calcium movement inside cardiac sarcoplasmic reticulum during release. Circ. Res. 2011, 108, 847–856. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  75. Sobie, E.A.; Lederer, W.J. Dynamic local changes in sarcoplasmic reticulum calcium: Physiological and pathophysiological roles. J. Mol. Cell. Cardiol. 2012, 52, 304–311. [Google Scholar] [CrossRef] [Green Version]
  76. Zima, A.V.; Picht, E.; Bers, D.M.; Blatter, L.A. Termination of cardiac Ca2+ sparks: Role of intra-SR [Ca2+], release flux, and intra-SR Ca2+ diffusion. Circ. Res. 2008, 103, e105–e115. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Brochet, D.X.P.; Yang, D.; Di Maio, A.; Lederer, W.J.; Franzini-Armstrong, C.; Cheng, H. Ca2+ blinks: Rapid nanoscopic store calcium signaling. Proc. Natl. Acad. Sci. USA 2005, 102, 3099–3104. [Google Scholar] [CrossRef] [Green Version]
  78. Terentyev, D.; Viatchenko-Karpinski, S.; Gyorke, I.; Volpe, P.; Williams, S.C.; Gyorke, S. Calsequestrin determines the functional size and stability of cardiac intracellular calcium stores: Mechanism for hereditary arrhythmia. Proc. Natl. Acad. Sci. USA 2003, 100, 11759–11764. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  79. Kubalova, Z.; Gyorke, I.; Terentyeva, R.; Viatchenko-Karpinski, S.; Terentyev, D.; Williams, S.C.; Gyorke, S. Modulation of cytosolic and intra-sarcoplasmic reticulum calcium waves by calsequestrin in rat cardiac myocytes. J. Physiol. 2004, 561, 515–524. [Google Scholar] [CrossRef] [PubMed]
  80. Michailova, A.; Del Principe, F.; Egger, M.; Niggli, E. Spatiotemporal Features of Ca2+ Buffering and Diffusion in Atrial Cardiac Myocytes with Inhibited Sarcoplasmic Reticulum. Biophys. J. 2002, 83, 3134–3151. [Google Scholar] [CrossRef] [Green Version]
  81. Harkins, A.B.; Kurebayashi, N.; Baylor, S.M. Resting myoplasmic free calcium in frog skeletal muscle fibers estimated with fluo-3. Biophys. J. 1993, 65, 865–881. [Google Scholar] [CrossRef] [Green Version]
  82. Baylor, S.M.; Hollingworth, S. Model of sarcomeric Ca2+ movements, including ATP Ca2+ binding and diffusion, during activation of frog skeletal muscle. J. Gen. Physiol. 1998, 112, 297–316. [Google Scholar] [CrossRef] [Green Version]
  83. Loughrey, C.M.; MacEachern, K.E.; Cooper, J.; Smith, G.L. Measurement of the dissociation constant of Fluo-3 for Ca2+ in isolated rabbit cardiomyocytes using Ca2+ wave characteristics. Cell Calcium 2003, 34, 1–9. [Google Scholar] [CrossRef]
  84. Williams, G.S.B.; Chikando, A.C.; Tuan, H.-T.M.; Sobie, E.A.; Lederer, W.J.; Jafri, M.S. Dynamics of Calcium Sparks and Calcium Leak in the Heart. Biophys. J. 2011, 101, 1287–1296. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  85. HDF-Group. Available online: http://www.hdfgroup.org/solutions/HDF5/ (accessed on 14 December 2021).
  86. Lukyanenko, V.; Viatchenko-Karpinski, S.; Smirnov, A.; Wiesner, T.F.; Györke, S. Dynamic regulation of sarcoplasmic reticulum Ca2+ content and release by luminal Ca2+-sensitive leak in rat ventricular myocytes. Biophys. J. 2001, 81, 785–798. [Google Scholar] [CrossRef] [Green Version]
  87. Qin, J.; Valle, G.; Nani, A.; Nori, A.; Rizzi, N.; Priori, S.G.; Volpe, P.; Fill, M. Luminal Ca2+ regulation of single cardiac ryanodine receptors: Insights provided by calsequestrin and its mutants. J. Gen. Physiol. 2008, 131, 325–334. [Google Scholar] [CrossRef] [Green Version]
  88. Gyorke, S.; Stevens, S.C.; Terentyev, D. Cardiac calsequestrin: Quest inside the SR. J. Physiol. 2009, 587, 3091–3094. [Google Scholar] [CrossRef]
  89. Brochet, D.X.P.; Lederer, W.J. Decomposition of a Calcium Spark in Cardiac Myocytes. Biophys. J. 2014, 104, 483a. [Google Scholar] [CrossRef] [Green Version]
  90. Macquaide, N.; Tuan, H.-T.M.; Hotta, J.-I.; Sempels, W.; Lenaerts, I.; Holemans, P.; Hofkens, J.; Jafri, M.S.; Willems, R.; Sipido, K.R. Functional consequences of RyR cluster fragmentation and redistribution in persistent atrial fibrillation. Cardiovasc. Res. 2014, 108, 387–398. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  91. Baddeley, D.; Jayasinghe, I.D.; Lam, L.; Rossberger, S.; Cannell, M.B.; Soeller, C. Optical single-channel resolution imaging of the ryanodine receptor distribution in rat cardiac myocytes. Proc. Natl. Acad. Sci. USA 2009, 106, 22275–22280. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  92. MacQuaide, N.; Ramay, H.R.; Sobie, E.A.; Smith, G.L. Differential sensitivity of Ca2+ wave and Ca2+ spark events to ruthenium red in isolated permeabilised rabbit cardiomyocytes. J. Physiol. 2010, 588, 4731–4742. [Google Scholar] [CrossRef]
  93. Parker, I.; Zang, W.-J.; Wier, W.G. Ca2+ sparks involving multiple Ca2+ release sites along Z-lines in rat heart cells. J. Physiol. 1996, 497, 31–38. [Google Scholar] [CrossRef] [Green Version]
  94. Blatter, L.A.; Huser, J.; Rios, E. Sarcoplasmic reticulum Ca2+ release flux underlying Ca2+ sparks in cardiac muscle. Proc. Natl. Acad. Sci. USA 1997, 94, 4176–4181. [Google Scholar] [CrossRef] [Green Version]
  95. Trafford, A.W.; Lipp, P.; O’Neill, S.C.; Niggli, E.; Eisner, D.A. Propagating calcium waves initiated by local caffeine application in rat ventricular myocytes. J. Physiol. 1995, 489, 319–326. [Google Scholar] [CrossRef] [Green Version]
  96. Smith, G.L.; O’Neill, S.C. A comparison of the effects of ATP and tetracaine on spontaneous Ca2+ release from rat permeabilised cardiac myocytes. J. Physiol. 2001, 534, 37–47. [Google Scholar] [CrossRef]
  97. Jaffe, L.F. The path of calcium in cytosolic calcium oscillations: A unifying hypothesis. Proc. Natl. Acad. Sci. USA 1991, 88, 9883–9887. [Google Scholar] [CrossRef] [Green Version]
  98. Swietach, P.; Spitzer, K.W.; Vaughan-Jones, R.D. Modeling calcium waves in cardiac myocytes: Importance of calcium diffusion. Front. Biosci. 2010, 15, 661–680. [Google Scholar] [CrossRef] [Green Version]
  99. Kort, A.A.; Capogrossi, M.C.; Lakatta, E.G. Frequency, amplitude, and propagation velocity of spontaneous Ca++-dependent contractile waves in intact adult rat cardiac muscle and isolated myocytes. Circ. Res. 1985, 57, 844–855. [Google Scholar] [CrossRef] [Green Version]
  100. Stuyvers, B.D.; Boyden, P.A.; Keurs, H.E.D.J.T. Calcium Waves. Physiological Relevance in Cardiac Function. Circ. Res. 2000, 86, 1016–1018. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  101. Kaneko, T.; Tanaka, H.; Oyamada, M.; Kawata, S.; Takamatsu, T. Three distinct types of Ca2+ waves in Langendorff-perfused rat heart revealed by real-time confocal microscopy. Circ. Res. 2000, 86, 1093–1099. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  102. Jafri, M.S. On the Roles of Ca2+ Diffusion, Ca2 + Buffers, and the Endoplasmic Reticulum in IP3-lnduced Ca2+ Waves. Biophys. J. 1995, 69, 2139–2153. [Google Scholar] [CrossRef] [Green Version]
  103. Rossi, F.M.; Kao, J.P. Nmoc-DBHQ, a new caged molecule for modulating sarcoplasmic/endoplasmic reticulum Ca2+ ATPase activity with light flashes. J. Biol. Chem. 1997, 272, 3266–3271. [Google Scholar] [CrossRef] [Green Version]
  104. Keller, M.; Kao, J.P.; Egger, M.; Niggli, E. Calcium waves driven by "sensitization" wave-fronts. Cardiovasc. Res. 2007, 74, 39–45. [Google Scholar] [CrossRef]
  105. Sobie, E.A.; Guatimosim, S.; Gómez-Viquez, L.; Song, L.-S.; Hartmann, H.; Jafri, S.M.; Lederer, W.J. The Ca2+ leak paradox and rogue ryanodine receptors: SR Ca2+ efflux theory and practice. Prog. Biophys. Mol. Biol. 2006, 90, 172–185. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  106. Shah, M.; Sikkel, M.B.; Desplantez, T.; Collins, T.P.; O’Gara, P.; Harding, S.E.; Lyon, A.R.; Macleod, K.T. Flecainide reduces wave frequency and mean spark amplitude in isolated rat ventricular cardiomyocytes. Cardiovasc. Res. 2012, 93, S90. [Google Scholar]
  107. Dedkova, E.N.; Blatter, L.A. Mitochondrial Ca2+ and the heart. Cell Calcium 2008, 44, 77–91. [Google Scholar] [CrossRef] [PubMed]
  108. Aon, M.A.; Cortassa, S.; Akar, F.G.; Brown, D.A.; Zhou, L.; O’Rourke, B. From mitochondrial dynamics to arrhythmias. Int. J. Biochem. Cell Biol. 2009, 41, 1940–1948. [Google Scholar] [CrossRef] [Green Version]
  109. O’Rourke, B.; Blatter, L.A. Mitochondrial Ca2+ uptake: Tortoise or hare? J. Mol. Cell. Cardiol. 2009, 46, 767–774. [Google Scholar] [CrossRef] [Green Version]
  110. Setterberg, I.E.; Le, C.; Frisk, M.; Perdreau-Dahl, H.; Li, J.; Louch, W.E. The Physiology and Pathophysiology of T-Tubules in the Heart. Front. Physiol. 2021, 12, 1193. [Google Scholar] [CrossRef] [PubMed]
  111. Gadeberg, H.C.; Bond, R.C.; Kong, C.H.T.; Chanoit, G.P.; Ascione, R.; Cannell, M.B.; James, A.F. Heterogeneity of T-Tubules in Pig Hearts. PLoS ONE 2016, 11, e0156862. [Google Scholar] [CrossRef] [Green Version]
  112. Ibrahim, M.; Gorelik, J.; Yacoub, M.H.; Terracciano, C.M. The structure and function of cardiac t-tubules in health and disease. Proceedings Biol. Sci. 2011, 278, 2714–2723. [Google Scholar] [CrossRef] [Green Version]
  113. Pérez-Treviño, P.; Pérez-Treviño, J.; Borja-Villa, C.; García, N.; Altamirano, J. Changes in T-Tubules and Sarcoplasmic Reticulum in Ventricular Myocytes in Early Cardiac Hypertrophy in a Pressure Overload Rat Model. Cell. Physiol. Biochem. 2015, 37, 1329–1344. [Google Scholar] [CrossRef] [PubMed]
  114. Schobesberger, S.; Wright, P.; Tokar, S.; Bhargava, A.; Mansfield, C.; Glukhov, A.V.; Poulet, C.; Buzuk, A.; Monszpart, A.; Sikkel, M.; et al. T-tubule remodelling disturbs localized β2-adrenergic signalling in rat ventricular myocytes during the progression of heart failure. Cardiovasc. Res. 2017, 113, 770–782. [Google Scholar] [CrossRef] [PubMed]
  115. Song, Z.; Liu, M.B.; Qu, Z. Transverse tubular network structures in the genesis of intracellular calcium alternans and triggered activity in cardiac cells. J. Mol. Cell. Cardiol. 2018, 114, 288–299. [Google Scholar] [CrossRef] [PubMed]
  116. Nivala, M.; Ko, C.Y.; Nivala, M.; Weiss, J.N.; Qu, Z. Criticality in Intracellular Calcium Signaling in Cardiac Myocytes. Biophys. J. 2012, 102, 2433–2442. [Google Scholar] [CrossRef] [Green Version]
  117. Colman, M.A.; Pinali, C.; Trafford, A.W.; Zhang, H.; Kitmitto, A. A computational model of spatio-temporal cardiac intracellular calcium handling with realistic structure and spatial flux distribution from sarcoplasmic reticulum and t-tubule reconstructions. PLoS Comput. Biol. 2017, 13, e1005714. [Google Scholar] [CrossRef] [Green Version]
  118. Hoang-Trong, T.M.; Lederer, W.J.; Jafri, M.S. Exploring SR Calcium and Cytosolic Calcium Wave Dynamics using a 3D Stochastic Myocyte Model. Biophys. J. 2014, 106, 320a. [Google Scholar] [CrossRef] [Green Version]
  119. Chen, X.; Feng, Y.; Huo, Y.; Tan, W. The Interplay of Rogue and Clustered Ryanodine Receptors Regulates Ca2+ Waves in Cardiac Myocytes. Front. Physiol. 2018, 9, 393. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  120. Cannell, M.B.; Kong, C.H.T.; Imtiaz, M.S.; Laver, D.R. Control of sarcoplasmic reticulum Ca2+ release by stochastic RyR gating within a 3D model of the cardiac dyad and importance of induction decay for CICR termination. Biophys. J. 2013, 104, 2149–2159. [Google Scholar] [CrossRef] [Green Version]
Figure 1. A schematic diagram of ventricular myocytes modeled as a rectangular solid.
Figure 1. A schematic diagram of ventricular myocytes modeled as a rectangular solid.
Membranes 11 00989 g001
Figure 2. The placement of calcium release sites (A) at one Z-depth, and (B) at one Z-disc. The inset in (A) shows the CRUs on two T-tubules at two adjacent Z-discs. The distribution of inter-CRU distance is derived based on the experimental data.
Figure 2. The placement of calcium release sites (A) at one Z-depth, and (B) at one Z-disc. The inset in (A) shows the CRUs on two T-tubules at two adjacent Z-discs. The distribution of inter-CRU distance is derived based on the experimental data.
Membranes 11 00989 g002
Figure 3. The schematic diagram shows a grid point as a cube of size 200 nm on each dimension in the 3-dimensional space. In this cube, we put a single dyad. Assuming that the width of the dyad is 300 nm, we model the efflux of calcium flowing into 4 adjacent grid points, and it is equally split into 4 parts. The grid location is given by the coordinate (i,j,k).
Figure 3. The schematic diagram shows a grid point as a cube of size 200 nm on each dimension in the 3-dimensional space. In this cube, we put a single dyad. Assuming that the width of the dyad is 300 nm, we model the efflux of calcium flowing into 4 adjacent grid points, and it is equally split into 4 parts. The grid location is given by the coordinate (i,j,k).
Membranes 11 00989 g003
Figure 4. Using the simulation data, a multitude of information can be extracted. Apart from a snapshot of a pseudo-line scan calcium transient during Ca2+ transient along the longitudinal direction, as shown in (A), we are able to show (B) the dynamics of free calcium content, (C) the dynamics of calcium-bound fluorescent, (D) the dynamics of total calcium content.
Figure 4. Using the simulation data, a multitude of information can be extracted. Apart from a snapshot of a pseudo-line scan calcium transient during Ca2+ transient along the longitudinal direction, as shown in (A), we are able to show (B) the dynamics of free calcium content, (C) the dynamics of calcium-bound fluorescent, (D) the dynamics of total calcium content.
Membranes 11 00989 g004aMembranes 11 00989 g004b
Figure 5. The probability of one CRU triggering the neighboring one at different distance and the delay. (A) The probability of triggering a Ca2+ spark under normal diastolic conditions ([Ca2+]myo = 0.096 µM, and [Ca2+]nsr = 1.02 mM). (B) The Ca2+ spark triggering delay under normal conditions. (C) The probability of triggering a Ca2+ spark under high cytosolic calcium ([Ca2+]myo = 0.4 µM, and [Ca2+]nsr = 1.02 mM). (D) The Ca2+ spark triggering delay under high cytosolic calcium.
Figure 5. The probability of one CRU triggering the neighboring one at different distance and the delay. (A) The probability of triggering a Ca2+ spark under normal diastolic conditions ([Ca2+]myo = 0.096 µM, and [Ca2+]nsr = 1.02 mM). (B) The Ca2+ spark triggering delay under normal conditions. (C) The probability of triggering a Ca2+ spark under high cytosolic calcium ([Ca2+]myo = 0.4 µM, and [Ca2+]nsr = 1.02 mM). (D) The Ca2+ spark triggering delay under high cytosolic calcium.
Membranes 11 00989 g005
Figure 6. The probability of one CRU triggering the neighboring one at different distances and the delay. (A) The probability of triggering a spark at the second CRU at different distances from an activated CRU under the high overload condition ([Ca2+] = 0.156 µM, and [Ca2+]nsr = 1.70 mM). (B) The delay in triggering the second CRU under that high overload condition. (C) The probability of triggering a spark at the second CRU at different distances from an activated CRU under the low overload condition ([Ca2+]myo = 0.156 µM, and [Ca2+]nsr = 1.30 mM). (D) The delay in triggering the second CRU under that low overload condition.
Figure 6. The probability of one CRU triggering the neighboring one at different distances and the delay. (A) The probability of triggering a spark at the second CRU at different distances from an activated CRU under the high overload condition ([Ca2+] = 0.156 µM, and [Ca2+]nsr = 1.70 mM). (B) The delay in triggering the second CRU under that high overload condition. (C) The probability of triggering a spark at the second CRU at different distances from an activated CRU under the low overload condition ([Ca2+]myo = 0.156 µM, and [Ca2+]nsr = 1.30 mM). (D) The delay in triggering the second CRU under that low overload condition.
Membranes 11 00989 g006aMembranes 11 00989 g006b
Figure 7. (A) A simulated calcium spark. (B) Free calcium shows the underlying structure of the release site (the delayed activation of the two satellite clusters are invisible under fluorescence profile. (C) The profile of a calcium spark giving FWHM = 1.85 um (each color represents the snapshot at different time points after the peak (e.g., bk = 0 means the black line at 0 ms is delayed). (D) The free calcium profile using back-calculation method agrees with experimental estimates, however, it underestimates the real free myoplasmic calcium amplitude.
Figure 7. (A) A simulated calcium spark. (B) Free calcium shows the underlying structure of the release site (the delayed activation of the two satellite clusters are invisible under fluorescence profile. (C) The profile of a calcium spark giving FWHM = 1.85 um (each color represents the snapshot at different time points after the peak (e.g., bk = 0 means the black line at 0 ms is delayed). (D) The free calcium profile using back-calculation method agrees with experimental estimates, however, it underestimates the real free myoplasmic calcium amplitude.
Membranes 11 00989 g007
Figure 8. The probability of one CRU triggering the neighboring CRU at (A) different distance and (B) the delay. The overload conditions ([Ca2+]myo = 0.156 µM, and [Ca2+]nsr = 1.30 mM) where each CRU has 3 satellite clusters of 10 RyR2s, each at distance 0.2 µm were used.
Figure 8. The probability of one CRU triggering the neighboring CRU at (A) different distance and (B) the delay. The overload conditions ([Ca2+]myo = 0.156 µM, and [Ca2+]nsr = 1.30 mM) where each CRU has 3 satellite clusters of 10 RyR2s, each at distance 0.2 µm were used.
Membranes 11 00989 g008
Figure 9. (A) Simulation geometry to study spark-induced waves (X = CRU location). Black X = the activated CRU, Blue X = the CRU to be activated by diffusing Ca2+. (B) The proposed simulation geometry with an intermediate cluster (Green X = intermediate RyR2 cluster).
Figure 9. (A) Simulation geometry to study spark-induced waves (X = CRU location). Black X = the activated CRU, Blue X = the CRU to be activated by diffusing Ca2+. (B) The proposed simulation geometry with an intermediate cluster (Green X = intermediate RyR2 cluster).
Membranes 11 00989 g009
Figure 10. (A) The Po,trigger of Ca2+ release from 9 activated CRUs on one Z-line on the CRU at different distances. (B) The time delay for the activation at different distances.
Figure 10. (A) The Po,trigger of Ca2+ release from 9 activated CRUs on one Z-line on the CRU at different distances. (B) The time delay for the activation at different distances.
Membranes 11 00989 g010
Figure 11. (A) The Po,trigger of Ca2+ release from 9 activated CRUs on one Z-line, with 1 intermediate RyR2 cluster in the middle, on the CRU at the next Z-line of different distance. (B) The time delay for the activation at different distances.
Figure 11. (A) The Po,trigger of Ca2+ release from 9 activated CRUs on one Z-line, with 1 intermediate RyR2 cluster in the middle, on the CRU at the next Z-line of different distance. (B) The time delay for the activation at different distances.
Membranes 11 00989 g011
Figure 12. (A) Calcium overload ([Ca]nsr = 1.7 mM, [Ca]i = 0.15 μM), this computational model of the rat ventricular myocyte can reproduce a repetitive sustained calcium wave which typically initiates at one end of the cell. The initiation site typically occurs where release sites are closer together or at a boundary. (B) Calcium waves in the rat ventricular myocyte under [Ca]o = 5 mM overload experimental conditions, the repetitive waves occur at a particular site for each cell. This suggests that there are more density of release sites surrounding the region that allows mass calcium release high enough to trigger the wave. Some waves can sustain to the next end, while some decay and stop in between, which suggests a stochastic nature of the waves (personal communication from Brian Hagen).
Figure 12. (A) Calcium overload ([Ca]nsr = 1.7 mM, [Ca]i = 0.15 μM), this computational model of the rat ventricular myocyte can reproduce a repetitive sustained calcium wave which typically initiates at one end of the cell. The initiation site typically occurs where release sites are closer together or at a boundary. (B) Calcium waves in the rat ventricular myocyte under [Ca]o = 5 mM overload experimental conditions, the repetitive waves occur at a particular site for each cell. This suggests that there are more density of release sites surrounding the region that allows mass calcium release high enough to trigger the wave. Some waves can sustain to the next end, while some decay and stop in between, which suggests a stochastic nature of the waves (personal communication from Brian Hagen).
Membranes 11 00989 g012
Figure 13. Given the initial [Ca]nsr = 1.7 mM, to derive the triggering of the wave, the simulation suggested that an overload of 1.5 mM is enough to trigger the repetitive calcium waves.
Figure 13. Given the initial [Ca]nsr = 1.7 mM, to derive the triggering of the wave, the simulation suggested that an overload of 1.5 mM is enough to trigger the repetitive calcium waves.
Membranes 11 00989 g013
Figure 14. Effect of SERCA on the probability of a Ca2+ triggering a Ca2+ wave. (A) uniform SERCA distribution; (B) 90% reduction in SERCA flux.
Figure 14. Effect of SERCA on the probability of a Ca2+ triggering a Ca2+ wave. (A) uniform SERCA distribution; (B) 90% reduction in SERCA flux.
Membranes 11 00989 g014
Table 1. Percentage of Ca2+ sparks resulting in a Ca2+ wave.
Table 1. Percentage of Ca2+ sparks resulting in a Ca2+ wave.
CRU Dist.0.60.81.0
RyR2 Sense
[Ca2+]sr,max = 1.13 mM93.22%1.27%0.8%
[Ca2+]sr,max = 1.3 mM95.48%1.55%1%
Po,trigger when [Ca2+]myo = 0.156 µM, and [Ca2+]sr = 1.3 mM. The latter case assumes lumenal Ca2+ sensitivity saturate at 1.13 mM.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Hoang-Trong, T.M.; Ullah, A.; Lederer, W.J.; Jafri, M.S. A Stochastic Spatiotemporal Model of Rat Ventricular Myocyte Calcium Dynamics Demonstrated Necessary Features for Calcium Wave Propagation. Membranes 2021, 11, 989. https://doi.org/10.3390/membranes11120989

AMA Style

Hoang-Trong TM, Ullah A, Lederer WJ, Jafri MS. A Stochastic Spatiotemporal Model of Rat Ventricular Myocyte Calcium Dynamics Demonstrated Necessary Features for Calcium Wave Propagation. Membranes. 2021; 11(12):989. https://doi.org/10.3390/membranes11120989

Chicago/Turabian Style

Hoang-Trong, Tuan Minh, Aman Ullah, William Jonathan Lederer, and Mohsin Saleet Jafri. 2021. "A Stochastic Spatiotemporal Model of Rat Ventricular Myocyte Calcium Dynamics Demonstrated Necessary Features for Calcium Wave Propagation" Membranes 11, no. 12: 989. https://doi.org/10.3390/membranes11120989

APA Style

Hoang-Trong, T. M., Ullah, A., Lederer, W. J., & Jafri, M. S. (2021). A Stochastic Spatiotemporal Model of Rat Ventricular Myocyte Calcium Dynamics Demonstrated Necessary Features for Calcium Wave Propagation. Membranes, 11(12), 989. https://doi.org/10.3390/membranes11120989

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop