## Abstract

The objective of the study is to investigate the propagation of Ca^{2+} waves in full-width cardiac myocytes and carry out sensitivity analysis to study the effects of various physiological parameters on global Ca^{2+} waves. Based on the anomalous subdiffusion of Ca^{2+} sparks, a mathematical model was proposed to characterize the Ca^{2+} waves. The computed results were in agreement with the experimental measurements using confocal microscopy. This model includes variables of current through the Ca^{2+} release unit (CRU; *I*_{CRU}), duration of current flow through CRU (*T*_{open}), Ca^{2+} sensitivity parameter (*K*), the longitudinal and transverse spatial separation of CRUs (*l _{x}* and

*l*, where

_{y}*x*denotes longitudinal direction (

*x*-axis) and

*y*denotes transverse direction (

*y*-axis)) and Ca

^{2+}diffusion coefficients (

*D*,

_{x}*D*). The spatio-temporal mechanism of the anomalous Ca

_{y}^{2+}sparks led to results that were different from those based on Fick's law. The major findings were reported as:

*I*

_{CRU}affected the dynamic properties of Ca

^{2+}waves more significantly than

*T*

_{open}; the effect of

*K*on the properties of Ca

^{2+}waves was negligible;

*l*affected the amplitude significantly, but

_{y}*l*affected the longitudinal velocity significantly; in turn, the limitation and significance of the study are discussed.

_{x}## 1. Introduction

The concentration of Ca^{2+} in the endoplasmic and sarcoplasmic reticulum (SR) of cardiac myocytes is two to three orders of magnitude higher than that in cytosol. The process involving Ca^{2+} release from the SR, cytosolic Ca^{2+} diffusion, reaction of Ca^{2+} with various Ca^{2+}-binding proteins and resequestration of Ca^{2+} into the SR is called ‘Ca^{2+} spark’ [1–3], which results in a Ca^{2+} transient [1,4,5], localized increase of cytoplasmic Ca^{2+} concentration. Although the calcium homeostasis is important for the contraction and relaxation of cardiac muscle, spontaneous Ca^{2+} sparks can trigger a propagating wave of high Ca^{2+} concentration in some pathological conditions. This ‘fire–diffuse–fire’ event is called ‘calcium wave’ [6]. Fabiato & Fabiato [7] observed the event of Ca^{2+} waves in cardiac muscle. The Ca^{2+} waves were also found in other cell types such as skeletal [8], medaka eggs [9] and astrocytes [10]. Because Ca^{2+} waves may induce some heart diseases, such as ventricular arrhythmias, ventricular fibrillation and heart failure [11,12], it is of particular importance to understand the mechanism of Ca^{2+} waves.

The Fick's law models for Ca^{2+} sparks showed the full-width at half maximum (FWHM) of approximately 1.0 μm in comparison with the measured FWHM of approximately 2.0 μm [13–17]. This was called ‘FWHM paradox’. Izu *et al*. [15] tried to increase the current through the ryanodine receptor (RyR) to obtain larger FWHM, but the spark amplitude also increased (approx. 10 times), which was far beyond the experimental measurements. Kong *et al*. [18] explained the FWHM paradox experimentally by using microscope blurring. On the other hand, Tan *et al*. [19–22] developed an anomalous subdiffusion model of Ca^{2+} sparks to explain FWHM paradox theoretically. Physiologically, the evidence for natural phenomena exhibiting anomalous diffusion had grown so compelling as to prompt punch lines such as ‘anomalous is normal’ [23]. Recently, anomalous diffusion processes in the cytoplasm were observed [24], which supported the anomalous subdiffusion of Ca^{2+} in cytoplasm.

A Ca^{2+} wave is formed by the propagation of Ca^{2+} sparks such that it may also obey anomalous subdiffusion. The previous mathematical models on Ca^{2+} waves were, however, based on Fick's law [12,25–28]. Backx *et al*. [25] presented a model of Ca^{2+} waves based on the uniform distribution of Ca^{2+} release sites throughout the cytoplasm. Anisotropic Ca^{2+} diffusion was studied by Girard *et al*. [26]. Keizer *et al*. [27] investigated Ca^{2+} waves under stochastic firing of Ca^{2+} release units (CRUs). Izu *et al*. [28] combined large CRU currents, stochastic firing of CRUs, asymmetric distribution of CRUs and anisotropic Ca^{2+} diffusion. Lu *et al*. [12] studied the effect of rogue RyRs on Ca^{2+} waves in ventricular myocytes with heart failure. By contrast, based on the anomalous subdiffusion of Ca^{2+} sparks, we have recently proposed a mathematical model for Ca^{2+} waves and predicted its propagation that agreed with the measurement in local cardiac myocytes [29]. However, the propagation of Ca^{2+} waves in full-width cardiac myocytes and the effects of various physiological parameters are still unknown.

In this work, based on anomalous subdiffusion of Ca^{2+} sparks, we developed a mathematical two-dimensional model for Ca^{2+} waves in full-width cardiac myocytes. We experimentally and theoretically investigated Ca^{2+} waves initiating from one transverse row of CRUs near the boundary or at the middle of the full-width cardiac myocytes. Moreover, we carried out sensitivity analysis of various physiological parameters, such as current through the CRU (*I*_{CRU}), duration of current flow through CRU (*T*_{open}), Ca^{2+} sensitivity parameter (*K*), spatial separation of CRUs (*l _{x}*,

*l*) and Ca

_{y}^{2+}diffusion coefficients (

*D*,

_{x}*D*). It was found that the anomalous diffusion mechanism of Ca

_{y}^{2+}sparks significantly affects the properties of Ca

