Abstract

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.

1. Introduction

The genesis of cardiac arrhythmias in the whole heart scale has been linked to dynamical instabilities at the cellular level [16]. 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 [68]. 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, 716]. 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 beat-to-beat dynamics of cardiac cells [3, 79, 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 [1724]. 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 period-doubling 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.

2. 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 and denote, respectively, the APD and peak cytoplasmic calcium for the th beat, function represents the restitution properties of APD and the Ca-to-APD effect on membrane potential, and function 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 as follows: where with the pacing period, is the Ca-to-APD coupling parameter, and in which , and are positive parameters chosen to simulate experimentally observed APD restitution curve.

However, was not explicitly given in [8], we now deduce an explicit expression of , 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 . 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 is exactly the same as in [8], and is expressed as therein, 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 ; that is, . 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 .

3. 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.

3.1. Dynamics of Calcium Cycling

Let the APD-to-Ca coupling coefficient ; 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 . Let be the fixed point of ; then which is equivalent to Notice that the functions 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 ; is unstable if ; and the fixed point is nonhyperbolic if . In particular, period-doubling bifurcation can take place if . We will consider the bifurcation of the fixed point by taking as a bifurcation parameter.

In the case where , one can derive the following equation from the expression (16): Equivalently, at , we have where We will further assume that Under the above assumption, (18) has two positive solutions where and are the expressions with respect to the parameters , and of the map .

Therefore, at the fixed points , . Plugging the formula of into (15), we have

For fixed parameters , and of  , we let be the values of that satisfy (22) (at , resp.). Moreover, if the following nondegeneracy conditions are satisfied (i); (ii), Then, by applying period-doubling bifurcation theorem in [28] to the map , we have the following theorem.

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 (, resp.). Moreover, the period-two 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 from (22). For given values of , and , we can give numerical computation to show that there exists the value of 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 () (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 period-doubling bifurcation takes place at the left fixed point , and as increases, an inverse period-doubling bifurcation takes place at the right fixed point . 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 . 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.

3.2. 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 , and then it satisfies Linearizing the unperturbed system (10) at , we obtain the Jacobian matrix as follows: where

Let and be the eigenvalues of , and then where Note that can be complex numbers, and for simplicity, let denote the eigenvalue with smaller modulus; hence, . 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, ; is repelling if both of the eigenvalues are greater than one in modulus; that is, ; is called a saddle provided that . 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.

We first pay our attention to the occurrence of period-doubling bifurcations for the fixed point , which corresponds to alternans of the coupled system (10). When we obtain that the fixed point has two eigenvalues and . 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 . Moreover, is invariant under the coupled system (10) when ; that is, if . One can derive a one-dimensional map of the couple system (10) with 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 .

Numerical simulations for alternans of the coupled system are shown in Figure 2. Figure 2(a) gives an example of discordant alternans in APD and , where the corresponding Jacobi matrix at the unstable fixed point is , with eigenvalues and , and the eigenvector for is , and the slope of is negative. While Figure 2(b) gives an example of concordant alternans in APD and , where the corresponding Jacobi matrix is , with eigenvalues and , and the eigenvector for is , and the slope of is positive. The numerical computations illustrate the above theoretical analysis, and they are performed by the software MATLAB 7.0.

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 , from which one can derive that , the eigenvalues are a pair of conjugate complex numbers Thus, the fixed point is locally stable if , and it is unstable if . The fixed point is nonhyperbolic if . 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) by . Rewrite as Then, the nondegeneracy conditions for Hopf bifurcation are (iH); (iiH), for .

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) and 3(b), where the total calcium is chosen as the bifurcation parameter, with other parameters fixed at given values. In Figure 3(a), a stable fixed point exists for the unperturbed system, where the fixed point , the Jacobi matrix , , and . While in Figure 3(b), as decreases, the fixed point becomes unstable, and an attractive closed invariant curve bifurcates from the fixed point, where , , , and . Such computations illustrate well the above theoretical analysis.

Under noise perturbations, that is, , we will pay our attention to the effect of the noise on the closed invariant curve. Let denote the invariant curve, and write for the omega-limit set of under all possible random iterations, where and denotes the corresponding map for system (10). Note that 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 (). While in Figure 3(d) when the noise intensity is enhanced (), 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.

4. 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).

It is well known that power spectra reflect the contribution of frequencies to the variance of , according to Parseval's energy theorem [33]. Let , and let be a time series of length , and then it has a well-defined discrete Fourier transform, , where , and is referred to as Fourier frequency. The power spectrum, defined as , is always used to study the frequency composition of time series [34].

We first study the effect of the noise on alternans by applying power spectra analysis to time series of subsystem (11). Let be a time series obtained according to subsystem (11); that is, , where is a realization of the uniform noise in . Then the power spectrum for 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 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 , and are given in expression (25).

Motivated by the work in [34], we obtain the following theorem.

Theorem 5. Let denote the power spectra of and in the linearized system (37), respectively; then, one can relate the power spectrum of the noise to as follows: Moreover, the quotient is an increasing (decreasing) function with respect to in the interval if and only if ().

Proof. Denote for simplicity, by applying discrete Fourier transform to both sides of (37); we obtain where is the Fourier coefficient with respect to time series and is the Fourier coefficient for the uniform noise. Solving (39) with matrix algebra, we have where is a unit matrix and . Thus, the power spectrum of is with the power spectrum of the noise .
Since , then it is easy to compute that Hence, where
Moreover, the quotient between the two power spectra is which is an increasing (decreasing) function with respect to in the interval if and only if (or ). 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 , and its expected power spectrum .

Remark 7. The quotient provides a quantitative characteristic to estimate the comparison between the power spectra of and . For instance, if , then the quotient between the power spectra of 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 above 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, with . Hence, the quotient between the power spectra of 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 ( 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.

5. 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 calcium-channel noise on dynamical behavior of EC coupling in paced cardiac cells. We first analyzed the subsystem with respect to peak calcium , where period-doubling bifurcation and inverse period-doubling can occur. For the coupled system describing the interplay between APD and peak calcium , the three distinct modes of instability, concordant alternans and discordant alternans, as well as quasiperiodic oscillations of APD and , can be found in the present model. Power spectrum analysis was performed to quantitatively estimate the effects of noise on the alternans and quasi-period oscillations. In fact, such frequency analysis has been widely applied to study the variability in cardiology [3638]. By calculating power spectra for time series obtained in the model, we conclude that the existence of high-frequency peak can be used to measure the alternans, while the existence of low-frequency can be applied to estimate the quasiperiodicity of APD and (see Figure 7). Since time series of APDs and in this study are obtained numerically from our model and not from real experimental data, these results here provide a theoretical viewpoint on the effects of calcium-channel noise on complex dynamics in cardiac cells. This may have important implications in future study for cardiac diseases.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors are very grateful to the anonymous referees for careful reading and helpful suggestions which led to an improvement of our original paper. This research is supported by the NSFC Grants (no. 10925102) and the Program of Shanghai Subject Chief Scientists (no. 10XD1406200).