Effects of Calcium-Channel Noise on Dynamics of Excitation-Contraction Coupling in Paced Cardiac Cells

We study a simple discrete model with the impact of calcium-channel noise on the beat-to-beat dynamics of cardiac cells. The effects of the noise are assessed by bifurcation analysis and power spectrum analysis, respectively. It is shown that this model can undergo period-doubling bifurcation and Hopf bifurcation if there are not random perturbations. Under random perturbations, the period-doubling bifurcations of the model can be observed, and the invariant curve from Hopf bifurcation is perturbed to an annulus on the plane and then becomes a totally disordered and randomly scattered region. By the power spectrum analysis, we find that the existence of high-frequency peak in the power spectra links to the period-doubling orbits, while the existence of low-frequency peak corresponds to quasiperiodic orbit.


Introduction
The genesis of cardiac arrhythmias in the whole heart scale has been linked to dynamical instabilities at the cellular level [1][2][3][4][5][6].Cardiac cells are excitable cells in which the physiological link between electrical excitation and the production of contraction underlies its excitation-contraction (EC) coupling.It is well established that the contraction of cardiac cells is triggered by a considerable increase in intracellular calcium concentration and the interplay between membrane voltage and calcium cycling forms the basis of EC coupling [6][7][8].To study the mechanisms of the coupling, a range of experimental and theoretical studies has been conducted concerning the scale from cells to tissues (see [1,3,6] for a recent review).Most of these works focused on the dynamical instabilities in cardiac cells, such as the alternations and quasiperiodic oscillations in action potential duration (APD) and calcium transient [4,5,[7][8][9][10][11][12][13][14][15][16].Among various mathematical models in the literatures, discrete dynamical systems (or iterated maps) are considered to be a simple but effective approach in understanding the beatto-beat dynamics of cardiac cells [3,[7][8][9]14].A seminal mathematical interpretation of APD alternans proposed by Nolasco and Dahlen [14] was an iterated map, under the basic assumption that the present APD depends on the diastolic interval of the previous beat.Shiferaw et al. [7] developed a physiologically detailed model of calcium cycling and analyzed its nonlinear dynamics by reducing the model to a discrete map.
In particular, nowadays there is an increasing interest to take into account the stochastic factors in heart dynamics [17][18][19][20][21][22][23][24].There are millions of calcium channels distributed in the membrane as well as intracellular space of cardiac cells, leading to random fluctuations in calcium cycling system [21].For instance, Tanskanen et al. [23] took into account the stochastic gating of L-type calcium channels in their model, and Shuai and Jung [24] studied the stochasticity of intracellular calcium release mechanism due to the random opening and closing of calcium channels.In addition, to study the effects of random fluctuations, Sato et al. [17] added Gaussian random variables to APD restitution in their iterated map analysis.Such stochastic factors have been correlated with the beat-to-beat repolarization variability in cardiac cells [20].
In this paper, by combining the models in [7,8], we reduce a discrete model on the interaction between APD and peak cytoplasmic calcium, where stochastic events of calcium channels are incorporated into the model by using uniform random numbers.We first study the nonlinear dynamics of the system by bifurcation analysis.Alternans via perioddoubling bifurcation and quasiperiodic oscillations via Hopf bifurcation can occur for this model.To better characterize the effects of noise on the alternans and quasiperiodic oscillations, we further perform power spectrum analysis.It is interesting to note that the existence of high-frequency peak in the power spectra corresponds to alternations, while the existence of low-frequency peak corresponds to quasiperiodicity.Moreover, as the noise intensity increases, the prominent peaks will gradually disappear, and power spectra tend to distribute more evenly in all frequencies.

Model Formulation
Qu et al. [8] proposed an iterated map with respect to APD and the calcium load in the sarcoplasm reticulum (SR), and some interesting dynamics such as alternans and quasiperiodicity of APD and calcium can be observed for that model.However, note that calcium concentration in the SR is difficult to measure experimentally [7], and calcium alternans is typically referred to the alternations in peak cytoplasmic calcium [11].Thus, in this study, we start with a general two-dimensional iterated map summarized by in [9], , where APD  and    denote, respectively, the APD and peak cytoplasmic calcium for the th beat, function  1 represents the restitution properties of APD and the Ca-to-APD effect on membrane potential, and function  2 describes the mechanisms of calcium cycling and the APD-to-Ca effect on calcium current.Qu et al. [8] came up with a concrete formula of  1 as follows: where   =  − APD  with  the pacing period,  is the Ca-to-APD coupling parameter, and in which  0 ,  0 ,  0 ,  1 ,  1 , and  1 are positive parameters chosen to simulate experimentally observed APD restitution curve.However,  2 was not explicitly given in [8], we now deduce an explicit expression of  2 , according to a detailed intracellular calcium cycling model in [7].Let    and    denote, respectively, the calcium of the SR and the calcium in the cytoplasm, at the beginning of the th beat.Hence, the total calcium in the cell, denoted by   , is the sum of calcium transient in the SR and cytoplasm; that is, where  is a tissue-dependent positive parameter and   is assumed to be a constant from beat to beat [7,8].Let   be the total calcium released from the SR to cytoplasm during the th beat, and then we have Shiferaw et al. [7] proposed a one-dimensional map with respect to   as follows: where  = () is a proportionality parameter depending on the pacing period, with 0 <  < 1. Substituting ( 3) and ( 4) to (5), we obtain Combining ( 6) with (3), we have Finally, we derive the iterated map with respect to   as follows: The notation   +1 is exactly the same as   in [8], and  +1 is expressed as with   , , and   chosen to fit experimental data.Here, (  ) represents the dependence of SR calcium release on the restitution properties, and  is the APD-to-Ca coupling coefficient.
Finally, following the analysis in [17], we add random perturbations to calcium cycling and obtain the following coupled map: The term   denotes the random perturbations due to calcium-channel noise [23,24].Notice that   ≤   under the relationships (3) and ( 4).Hence, it is natural to assume that   is bounded.Moreover, to understand the effects of gating noise from massive calcium channels, we further assume that   is uniformly distributed in the interval [−1, 1]; that is,   ∼ (−1, 1).And uniform random numbers are often adopted to represent random perturbations in stochastic modeling [25,26].Here,  denotes the intensity of the noise.The value of  is determined by the strength of random perturbations [25], and  is supposed to be relatively small compared with the values of the peak cytoplasmic calcium    .

Dynamics of the Model
In this section, we will study the dynamical behavior of system (10) based on bifurcation analysis.The coupled system consists of two subsystems, one is with respect to the mechanism of APD restitution, and the other is about the mechanism of intracellular calcium cycling.The dynamical instability due to the shape of APD restitution curve has been analyzed by Qu et al. [8].And we will focus on the dynamical instability due to intracellular calcium cycling in the following subsection.

Dynamics of Calcium Cycling.
Let the APD-to-Ca coupling coefficient  = 0; we have a nonlinear iterated map only with respect to peak cytoplasmic calcium where We first analyze the bifurcation behavior of the unperturbed system (11); that is, set  = 0. Let  * be the fixed point of ; then which is equivalent to Notice that the functions  (  −((1−)/) * −  )/ and (( +   −   ) * −     )/(  −  * ) are both monotonic with respect to  * in the first quadrant, and it is easy to verify that there exists a unique solution for (14).Further, the slope of the map  at  * is Plugging ( 14) into the above equation, we obtain According to the local stability theory in [27], the fixed point  * is locally asymptotically stable if In particular, period-doubling bifurcation can take place if   ( * ) = −1.We will consider the bifurcation of the fixed point  * by taking   as a bifurcation parameter.
In the case where   ( * ) = −1, one can derive the following equation from the expression (16): Equivalently, at   ( * ) = −1, we have where We will further assume that Under the above assumption, (18) has two positive solutions where  * 1 <  * 2 and  *  ( = 1, 2) are the expressions with respect to the parameters , ,   ,   , , and   of the map .

Theorem 1.
Under assumption (20) and the above nondegeneracy conditions, the map  undergoes a period-doubling bifurcation at the fixed point  *  when parameter   varies in a small neighborhood of c ( = 1, 2, resp.).Moreover, the periodtwo orbit bifurcated from  *  is locally asymptotically stable if the sign in condition (i) is positive.Remark 2. (i) It is difficult to give an the explicit expression of c ( = 1, 2) from (22).For given values of , ,   , , and   , we can give numerical computation to show that there exists the value of c such that the period-doubling bifurcation occurs.
(ii) Due to the uniqueness of fixed points of map , we can see that the period-doubling bifurcation at the other point is inverse bifurcation if a period-doubling bifurcation takes place at one of  *  ( = 1, 2) (see [29]).
The bifurcation diagrams of the map  are shown in Figures 1(a) and 1(b), and in (a), the total calcium   is chosen as the bifurcation parameter with the other parameters fixed, while in (b), the uptake proportionality coefficient  is chosen as the bifurcation parameter with the other parameters fixed.It is easy to check that the parameter values here satisfy assumption (20).From Figure 1(a), one can see that a perioddoubling bifurcation takes place at the left fixed point  * 1 , and as   increases, an inverse period-doubling bifurcation takes place at the right fixed point  * 2 .This numerical result illustrates our theoretical analysis.In addition, the bifurcation diagram of Figure 1(a) coincides well with the plot in [8] (see Figure 5(f) therein), although we used a simplified uptake function in calcium cycling.In Figure 1(b), there are also a period-doubling bifurcation and an inverse period-doubling bifurcation for fixed points, and moreover, the map  can undergo further period-doubling bifurcations for periodic orbits as  varies.
We now turn to system (11) with noise perturbations; that is, set  ̸ = 0.It is known that random diffeomorphisms with bounded absolutely continuous noise possess a finite number of stationary measures, which are supported on invariant sets [25].For the perturbed system (11), since the noise is uniformly distributed, we can see that system (11) with small  has a stationary measure with support near  * .More precisely, let   ( * ) denote the basin of attraction for the fixed point  * , and write for the omega-limit set of  * under all possible random iterations, where random perturbations take values in Note that   is an invariant set for system (11) [25].If the noise intensity  is sufficiently small, then   ⊂   ( * ), and   depends continuously on  with respect to Hausdorff metric [25].Numerical simulations for the perturbed system (11) are given in Figures 1(c) and 1(d), from which one can see that period-doubling bifurcation and its inverse bifurcation still exist in the sense of perturbations, and time series fluctuate in bounded area.

Dynamics of the Coupled System.
In this subsection, we will analyze the dynamics of the coupled system concerning the interactions between membrane voltage and intracellular calcium cycling.Let ( * ,  * ) be the fixed point of system (10) when  = 0, and then it satisfies Linearizing the unperturbed system (10) at ( * ,  * ), we obtain the Jacobian matrix as follows: where Let  1 and  2 be the eigenvalues of ( * ,  * ), and then where Note that  We first pay our attention to the occurrence of perioddoubling bifurcations for the fixed point ( * ,  * ), which corresponds to alternans of the coupled system (10).When (tr ) 2 − 4 det  > 0, −2 < tr  < 0, we obtain that the fixed point ( * ,  * ) has two eigenvalues  1 = −1 and | 2 | < 1.Thus, the fixed point is nonhyperbolic, and the period-doubling bifurcation may take place.By the center manifold theory, we can study the dynamics of the coupled system (10) near the fixed point ( * ,  * ) from the map restricted to the center manifold [27,30], In the neighborhood of ( * ,  * ),   ( * ,  * ) is a smooth curve on the plane which passes the point ( * ,  * ) and is tangent to the eigenvector for the eigenvalue  1 = −1.Moreover,   ( * ,  * ) is invariant under the coupled system (10) when  = 0; that is, One can derive a one-dimensional map of the couple system (10) with  = 0 on the center manifold   ( * ,  * ).
Therefore, applying period-doubling bifurcation theorem in [28] to this derived one-dimensional map, we have the following theorem.Theorem 3. Suppose that assumption (29) and the nondegeneracy conditions of the system restricted on the center manifold   ( * ,  * ) hold.Then, the unperturbed system (10) undergoes a period-doubling bifurcation at ( * ,  * ).
From the viewpoint of cardiology, the long (short) APD corresponds to a small (large) calcium transient for discordant alternans in Figure 2(a), while the long (short) APD corresponds to a large (small) calcium transient for concordant alternans in Figure 2(c).Such concordant and discordant alternans were also observed in the previous studies [8,10].Under random perturbations, as is shown in Figures 2(b) and 2(d), the noise due to calcium channels can prominently affect the action potential properties, leading to fluctuations in APD sequences.
In the following, we will focus on the occurrence of Hopf bifurcation for the coupled map, where the closed invariant curve bifurcated through Hopf corresponds to quasiperiodic oscillations in APD and   [10], which has also been observed experimentally [31,32].In the case when (tr ) 2 − 4 det  < 0, from which one can derive that det  > 0, the eigenvalues  1,2 are a pair of conjugate complex numbers Thus, the fixed point ( * ,  * ) is locally stable if det  < 1, and it is unstable if det  > 1.The fixed point is nonhyperbolic if det = 1.Therefore, we obtain the Hopf bifurcation condition as follows: We now give the nondegeneracy conditions for Hopf bifurcation.Choose   as a bifurcation parameter, let the other parameters be fixed, and denote the solution of (32)   Then, the nondegeneracy conditions for Hopf bifurcation are Applying the Neimark-Sacker bifurcation theorem in [28] to the unperturbed system (10), we obtain the following theorem.
Theorem 4. Suppose that assumption (33) and conditions (iH) and (iiH) hold, then the unperturbed system (10) can undergo a Hopf bifurcation at fixed point ( * ,  * ) when the parameter   varies in the small neighborhood of ĉ .
Numerical simulations for Hopf bifurcation are shown in Figures 3(a Under noise perturbations, that is,  ̸ = 0, we will pay our attention to the effect of the noise on the closed invariant curve.Let  0 denote the invariant curve, and write for the omega-limit set of  0 under all possible random iterations, where Δ  = [−, ] and  denotes the corresponding map for system (10).Note that   ( 0 ) is an invariant set [25].Numerical results with the impact of noise are shown in Figures 3(c) and 3(d).Figure 3(c) gives an invariant annulus on the plane when the noise intensity is sufficiently small ( = 0.02).While in Figure 3(d) when the noise intensity is enhanced ( = 0.1), the invariant annulus disappears, leading to a totally disordered and randomly scattered region on the plane.From Figures 3(c) and 3(d), one can see that the quasiperiodicity of APD and   weakens as the noise intensity increases, and the corresponding time series fluctuate randomly on the plane.

Power Spectrum Analysis
We have analyzed the dynamics of the randomly perturbed system (10) through bifurcation approach.As is shown above, considerable fluctuations always appear in the time series under random perturbations.To better characterize the effects of noise on the alternans and quasiperiodic oscillations, we will further perform power spectrum analysis for system (10).
We first study the effect of the noise on alternans by applying power spectra analysis to time series of subsystem (11).Let {   :  = 0, 1, 2, . . .,  − 1} be a time series obtained according to subsystem (11); that is,  It is clear that (  ) depends continuously on the noise intensity  as well as the parameters , ,   ,   , , and   in .Such continuity on parameters holds true for the power spectra of the coupled system.Hence, in the following numerical simulations, three values of  are chosen to illustrate the influence of different noise intensites.
The case when the unperturbed subsystem (11) has a stable fixed point is shown in Figure 4, and the case when the system has a period-two cycle as the steady state is given in Figure 5. From Figure 5(a), it is clear to see that time series for c  exhibit prominent fluctuations under random perturbations and the amplitudes of fluctuations become larger as the noise intensity increases.The power spectra for the corresponding time series in Figure 5(b) confirm this fact, the stronger the noise is, the higher the power spectra at the same frequencies.And the spectra tend to distribute evenly in all frequencies with the effects of noise.Moreover, comparing (b) in Figures 4 and 5, we can see that there exists a high-frequency peak in the power spectra for the period-two cycle, and such peaks gradually disappear as the noise intensity increases.This phenomenon reveals that the existence of high-frequency peak in the power spectra links to alternations, which could be used to estimate alternans experimentally.
We turn to the effect of the noise on quasiperiodic oscillations for the coupled system (10).We will start from the case when the unperturbed system (10) has a stable fixed point ( * ,  * ), where the orbits of perturbed system can stay in the neighborhood of ( * ,  * ) for  sufficiently small.Hence, we use the following linearized version of system (10) to approximate its dynamics [34] ( where   =   , e 2 = (0, 1)  , and ( * ,  * ) are given in expression (25).Motivated by the work in [34], we obtain the following theorem.
Theorem 5. Let  1,2 (  ) denote the power spectra of APD  and    in the linearized system (37), respectively; then, one can relate the power spectrum of the noise   (  ) to  1,2 (  ) as follows: where (  ) is the Fourier coefficient with respect to time series {  :  = 0, 1, . . .,  − 1} and (  ) is the Fourier coefficient for the uniform noise.Solving (39) with matrix algebra, we have where  is a unit matrix and (  ) =  i2   − ( * ,  * ).Thus, the power spectrum of   is with   (  ) the power spectrum of the noise .
Since e 2 = (0, 1)  , then it is easy to compute that (  ) where Moreover, the quotient between the two power spectra is which is an increasing (decreasing) function with respect to   in the interval 0 ≤   ≤ 0.5 if and only if  11 < 0 (or  11 > 0).The proof is finished.Remark 6.Since the successive noise terms   are independent and identically distributed, then the noise is white and has equal parts of all frequencies [34].Note that   ∼ (−, ), so the variance of   is (  ) =  2 /3, and its expected power spectrum   (  ) = () =  2 /3.Remark 7. The quotient  1 (  )/ 2 (  ) provides a quantitative characteristic to estimate the comparison between the power spectra of APD  and    .For instance, if  11 < 0, then the quotient between the power spectra of APD  and    is greater in high frequencies than that in low frequencies.
The case when the unperturbed system (10) has a stable fixed point is displayed in Figure 6, and the case when there exists an invariant closed curve is given in Figure 7.To facilitate comparisons, here the spectra are normalized to unit variance [34].Firstly, according to the theoretical analysis, in the first case, the power spectra properties of system (10) can be predicted by studying the linearized version.In fact, the value of the Jacobi matrix has been given in Section 3.2; that is, ( * ,  * ) = [−1.0100,0.0038; −3.5759, −0.9333] with  11 = −1.0100< 0. Hence, the quotient between the power spectra of APD  and    increases as the frequency increases, which is obvious in Figures 6(c) and 6(d).Such property about the quotient still holds true in Figures 7(c) and 7(d).Secondly, there are obvious high-frequency peaks for power spectra in both figures, which indicates that the corresponding time series converge to the steady state in oscillations.Finally, it turns out in our simulations that the power spectra in Figure 7(a) have prominent low-frequency peaks, compared with Figure 6(a).With a sufficiently small noise intensity ( = 0.02 in Figure 7(b)), the low-frequency peaks remain in existence but become less prominent.We conclude that the existence of low-frequency peak corresponds to quasiperiodicity.

Discussion
In a clinical setting, alternans in the electrocardiogram morphology is often referred to as a risk of sudden cardiac death [35].In this paper, we study the effects of calciumchannel noise on dynamical behavior of EC coupling in paced

Figure 4 :
Figure 4: (a) Time series of system (11) for different noise intensities  = 0,  = 0.01, and  = 0.1, respectively, where the total calcium   = 1.9 and the other parameters are the same as Figure 1(a).(b) Power spectra corresponding to the time series in (a).

Figure 5 :
Figure 5: (a) Time series of system (11) for different noise intensities  = 0,  = 0.05, and  = 0.5, respectively, where the total calcium   = 1.4 and the other parameters are the same as Figure 1(a).(b) Power spectra corresponding to the time series in (a).
[27]can be complex numbers, and for simplicity, let  1 denote the eigenvalue with smaller modulus; hence, | 1 | ≤ | 2 |.According to the stability theory for nonlinear maps in[27], the fixed point ( * ,  * ) is locally asymptotically stable if both of the eigenvalues are less than one in modulus; that is, | 1,2 | < 1; ( * ,  * ) is repelling if both of the eigenvalues are greater than one in modulus; that is, | 1,2 | > 1; ( * ,  * ) is called a saddle provided that | 1 | < 1 < | 2 |.It is clear that ( * ,  * ) is unstable in the two latter cases.Further, the fixed point ( * ,  * ) is nonhyperbolic if there exist eigenvalues with unit modulus, and corresponding bifurcations may occur for the system.