^{2+}waves, given the difference between the present model and the Fick's law model.

## 2. Material and methods

### 2.1. Experimental methods

In order to record calcium waves in full-width cardiac myocytes, ventricular cardiac myocytes were isolated from adult Sprague–Dawley (rats two to three months old, weight 225–300 g) using standard enzymatic techniques, as described previously [30]. Freshly isolated single cells were stored in Tyrode's solution containing (in millimolar) 135 NaCl, 1 CaCl_{2}, 4 KCl, 1 MgCl_{2}, 10 glucose, 10 HEPES and pH 7.4. Before imaging, isolated myocytes were loaded with the dye Fluo-4-AM for 5 min, then washed twice with normal Tyrode's solution. Two-dimensional imaging used a Zeiss LSM-5Live fast confocal microscope equipped with a 40 × 1.3 NA oil immersion objective at a sample rate of 100 frames s^{−1}. The myocytes were excited at 488 nm. All experiments were performed at room temperature (23–25°C).

### 2.2. Anomalous diffusion model for calcium waves

The size of a normal cardiac myocyte is approximately 15 × 30 × 140 μm. Rapid three-dimensional confocal imaging of calcium waves shows that there is little variation of the wave in the *z*-direction [31]. Similar to Izu *et al*. [28], we eliminate the gradients along the *z*-direction because of the small thickness and establish line sources at the *z*-direction to simplify the three-dimensional problem to a two-dimensional model of a cardiac myocyte. The regular intervals of CRUs are *l _{x}* along longitudinal direction (

*x*-axis) and

*l*along transverse direction (

_{y}*y*-axis). The governing equation for Ca

^{2+}waves based on the anomalous subdiffusion model of Ca

^{2+}sparks can be expressed as 2.1where [Ca

^{2+}] is the free Ca

^{2+}concentration;

*D*and

_{x}*D*are diffusion coefficients, for anisotropic diffusion

_{y}*D*= 0.30 μm

_{x}^{2}ms

^{–1}and

*D*= 0.15 μm

_{y}^{2}ms

^{–1}[28,32];

*J*

_{dye}and

*J*

_{buffer}are fluxes due to Ca

^{2+}fluorescent indicator dye and endogenous stationary buffers;

*J*

_{pump}is pumping rate of SR Ca

^{2+}–ATPase and

*J*

_{leak}is a SR leak that balances

*J*

_{pump}. The expressions of

*J*

_{dye},

*J*

_{buffer},

*J*

_{pump}and

*J*

_{leak}are 2.2 2.3 2.4 2.5 2.6 2.7where

*n*identifies each of the buffer species. [

*F*]

*and [*

_{T}*B*]

_{n}*represent total concentration of the indicator and buffers, respectively. [CaF] and [CaB*

_{T}*] are the concentrations of the Ca*

_{n}^{2+}-bound complexes. are reaction kinetic parameters.

*K*

_{pump}is the affinity constant for SR pumps,

*m*is the Hill constant and is the maximum rate. SR leak is used to balance

*J*

_{pump}in the resting state.

*J*_{CRU} is the flux of Ca^{2+} release from CRU, the expression of which is the same as that of Izu *et al*. [28],
2.8where *σ* = 0.64 *I*_{CRU}/2*F* is a molar flux of a clustered RyR channel (*I*_{CRU} is the current through the CRU and *F* is the Faraday's constant), *δ* is the Dirac delta function, *S* is a stochastic function that controls the firing of the CRU, and *T*_{open} the firing time. Within a time interval *Δ**t*, the probability that the CRU fires is *P*(*C*(*x*, *y*, *t*), *K*, *n*) · *Δ**t*, where *P*(*C*(*x*, *y*, *t*), *K*, *n*) = *P*_{max}[Ca^{2+}]* ^{n}*/(

*K*+ [Ca

^{n}^{2+}]

*) with*

^{n}*P*

_{max}= 0.3 per CRU ms

^{−1}the maximum probability of Ca

^{2+}spark occurrence,

*K*is the Ca

^{2+}sensitivity parameter and

*n*= 1.6 the Hill coefficient.

The anomalous diffusion model is used in equation (2.1) where are the Riemann–Liouville operators defined as
2.9
2.10in which *k* is an integer with *β* < *k* < *β* + 1 (1 < *β* < 3). 1 < *β* < 2 refers to anomalous space superdiffusion; *β* = 2 corresponds to Fick's law; 2 < *β* < 3 represents anomalous space subdiffusion. According to Tan *et al*. [19], Ca^{2+} sparks follow the anomalous subdiffusion of *β* = 2.25, so this value is also used here, because Ca^{2+} waves are the propagation of Ca^{2+} sparks.

### 2.3. Numerical methods

The simulation is performed on a rectangular region with size of 100 × 20 μm, which is meshed with a uniform grid size of 0.1 × 0.1 μm. The time-step size is 0.005 ms. In order to site a CRU at the middle of the cardiac myocyte (the point (50.0, 10.0)), the position of the CRU at the left corner is sited at the point (2.0, 0.4).

For the fractional differential term, the right-shifted Grünwald formula is used to make a finite difference approximation [33],
2.11
2.12where *N _{x}* and

*N*are positive integers, and

_{y}*h*=

_{x}*x*/

*N*,

_{x}*h*=

_{y}*y*/

*N*.

_{y}*Γ*denotes the Gamma function. The shifted Grünwald approximation for fractional order derivative is unconditionally stable [33].

