Modeling Calcium Wave Based on Anomalous Subdiffusion of Calcium Sparks in Cardiac Myocytes

sparks and waves play important roles in calcium release and calcium propagation during the excitation-contraction (EC) coupling process in cardiac myocytes. Although the classical Fick’s law is widely used to model sparks and waves in cardiac myocytes, it fails to reasonably explain the full-width at half maximum(FWHM) paradox. However, the anomalous subdiffusion model successfully reproduces sparks of experimental results. In this paper, in the light of anomalous subdiffusion of sparks, we develop a mathematical model of calcium wave in cardiac myocytes by using stochastic release of release units (CRUs). Our model successfully reproduces calcium waves with physiological parameters. The results reveal how concentration waves propagate from an initial firing of one CRU at a corner or in the middle of considered region, answer how large in magnitude of an anomalous spark can induce a wave. With physiological currents (2pA) through CRUs, it is shown that an initial firing of four adjacent CRUs can form a wave. Furthermore, the phenomenon of calcium waves collision is also investigated.

Introduction Nomenclature x,y spatial coordinates, mm t time, ms b fractional order of the spatial derivative. D Cx , D Cy Ca 2z diffusion coefficients along x-axis and y-axis, open time of CRU, ms S stochastic switching function equaling either 0 or 1 K Ca 2z sensitivity parameter, mM P probability of Ca 2z spark occurrence,/ calcium release unit/ms P max maximum probability of Ca 2z spark occurrence,/calcium release unit/ms v x wave velocity along x-axis, mm=s v y wave velocity along y-axis, mm=s In the endoplasmic or sarcoplasmic reticulum(SR) of cardiac cells, there stores plenty of Ca 2z , the concentration of which is 2-3 orders of magnitude greater than that in the cytosol. During the excitation-contraction(EC) coupling process, triggered by L-type Ca 2z channels, Ca 2z is released from SR through ryanodine receptors(RyRs) on the z-lines [1][2][3][4], where RyR is one kind of Ca 2z release units(CRUs). This event is called ''Ca 2z spark''. Ca 2z -induced Ca 2z release(CICR) makes RyRs fire in succession such that Ca 2z concentration rises [1,5], the process of which is called calcium transient. Physiologically, calcium homeostasis is important for the contraction and relaxation of the heart muscle. However, in some pathological conditions, spontaneous propagating wave of Ca 2z may occur, which is called ''calcium wave''. The occurrence of calcium wave can affect the heart's normal function, and may induce some disease, such as ventricular arrhythmias [6].
The model of Ca 2z spark using Fick's Law failed to reproduce the full-width at half maximum(FWHM) of experimental results for Ca 2z sparks. Simulated results for Ca 2z spark based on Fick's Law presented a lower FWHM (,1.0 mm), which was only half the width of experimental result (,2.0 mm). Izu et al. [7] tried to increase the current through RyR to get larger FWHM, however, the spark amplitude also increased (,10 times), which is far beyond experimental results and physiological conditions. In contrast, the results obtained with the anomalous subdiffusion model of Ca 2z spark were found to be in close agreement with the experimental ones so that the ''FWHM Paradox'' was successfully explained [8][9][10]. Therefore, it is confirmed that diffusion of Ca 2z in cytoplasm obeys no longer Fick's Law, but the anomalous subdiffusion.
A Ca 2z wave is formed from propagation of Ca 2z sparks. According to the results for Ca 2z sparks, Ca 2z wave should also obey the anomalous subdiffusion. However, all previous work on Ca 2z waves were based on Fick's Law [11][12][13][14]. Anisotropic Ca 2z diffusion was studied by Girard et al. [11]. Keizer and Smith [12] investigated Ca 2z waves under stochastic firing of CRUs. Izu et al. [13] combined large CRU currents [7], stochastic firing of CRUs, asymmetric distribution of CRUs and anisotropic Ca 2z diffusion to investigate the propagation of Ca 2z waves. Lu et al. [14] studied the effect of rogue RyRs on Ca 2z waves in ventricular myocytes with heart failure.
In this work, we develop a mathematical 2D model based on anomalous subdiffusion of Ca 2z sparks. The anomalous subdiffusion model is used to study Ca 2z waves propagation from an initial firing of one CRU at a corner or in the middle of the considered region. We reproduce wave velocities of experimental results using a small current through CRUs which is close to the physiological conditions. The phenomenon of calcium waves collision is also investigated. With physiological Ca 2z currents(2pA) through CRUs, an initial firing of four adjacent CRUs is shown to form a Ca 2z wave. Furthermore, study on how the system becomes unstable is also performed by changing the transverse distance of CRUs. Figure 1 shows a 2-dimensional schematic of a cardiac myocyte(establishing line resources along z-axis [13]) which contains plenty of CRUs. The regular intervals of CRUs are l x along x-axis and l y along y-axis. The governing equation for Ca 2z waves based on the anomalous subdiffusion model can be expressed as

