Original Research ARTICLE
Computational Model of Calcium Signaling in Cardiac Atrial Cells at the Submicron Scale
- Departament de Física, Universitat Politècnica de Catalunya, Barcelona, Spain
In cardiac cells, calcium is the mediator of excitation-contraction coupling. Dysfunctions in calcium handling have been identified as the origin of some cardiac arrhythmias. In the particular case of atrial myocytes, recent available experimental data has found links between these dysfunctions and structural changes in the calcium handling machinery (ryanodine cluster size and distribution, t-tubular network, etc). To address this issue, we have developed a computational model of an atrial myocyte that takes into account the detailed intracellular structure. The homogenized macroscopic behavior is described with a two-concentration field model, using effective diffusion coefficients of calcium in the sarcoplasmic reticulum (SR) and in the cytoplasm. The model reproduces the right calcium transients and dependence with pacing frequency. Under basal conditions, the calcium rise is mostly restricted to the periphery of the cell, with a large concentration ratio between the periphery and the interior. We have then studied the dependence of the speed of the calcium wave on cytosolic and SR diffusion coefficients, finding an almost linear relation with the former, in agreement with a diffusive and fire mechanism of propagation, and little dependence on the latter. Finally, we have studied the effect of a change in RyR cluster microstructure. We find that, under resting conditions, the spark frequency decreases slightly with RyR cluster spatial dispersion, but markedly increases when the RyRs are distributed in clusters of larger size, stressing the importance of RyR cluster organization to understand atrial arrhythmias, as recent experimental results suggest (Macquaide et al., 2015).
Calcium is one of the most important intracellular messengers, and thus the mechanisms that control the intracellular free calcium concentration are of fundamental physiological importance (Berridge, 1997). For instance, Ca2+ takes part in oocyte activation at fertilization (Poenie et al., 1985), axonal growth (Bixby and Harris, 1991), cell migration (Huttenlocher et al., 1997), gene expression (Bading et al., 1993), formation of nodules in plant root hairs (Ehrhardt et al., 1996), development of muscle (Ford and Podolsky, 1972), release of cytokines from epithelial cells (Kaufman and Roizman, 1989), cell death (Schanne et al., 1979; Farber, 1981), and excitation-contraction coupling in muscle cells (Fabiato and Fabiato, 1979).
In cardiac cells, calcium dysregulation has been related to the appearance of arrhythmias and sudden cardiac death. A life-threatening arrhythmia, fibrillation, results when an electrical wavebreak induces reentry and triggers a cascade of new wavebreaks. Ventricular fibrillation (VF) is the most common cause of sudden death, whereas atrial fibrillation (AF), the most prevalent clinical arrhythmia, accounts for nearly one third of strokes in the elderly (Weiss et al., 2005). Clinically, AF duplicates the mortality rate and increases the risk of ictus (in which poor blood flow to the brain results in cell death) 5-fold. In spite of this the treatment of AF remains deficient or inefficient, because of the incomplete knowledge of the complex pathophysiology of this disease. Often, AF has been linked to a dysregulation in the dynamics of intracellular calcium, thus the importance of a good knowledge of calcium handling dynamics in the cell. On the other hand, in the last ten years, the refinement of the experimental techniques, such as STED and dSTORM (Hell and Wichmann, 1994; Izu et al., 2006; Soeller and Baddeley, 2013) has provided, for instance, a link between the calcium handling microstructure and the occurrence of cardiac diseases, as AF (Macquaide et al., 2015), prompting the quest for more detailed models of calcium handling, able to mechanistically explain this relation.
Inside cardiac cells, most intracellular calcium is stored in a complex structure called sarcoplasmic reticulum (SR), see Figure 1. Ca2+ is released from this internal network via the Ryanodine Receptors (RyR, Franzini-Armstrong and Protasi, 1997) when a threshold calcium concentration in the cytoplasm is achieved. This happens due to a small influx of calcium through the L-type calcium channels (LCC) during the cardiac action potential. This current triggers calcium release from the SR by activating the RyRs. RyRs open and close collectively in clusters forming functional units known as Calcium Release Units (CaRU), which are often confronted to a cluster of LCCs. In each CaRU the number of RyR and LCC is small (of the order of 10–100 of the former and 5–10 of the latter), thus its dynamics is intrinsically stochastic. CaRUs are distributed inside the cell, resulting in random and discrete Ca2+ release events, known as Ca2+ sparks (Cheng et al., 1994). A Ca2+ spark has been considered as the unitary dynamical element which produces the cellular Ca2+ dynamics, such as Ca2+ waves and oscillations (Falcke, 2003). The (seemingly deterministic) global calcium signal appears from the coordination of several tens of thousands of these CaRUs.
Figure 1. Basic components of the CICR process. Calcium enters through the LCCs, stimulating release from the RyRs, that then is reuptaken into the SR by SERCA and taken out of the cell by the sodium-calcium exchanger.
After the excitation process, Ca2+ removal allows relaxation of the cardiac muscle. This requires Ca2+ transport out of the cytoplasm by several pathways. The concentration in the SR is recovered by the active pumping of calcium from the cytoplasm to the SR carried out by the Sarcoplasmic Reticulum Ca2+-ATPase (SERCA). Moreover, the Na-Ca exchanger pumps Ca2+ out of the cell. The whole process described is called calcium-induced calcium release (CICR, Berridge, 1993; Clapham, 1995). CaRUs not just couple SR and cytoplasm Ca2+ concentrations via Ca2+ release but they are also correlated due to the Ca2+ diffusion in both domains. Therefore, the behavior of a single CaRU depends on the behaviors of the neighboring CaRUs.
Even though the same mechanism (CICR) triggers the transient elevation of Ca2+ in both ventricular and atrial myocytes, there are substantial differences in the intracellular structures. The absence of transversal tubules (t-tubules) in atrial myocytes produces inhomogeneous spatio-temporal calcium patterns when the CICR occurs. In particular, the excitation starts at the cell membrane and then propagates inward, resulting in a delay in activation time between the subsarcolemma and the cell interior. This is a key difference between atrial and ventricular cells. In the latter, the opening of LCC channels along the t-tubules triggers the release of calcium from the SR, resulting in a homogenized calcium pattern. In the former, this trigger is due to the inward wave.
Detailed models of calcium handling have been first developed for ventricular cells, including the stochastic modeling of each individual CaRU, coupled then by diffusion. In this framework, each CaRU is typically divided into different subcompartments, in which the calcium concentration is assumed to be homogeneous (Restrepo et al., 2008; Rovetti et al., 2010), although some recent models consider also calcium diffusion within the CaRU (Nivala M. et al., 2012). These models have been very successful in reproducing several calcium dysfunctions, such as calcium alternans (Restrepo et al., 2008; Rovetti et al., 2010; Alvarez-Lacalle et al., 2015) or spontaneous calcium release induced delayed afterdepolarizations (Song et al., 2015). Current advances in microscopy have allowed the development of very detailed models of calcium release at the level of the CaRU, including realistic shape of the SR, the RyR cluster, myofibrils and the mitochondria (Kekenes-Huskey et al., 2012; Hatano et al., 2013; Hake et al., 2014; Rajagopal et al., 2015).
Modeling is less developed for the case of atrial cells (Heijman et al., 2016). Common pool models, in which calcium concentration is considered to be homogeneous in each of several compartments (SR, cytosol, dyadic space, etc) have been developed for rabbit (Lindblad et al., 1996), dog (Ramirez et al., 2000), mouse (Davies et al., 2014), and human (Courtemanche et al., 1998; Nygren et al., 1998; Grandi et al., 2011; Lugo et al., 2014). One of the first models that took into account inward wave propagation was by Koivumäki et al. (2011), where the bulk cytoplasm and SR spaces were divided into several compartments, being thus a one-dimensional model, allowing for centripetal but not lateral diffusion. A similar model was also used by Li et al. (2012), showing the presence of alternans. A model allowing for both centripetal and lateral diffusion, as well as stochastic RyR gating was developed by Voigt et al. (2013), in order to study the mechanisms of after-depolarizations and triggered activity in paroxysmal atrial fibrillation. Calcium wave initiation and propagation has been considered by Thul et al. (2012) in a three-dimensional geometry, assuming a diffuse and fire model for calcium release. Finally, Macquaide et al. (2015) developed a detailed three-dimensional bidomain model of calcium propagation to study intra-CaRU cluster interactions, supporting the idea that cluster fragmentation and redistribution sustains atrial fibrillation through the enhancement of calcium release.
Still, there are several open questions regarding CICR in atrial cells. To name some: (1) the role of buffers, RyR sensitivity and the level of cytosolic calcium in calcium wave propagation; (2) the effect of the RyR cluster spatial structure and size distribution; (3) the role of t-tubules (if present). Most subcellular ventricular and atrial models (Restrepo et al., 2008; Rovetti et al., 2010; Voigt et al., 2013) consider the cell divided into several thousands of functional units (CaRUs). Each CaRU is then divided into different compartments, replicating at the subcellular scale the structure of common pool models. Despite the success of such models to replicate calcium transients and spark characteristics, they are not well-suited to study the effects of changes in the microstructure (position of the RyR clusters, inhomogeneities, etc). Rather, to study the effect of RyR cluster distribution on wave propagation, continuum models of calcium diffusion with point release sites have been considered, although often with simplified release dynamics (Izu et al., 2006; Thul et al., 2012, 2015; Øyehaug et al., 2013). On the other hand, very detailed models at the level of the CaRU (Hake et al., 2014) are very demanding computationally, and typically not well-suited to study effects that require of long simulation times, as calcium homeostasis or spark rates. With that in mind, we present a subcellular calcium atrial model where the homogenized local behavior is described with a two-concentration field model, using effective diffusion coefficients of calcium in the SR and in the cytoplasm, with stochastic gating of the RyRs and LCCs. This model follows the spirit of earlier bidomain models (Jafri and Keizer, 1995; Keener and Sneyd, 1998), defining at each point in space cytosolic and SR calcium concentrations, with given volume fractions (Keizer and De Young, 1992). The model presents some important characteristics: (1) a very fine discretization, making it possible to describe (even if coarsely) the RyR cluster structure; (2) incorporation of the cell structure with distinction between z-lines and normal cytosol in terms of the volume ratio of SR and cytosolic volumes, diffusion constants and presence of buffers; (3) freedom to set the center of the RyR clusters arbitrarily, that do not need to be disposed in an homogeneous regular grid. In this paper, we focus on the effect of CaRU spatial structure and distribution, and find that a more disordered distribution of the CaRUs presents a lower frequency of sparks in resting conditions. On the contrary, when the spatial distribution is maintained constant, but the RyRs are distributed in a smaller number of larger CaRUs (so the total number of RyRs remains constant), the spark frequency increases, in accordance with experimental results in cells presenting AF (Macquaide et al., 2015).
Our computational model performs single cell simulations and is based on homogenization (Goel et al., 2006). Although it is well-known that the SR forms a branching network (largely interconnected), with an interior that is distinct from the cell cytoplasm, this fact has largely been ignored, with most models making the a priori assumption that a Ca2+ concentration for both the SR and the cytoplasm can be defined at each point in space. So that, the cytoplasm and the SR are assumed to coexist at every point in space. For this reason, a fraction of each volume is occupied by the cytoplasm (vi) and the complementary fraction by the SR (vsr), given that vi + vsr = 1.
We define ci, csr, and cbi as the concentration of calcium in the cytoplasm, the SR, and the concentration of calcium bound to buffers. This description assumes that there exist effective diffusion coefficients Di = Di(vi) and Dsr = Dsr(vsr) that, in an average sense, incorporate the effect of that complex geometry. Although in principle these coefficients could be calculated knowing the SR structure (Goel et al., 2006), we will take the functional forms used in Goel et al. (2006). Since both fractions, vi and vsr, vary in different parts of the cell, it implies that both diffusion coefficients are functions of the position, Di = Di(r) and Dsr = Dsr(r). In our simulations we take the values Di ~ 250 μm2/s and Dsr ~ 90 μm2/s, that are within the upper range considered in the literature (Louch et al., 2010; Bers and Shannon, 2013).
The cardiac cell is modeled as a two dimensional domain with Lx = 100 μm and Ly = 15 μm. The spatial grid belongs to the submicron scale and it is defined as dx = dy = 0.1μm. There are points of the grid with and without RyRs. A typical RyR has a size of 30 x 30 nm. The RyRs are transmembrane proteins located at the surface of the SR, so they form a 2D grid. Thus in each of our grid points we locate a maximum of 10 RyRs.
A collection of grid points presenting RyRs form a cluster, i.e., a CaRU. In atrial cells, CaRUs are arranged periodically in the longitudinal and transversal directions, with some—seemingly Gaussian—dispersion (Chen-Izu et al., 2006). In our model, we place the centers of the clusters on the perimeter following an exact periodic distribution with a period μm (see Figure 2). In front of all these exterior CaRUs there are LCC groups. Inside the cell, CaRUs are placed following a Gaussian distribution centered at the z-lines and with a fixed dispersion σ. We take σ = 0.4 μm as standard value. The average distance between CaRUs is Tx = 1.6 μm and Ty = 0.5 μm. Experimental data shows that the SR domain coincides with these z-lines (Soeller et al., 2007). In this sense, we identify the z-lines with periodic narrow strips (0.3 μm width) with a predefined period (Tx). Let be Ωc the sarcomere domain, that is, the zone between z-lines and let be Ωsr the zone contained in z-lines and all the contour (∂Ω). Notice that Ωc∩Ωsr = ∅. Besides, we consider the presence of Ca2+ buffers: troponin (TnC), Calmodulin (CaM), and SR Ca-binding sites. The TnC buffer affects the cytoplasmic concentration of calcium in the Ωc domain. The other buffers, calmodulin and SR, affect also ci but in all the cell, Ωc∪Ωsr. We assume that all the buffers are immobile.
Figure 2. Bottom left corner of the whole cell. Circles (red) represent points with RyRs, black crosses are LCC groups, green stripes indicate the z-line region (domain Ωsr), and small dots the sarcomere (Ωc).
Because of the homogenization coarse grain, we define ci(r, t) (free calcium concentration), csr(r, t) (calcium concentration in the SR), and cbi(r, t) (calcium attached to buffers: TnC, CaM, and SR buffer) in all points. Therefore, we state the problem with the following set of partial differential equations (PDEs).
where Ji and Jsr are the fluxes into the cytosol and the SR spaces, respectively, Jbi accounts for the binding of free calcium to the different buffers. In order to relate the fluxes between the cytoplasm and the SR we have multiplied by the volume fraction vi/vsr, that depends on r, that is, on the domain Ωc and Ωsr. In addition, each point could have different components (RyR or not, LCC or not) and could belong to the membrane or not. The fluxes that may contribute to the total flux into the cytosol Ji are the SR release flux Jrel, the SERCA pump Jup, the L-type calcium flux JCaL and the sodium-calcium exchanger flux JNaCa. The release flux Jrel carries Ca2+ ions from the SR to the cytoplasm through the RyRs. Thus, it exists only on those points that have a CaRU, indicated by a red dot in Figure 2. Jup pumps calcium from the cytoplasm to the SR and it is present in all cell domain (Ωc∪Ωsr). The sum of these two fluxes (when appropriate) constitute the total flux from the SR. Then, JCaL, the inward L-type calcium flux, depends on the LCC clusters, so that it will act on those points that contain this channel, indicated by a cross in Figure 2. Indeed, LCCs appear only in some points of the cell membrane, ∂Ω (those that also have a CaRU). Finally, the NaCa exchanger, JNaCa, acts along all the perimeter ∂Ω.
A detailed description of all the fluxes can be found in the Supplementary Material. Below we present some details of the release and L-type calcium fluxes.
2.1. Release Flux
As shown in Figure 2, we consider each CaRU formed by several grid points containing RyRs. As standard for a CaRU, we consider one containing 36 RyRs, divided equally among 4 grid points, each one containing 9 RyRs. We will change this configuration in section 3.4 to consider larger CaRUs, maintaining fixed the total number of RyRs in the cell. This resembles the situation found in cells presenting AF (Macquaide et al., 2015).
Following Stern et al. (1999) each RyR can be in one of four different states: close C, open O, and two inactivated states I1, I2 (Figure 3). Calcium release from the SR to the cytoplasm is taken to be proportional to the concentration difference and the number of RyR in the open state, ORyR,
This flux is only present in those points that present RyRs (highlighted in red in Figure 2).
2.2. L-Type Calcium Flux
The inward current of calcium from the extracellular medium toward each CaRU is dependent on the number of LCC channels in the open state OLCC, the voltage, and the local calcium concentration in these points, which are close to the membrane, according to
where z = VF/(RT) and zm = 0.341zF. The current ICaL is converted to the flux JCaL, with units of μ M/ms, through:
where vmyo is the volume of the cytosol.
We have used the LCC model described in Mahajan et al. (2008) with some changes in the parameters as in Alvarez-Lacalle et al. (2015). We consider the presence of 5 LCC channels in each CaRU (located all in the same grid point) with five possible states (Figure 3): two closed states (C1 and C2), two inactivated states (I1 and I2) and one open state (O). The stochastic dynamics of the transitions is implemented using a time-adaptive Gillespie's method (Nivala J. et al., 2012). The transition rates aij are described in the Supplementary Material.
2.3. Other Fluxes
There are extra fluxes that appear on the model. The Na-Ca exchanger and the SERCA pump are both explained in the Supplementary Material.
3.1. Calcium Handling Characteristics
The calcium trace results from the sum of calcium at different sites. Since in our model the volume fraction changes from site to site, we have to define the average calcium concentrations as:
Figure 4 shows typical traces during one beat of the average calcium over all the cell in both domains: cytoplasm and SR. The calcium peak, of ~700–800 nM, agrees well with experimental observations (Mackenzie et al., 2004). The calcium concentration in the SR, though, is larger than observed in experiments due to the lack of the SR buffer calsequestrin (CSQN) in our model. We also show in Figure 4 the four cytoplasmic fluxes, corresponding to the sodium-calcium exchanger (JNaCa), the L-type calcium flux (JCaL), SR release (Jrel), and SERCA (Jup). Due to the small number of LCC channels, the L-type calcium flux is particularly stochastic.
Figure 4. (Top) Temporal profiles of [Ca2+] for both domains, cytoplasm and SR. (Bottom left) Inward and outward cell currents of JNaCa (dashed line) and JCaL (solid line). (Bottom right) SERCA pump flux (solid line) and release flux (dashed line) in the cytoplasm at a pacing period of 800 ms.
Depending on the pacing period, the model shows different behaviors. We have quantified this effect by calculating the calcium peak and the calcium diastolic level in the cytoplasm and in the SR domain (Figure 5). To assure that the system is close to the steady state, we have paced the cell for 50 s at each pacing period, and then taken the average over the next 20 stimulations. As the pacing period decreases, the cytosolic calcium peak increases moderately, up to a pacing period period of ~200-300 ms, beyond which it decreases, due to the decrease in SR calcium content and fractional release. This behavior agrees qualitatively with the observed change in the contractile force as a function of pacing period observed in atrial cells (Maier et al., 2000; Schotten et al., 2002), that shows a peak at a period of ~500 ms, beyond which it decreases.
Figure 5. Ca2+ peak and Ca2+ basal level as function of the pacing period (Ts) in both domains: cytoplasm and SR. Each point has been averaged over 10 beats in steady state.
3.2. Inward Calcium Wave Propagation
In order to compare the spatial heterogeneities within the cell, we have considered longitudinal sections at the central and peripheral regions, averaged over a 1 μm width. The complete CaRU distribution is shown in Figure 6, where the longitudinal sections are plotted in green and blue. Simulations suggest strong differences between calcium levels in the subsarcolemmal space and the center of the cell (see Figure 6), as well as a delay between release at the peripheral and central regions.
Figure 6. (Top) Spatial distribution of the CaRU. Black dots show the position of the CaRUs. The green and blue lines indicate longitudinal lines at the cell center and periphery. (Down) Local Ca2+ measured in the subsarcolemmal space (blue dashed lines) and the center of the cell (solid lines) for both cytoplasm and SR domains at a pacing period of 800 ms. All traces have been averaged over a longitudinal section of 1 μm width.
The spatio-temporal and local correlation between ci and csr calcium is shown in the line scan profiles on Figure 7. The four profiles correspond to the same beat. In the subsarcolemmal region the presence of LCCs and CaRUs results in an important release activity causing a relevant SR depletion. On the other hand, in the central region, calcium does not penetrate, and the local activity is scarce. Still, there is a depletion of the SR content (visible also in Figure 6) due, not so much to release, almost negligible at the central region, but to diffusion of SR calcium to the periphery.
Figure 7. Longitudinal line scan during a single beat on the subsarcolemmal and the central region at a pacing period of 800 ms for cytosolic (Top) and SR (Bottom) calcium. The colorbar corresponds to calcium concentration in μM.
The spatio-temporal Ca2+ dynamics in the cytoplasm allows us to clearly understand how the standard inward wave propagation occurs. Figure 8 shows spatial profiles at different times during a single beat. Under normal conditions, the calcium wave starts on the cell membrane and propagates to the center but this propagation does not reach the central region. This situation is observed more clearly averaging the calcium concentration over the longitudinal direction, so we can observe the average inward propagation of the calcium wave (Figure 9) Typically, the inward wave propagates 4 or 5 μm in the transversal direction. From the figure, we can estimate an inward wave velocity of roughly 150 μm/s, that agrees well with typical observed calcium wave velocities of ~100 μm/s (Izu et al., 2013).
Figure 8. Inward wave propagation at a pacing period of 800 ms. The colorbar corresponds to calcium concentration in μM.
Figure 9. Line scans averaged over the longitudinal direction at pacing periods of 500 and 800ms. The colorbar corresponds to calcium concentration in μM.
Intracellular waves are Ca2+ release events that propagate across the cell at a constant velocity. To have a better control of the calcium wave and be able to study its speed and dependence on different parameters, we have created a new geometry with 10 equidistant z-lines (see Figure S2, in Supplementary Material). The distance between two z-lines is 1.5 μm. Initially, cytosolic calcium at the first z-line is increased and then the system let to evolve without being forced. The wave front is monitored and the wave front velocity calculated. This way we determine the wave velocity as a function of different parameters. The typical wave velocity is of the order of 200−300 μm/s, that agrees well with a diffusive process within z-lines, that would give a speed of v ~ 2D/d ~ 2·200 μm2s−1/1.5 μm ~ 260 μm/s. This velocity increases slightly with the calcium SR load (Figure 10). The dependence on intracellular calcium diffusion Di is roughly linear, as one would expect in saltatory dynamics (Dawson et al., 1999). The dependence on SR calcium diffusion Dsr, on the other hand, is not so pronounced.
Figure 10. Propagation velocity as function of the diffusion constants Di and Dsr for different values of the SR calcium load. As standard values we take Di = 0.25 μm2/ms and Dsr = 0.09 μm2/ms in the zone between z-lines.
3.3. Effect of the Cell Structure
It is interesting to compare also the calcium dynamics at the z-lines and in the space within z-lines, where no CaRUs are present. To this end we have performed transversal section measurements, as shown in Figure 11, in a situation when the cell is at rest, without external stimulation. The sarcomere measurement corresponds to the space between the z-lines. It is important to notice that, because of the proximity between the measurements, there exists a correlation between the resulting profiles. For instance, at t = 0.6 s there is a spark that starts in the first z-line, it propagates to the sarcomere region and, then, to the second z-line. The second thing to notice is that, due to the presence of random Ca2+ releases associated to the position of the RyRs, at the z-lines the calcium trace is more stochastic (Figure 11).
Figure 11. (Top) Spatial distribution of the CaRUs. Black dots show the position of the CaRUs. The blue and red lines indicate two transversal lines within neighboring z-lines. The green line a transversal line within those two z-lines. (Down) Local Ca2+ measured in the three transversal lines in post-rest potential conditions. All traces have been averaged over a transversal section of 0.3 μm width.
3.4. RyR Cluster Structure and Distribution
Calcium sparks are the basic calcium release events. A good understanding of their characteristics (size, amplitude, and frequency) is thus very important to properly characterize the process of CICR. Due to the fine discretization of our model, we can observe their detailed spatio-temporal profiles. Figure 12 shows, for instance, the standard time evolution of a spark. We have also studied the effect of the microstructure in the frequency of sparks. We have modified the microstructure, changing the size and distribution of the CaRUs. This is particularly important since it has been observed that the RyR distribution changes in a particular way under conditions of AF (Macquaide et al., 2015). We have then calculated the spark frequency under resting conditions for different configurations defined by the Gaussian distribution of position sites and size of the CaRUs. For the standard size of the CaRU (36 RyRs divided equally among 4 grid points, each one containing 9 RyRs), we have considered three values of the dispersion in the Gaussian distribution, the standard value of σ = 0.4 μm, and two cases with larger dispersion of σ = 2 and 3.6 μm, see Figure 13. Besides, we also consider the effect of a change in the CaRU size, considering CaRUs with 54 RyRs, divided equally among 6 grid points, each one containing 9 RyRs (Figure 13 down right), but maintaining fixed the total number of RyRs in the cell. Since the total number of RyRs is the same, this means that there are larger, but less CaRUs in the cell.
Figure 13. Different RyR distributions considered in the text. The bottom left corner of the whole cell is displayed. The red dots represent grid points with RyRs. First: control configuration, with standard dispersion in the position of the RyRs σ = 0.4 μm. Second: increased the dispersion of σ to 2 μm. Third: σ = 3.6 μm. Fourth: New structure configuration, where each grid point contains 9 RyRs and a CaRU is formed by 6 grid points, so that, each CaRU represents 54 RyRs. The total number of RyRs remains constant but now, in the new structure, they are more grouped, that is, the CaRUs are bigger. The dispersion is the same as in the standard case: σ = 0.4 μm.
In Figure 14A the average mean CaRU size is shown. To define the size of a CaRU, we follow the results by Macquaide et al. (2015), that showed using a computational model that clusters closer than 150 nm triggered together functionally as a single cluster. We thus assume that a group of clusters belong to the same CaRU if they are separated, at most, by 0.15 μm edge to edge from each other. By definition, there is a random component in the structure, so that, there is a non-zero dispersion on the distance with respect to the z-line (Figure 14C). Finally, the mean nearest neighbor distance is shown in Figure 14B. In the new structure, since the RyRs are more grouped, the total number of CaRUs is smaller, such that the mean distance among them increases.
Figure 14. (A) Mean area occupied by a CaRU. (B) Nearest-neighbor distance among CaRUs, measured from the center of each CaRU. (C) Longitudinal dispersion, measured as the mean distance between the RyRs and the center of the z-line. All calculations have been averaged over 35 configurations.
The spark frequency is one of the most important indicators of the stochastic activity during a post-rest potential period. The total number of spontaneous calcium sparks has been recorded. In order to consider the sparks, we have counted all the calcium release events greater than a certain spatial radial threshold of 1.6 μm. As shown in Figure 15, we observe that the frequency decreases with the longitudinal dispersion meaning that the cluster-cluster communication plays an important role in stochastic activity. On the other hand, when the CaRUs are bigger (structure configuration) this probability increases.
Figure 15. Spark frequency for the four different configurations. The spark frequency decreases with nearest-neighbor distance but increases with the size of the CaRU.
4. Discussion and Conclusions
In this work we present a novel model to fully simulate a 2D longitudinal plane of a cardiac atrial cell. By modeling the intracellular calcium dynamics and solving the spatio-temporal Ca2+ reaction-diffusion equations, both local and global behavior have been recorded. Because of the high spatial resolution, the model allows us to study in detail the dynamics in the surroundings of the CaRUs, that is, the local calcium concentrations and the spark activity. It is also well-suited to study the effect of changes in the spatial distribution and form of CaRUs. In this regard, important differences have been noticed using different spatial configurations of RyRs, showing that the resulting calcium dynamics is highly dependent on the spatial distribution. We observe, for instance, a decrease in spark frequency with CaRU spatial dispersion (Figure 15). This is correlated with an increase in CaRU nearest-neighbor distance, suggesting that cooperativity among local release events at nearby CaRUs could play an important role in the generation of sparks. In fact, sparks (or macrosparks) encompassing several CaRUs have also been observed experimentally (Kockskämper et al., 2001). When the same total number of RyRs of the cell are distributed in larger CaRUs, we observe an increase in spark frequency (Figure 15). These larger CaRUs are obtained increasing their size, but maintaining the density of RyRs, similar to what is observed in atrial cells presenting AF (Macquaide et al., 2015). Due to the larger number of RyRs per CaRU, an increase in the spark probability of each individual CaRU is expected, but, since the number of CaRUs is decreased, the increase in spark frequency per CaRU must be non-linear with size. In this case, the increase in nearest-neighbor distance (Figure 14) does not result in a decrease in spark rate. However, it should be noted that we have measured the nearest-neighbor distance from the centers of the CaRUs, but since their size is larger, the distance between the edges of the CaRUs could actually be similar, or smaller. A more detailed study of the influence of the CaRU structure in spark frequency, the appearance of macrosparks, and the transition to waves is deferred to a future study.
In our simulations, the opening of the L-type calcium channels induce a calcium increase in the periphery of the cell, that hardly reaches the interior (Figure 6). This result agrees with what is observed in atrial cells without t-tubules where, under basal conditions, calcium signals are restricted to subsarcolemmal regions (Mackenzie et al., 2004). The observed values of the calcium transients, slightly over 1 μM at the periphery and 0.2–0.3 μM at the center are similar to what we obtain in our simulations (Figure 6). As observed in previous works (Dawson et al., 1999), the speed of the wave varies roughly linearly with the diffusion coefficient in the cytosol. Although we have not changed the distance among z-planes in this work, the fact that this linear relation continues at low values of the diffusion constant, seems to indicate the lack of a threshold for propagation in the model. This contrasts with results by Izu et al. (2006) and Hoang-Trong et al. (2015) that found a strong dependence on wave initiation with the distance among CaRUs. However, one should notice that, for the calculation of wave propagation (Figure 10), we have considered an increased strength of the L-type calcium current, probably pushing the threshold to smaller values of the diffusion constant than we have considered. In our simulations, calcium waves are obtained at transients larger than in experiments and with diffusion and RyR cluster spacing in the upper and lower ranges, respectively, allowing for diffuse and fire behavior. However, one should mention the difficulty to reconciliate calcium wave propagation at low calcium concentrations with a stochastic description of the RyR cluster (Izu et al., 2013). Recent studies suggest also an important role of SR calcium diffusion for the propagation of the calcium wave, through junctional SR calcium depletion and sensitization of the RyRs (Keller et al., 2007). We find that the dependence with diffusion in the SR is not so pronounced (Figure 10), seemingly ruling out an important role of SR Ca diffusion in wave propagation. This effect, however, should be studied in more detail, as well as the role of buffers, RyR sensitivity and the level of cytosolic calcium in calcium wave propagation.
The present study presents several limitations. To cite some, that we consider a two-dimensional geometry, a voltage clamp protocol, isotropic diffusion, or immobile buffers. The main reason to use a two-dimensional geometry was computational cost. A generalization to three-dimensions is straight-forward and we are implementing it to study some of the questions posed in this article in more detail. The dynamics of transmembrane voltage can be also readily incorporated into the model, and could be used to study the arrhythmogenic effect of spontaneous calcium release events, for instance. The correct characterization of calcium diffusion in the cell represents a harder challenge. We have considered typical values of the diffusion coefficients in the cytosol and SR and assumed that they depend linearly on the cytosolic/SR volume fractions, as suggested by homogenization (Goel et al., 2006). However, diffusion (particularly in the SR) is most likely to be anisotropic, and this could importantly affect wave characteristics. A better knowledge of the SR microstructure could help to estimate these diffusion coefficients and give a better representation of calcium wave propagation. An important addition would be to incorporate a (partial) network of t-tubules. The presence of t-tubules in atrial cells has been found to depend on the species, and there is evidence of their presence in large mammals' atrial cells (Richards et al., 2011). Besides transversal, axial tubules have also been found to contribute to rapid activation of the atrial cell (Brandenburg et al., 2016). In heart disease, including human heart failure (HF), there is extensive remodeling, resulting in loss and disorganization of t-tubules (Dibb et al., 2013). Besides, there are other important factors that may affect calcium transients. For instance, the presence of IP3R may affect the form of Ca2+ sparks, leading to a difference between calcium handling at the peripheral and central regions (Mackenzie et al., 2004; Kim et al., 2010). Another important effect, not included in our model, is the presence of mitochondria. In venticular myocytes, there is evidence suggesting that the mitochondrial outer membrane is linked to t-tubules (Hayashi et al., 2009). Models of excitation-contraction coupling, including mitochondrial calcium handling have been developed for ventricular myocytes (Cortassa et al., 2003, 2006; Matsuoka et al., 2004; Maack and O'rourke, 2007; Hatano et al., 2011) and used, to study, for instance, the influence of the distance between mitochondria and Ca2+ release sites (Hatano et al., 2013). In the atria, the mitochondria has been suggested to act as a buffer that prevents inward calcium propagation (Mackenzie et al., 2004). The effect of these structural factors on wave propagation is an important matter for future work.
MM assisted with the design of the work, performed the simulations, analyzed the data and prepared the manuscript. BE obtained funds, designed the work and revised and edited the manuscript. All authors approved the final manuscript.
Financial support was provided by the Fundació La Marató de TV3, under grant number 20151110, and from the Spanish Ministerio de Economía y Competitividad (MINECO) under grant number SAF2014-58286-C2-2-R.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank L. Hove-Madsen and E. Alvarez-Lacalle for fruitful discussions.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01760/full#supplementary-material
Alvarez-Lacalle, E., Echebarria, B., Spalding, J., and Shiferaw, Y. (2015). Calcium alternans is due to an order-disorder phase transition in cardiac cells. Phys. Rev. Lett. 114:108101. doi: 10.1103/PhysRevLett.114.108101
Brandenburg, S., Kohl, T., Williams, G. S., Gusev, K., Wagner, E., Rog-Zielinska, E. A., et al. (2016). Axial tubule junctions control rapid calcium signaling in atria. J. Clin. Invest. 126:3999. doi: 10.1172/JCI88241
Chen-Izu, Y., McCulle, S. L., Ward, C. W., Soeller, C., Allen, B. M., Rabang, C., et al. (2006). Three-dimensional distribution of ryanodine receptor clusters in cardiac myocytes. Biophys. J. 91, 1–13. doi: 10.1529/biophysj.105.077180
Cortassa, S., Aon, M. A., Marbán, E., Winslow, R. L., and O'Rourke, B. (2003). An integrated model of cardiac mitochondrial energy metabolism and calcium dynamics. Biophys. J. 84, 2734–2755. doi: 10.1016/S0006-3495(03)75079-6
Cortassa, S., Aon, M. A., O'Rourke, B., Jacques, R., Tseng, H. J., Marbán, E., et al. (2006). A computational model integrating electrophysiology, contraction, and mitochondrial bioenergetics in the ventricular myocyte. Biophys. J. 91, 1564–1589. doi: 10.1529/biophysj.105.076174
Courtemanche, M., Ramirez, R. J., and Nattel, S. (1998). Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am. J. Physiol. Heart Circ. Physiol. 275, H301–H321. doi: 10.1152/ajpheart.1998.275.1.H301
Davies, L., Jin, J., Shen, W., Tsui, H., Shi, Y., Wang, Y., et al. (2014). Mkk4 is a negative regulator of the transforming growth factor beta 1 signaling associated with atrial remodeling and arrhythmogenesis with age. J. Am. Heart Assoc. 3:e000340. doi: 10.1161/JAHA.113.000340
Dibb, K. M., Clarke, J. D., Eisner, D. A., Richards, M. A., and Trafford, A. W. (2013). A functional role for transverse (t-) tubules in the atria. J. Mol. Cell. Cardiol. 58, 84–91. doi: 10.1016/j.yjmcc.2012.11.001
Franzini-Armstrong, C., and Protasi, F. (1997). Ryanodine receptors of striated muscles: a complex channel capable of multiple interactions. Physiol. Rev. 77, 699–729. doi: 10.1152/physrev.19184.108.40.2069
Grandi, E., Pandit, S. V., Voigt, N., Workman, A. J., Dobrev, D., Jalife, J., et al. (2011). Human atrial action potential and ca2+ model: sinus rhythm and chronic atrial fibrillation. Circ. Res. 109, 1055–1066. doi: 10.1161/CIRCRESAHA.111.253955
Hatano, A., Okada, J., Washio, T., Hisada, T., and Sugiura, S. (2011). A three-dimensional simulation model of cardiomyocyte integrating excitation-contraction coupling and metabolism. Biophys. J. 101, 2601–2610. doi: 10.1016/j.bpj.2011.10.020
Hatano, A., Okada, J.-I., Washio, T., Hisada, T., and Sugiura, S. (2013). Mitochondrial colocalization with ca2+ release sites is crucial to cardiac metabolism. Biophys. J. 104, 496–504. doi: 10.1016/j.bpj.2012.12.004
Hayashi, T., Martone, M. E., Yu, Z., Thor, A., Doi, M., Holst, M. J., et al. (2009). Three-dimensional electron microscopy reveals new details of membrane systems for ca2+ signaling in the heart. J. Cell Sci. 122, 1005–1013. doi: 10.1242/jcs.028175
Heijman, J., Erfanian Abdoust, P., Voigt, N., Nattel, S., and Dobrev, D. (2016). Computational models of atrial cellular electrophysiology and calcium handling, and their role in atrial fibrillation. J. Physiol. 594, 537–553. doi: 10.1113/JP271404
Hell, S. W., and Wichmann, J. (1994). Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy. Opt. Lett. 19, 780–782. doi: 10.1364/OL.19.000780
Huttenlocher, A., Palecek, S. P., Lu, Q., Zhang, W., Mellgren, R. L., Lauffenburger, D. A., et al. (1997). Regulation of cell migration by the calcium-dependent protease calpain. J. Biol. Chem. 272, 32719–32722. doi: 10.1074/jbc.272.52.32719
Izu, L. T., Means, S. A., Shadid, J. N., Chen-Izu, Y., and Balke, C. W. (2006). Interplay of ryanodine receptor distribution and calcium dynamics. Biophys. J. 91, 95–112. doi: 10.1529/biophysj.105.077214
Jafri, M. S., and Keizer, J. (1995). On the roles of ca2+ diffusion, ca2+ buffers, and the endoplasmic reticulum in ip3-induced ca2+ waves. Biophys. J. 69, 2139–2153. doi: 10.1016/S0006-3495(95)80088-3
Kekenes-Huskey, P. M., Cheng, Y., Hake, J. E., Sachse, F. B., Bridge, J. H., Holst, M. J., et al. (2012). Modeling effects of l-type ca2+ current and na+-ca2+ exchanger on ca2+ trigger flux in rabbit myocytes with realistic t-tubule geometries. Front. Physiol. 3:351. doi: 10.3389/fphys.2012.00351
Kim, J.-C., Son, M.-J., Subedi, K. P., Li, Y., Ahn, J. R., and Woo, S.-H. (2010). Atrial local ca2+ signaling and inositol 1, 4, 5-trisphosphate receptors. Prog. Biophys. Mol. Biol. 103, 59–70. doi: 10.1016/j.pbiomolbio.2010.02.002
Kockskämper, J., Sheehan, K. A., Bare, D. J., Lipsius, S. L., Mignery, G. A., and Blatter, L. A. (2001). Activation and propagation of ca 2+ release during excitation-contraction coupling in atrial myocytes. Biophys. J. 81, 2590–2605. doi: 10.1016/S0006-3495(01)75903-6
Koivumäki, J. T., Korhonen, T., and Tavi, P. (2011). Impact of sarcoplasmic reticulum calcium release on calcium dynamics and action potential morphology in human atrial myocytes: a computational study. PLoS Comput. Biol. 7:e1001067. doi: 10.1371/journal.pcbi.1001067
Li, Q., O'Neill, S. C., Tao, T., Li, Y., Eisner, D., and Zhang, H. (2012). Mechanisms by which cytoplasmic calcium wave propagation and alternans are generated in cardiac atrial myocytes lacking t-tubules-insights from a simulation study. Biophys. J. 102, 1471–1482. doi: 10.1016/j.bpj.2012.03.007
Lindblad, D. S., Murphey, C. R., Clark, J. W., and Giles, W. R. (1996). A model of the action potential and underlying membrane currents in a rabbit atrial cell. Am. J. Physiol. Heart Circ. Physiol. 271, H1666–H1696. doi: 10.1152/ajpheart.1996.271.4.H1666
Louch, W. E., Hake, J., Jølle, G. F., Mørk, H. K., Sjaastad, I., Lines, G. T., et al. (2010). Control of ca2+ release by action potential configuration in normal and failing murine cardiomyocytes. Biophys. J. 99, 1377–1386. doi: 10.1016/j.bpj.2010.06.055
Lugo, C. A., Cantalapiedra, I. R., Peñaranda, A., Hove-Madsen, L., and Echebarria, B. (2014). Are sr ca content fluctuations or sr refractoriness the key to atrial cardiac alternans?: insights from a human atrial model. Am. J. Physiol. Heart Circ. Physiol. 306, H1540–H1552. doi: 10.1152/ajpheart.00515.2013
Mackenzie, L., Roderick, H. L., Berridge, M. J., Conway, S. J., and Bootman, M. D. (2004). The spatial pattern of atrial cardiomyocyte calcium signalling modulates contraction. J. Cell Sci. 117, 6327–6337. doi: 10.1242/jcs.01559
Macquaide, N., Tuan, H.-T., Hotta, J., Sempels, W., Lenaerts, I., Holemans, P., et al. (2015). Ryanodine receptor cluster fragmentation and redistribution in persistent atrial fibrillation enhance calcium release. Cardiovasc. Res. 108, 387–398. doi: 10.1093/cvr/cvv231
Mahajan, A., Shiferaw, Y., Sato, D., Baher, A., Olcese, R., Xie, L.-H., et al. (2008). A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophys. J. 94, 392–410. doi: 10.1529/biophysj.106.98160
Maier, L. S., Barckhausen, P., Weisser, J., Aleksic, I., Baryalei, M., and Pieske, B. (2000). Ca2+ handling in isolated human atrial myocardium. Am. J. Physiol. Heart Circ. Physiol. 279, H952–H958. doi: 10.1152/ajpheart.2000.279.3.H952
Matsuoka, S., Sarai, N., Jo, H., and Noma, A. (2004). Simulation of atp metabolism in cardiac excitation–contraction coupling. Prog. Biophys. Mol. Biol. 85, 279–299. doi: 10.1016/j.pbiomolbio.2004.01.006
Nivala, J., Knowles, P., Dotro, G., García, J., and Wallace, S. (2012). Clogging in subsurface-flow treatment wetlands: measurement, modeling and management. Water Res. 46, 1625–1640. doi: 10.1016/j.watres.2011.12.051
Nivala, M., de Lange, E., Rovetti, R., and Qu, Z. (2012). Computational modeling and numerical methods for spatiotemporal calcium cycling in ventricular myocytes. Front. Physiol. 3:114. doi: 10.3389/fphys.2012.00114
Nygren, A., Fiset, C., Firek, L., Clark, J. W., Lindblad, D. S., Clark, R. B., et al. (1998). Mathematical model of an adult human atrial cell: the role of k+ currents in repolarization. Circ. Res. 82, 63–81.
Øyehaug, L., Loose, K. Ø., Jølle, G. F., Røe, Å. T., Sjaastad, I., Christensen, G., et al. (2013). Synchrony of cardiomyocyte ca2+ release is controlled by t-tubule organization, sr ca2+ content, and ryanodine receptor ca2+ sensitivity. Biophys. J. 104, 1685–1697. doi: 10.1016/j.bpj.2013.03.022
Rajagopal, V., Bass, G., Walker, C. G., Crossman, D. J., Petzer, A., Hickey, A., et al. (2015). Examination of the effects of heterogeneous organization of ryr clusters, myofibrils and mitochondria on ca2+ release patterns in cardiomyocytes. PLoS Comput. Biol. 11:e1004417. doi: 10.1371/journal.pcbi.1004417
Ramirez, R. J., Nattel, S., and Courtemanche, M. (2000). Mathematical analysis of canine atrial action potentials: rate, regional factors, and electrical remodeling. Am. J. Physiol. Heart Circ. Physiol. 279, H1767–H1785. doi: 10.1152/ajpheart.2000.279.4.H1767
Richards, M. A., Clarke, J. D., Saravanan, P., Voigt, N., Dobrev, D., Eisner, D. A., et al. (2011). Transverse tubules are a common feature in large mammalian atrial myocytes including human. Am. J. Physiol. Heart Circ. Physiol. 301, H1996–H2005. doi: 10.1152/ajpheart.00284.2011
Rovetti, R., Cui, X., Garfinkel, A., Weiss, J. N., and Qu, Z. (2010). Spark-induced sparks as a mechanism of intracellular calcium alternans in cardiac myocytes. Circ. Res. 106, 1582–1591. doi: 10.1161/CIRCRESAHA.109.213975
Schotten, U., Greiser, M., Benke, D., Buerkel, K., Ehrenteidt, B., Stellbrink, C., et al. (2002). Atrial fibrillation-induced atrial contractile dysfunction: a tachycardiomyopathy of a different sort. Cardiovasc. Res. 53, 192–201. doi: 10.1016/S0008-6363(01)00453-9
Soeller, C., Crossman, D., Gilbert, R., and Cannell, M. B. (2007). Analysis of ryanodine receptor clusters in rat and human cardiac myocytes. Proc. Natl. Acad. Sci. U.S.A. 104, 14958–14963. doi: 10.1073/pnas.0703016104
Song, Z., Ko, C. Y., Nivala, M., Weiss, J. N., and Qu, Z. (2015). Calcium-voltage coupling in the genesis of early and delayed afterdepolarizations in cardiac myocytes. Biophys. J. 108, 1908–1921. doi: 10.1016/j.bpj.2015.03.011
Stern, M. D., Song, L.-S., Cheng, H., Sham, J. S., Yang, H. T., Boheler, K. R., et al. (1999). Local control models of cardiac excitation–contraction coupling: a possible role for allosteric interactions between ryanodine receptors. J. Gen. Physiol. 113, 469–489. doi: 10.1085/jgp.113.3.469
Thul, R., Coombes, S., Roderick, H. L., and Bootman, M. D. (2012). Subcellular calcium dynamics in a whole-cell model of an atrial myocyte. Proc. Natl. Acad. Sci. U.S.A. 109, 2150–2155. doi: 10.1073/pnas.1115855109
Thul, R., Rietdorf, K., Bootman, M. D., and Coombes, S. (2015). Unifying principles of calcium wave propagation-insights from a three-dimensional model for atrial myocytes. Biochim. Biophys. Acta, 1853, 2131–2143. doi: 10.1016/j.bbamcr.2015.02.019
Voigt, N., Heijman, J., Wang, Q., Chiang, D. Y., Li, N., Karck, M., et al. (2013). Cellular and molecular mechanisms of atrial arrhythmogenesis in patients with paroxysmal atrial fibrillation. Circulation 129, 145–156. doi: 10.1161/CIRCULATIONAHA.113.006641
Keywords: calcium modeling, atrial cells, local calcium signaling, calcium release unit, ryanodine receptor
Citation: Marchena M and Echebarria B (2018) Computational Model of Calcium Signaling in Cardiac Atrial Cells at the Submicron Scale. Front. Physiol. 9:1760. doi: 10.3389/fphys.2018.01760
Received: 05 December 2018; Accepted: 21 November 2018;
Published: 10 December 2018.
Edited by:Simonetta Filippi, Università Campus Bio-Medico, Italy
Reviewed by:Sonia Cortassa, National Institutes of Health (NIH), United States
Sanjay Ram Kharche, University of Western Ontario, Canada
Copyright © 2018 Marchena and Echebarria. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Blas Echebarria, email@example.com