Considering simple impermeability of the cell boundary to the diffusing ions, reflecting boundary conditions *∂C*/*∂n*|_{boundary} = 0 are taken on all edges [12]. The computation time is 500–1000 ms. Therefore, we assume that a CRU will not reopen after firing and closing (a long refractory period), then the effect of the reopened CRUs can be ignored.

Standard values of parameters used in the present study are listed in tables 1 and 2 [13]. *I*_{CRU}, *T*_{open}, *D _{x}*,

*D*,

_{y}*l*,

_{x}*l*and

_{y}*K*are variable parameters whose effects on the results will be discussed. Their default values are

*T*

_{open}= 10 ms,

*D*= 0.30 μm

_{x}^{2}ms

^{–1},

*D*= 0.15 μm

_{y}^{2}ms

^{–1},

*l*= 2.0 μm,

_{x}*l*= 0.8 μm and

_{y}*K*= 15 μM. The default value of

*I*

_{CRU}will be determined later.

## 3. Results

### 3.1. Experimental and simulated calcium waves

Because the spatial separation of CRUs along the transverse direction (*l _{y}*) is shorter than that along the longitudinal direction (

*l*) in full-width cardiac myocytes, one or more pathological Ca

_{x}^{2+}sparks can rapidly trigger the CRUs along the transverse direction. Figure 1 shows that the initial sparks are always located on a transverse line of a cardiac myocyte. Accordingly, simulations only consider the longitudinal propagation (wave velocity) of Ca

^{2+}waves and the amplitudes. The firing of CRUs in a transverse row is chosen to be the initial condition in our simulations to reduce the effects of the transverse propagation and the transverse reflecting boundaries. Here, we will study how large one transverse row of anomalous Ca

^{2+}sparks (also to determine the default value of

*I*

_{CRU}) can induce a Ca

^{2+}wave with experimental wave velocity and how the Ca

^{2+}waves propagate from the boundary or the middle of the full-width cardiac myocytes. In addition, collision of full-width Ca

^{2+}waves will also be presented. Because the experimental Ca

^{2+}waves are always spontaneous and stochastic, it is difficult to observe the Ca

^{2+}waves which exactly initiate from the boundaries or the middle of a cardiac myocyte. The initiation of experimental Ca

^{2+}waves is the firing of the CRUs near the boundaries or the middle of a cardiac myocyte.

#### 3.1.1. Calcium waves initiating from the boundary

Figure 1 shows experimental two-dimensional images of a Ca^{2+} wave initiating from the boundary of a cardiac myocyte. At the beginning, the CRUs near the left boundary of the cardiac myocyte fire in rapid succession. In a short time, a high Ca^{2+} concentration region is formed. The transverse propagation of Ca^{2+} waves in full-width cardiac myocytes cannot be observed (just propagate for a short time), the wavefronts are always a transverse line. From snapshots 1 to 8, the Ca^{2+} wave propagates from the left boundary to the right boundary with a constant longitudinal velocity of *v _{x}* = 115 μm s

^{−1}, and the amplitudes of each CRUs are also nearly the same. The results are similar to those in previous reports [34,35]. In snapshot 8, the Ca

^{2+}wave reaches the right boundary of the cardiac myocyte.

As shown in figure 2, when *I*_{CRU} = 4 pA, the firing of CRUs in a transverse row near the boundary can trigger a longitudinal propagating Ca^{2+} wave. Its longitudinal wave velocity (*v _{x}* = 120 μm s

^{−1}) is in good agreement with the experimental result. Other parameters were taken as their default values (

*T*

_{open}= 10 ms,

*D*= 0.30 μm

_{x}^{2}ms

^{–1}and

*D*= 0.15 μm

_{y}^{2}ms

^{–1},

*l*= 20 μm,

_{x}*l*= 0.8 μm and

_{y}*K*= 15 μM). Initial calcium sparks are formed by 10 ms opening of all the CRUs along the transverse line of

*x*= 20 μm. Snapshots are taken at 10, 170, 330, 490, 650 and 810 ms (from 1 to 6). In the three-dimensional images,

*x*-axis and

*y*-axis denote the longitudinal direction and the transverse direction of the cardiac myocyte, respectively, and

*z*-axis denotes the Ca

^{2+}concentration. Through the propagation of the Ca

^{2+}wave, the wavefront nearly maintains a transverse line. However, the stochastic firing of CRUs leads to the peak Ca

^{2+}concentration value in the same transverse row in a short succession. Hence, there are different Ca

^{2+}concentration peaks near the wavefronts. Because the mean waiting time (MWT) before firing of the CRU apart

*l*through longitudinal direction and its standard deviation (s.d.) is small (16.67 ± 4.88 ms), Ca

_{x}^{2+}concentration peaks cannot be observed far away from the wavefronts. Here, the MWT and s.d. are smaller than that based on Fick's law. The physical reason is that the high Ca

^{2+}concentration may stay in a larger region around the firing CRU (for one Ca

^{2+}spark, FWHM is 2.0 μm for subdiffusion and 1.0 μm for Fick's diffusion), the firing probabilities of adjacent CRUs become higher, then the MWT and s.d. becomes smaller, it will be easier for the calcium wave to occur and sustain. The subdiffusive model can reproduce a global Ca

^{2+}wave with the experimental velocity using a reasonable CRU current (4 pA). Fick's law model needs a large CRU current (20 pA) to trigger a normal Ca

^{2+}wave [28]. The CRU current of 20 pA is hard to obtain ([36], figure2

*c*). Therefore, the propagation of Ca

^{2+}waves in whole cardiac myocytes may obey the anomalous subdiffusion mechanism such as the single Ca

^{2+}sparks.

When a transverse row of CRUs fires together initially, the Ca^{2+} sparks affect each other, then a high Ca^{2+} concentration band appears. The Ca^{2+} concentration at the next CRUs through the longitudinal direction increases. Therefore, the initiation of Ca^{2+} sparks triggered by a transverse row of CRUs needs a smaller *I*_{CRU} to form a Ca^{2+} wave with the experimental longitudinal velocity. Four picoamps is chosen to be the default value of *I*_{CRU} in the following work if not mentioned otherwise. At *t* = 810 ms, the Ca^{2+} wave propagates to the right boundary of the cardiac myocyte, which is similar to our experimental result (approx. 800 ms).

Figure 3 shows time courses of Ca^{2+} concentration in the Ca^{2+} wave in figure 2. The locations are a (50.0, 10.0), b (50.0, 11.6), c (50.0, 13.2), d (48.0, 10.0) and e (46.0, 10.0). The firing time and the amplitudes of transverse adjacent CRUs at points a, b and c are very close each other. The spatial interval between CRUs at a, d and e is *l _{x}*. The firing time interval between the CRUs at a and d is 13.92 ms, and that between the CRUs at d and e is 17.60 ms, which are both in the range of MWT and s.d. (16.67 ± 4.88 ms). It is shown that the longitudinal velocity of the propagating Ca

^{2+}wave is almost a constant and all the Ca

^{2+}sparks are almost the same. Moreover, the amplitudes of all the CRUs are close to the average amplitude of the Ca

^{2+}wave (80 μM). It is shown that the amplitude of the propagating Ca

^{2+}wave is also almost a constant.

#### 3.1.2. Calcium waves initiating from the middle

It is experimentally and numerically found that Ca^{2+} waves occur more readily at the boundaries of cardiac myocytes [29,34,37,38]. The physical reason is that the impermeable boundaries of cardiac myocytes can increase the initial local Ca^{2+} concentration in a small region, then the firing probability of adjacent CRUs will rise. However, the initial conditions in our simulations were firing of all the CRUs in a transverse row. Do the reflecting boundary conditions still affect the longitudinal velocity of Ca^{2+} waves effectively?

Figure 4 shows experimental images of a Ca^{2+} wave initiating near the middle of a cardiac myocyte. In snapshot 6, the Ca^{2+} wave reaches the right boundary, but the left wave keeps on propagating. In snapshot 8, the Ca^{2+} wave reaches the left boundary. The longitudinal wave velocity is *v _{x}* = 117 μm s

^{−1}which is close to that of the wave starting from the boundary. As shown in figure 5, when

*I*

_{CRU}= 4 pA, the firing of CRUs along the transverse line of

*x*= 50.0 μm can trigger a Ca

^{2+}wave with

*v*= 118 μm s

_{x}^{−1}. Snapshots are taken at 10, 90, 170, 250, 330 and 410 ms (from 1 to 6). The average amplitude of the Ca

^{2+}wave is 78 μM. The longitudinal velocity and average amplitude of the Ca

^{2+}wave are only a little smaller than that starting from the boundary. The effect of the reflecting boundaries is very weak, because the initial firing of the transverse CRUs raises the local Ca

^{2+}concentration considerably. Thus, when physiological parameters are determined, the velocities and amplitudes of full-width Ca

^{2+}waves are also determined.

#### 3.1.3. Collision of full-width calcium waves

Figure 6 shows collision of two full-width Ca^{2+} waves. They are triggered initially by 10 ms opening of CRUs along the transverse lines of *x* = 20 and *x* = 98.0 μm. The two opposite Ca^{2+} waves travel at the longitudinal velocities of *v _{x}* = 115 μm s

^{−1}from

*x*= 2.0 μm and

*v*= 116 μm s

_{x}^{−1}from

*x*= 98.0 μm. The two waves collide at the centreline of

*x*= 50.0 μm at

*t*= 420 ms, and then annihilate, which is consistent with previous experimental results [38–40]. Time courses of Ca

^{2+}concentration at locations a (48.0, 10.0), b (50.0, 10.0) and c (50.0, 11.6) are displayed in figure 7. Because Ca

^{2+}waves become global and fully developed, the amplitude of the combined wave (79 μM) is similar to that of the two waves before collision (78 μM). This is because the local amplitude of a Ca

^{2+}wave is related to the intensity of Ca

^{2+}sparks and the effect of adjacent sparks. There are no differences between the time course or magnitude of the Ca

^{2+}wave during its propagation and at the point of collision.

### 3.2. Physiological parameter sensitivity analysis

Physiological parameter sensitivity analysis of Ca^{2+} waves based on the Fick's diffusion has been studied by many researchers [26–28]. However, there is no parameter sensitivity analysis of Ca^{2+} waves based on the anomalous subdiffusion. Here, based on anomalous subdiffusion of Ca^{2+} sparks, we perform sensitivity analysis of several parameters to identify the possible determinants of full-width Ca^{2+} waves. The default values are *T*_{open} = 10 ms, *D _{x}* = 0.30 μm

^{2}ms

^{–1}and

*D*= 0.15 μm

_{y}^{2}ms

^{–1},

*l*= 20 μm,

_{x}*l*= 0.8 μm,

_{y}*K*= 15 μM and

*I*

_{CRU}= 4 pA. The initial conditions are all the CRUs along the transverse line of

*x*= 2.0 μm open at the time of

*T*

_{open}.

#### 3.2.1. Effect of *I*_{CRU}

The total amount of Ca^{2+} release per CRU is determined by the current through CRU (*I*_{CRU}) and duration of current flow through CRU (*T*_{open}). When *T*_{open} is invariable, *I*_{CRU} determines the spatio-temporal property of each Ca^{2+} spark. Thus, how *I*_{CRU} affects the properties of Ca^{2+} waves is very important to the stability of cardiac myocytes.

Figure 8 shows the effect of *I*_{CRU} on the longitudinal velocity of the Ca^{2+} wave. *I*_{CRU} changes from 2 to 5.5 pA. For subdiffusive model, *I*_{CRU} = 4 pA is a inflection point. When *I*_{CRU} increases from 2 to 4 pA, the curve is a concave function and its slope increases with *I*_{CRU}. *v _{x}* increases rapidly with

*I*

_{CRU}near the value of 2 pA, but slowly near the value of 4 pA. When

*I*

_{CRU}increases from 4 to 5.5 pA, the curve is a convex function and its slope decreases with

*I*

_{CRU},

*v*increases rapidly with

_{x}*I*

_{CRU}again when

*I*

_{CRU}is beyond 4 pA. Physiologically, 2 pA is the normal physiological current through CRUs, whereas 4 pA is the release current of a Ca

^{2+}spark [1,36,41,42]. When

*I*

_{CRU}takes the normal physiological current through CRUs (2 pA), it can only form a weak Ca

^{2+}wave with

*v*= 25 μm s

_{x}^{−1}. The intensity of Ca

^{2+}release per CRU is so small that the longitudinal velocity of the Ca

^{2+}wave is rarely low.

*I*

_{CRU}becomes the determinant of

*v*, so

_{x}*v*increases rapidly with

_{x}*I*

_{CRU}near 2 pA. Accompanied by the enlargement of the spatial width, the duration of one subdiffusive Ca

^{2+}spark is shortened [19]. When

*I*

_{CRU}is near 4 pA, the intensity of Ca

^{2+}release is not so large that the short duration of the subdiffusive spark becomes the most significant effect on the increase of adjacent CRU's firing probabilities. Then,

*v*increases slowly with

_{x}*I*

_{CRU}. When the intensity of Ca

^{2+}release is large enough, it connects with the large spatial width of the subdiffusive spark, making

*v*increase rapidly again. Close to the

_{x}*I*

_{CRU}value of 4 pA, the variation of

*v*is nearly robust with

_{x}*I*

_{CRU}, Ca

^{2+}sparks and Ca

^{2+}waves are relatively stable with the increase of

*I*

_{CRU}, which may be the physical reason why 4 pA is always the release current of a Ca

^{2+}spark [1,36,41,42]. However, different from our results,

*v*is linear with

_{x}*I*

_{CRU}in Fick's model.

The amplitudes are almost linear for the subdiffusive model and Fick's model. The diffusion model of Ca^{2+} sparks is not the determinant of the amplitude in Ca^{2+} waves. The spatio-temporal properties of anomalous sparks do not affect the amplitude directly. The determinant of the amplitude is the total amount of Ca^{2+} release per CRU (here is *I*_{CRU}), and the effect is almost linear. When *I*_{CRU} = 5 pA, the average amplitude of the Ca^{2+} wave triggered by several sparks is 115 μM, which is close to that of the local Ca^{2+} wave triggered by one single spark (110 μM [29]). Therefore, based on anomalous diffusion of Ca^{2+} sparks, initial firing number of CRUs barely affects the average amplitude of Ca^{2+} waves. However, the longitudinal velocity *v _{x}* increases obviously with the initial firing number of CRUs. When

*I*

_{CRU}= 5 pA, for one transverse row of sparks

*v*= 160 μm s

_{x}^{−1}, and for one single spark

*v*= 96 μm s

_{x}^{−1}[29]. The physical reason may be that, the synergistic effect of CRUs induces a large advancement of the firing probability of adjacent CRUs, the wave velocity is larger than that triggered by one single spark. In general, the synergistic effect of CRUs only affects the wave velocity but not the amplitude, and the simple summation of Ca

^{2+}sparks does not exist.

#### 3.2.2. Effect of *T*_{open}

Another determinant to the total amount of Ca^{2+} release per CRU is duration of current flow through CRU (*T*_{open}). Thus, *T*_{open} directly determines the intensity of a Ca^{2+} spark, then determines the amplitude of the Ca^{2+} wave. On the other hand, when *T*_{open} is small, the duration of high Ca^{2+} concentration at the adjacent CRUs is very short, so the firing probability will be low. It may cause the small wave velocity.

We study the propagation of full-width Ca^{2+} waves with different values of *T*_{open} (*T*_{open} = 5, *T*_{open} = 8 and *T*_{open} = 15 ms). When *T*_{open} reduces to one half of the standard value *T*_{open} = 10 ms, the full duration at half maximum (FDHM), the FWHM and the amplitude of each Ca^{2+} spark become smaller obviously. The MWT and s.d. for this Ca^{2+} wave is 128 ± 90 ms. A slow propagating Ca^{2+} wave is formed, the longitudinal velocity (15.6 μm s^{−1}) and amplitude (30 μM) are both very small. When *T*_{open} = 8 ms, the Ca^{2+} wave is close to the standard one (*v _{x}* = 75 μm s

^{−1}and amplitude is 65 μM). The Ca

^{2+}wave with a longer duration of sparks (

*T*

_{open}= 15 ms) propagates faster (

*v*= 160 μm s

_{x}^{−1}), and its amplitude also rises obviously (126 μM).

Similar to the effect of *I*_{CRU}, the average amplitude almost linearly increases with *T*_{open}. It also shows that the diffusion model of Ca^{2+} sparks is not the determinant of the amplitude in Ca^{2+} waves. However, the increase of *v _{x}* has an obvious inflection point at

*T*

_{open}= 10 ms. The physical reason is that although each subdiffusive spark reaches a wide region, the duration is very short. When

*T*

_{open}becomes smaller, the duration of each spark will be shorter, and the firing probability of adjacent CRUs will be lower, then the longitudinal velocity will be smaller. So near the small value of

*T*

_{open}, the duration of each spark is the determinant of

*v*. When

_{x}*T*

_{open}is close to the normal duration of current flow through CRU (10 ms), the spatial property of the subdiffusive spark determines the increase of

*v*, and the increase of

_{x}*v*is relatively stable. With a larger

_{x}*T*

_{open}, the anomalous diffusion model is not the determinant of

*v*any more under the restrain function of the other physiological parameters, such as

_{x}*I*

_{CRU}and

*l*. Then

_{x}*v*will increase at a slow and stable rate.

_{x}#### 3.2.3. Effect of total calcium release

The total amount of Ca^{2+} release per CRU can be written as *I*_{CRU} × *T*_{open}, which refers to the intensity of each Ca^{2+} spark. Based on Fick's diffusion of Ca^{2+} sparks, *I*_{CRU} × *T*_{open} is found to be a strong determinant of wave properties, such as the sharpness of the wavefront [28]. Whether it is also the determinant of wave properties based on anomalous subdiffusion of Ca^{2+} sparks requires further study.

The full-width calcium wave with the parameters of *I*_{CRU} = 5 pA and *T*_{open} = 8 ms is studied. Table 3 shows the longitudinal velocities and the average amplitudes of two cases with the same total calcium release. The values of (*I*_{CRU}, *T*_{open}) are (4, 10) and (5, 8), respectively. *v _{x}* for (5, 8) is 5% larger than that for (4, 10), and the amplitude for (5, 8) is 14% larger than that for (4, 10). Their longitudinal velocities are very close, but the amplitudes have a little difference. For the morphology of a single Ca

^{2+}spark, the current

*I*

_{CRU}denotes the amplitude of the spark, and

*T*

_{open}denotes the duration of the spark. The spatial width of anomalous Ca

^{2+}spark affects the properties of Ca

^{2+}waves significantly. Therefore, the amplitude of one spark affects the amplitude and velocity of the wave more significantly than the duration. (Here, the precondition is that

*T*

_{open}is close to the standard value, if it is too small, then the effect will be more significant.) However, the enlargements of 5% for

*v*and 14% for amplitude are not notable. It can be confirmed that the total amount of Ca

_{x}^{2+}release per CRU

*I*

_{CRU}×

*T*

_{open}is the determinant of wave properties based on anomalous subdiffusion of Ca

^{2+}sparks.

Physiologically, *I*_{CRU} and *T*_{open} are related to release intensity of RyRs. Being a protein, release intensity of RyRs can be controlled by some diseases or by some drugs. Therefore, our results about the sensitivities of *I*_{CRU} and *T*_{open} based on the anomalous subdiffusion may be used to predict the pathological transformation of cardiac myocytes.

#### 3.2.4. Effect of Ca^{2+} release unit calcium sensitivity

CRU Ca^{2+} sensitivity parameter *K* denotes the sensitivity of high Ca^{2+} concentration to the firing probability of the CRU. Physically, when *K* becomes smaller, the CRU will be easier to fire under high Ca^{2+} concentration, then the longitudinal velocity will be larger. Here, based on the anomalous subdiffusion of Ca^{2+}, the firing probability of the adjacent CRU is already high. In that case, is the effect of CRU Ca^{2+} sensitivity on the properties of Ca^{2+} waves also obvious?

The quantitative analysis of the effect of *K* on the longitudinal velocities and the amplitudes for both subdiffusive model and Fick's model is performed in figure 9. The result with *K* = 15 μM is chosen to be the standard value of 1, and the normalized response is showed here. Based on the subdiffusive model, for *K* = 7.5 μM which reduces to one half of the standard value *K* = 15 μM, the MWT and s.d. of the Ca^{2+} wave is 16.13 ± 4.67 ms. The MWT is longer than that of *K* = 15 μM (16.67 ms), and the s.d. is smaller than that of *K* = 15 μM (4.88 ms). The firing probability of the adjacent CRU is enlarged. Although the longitudinal velocity decreases with *K*, the change is not obvious. *v _{x}* for

*K*= 7.5 μM is only 2.5% larger than that for

*K*= 15 μM. However, for Fick's model,

*v*decreases with

_{x}*K*significantly.

*v*for

_{x}*K*= 7.5 μM is 38% larger than that for

*K*= 15 μM, which agrees with the result performed by Izu

*et al*. [28] (figure 7, curve a). Based on the subdiffusive model, the spatial width of one Ca

^{2+}spark is much larger than that based on Fick's law such that the firing probability of the adjacent CRU is much larger. Under a high firing probability, when

*K*changes, the firing probability will not change obviously.

On the other hand, the average amplitude of the Ca^{2+} wave barely changes with *K*. That is to say, the amplitude of the full-width Ca^{2+} wave is determined by the total amount of Ca^{2+} release per CRU but not the firing probability of the adjacent CRU.

In some extreme conditions, for example, *K* decreases to near 0.1 μM or increases to ∞, the firing probability will apparently change. Then, the Ca^{2+} wave will become so easy to propagate (*v _{x}* becomes so large) or the Ca

^{2+}wave cannot propagate (the firing probability tends to 0). These conditions can be triggered by some drugs or some diseases.

Physiologically, *K* can be affected by the Mg^{2+} concentration, drugs and SR Ca^{2+} loading [43,44]. According to our results, near the standard value the effect of *K* on the anomalous subdiffusive wave is very weak. In other words, when the abnormal conditions are not excessive, the cardiac myocytes are relatively stable.

#### 3.2.5. Effect of longitudinal and transverse lattice spacing

All the subdiffusive Ca^{2+} sparks have almost the same spatial width (FWHM) and temporal duration (FDHM). Intuitively, when the longitudinal or transverse lattice spacing of CRUs becomes smaller, the effect of one spark on adjacent CRUs will be more significant, higher Ca^{2+} concentrations will reach the adjacent CRUs and increase their firing probabilities. Then, the local velocities of the Ca^{2+} wave will increase. For the local amplitude of each CRU, the effect of the adjacent sparks will induce a high local Ca^{2+} concentration and enlarge the local amplitude. The longitudinal lattice spacing of CRUs (*l _{x}*) affects the longitudinal velocity (

*v*) directly, and the transverse lattice spacing of CRUs (

_{x}*l*) affects the transverse velocity (

_{y}*v*) directly. Here, we consider only

_{y}*v*. There is a question, does the change of

_{x}*l*affect

_{y}*v*? If yes, is the effect weak or considerable?

_{x}The full-width Ca^{2+} waves with different values of *l _{x}* and

*l*are studied, in which (

_{y}*l*= 1.0 μm,

_{x}*l*= 0.8 μm), (

_{y}*l*= 2.0 μm,

_{x}*l*= 0.4 μm) and (

_{y}*l*= 2.0 μm,

_{x}*l*= 1.2 μm). Table 4 presents the effect of (

_{y}*l*,

_{x}*l*) on the longitudinal velocity and average amplitude. The standard

_{y}*v*and amplitude (

_{x}*l*= 2.0 μm,

_{x}*l*= 0.8 μm) are shown for comparison.

_{y}When the intervals of CRUs along the longitudinal direction reduce to one half of the standard value *l _{x}* = 2.0 μm, the longitudinal velocity increases obviously, which becomes almost two times as standard value (

*v*for

_{x}*l*= 2.0 μm is 120 μm s

_{x}^{−1}and for

*l*= 1.0 μm is 220 μm s

_{x}^{−1}). The average amplitude also increases obviously (amplitude for

*l*= 2.0 μm is 80 μM and for

_{x}*l*= 1.0 μm is 142 μM). The enlargement of

_{x}*v*is easy to understand as explained before. The enlargement of average amplitude is due to the FWHM of anomalous sparks (approx. 2.0 μm). When the longitudinal intervals of CRUs reduce to 1.0 μm, the local amplitude of each CRU will be affected by several adjacent sparks, the number of which will be more than the situation of

_{x}*l*= 2.0 μm. Then, the average amplitude of the Ca

_{x}^{2+}wave will increase significantly.

When *l _{x}* = 2.0 μm and

*l*= 0.4 μm, the longitudinal velocities (

_{y}*v*= 217 μm s

_{x}^{−1}) are almost the same as that under the condition (

*l*= 1.0 μm and

_{x}*l*= 0.8 μm). The decrease of

_{y}*l*induces the superposition of more transverse sparks, the Ca

_{y}^{2+}concentration at the location of adjacent CRU will be higher, the duration of which will be longer. Then, the firing probability of adjacent CRU will be larger. For the average amplitude, the decrease of

*l*induces more obvious increase of the amplitude (160 μM for (2.0, 0.4) and 142 μM for (1.0, 0.8)). This conclusion is related to the physiological structure of cardiac myocytes. The transverse lattice spacing of CRUs is physiologically smaller than the longitudinal spacing. When the interval of CRUs is smaller, the effect of adjacent sparks on one site of CRU will be more significantly and then the average amplitude will be much larger.

_{y}When *l _{y}* increases to 1.2 μm, the properties of the Ca

^{2+}wave have obvious change. The connection of transverse CRUs is weakened because of the larger

*l*, the firing probability of adjacent CRU descends considerably, and

_{y}*v*decreases significantly. On the other hand, owing to the relative independence of each Ca

_{x}^{2+}spark, the effect of adjacent sparks on the local amplitude at each site of CRU is very weak. Consequently, the average amplitude also decreases.

Physiologically, if the longitudinal or transverse interval of CRUs becomes smaller owing to some reasons, such as cardiac myocytes deformation, then Ca^{2+} waves will easily occur with large amplitude, which then causes cardiac myocytes to be unstable. By contrast, if the lattice spacing of CRUs can be enlarged by some method, then we can restrain the propagation of Ca^{2+} waves and cardiac myocytes will be stable.

#### 3.2.6. Effect of diffusion coefficient

In the above work, the diffusion of Ca^{2+} is anisotropic. The diffusion coefficient along the *x*-axis (*D _{x}* = 0.30 μm

^{2}ms

^{–1}) is two times that along the

*y*-axis (

*D*= 0.15 μm

_{y}^{2}ms

^{–1}). It corresponds to the physiological structure of cardiac myocytes. Under some pathological conditions, such as deformation of the cell or change in protein conformation, the diffusion coefficient may change. Here, we investigate the Ca

^{2+}wave based on the isotropic diffusion of Ca

^{2+}.

The full-width Ca^{2+} wave based on the isotropic diffusion of Ca^{2+} is studied. The parameters are *D _{x}* = 0.25 μm

^{2}ms

^{–1}and

*D*= 0.25 μm

_{y}^{2}ms

^{–1}. Table 5 gives the longitudinal velocities and average amplitudes of the two cases (anisotropic diffusion and isotropic diffusion). The longitudinal velocity decreases with the decrease in the longitudinal diffusion coefficient, whereas the average amplitude increases with the increase in the transverse diffusion coefficient. This conclusion is similar to the effect of longitudinal and transverse lattice spacing. Because the transverse interval of CRUs is small, the increase of

*D*enlarges the amplitude significantly. Owing to the large longitudinal interval of CRUs, the decrease of

_{y}*D*diminishes the firing probability of adjacent CRU considerably. Although the increase of

_{x}*D*also enlarges

_{y}*v*, the effect is not as obvious as

_{x}*D*.

_{x}## 4. Conclusion

Propagation of Ca^{2+} waves initiating from Ca^{2+} sparks in cardiac myocytes is related to ventricular arrhythmias, fibrillation or heart failure. It is important for our understanding of Ca^{2+} waves. Here, we carried out experimental and numerical studies to investigate Ca^{2+} waves in full-width cardiac myocytes. Based on anomalous subdiffusion of Ca^{2+} sparks, we proposed a mathematical model for Ca^{2+} waves. There was an good agreement between the experimental measurements and numerical predictions. It was found that one transverse row of anomalous Ca^{2+} sparks can trigger a Ca^{2+} wave with the experimental longitudinal velocity with 4 pA current through the CRU. The 4 pA current is more close to the physiological current than that in the Fick's law model. For the initial transverse row of Ca^{2+} sparks, the effect of the reflecting boundaries on the properties of full-width Ca^{2+} waves is very weak. The simulation of full-width Ca^{2+} waves collision revealed that the simple summation of the waves did not exist at collision. Moreover, we performed sensitivity analysis of physiological parameters for calcium waves in full-width cardiac myocytes.

The total amount of Ca^{2+} release per CRU can be written as *I*_{CRU} × *T*_{open}, which refers to the intensity of each Ca^{2+} spark. The determinant of *v _{x}* is

*I*

_{CRU};

*v*increases rapidly with

_{x}*I*

_{CRU}near the value of 2 pA; near the value of 4 pA,

*v*increases slowly with

_{x}*I*

_{CRU}owing to the short duration of the subdiffusive spark;

*v*increases rapidly again with

_{x}*I*

_{CRU}because of the large spatial width of the spark. On the other hand, the average amplitude is almost linear with

*I*

_{CRU}. The variation of

*v*with

_{x}*T*

_{open}has an inflection point at

*T*

_{open}= 10 ms which is also related to the spatio-temporal properties of the subdiffusive sparks, and the average amplitude is almost linear. The total amount of Ca

^{2+}release per CRU

*I*

_{CRU}×

*T*

_{open}is the determinant of wave properties based on anomalous subdiffusion of Ca

^{2+}sparks.

Based on the anomalous subdiffusion of Ca^{2+}, the firing probability of the adjacent CRU is very high. Hence, the effect of CRU calcium sensitivity *K* on the properties of Ca^{2+} waves is very weak, which is different from the predictions of the Fick's law model. The longitudinal and transverse lattice spacing of CRUs are both important to Ca^{2+} waves. The decrease of *l _{y}* can increase

*v*considerably, and

_{x}*l*affects the amplitude of Ca

_{y}^{2+}wave more significantly than

*l*owing to the physiological structure of cardiac myocytes. In the same way, the increase of

_{x}*D*enlarges the amplitude significantly, but the decrease of

_{y}*D*reduces

_{x}*v*significantly.

_{x}This study shows some implications to fibrilation and arrhythmias in cardiac myocytes. Ca^{2+} wave is related to pathologies such as fibrillation and arrhythmias [45,46]. Ca^{2+} alternans [47] and delayed after depolarization [48] triggered by Ca^{2+} waves are potentially arrhythmogenic abnormalities of cardiac Ca^{2+} signalling. In general, the amplitudes and wave velocities of Ca^{2+} waves based on the subdiffusive model are larger than those based on the Fick's law model under the same conditions. The abnormal Ca^{2+} signal caused the abnormal action potential duration and contraction amplitude. Fibrilation and arrhythmias may readily occur. In addition, the amplitudes and wave velocities of Ca^{2+} waves alter considerably with the variation of *I*_{CRU}, *T*_{open}, *l _{x}* and

*l*. This may significantly alter Ca

_{y}^{2+}alternans in whole cells, and then induce fibrillation and arrhythmias. By contrast, cardiac myocytes would be relatively stable when

*K*changes.

One limitation of the model is that *I*_{CRU} and *T*_{open} are fixed for all CRUs in a particular simulation. Because the wave properties are sensitive to *I*_{CRU} and *T*_{open}, some calcium waves may terminate after a short propagation owing to the decreases of *I*_{CRU} and *T*_{open} at propagation. We will carry out more experiments and simulations to investigate how the variations of *I*_{CRU} and *T*_{open} at propagation affect calcium waves in future studies.

Based on anomalous subdiffusion of Ca^{2+} sparks, the effects of physiological parameters qualitatively coincide with that based on Fick's diffusion. However, there exist quantitative differences between the two models. These numerical results need experimental confirmation in further studies before they can be used to predict the pathological mechanism of cardiac myocytes and prevent Ca^{2+} waves from propagating to the whole cardiac myocyte.

All rat experiments were carried out in strict accordance with the regulation of Research Ethics Committee of State Key Laboratory of Biomembrane and Membrane Biotechnology, Peking University.

## Funding statement

This work was supported by the National Natural Science Foundation of China (grant no. 11272014), and National Key Basic Research Programme of China (grant no. 2013CB531200).

- Received October 12, 2013.
- Accepted November 21, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.