Anomalous Diffusion Model for Calcium Waves
where ½Ca 2z is the free Ca 2z concentration; D x and D y are diffusion coefficients for anisotropic diffusion with D x~0 :30mm 2 ms {1 and D y~0 :15mm 2 ms {1 [15]; J dye and J buffer are fluxes due to Ca 2z fluorescent indicator dye and endogenous stationary buffers; J pump is pumping rate of SR Ca 2z -ATPase, and J leak is a SR leak that is to balance J pump ; The expressions of J dye , J buffer , J pump , and J leak are L½CaB n Lt~k where n identifies each of buffer species. ½F T and ½B n T represent total concentration of the indicator and buffers, respectively. ½CaF and ½CaB n are concentration of the Ca 2z -bound complexes. k z F , k { F , k z n and k { n are reaction kinetics. K pump is the affinity constant for SR pumps, m the Hill constant, and V max pump the maximum rate. SR leak is used to balance J pump in resting state.
J CRU is the flux of Ca 2z release from CRU, the expression of which is the same as that of Izu et al. [13], where s~0:64I CRU =2F is a molar flux of a clustered RyR channel(I CRU is current through the CRU and F is Faraday's constant), and d the Dirac delta function, S a stochastic function which controls the firing of the CRU, and T open the firing time.
Within a time interval Dt, the probability that the CRU fires is P(C(x,y,t),K,n) : Dt, where P(C(x,y,t),K,n)~P max ½Ca 2z n =(K n z½Ca 2z n ) with P max~0 :3=CRU=ms the maximum probability of Ca 2z spark occurrence, K~15mM Ca 2z the sensitivity parameter and n~1:6 the Hill coefficient.
The anomalous space diffusion is model used in Eq.(1), where L b Lx b is the Riemann-Liouville operators which is defined as where k is an integer with bvkvbz1 (1vbv3). In Eq.(1), when 1vbv2, anomalous space superdiffusion occurs, while anomalous space subdiffusion occurs when 2vbv3. Particularly, our model reduces to Fick's Law when b~2. According to Li et al. [10], calcium sparks follow the anomalous space subdiffusion of b~2:25, so we only consider the space subdiffusion in the following sections.

Numerical Methods
Our simulation is performed on a rectangular region with size of 20mm|20mm, which is meshed with a uniform grid size of 0:1mm|0:1mm. The time-step size is 0.005ms. For the fractional differential term, we used the right-shifted Gr€ u unwald formula to make a finite difference approximation [16].
where N x and N y are positive integers, and h x~x =N x , h y~y =N y . C denotes the gamma function. The shifted Gr€ u unwald approximation for fractional order derivative has been shown to be unconditionally stable [16]. Considering simple impermeability of the cell boundary to the diffusing ions, reflecting boundary conditions LC=LnD boundary~0 are taken on all edges [14]. The scale of our computation time is 200-500ms so that a CRU would not reopen after firing and closing.
Standard values of parameters used in the current study are listed in Table 1 and Table 2. b, I CRU and l x are changeable parameters whose effects on the results will be investigated.

Results and Discussion
Modeling a Ca 2z wave from a Single Ca 2z Spark Ca 2z waves have been shown to be initiated and sustained by Ca 2z sparks [17]. Under pathological conditions, Ca 2z sparks fire spontaneously and stochastically, so whether a single spark can trigger a Ca 2z wave is important to the stability of cardiac myocytes [4]. Simulations based on Fick's Law reveal that large currents through CRUs and high calcium concentrations are needed to trigger a Ca 2z wave [13]. In this work, based on the anomalous subdiffusion model, we find the current which can trigger a normal Ca 2z wave initiating from a single Ca 2z spark at the corner of considered region.
According to Li et al. [10], calcium sparks follow the anomalous space subdiffusion of b~2:25, so this subdiffusion order is also taken in our simulation. In our model, initial source is a 10ms opening of one CRU, the longitudinal intervals are l x~2 :0mm. When we take I CRU~5 pA, the longitudinal wave velocity (v x~9 6mm=s initiating from the corner) is in good agreement with the experimental result(&100mm=s [17]). Figure 2 shows Ca 2z waves propagating on a discrete rectangle lattice initiating from a 10ms opening of the CRU at point (2,0.8). The snapshots are the Ca 2z concentration distribution at 10, 30, 50, 70, 90, 110, 130, 150, 170, 190 and 200ms (left to right, top to bottom). From image to image, we can see CRUs fire stochastically by turns while Ca 2z concentration wave propagates to the points of CRUs. At the beginning, CRUs fire one by one, and the amplitude(maximum of Ca 2z concentration in a region) is not very large; but with the increase of time, some of CRUs fire simultaneously in a short time so that Ca 2z sparks influence each other, and the amplitudes of Ca 2z sparks become larger and larger. For example, the CRU at (2,1.6) fires at t~16ms, the CRU at (2,2.4) fires at t~31ms; at t~67ms, the CRUs at (4,0.8) and (4,2.4) fires simultaneously, in a short time interval, at t~68ms the CRU at (2,4.8) fires. In image 10, sparks occur at (18,4.0), (18,4.8), (18,5.6), (18,6.4), (18,7.2), (18,8.0), (18,8.8) and (18,9.6) in rapid succession, which may trigger a calcium transient. Image 11 shows that the boundary limits the propagation of Ca 2z wave, but in an actual cardiac myocyte with size 100mm|20mm, calcium transient will be observed. In addition, though the space intervals of CRUs along y-axis are more compact than along x-axis, transverse wave velocity v y~8 2mm=s is smaller than that along xaxis v x (because the diffusion coefficient along y-axis D Cy~0 :15mm 2: ms {1 is smaller than that along x-axis D Cx~0 :30mm 2: ms {1 ).
Our numerical results show that a single Ca 2z spark can trigger a normal Ca 2z wave under the pathological condition of Because of the large FWHM for one spark due to anomalous subdiffusion, one firing CRU will trigger a Ca 2z wave. Then we prohibit the event of another spontaneous spark so that it will not affect the initial Ca 2z wave. In other words, when local Ca 2z concentration is larger than resting Ca 2z concentration, Ca 2z sparks may occur. So the ''wall'' of high Ca 2z concentration spreads from the left corner to the top of cardiac myocyte. Image 12 shows the sequence of CRU firing along y~2:4mm(the longitudinal linescan). The horizontal axis denotes time t(from left to right,200ms), and the vertical axis denotes spatial coordinate x(20mm). Except for initial two sparks, CRUs fire at nearly regular intervals, and the ''wall'' of high Ca 2z concentration is nearly a straight line. Here, from the initial spark to the second spark, it takes more time than those of the subsequent sparks, which is different from the results by Izu et al. [13](in their simulation, the transverse linescan is adopted, but the qualitative profile must be the same). It is because the subsequent sparks are triggered by two or more adjacent sparks, and the longitudinal wave velocity approaches to a constant value, but it takes more time to trigger the next CRUs from the initial signal spark.

Effect of the Anomalous Subdiffusion Order
The anomalous diffusion order b, which determines the diffusion mode of Ca 2z waves, was shown to affect the wave velocity considerably in last subsection(comparing with Fick's Law). When bw2, a large value of b means a wild spread of initial concentration, but with the increase of time, the remanent concentration at initial point will be smaller due to the wild spread of calcium concentration. So the anomalous diffusion order b affects not only the wave velocity, but also the amplitude of each CRU, and further the average amplitude of a Ca 2z wave. Here, b is taken to be 2.00, 2.05, 2.15 and 2.25 in order to figure out whether the amplitudes of Ca 2z waves will change obviously with the variation of wave velocities. The initial condition is still a 10ms opening of the CRU at the point (2,0.8). Table 3 presents the effect of anomalous fractional order b on the longitudinal wave velocity v x and the average amplitude, respectively. Comparing with the results based on Fick's law, velocities of Ca 2z waves increase considerably when anomalous subdiffusion order b becomes bigger. For b~2:25, the wave velocity(96mm=s) is almost twice as big as that based on Fick's law(b~2:00, 57mm=s). It is because FWHM along x-axis for b~2:25 (&2:2mm) is almost twice as large as that for b~2:00(&1:1mm). Here, in Ca 2z waves, FWHM for one CRU is affected by adjacent sparks, so it is a little bigger than FWHM of a single Ca 2z spark(1:97mm [10]). In contrast to wave velocity, the variation of amplitudes is not very considerable. For b~2:25, amplitude is 81% as that for b~2:00. The physical reason is that although for b~2:25, FWHM along x-axis is twice as big as that for b~2:00, the full duration at half maximum(FDHM) along xaxis for one spark still has a obvious decrease. So when the total release of Ca 2z concentration is almost the same, under the expansion of spatial affection and the decrease of temporal continuity, amplitude of Ca 2z waves does not decreases obviously. In addition, wave velocities and amplitudes do not vary linearly with b. When b is larger, the effect of subdiffusion on Ca 2z waves is greater. It is because when the variation of b is small, the other parameters, such as the speed of diffusion D Cx and the release strength of sparks I CRU , play important roles in Ca 2z waves. When b becomes bigger, replacing the primary position of the diffusion speed and release strength, diffusion mode affects Ca 2z waves significantly(wave velocities).

Effect of Initial Location
Propagation of a Ca 2z wave from a corner of the cardiac myocyte has been studied. It is found that the reflecting boundaries increase the amplitude of the initial spark, then further promote the propagation of the Ca 2z wave. In order to figure out the boundaries determine the propagation of the Ca 2z wave or just affect the wave velocity, we change the location of the initial Ca 2z spark and study how the reflecting boundaries affect Ca 2z waves. In general, the process in which more than two sparks firing together, then several Ca 2z waves propagating, meeting and dissipating is very common in cardiac myocytes. This event is called Ca 2z waves collision, and it was observed in experiments [17]. We will discuss in the following the interaction of several Ca 2z waves.
As shown in Fig. 3, it takes only 120ms for a Ca 2z wave initiating from a 10ms opening of the CRU at a middle point (10,9.6) to propagates to the left and the bottom boundaries. Triggering from the middle of the region, the ''walking distance'' of a Ca 2z wave becomes shorter, and it will spread more quickly to the boundary. Due to the shorter ''walking distance'', less sparks will occur simultaneously, and it will not make the Ca 2z wave develop sufficiently; but trigger from the corner, while the Ca 2z wave spreads wildly, large amount of sparks will fire together in a small region (Figure 2, right corner of Image 10). Comparing with the Ca 2z wave initiating from a corner without the effect of the boundaries, the initial concentration of the region will be smaller, then the probability of CRUs firing will be lower, so the events of Ca 2z sparks are more stochastic and irregular. Affected by the absence of the reflecting boundaries and the shorter ''walking distance'', longitudinal wave velocity v x reduces to 84mm=s, and v y reduces to 69mm=s. So Ca 2z waves are easier to occur at the boundaries of cardiac myocytes, whcih can be compared with the experimental results [17], [18], [19]. Initiation of Free Ca 2z waves [18] and spontaneous Ca 2z waves [19] is kinetically favored near the boundaries, and the waves initialing from the boundaries are also easier to propagate. In ref. [17], though the results are obtained by line-scan, initiation of the waves is always near the endpoints of the line, and the waves are always triggered near the boundaries of cardiac myocytes. However, the amplitude is smaller than that of the wave from the corner though the change is not obvious. It is because the initial condition is an opening of only one spark, the reflecting effects of the boundaries are not sufficiently obvious. So the reflecting boundaries can increase the propagating probability of Ca 2z waves, though they are not the crucial factor for the propagation of Ca 2z waves. In contrast, the anomalous subdiffusion mode of Ca 2z concentration is the decisive factor for whether the Ca 2z wave can be formed by a single Ca 2z spark. Figure 4 presents the event of two Ca 2z waves collision(3D images). The initial condition is 10ms opening of CRUs at points (2,9.6) and (18,9.6), and they will form two Ca 2z waves. Two Ca 2z waves will meet at the middle as shown in Image 4, at t~120ms. Several sparks fire at the same time, and local Ca 2z concentration reaches a peak value. With increasing time, Ca 2z concentration will return to a lower value under the effect of buffers and pump, and no sparks will occur because of the CRUs' ''refractory period''. The 2D linescan image shows the process of two Ca 2z waves collision and vanishment (Fig. 8 in [17]), but it cannot show how the propagating direction changes. When Ca 2z concentration reaches a peak value in the middle, CRUs along xaxis(y~9:6mm) are closed, but CRUs along y-axis(x~10mm) have never been opened before. Therefore, Ca 2z waves can propagate along the line of x~10mm. Finally, all CRUs are closed and will not reopen.

Modeling Ca 2z Waves Under a Physiological Current
To reproduce the feature of calcium waves found in experiments(primary result is wave velocity), a large current through CRU has been used in the former subsection. However, physiological value of I CRU is about 2pA. So in the following discussion, I CRU~2 pA is adopted to study how many adjacent normal Ca 2z sparks can trigger a Ca 2z wave and find out the longitudinal interval of CRUs which could make a single Ca 2z  spark trigger a normal Ca 2z wave. Because of the small value of I CRU , wave velocity and amplitude will be smaller. In order to make Ca 2z waves spread all over the region, the computation time is prolonged to 500ms. To diminish the effect of the reflecting boundaries, the initial location is chosen at the middle of the region. Figure 5a shows whether a Ca 2z wave can be triggered by one spark at the middle of the region for l x~2 :0mm. When t~490ms, the wave only spreads through half the region, and the events of Ca 2z sparks are almost isolate. Under physiological conditions, even considering anomalous subdiffusion, high value of Ca 2z concentration may stay in a larger region around the firing CRU. For a small current, the amplitude of one spark is still small, the Ca 2z wave cannot propagate to the whole region, so the cardiac myocyte is stable when a normal spontaneous spark occurs.
The initial number of firing sparks is changed to study how many adjacent normal Ca 2z sparks can trigger a Ca 2z wave. The result is shown in Fig. 5b. It can be seen that four CRUs firing simultaneously at the middle will form a ''weak'' Ca 2z wave in the region. At t~490ms, the wave reaches the top, bottom and right boundaries, and several sparks can be found at the same time. However, three adjacent normal Ca 2z sparks can only form a local Ca 2z wave. So with the computation time of 500ms, for b~2:25, I CRU~2 :0pA and l x~2 :0mm, four adjacent CRUs firing together is the critical initial condition to trigger a Ca 2z wave. However, wave velocity and amplitude here is very small(v x &20mm=s), and the Ca 2z concentration of the whole region is much smaller than that in Figs. 2, 3, 4.
In Figure 6a, the longitudinal interval is changed. For the case of l x~1 :0mm, it takes only 110ms for the wave to reach the left boundary. With the simultaneous firing of several CRUs, an obvious Ca 2z concentration wave is observed. Although the amplitude is smaller, longitudinal wave velocity(v x~9 4mm=s) is comparative with the case of I CRU~5 pA, l x~2 :0mm. The physical reason for such a significant change which happens by changing l x~2 :0mm to l x~1 :0mm is that FWHM for I CRU~2 :0pA is about 2:0mm, and if the interval between two CRUs reduces to 1:0mm, the half maximum value of a spark can ''reach'' adjacent CRUs easily. In addition, less interval makes more CRUs fire together, and the wave will be easier to propagate.
Our results have revealed that two factors(l x and number of firing CRUs) can both make a Ca 2z wave propagate. But which the effect is more significant? Figure 6 shows when l x~1 :0mm, three Ca 2z waves trigger from one, four, and nine initial adjacent Ca 2z sparks, respectively. From 6a to 6c, both longitudinal wave velocity and amplitude become larger (v x~9 4,100,113mm=s, amplitudes are 51,65,88mM), but the difference is not obvious as that between Figs. 5a and 6a. So the longitudinal interval of CRUs affects Ca 2z waves more significantly than the number of firing CRUs. If the longitudinal interval of CRUs becomes smaller due to some reasons, such as cardiac myocytes deformation, Ca 2z waves will easily occur, then cardiac myocytes will be unstable.

Conclusion
In this work, we present a mathematical model based on anomalous subdiffusion of Ca 2z concentration in the process of Ca 2z wave triggered by Ca 2z sparks. Ca 2z waves propagating from an initial firing of one single CRU at a corner or in the middle of a 2D rectangular region is numerically simulated. Our results can reproduce wave velocities of experimental results using a small current. We show that Ca 2z waves can be triggered by one single Ca 2z spark under a small CRU current(5pA). When anomalous subdiffusion order b becomes bigger, velocities of Ca 2z waves increase obviously, but the variation of amplitude is not very considerable. The phenomenon of calcium waves collision is also simulated. Under physiological Ca 2z currents(2pA) through CRUs, an initial firing of four adjacent CRUs is  shown to form a Ca 2z wave. When l x~2 :0mm, an isolated spark cannot trigger a Ca 2z wave, so the system is stable under physiological condition. Then the longitudinal interval of CRUs is changed to study how the system becomes unstable and how an obvious Ca 2z wave is formed. Our work is based on a more realistic diffusion model of Ca 2z sparks with the parameters close to physiological values. The simulation results may be useful in further studies about Ca 2z waves.  (10,9.6), snapshots was taken at 10, 60 and 110ms. (b)Snapshots of Ca 2z waves initiating from a 10ms opening of the CRUs at the point (9,9.6), (10,9.6), (9,10.4) and (10,10.4), snapshots are taken at 10, 50 and 90ms. (c)Snapshots of Ca 2z waves initiating from 10ms opening of the CRUs at the point (9,9.6), (10,9.6), (11,9.6), (9,10.4), (10,10.4), (11,10.4), (9,11.2), (10,11.2), (11,11.2),snapshots are taken at 10, 40 and 70ms. The value of parameters for a, b and c are b~2:25, I CRU~2 pA and l x~1 :0mm. doi:10.1371/journal.pone.0057093.g006