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

Authors: Jiying Ma; Dongmei Xiao
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

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


---

## Body

## 1. Introduction

The genesis of cardiac arrhythmias in the whole heart scale has been linked to dynamical instabilities at the cellular level [1–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–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–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 beat-to-beat dynamics of cardiac cells [3, 7–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–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 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], APDn+1=F1(APDn,cn+1p), cn+1p=F2(APDn,cnp), where APDn and cnp denote, respectively, the APD and peak cytoplasmic calcium for the nth beat, function F1 represents the restitution properties of APD and the Ca-to-APD effect on membrane potential, and function F2 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 F1 as follows:
(1)F1(APDn,cn+1p)=f(Dn)1-γcn+1p,
where Dn=T-APDn with T the pacing period, γ is the Ca-to-APD coupling parameter, and
(2)f(Dn)=A0[1-1(1+e(Dn-D0)/τ0)]+A1e-(Dn-D1)2/τ1,
in which A0,D0,τ0,A1,D1, and τ1 are positive parameters chosen to simulate experimentally observed APD restitution curve.However,F2 was not explicitly given in [8], we now deduce an explicit expression of F2, according to a detailed intracellular calcium cycling model in [7]. Let cnj and cni denote, respectively, the calcium of the SR and the calcium in the cytoplasm, at the beginning of the nth beat. Hence, the total calcium in the cell, denoted by ct, is the sum of calcium transient in the SR and cytoplasm; that is,
(3)ct=cniα+cnj,
where α is a tissue-dependent positive parameter and ct is assumed to be a constant from beat to beat [7, 8]. Let Rn be the total calcium released from the SR to cytoplasm during the nth beat, and then we have
(4)cnp=cni+αRn.
Shiferaw et al. [7] proposed a one-dimensional map with respect to cj as follows:
(5)cn+1j=wct+(1-w)(cnj-Rn),
where w=w(T) is a proportionality parameter depending on the pacing period, with 0<w<1. Substituting (3) and (4) to (5), we obtain
(6)cn+1j=wct+(1-w)(ct-cniα-Rn)=ct-(1-w)cnpα.
Combining (6) with (3), we have
(7)cn+1i=αct-αcn+1j=(1-w)cnp.
Finally, we derive the iterated map with respect to cp as follows:
(8)cn+1p=cn+1i+αRn+1=(1-w)cnp+αRn+1.
The notation cn+1j is exactly the same as ln in [8], and Rn+1 is expressed as Rn+1=q(Dn)g(cn+1j) therein,
(9)q(Dn)=1-σe-Dn/τq,g(x)=x[1-(1-αg)(1+e(x-lc)/β)],
with αg,β, and lc chosen to fit experimental data. Here, q(Dn) 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:
(10)APDn+1=F1(APDn,cn+1p)=f(T-APDn)(1-γcn+1p),cn+1p=F2(APDn,cnp)=(1-w)cnp+αq(T-APDn)g(ct-(1-w)cnpα)+ϵξn.
The term ϵξn denotes the random perturbations due to calcium-channel noise [23, 24]. Notice that cp≤αct under the relationships (3) and (4). Hence, it is natural to assume that ξn is bounded. Moreover, to understand the effects of gating noise from massive calcium channels, we further assume that ξn is uniformly distributed in the interval [-1,1]; that is, ξn~U(-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 cnp.


## 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σ=0; we have a nonlinear iterated map only with respect to peak cytoplasmic calcium
(11)cn+1p=G(cnp)+ϵξn:=G~(cnp,ϵξn),
where
(12)G(cnp)=(1-w)cnp+αg(ct-(1-w)cnpα)=(1-w)cnp+α(ct-1-wαcnp)×[1-(1-αg)(1+e(ct-((1-w)/α)cnp-lc)/β)].We first analyze the bifurcation behavior of the unperturbed system (11); that is, set ϵ=0. Let c* be the fixed point of G; then(13)(1-w)c*+α(ct-1-wαc*)×[1-1-αg1+e(ct-((1-w)/α)c*-lc)/β]=c*,
which is equivalent to
(14)e(ct-((1-w)/α)c*-lc)/β=(w+αg-wαg)c*-ααgctαct-c*.
Notice that the functions e(ct-((1-w)/α)c*-lc)/β and ((w+αg-wαg)c*-ααgct)/(αct-c*) are both monotonic with respect to c* in the first quadrant, and it is easy to verify that there exists a unique solution for (14). Further, the slope of the map G at c* is
(15)G′(c*)=1-w-(1-w)[1-1-αg1+e(ct-((1-w)/α)c*-lc)/β]-(1-w)(1-αg)αβ(αct-(1-w)c*)×e(ct-((1-w)/α)c*-lc)/β[1+e(ct-((1-w)/α)c*-lc)/β]2.
Plugging (14) into the above equation, we obtain
(16)G′(c*)=1-w-w(1-w)c*αct-(1-w)c*-1-wαβ(1-αg)×(αct-c*)[(w+αg-wαg)c*-ααgct]αct-(1-w)c*.According to the local stability theory in [27], the fixed point c* is locally asymptotically stable if |G′(c*)|<1; c* is unstable if |G′(c*)|>1; and the fixed point is nonhyperbolic if |G′(c*)|=1. In particular, period-doubling bifurcation can take place if G′(c*)=-1. We will consider the bifurcation of the fixed point c* by taking ct as a bifurcation parameter.In the case whereG′(c*)=-1, one can derive the following equation from the expression (16):
(17)1-w-w(1-w)c*αct-(1-w)c*-1-wαβ(1-αg)×(αct-c*)[(w+αg-wαg)c*-ααgct]αct-(1-w)c*=-1.
Equivalently, at G′(c*)=-1, we have
(18)C2c*2+C1c*+C0=0,
where
(19)C2=(1-w)(w+αg-wαg)αβ(1-αg),C1=-2(1-w)-(1-w)(w+2αg-wαg)αctαβ(1-αg),C0=(2-w)αct+(1-w)αgα2ct2αβ(1-αg).
We will further assume that
(20)C2>0,C1<0,C0>0,C12-4C0C1>0.
Under the above assumption, (18) has two positive solutions
(21)c1*=-C1-C12-4C0C12C2,c2*=-C1+C12-4C0C12C2,
where c1*<c2* and ci*(i=1,2) are the expressions with respect to the parameters w,α,ct,αg,β, and lc of the map G.Therefore, at the fixed pointsci*(i=1,2), G′(ci*)=-1. Plugging the formula of ci* into (15), we have
(22)1-w-(1-w)[1-1-αg1+e(ct-((1-w)/α)ci*-lc)/β]-(1-w)(1-αg)(αct-(1-w)ci*)e(ct-((1-w)/α)ci*-lc)/βαβ[1+e(ct-((1-w)/α)ci*-lc)/β]2=-1.For fixed parametersw,α,αg,β, and lc of  G, we let c~ti be the values of ct that satisfy (22) (at ci*,i=1,2, resp.). Moreover, if the following nondegeneracy conditions are satisfied  (i)
[(1/2)(G′′(ci*))2+(1/3)G′′′(ci*)]∣ct=c~ti≠0;(ii)
(∂2G(ci*,c~ti)/∂cnp∂ct)≠0, Then, by applying period-doubling bifurcation theorem in [28] to the map G, we have the following theorem.Theorem 1.
Under assumption (20) and the above nondegeneracy conditions, the map G undergoes a period-doubling bifurcation at the fixed point ci* when parameter ct varies in a small neighborhood of c~ti (i=1,2, resp.). Moreover, the period-two orbit bifurcated from ci* is locally asymptotically stable if the sign in condition (i) is positive.Remark 2.
(i) It is difficult to give an the explicit expression ofc~ti(i=1,2) from (22). For given values of w,α,αg,β, and lc, we can give numerical computation to show that there exists the value of c~ti such that the period-doubling bifurcation occurs.
(ii) Due to the uniqueness of fixed points of mapG, 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 ci*(i=1,2) (see [29]).The bifurcation diagrams of the mapG are shown in Figures 1(a) and 1(b), and in (a), the total calcium ct is chosen as the bifurcation parameter with the other parameters fixed, while in (b), the uptake proportionality coefficient w 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 c1*, and as ct increases, an inverse period-doubling bifurcation takes place at the right fixed point c2*. 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 G can undergo further period-doubling bifurcations for periodic orbits as w varies.Bifurcation diagrams of system (11) with respect to peak calcium. (a) Peak calcium cp versus total calcium ct,ϵ=0, with g(x)=x[1-(1-0.05)/(1+e(x-0.9)/0.06)], w=0.6, and α=2. (b) Peak calcium cp versus the proportionality parameter w,ϵ=0, with the same α and g in (a), and ct=2.1. (c) ϵ=0.2, under the other conditions same as in (a). (d) ϵ=0.05, under the other conditions same as in (b).
(a)
(b)
(c)
(d)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 c*. More precisely, let Ws(c*) denote the basin of attraction for the fixed point c*, and write
(23)Wϵ=⋂n≥0‍⋃i≥n‍G~i(c*,ΔϵN)¯
for the omega-limit set of c* under all possible random iterations, where random perturbations take values in Δϵ=[-ϵ,ϵ]. Note that Wϵ is an invariant set for system (11) [25]. If the noise intensity ϵ is sufficiently small, then Wϵ⊂Ws(c*), and Wϵ 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(a*,c*) be the fixed point of system (10) when ϵ=0, and then it satisfies
(24)f(T-a*)(1-γc*)=a*,(1-w)c*+αq(T-a*)g(ct-(1-w)c*α)=c*.
Linearizing the unperturbed system (10) at (a*,c*), we obtain the Jacobian matrix as follows:
(25)J(a*,c*)=(-γαq′gf(1-γc*)2-f′1-γc*γ(1-w)(1-qg′)f(1-γc*)2-αq′g(1-w)(1-qg′)):=(J11J12J21J22),
where
(26)f:=f(T-a*)=A0[1-1(1+e(T-a*-D0)/τ0)]+A1e-(T-a*-D1)2/τ1,f′:=f′(T-a*)=A0e(T-a*-D0)/τ0τ0[1+e(T-a*-D0)/τ0]2-2A1(T-a*-D1)τ1e-(T-a*-D1)2/τ1,q:=q(T-a*)=1-σe-(T-a*)/τq,q′:=q′(T-a*)=στqe-(T-a*)/τq,g:=g(ct-1-wαc*)=(ct-1-wαc*)[1-1-αg1+e(ct-((1-w)/α)c*-lc)/β],g′:=g′(ct-1-wαc*)=1-1-αg1+e(ct-((1-w)/α)c*-lc)/β-(1-αg)(ct-((1-w)/α)c*)e(ct-((1-w)/α)c*-lc)/ββ[1+e(ct-((1-w)/α)c*-lc)/β]2.Letλ1 and λ2 be the eigenvalues of J(a*,c*), and then
(27)λ1,2=trJ±(trJ)2-4detJ2,
where
(28)trJ≔J11+J22=(1-w)(1-qg′)-γαq′gf(1-γc*)2-f′1-γc*,detJ:=J11J22-J12J21=-(1-w)f′(1-qg′)1-γc*.
Note that λ1,2 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 (a*,c*) is locally asymptotically stable if both of the eigenvalues are less than one in modulus; that is, |λ1,2|<1; (a*,c*) is repelling if both of the eigenvalues are greater than one in modulus; that is, |λ1,2|>1; (a*,c*) is called a saddle provided that |λ1|<1<|λ2|. It is clear that (a*,c*) is unstable in the two latter cases. Further, the fixed point (a*,c*) 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(a*,c*), which corresponds to alternans of the coupled system (10). When
(29)(trJ)2-4detJ>0,-2<trJ<0,trJ+detJ=-1,
we obtain that the fixed point (a*,c*) 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 (a*,c*) from the map restricted to the center manifold [27, 30],
(30)Wc(a*,c*)={(x,y)∈R2∣y=h(x),h(a*)=c*}.
In the neighborhood of (a*,c*), Wc(a*,c*) is a smooth curve on the plane which passes the point (a*,c*) and is tangent to the eigenvector for the eigenvalue λ1=-1. Moreover, Wc(a*,c*) is invariant under the coupled system (10) when ϵ=0; that is, cn+1p=h(APDn+1) if (APDn,cnp)∈Wc(a*,c*). One can derive a one-dimensional map of the couple system (10) with ϵ=0 on the center manifold Wc(a*,c*).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 Wc(a*,c*) hold. Then, the unperturbed system (10) undergoes a period-doubling bifurcation at (a*,c*).Numerical simulations for alternans of the coupled system are shown in Figure2. Figure 2(a) gives an example of discordant alternans in APD and cp, where the corresponding Jacobi matrix at the unstable fixed point is J(a*,c*)=[-0.2411,0.0085;-0.0496,-1.6841], with eigenvalues λ1=-0.2413 and λ2=-1.6838, and the eigenvector for λ2 is v=(-0.0059,1), and the slope of v is negative. While Figure 2(b) gives an example of concordant alternans in APD and cp, where the corresponding Jacobi matrix is J(a*,c*)=[-0.5310,-0.0107;-0.2108,-1.6754], with eigenvalues λ1=-0.5290 and λ2=-1.6774, and the eigenvector for λ2 is v=(0.0094,1), and the slope of v is positive. The numerical computations illustrate the above theoretical analysis, and they are performed by the software MATLAB 7.0.Plots of peak calciumcp versus APD for the coupled system (10), with the same function g in Figure 1, and f(Dn)=0.2[1-1/(1+e(Dn-0.02)/0.04)]+0.1e-(Dn-0.06)2/0.8, and  q(Dn)=1-0.5e-Dn/0.02 at the parameters T=0.4,α=2,w=0.4, and ct=2. (a) ϵ=0, and γ=-0.02. (b) ϵ=0.5, under the other conditions same as in (a). (c) ϵ=0 and γ=0.02. (d) ϵ=0.5, under the other conditions same as in (c).
(a)
(b)
(c)
(d)From the viewpoint of cardiology, the long (short) APD corresponds to a small (large) calcium transient for discordant alternans in Figure2(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 andcp [10], which has also been observed experimentally [31, 32]. In the case when (trJ)2-4detJ<0, from which one can derive that detJ>0, the eigenvalues λ1,2 are a pair of conjugate complex numbers
(31)λ1,2=trJ2±i4detJ-(trJ)22,|λ1|=|λ2|=detJ=-(1-w)f′(1-qg′)1-γc*.
Thus, the fixed point (a*,c*) is locally stable if detJ<1, and it is unstable if detJ>1. The fixed point is nonhyperbolic if det=1. Therefore, we obtain the Hopf bifurcation condition as follows:
(32)-(1-w)f′(1-qg′)1-γc*=1,(33)(trJ)2-4detJ=[(1-w)(1-qg′)-γαq′gf(1-γc*)2-f′1-γc*]2-4<0.We now give the nondegeneracy conditions for Hopf bifurcation. Choosect as a bifurcation parameter, let the other parameters be fixed, and denote the solution of (32) by c^t. Rewrite λ1,2 as
(34)λ1,2=detJe±iφ(ct),φ(ct)=arccostrJ2detJ.
Then, the nondegeneracy conditions for Hopf bifurcation are  (iH)
(∂detJ/∂ct)|ct=c^t≠0;(iiH)
eikφ(c^t)≠1, for k=1,2,3,4.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 (a*,c*) when the parameter ct varies in the small neighborhood of c^t.Numerical simulations for Hopf bifurcation are shown in Figures3(a) and 3(b), where the total calcium ct 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 (a*,c*)=(0.2203,3.4375), the Jacobi matrix J(a*,c*)=[-1.0100,0.0038;-3.5759,-0.9333], λ1,2=-0.9717±i0.1109, and detJ=0.9564<1. While in  Figure 3(b), as ct decreases, the fixed point becomes unstable, and an attractive closed invariant curve bifurcates from the fixed point, where (a*,c*)=(0.2206,3.2914), J(a*,c*)=[-1.0159,0.0043;-3.4776,-1.0505], λ1,2=-1.0332±i0.1218, and detJ=1.0823>1. Such computations illustrate well the above theoretical analysis.Hopf bifurcation for the coupled system (10), plots of peak calcium cp versus APD, where T=0.27,w=0.478, and the other parameters are the same as that in Figure 2. (a) ϵ=0, and ct=1.9. (b) ϵ=0, and ct=1.85. (c) ϵ=0.02, under the other conditions same as in (b). (d) ϵ=0.1, under the other conditions same as in (b).
(a)
(b)
(c)
(d)Under noise perturbations, that is,ϵ≠0, we will pay our attention to the effect of the noise on the closed invariant curve. Let S0 denote the invariant curve, and write
(35)Wϵ(S0)=⋂n≥0‍⋃i≥n‍Fi(S0,ΔϵN)¯
for the omega-limit set of S0 under all possible random iterations, where Δϵ=[-ϵ,ϵ] and F denotes the corresponding map for system (10). Note that Wϵ(S0) 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 cp weakens as the noise intensity increases, and the corresponding time series fluctuate randomly on the plane.


## 3.1. 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
(11)cn+1p=G(cnp)+ϵξn:=G~(cnp,ϵξn),
where
(12)G(cnp)=(1-w)cnp+αg(ct-(1-w)cnpα)=(1-w)cnp+α(ct-1-wαcnp)×[1-(1-αg)(1+e(ct-((1-w)/α)cnp-lc)/β)].We first analyze the bifurcation behavior of the unperturbed system (11); that is, set ϵ=0. Let c* be the fixed point of G; then(13)(1-w)c*+α(ct-1-wαc*)×[1-1-αg1+e(ct-((1-w)/α)c*-lc)/β]=c*,
which is equivalent to
(14)e(ct-((1-w)/α)c*-lc)/β=(w+αg-wαg)c*-ααgctαct-c*.
Notice that the functions e(ct-((1-w)/α)c*-lc)/β and ((w+αg-wαg)c*-ααgct)/(αct-c*) are both monotonic with respect to c* in the first quadrant, and it is easy to verify that there exists a unique solution for (14). Further, the slope of the map G at c* is
(15)G′(c*)=1-w-(1-w)[1-1-αg1+e(ct-((1-w)/α)c*-lc)/β]-(1-w)(1-αg)αβ(αct-(1-w)c*)×e(ct-((1-w)/α)c*-lc)/β[1+e(ct-((1-w)/α)c*-lc)/β]2.
Plugging (14) into the above equation, we obtain
(16)G′(c*)=1-w-w(1-w)c*αct-(1-w)c*-1-wαβ(1-αg)×(αct-c*)[(w+αg-wαg)c*-ααgct]αct-(1-w)c*.According to the local stability theory in [27], the fixed point c* is locally asymptotically stable if |G′(c*)|<1; c* is unstable if |G′(c*)|>1; and the fixed point is nonhyperbolic if |G′(c*)|=1. In particular, period-doubling bifurcation can take place if G′(c*)=-1. We will consider the bifurcation of the fixed point c* by taking ct as a bifurcation parameter.In the case whereG′(c*)=-1, one can derive the following equation from the expression (16):
(17)1-w-w(1-w)c*αct-(1-w)c*-1-wαβ(1-αg)×(αct-c*)[(w+αg-wαg)c*-ααgct]αct-(1-w)c*=-1.
Equivalently, at G′(c*)=-1, we have
(18)C2c*2+C1c*+C0=0,
where
(19)C2=(1-w)(w+αg-wαg)αβ(1-αg),C1=-2(1-w)-(1-w)(w+2αg-wαg)αctαβ(1-αg),C0=(2-w)αct+(1-w)αgα2ct2αβ(1-αg).
We will further assume that
(20)C2>0,C1<0,C0>0,C12-4C0C1>0.
Under the above assumption, (18) has two positive solutions
(21)c1*=-C1-C12-4C0C12C2,c2*=-C1+C12-4C0C12C2,
where c1*<c2* and ci*(i=1,2) are the expressions with respect to the parameters w,α,ct,αg,β, and lc of the map G.Therefore, at the fixed pointsci*(i=1,2), G′(ci*)=-1. Plugging the formula of ci* into (15), we have
(22)1-w-(1-w)[1-1-αg1+e(ct-((1-w)/α)ci*-lc)/β]-(1-w)(1-αg)(αct-(1-w)ci*)e(ct-((1-w)/α)ci*-lc)/βαβ[1+e(ct-((1-w)/α)ci*-lc)/β]2=-1.For fixed parametersw,α,αg,β, and lc of  G, we let c~ti be the values of ct that satisfy (22) (at ci*,i=1,2, resp.). Moreover, if the following nondegeneracy conditions are satisfied  (i)
[(1/2)(G′′(ci*))2+(1/3)G′′′(ci*)]∣ct=c~ti≠0;(ii)
(∂2G(ci*,c~ti)/∂cnp∂ct)≠0, Then, by applying period-doubling bifurcation theorem in [28] to the map G, we have the following theorem.Theorem 1.
Under assumption (20) and the above nondegeneracy conditions, the map G undergoes a period-doubling bifurcation at the fixed point ci* when parameter ct varies in a small neighborhood of c~ti (i=1,2, resp.). Moreover, the period-two orbit bifurcated from ci* is locally asymptotically stable if the sign in condition (i) is positive.Remark 2.
(i) It is difficult to give an the explicit expression ofc~ti(i=1,2) from (22). For given values of w,α,αg,β, and lc, we can give numerical computation to show that there exists the value of c~ti such that the period-doubling bifurcation occurs.
(ii) Due to the uniqueness of fixed points of mapG, 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 ci*(i=1,2) (see [29]).The bifurcation diagrams of the mapG are shown in Figures 1(a) and 1(b), and in (a), the total calcium ct is chosen as the bifurcation parameter with the other parameters fixed, while in (b), the uptake proportionality coefficient w 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 c1*, and as ct increases, an inverse period-doubling bifurcation takes place at the right fixed point c2*. 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 G can undergo further period-doubling bifurcations for periodic orbits as w varies.Bifurcation diagrams of system (11) with respect to peak calcium. (a) Peak calcium cp versus total calcium ct,ϵ=0, with g(x)=x[1-(1-0.05)/(1+e(x-0.9)/0.06)], w=0.6, and α=2. (b) Peak calcium cp versus the proportionality parameter w,ϵ=0, with the same α and g in (a), and ct=2.1. (c) ϵ=0.2, under the other conditions same as in (a). (d) ϵ=0.05, under the other conditions same as in (b).
(a)
(b)
(c)
(d)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 c*. More precisely, let Ws(c*) denote the basin of attraction for the fixed point c*, and write
(23)Wϵ=⋂n≥0‍⋃i≥n‍G~i(c*,ΔϵN)¯
for the omega-limit set of c* under all possible random iterations, where random perturbations take values in Δϵ=[-ϵ,ϵ]. Note that Wϵ is an invariant set for system (11) [25]. If the noise intensity ϵ is sufficiently small, then Wϵ⊂Ws(c*), and Wϵ 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(a*,c*) be the fixed point of system (10) when ϵ=0, and then it satisfies
(24)f(T-a*)(1-γc*)=a*,(1-w)c*+αq(T-a*)g(ct-(1-w)c*α)=c*.
Linearizing the unperturbed system (10) at (a*,c*), we obtain the Jacobian matrix as follows:
(25)J(a*,c*)=(-γαq′gf(1-γc*)2-f′1-γc*γ(1-w)(1-qg′)f(1-γc*)2-αq′g(1-w)(1-qg′)):=(J11J12J21J22),
where
(26)f:=f(T-a*)=A0[1-1(1+e(T-a*-D0)/τ0)]+A1e-(T-a*-D1)2/τ1,f′:=f′(T-a*)=A0e(T-a*-D0)/τ0τ0[1+e(T-a*-D0)/τ0]2-2A1(T-a*-D1)τ1e-(T-a*-D1)2/τ1,q:=q(T-a*)=1-σe-(T-a*)/τq,q′:=q′(T-a*)=στqe-(T-a*)/τq,g:=g(ct-1-wαc*)=(ct-1-wαc*)[1-1-αg1+e(ct-((1-w)/α)c*-lc)/β],g′:=g′(ct-1-wαc*)=1-1-αg1+e(ct-((1-w)/α)c*-lc)/β-(1-αg)(ct-((1-w)/α)c*)e(ct-((1-w)/α)c*-lc)/ββ[1+e(ct-((1-w)/α)c*-lc)/β]2.Letλ1 and λ2 be the eigenvalues of J(a*,c*), and then
(27)λ1,2=trJ±(trJ)2-4detJ2,
where
(28)trJ≔J11+J22=(1-w)(1-qg′)-γαq′gf(1-γc*)2-f′1-γc*,detJ:=J11J22-J12J21=-(1-w)f′(1-qg′)1-γc*.
Note that λ1,2 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 (a*,c*) is locally asymptotically stable if both of the eigenvalues are less than one in modulus; that is, |λ1,2|<1; (a*,c*) is repelling if both of the eigenvalues are greater than one in modulus; that is, |λ1,2|>1; (a*,c*) is called a saddle provided that |λ1|<1<|λ2|. It is clear that (a*,c*) is unstable in the two latter cases. Further, the fixed point (a*,c*) 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(a*,c*), which corresponds to alternans of the coupled system (10). When
(29)(trJ)2-4detJ>0,-2<trJ<0,trJ+detJ=-1,
we obtain that the fixed point (a*,c*) 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 (a*,c*) from the map restricted to the center manifold [27, 30],
(30)Wc(a*,c*)={(x,y)∈R2∣y=h(x),h(a*)=c*}.
In the neighborhood of (a*,c*), Wc(a*,c*) is a smooth curve on the plane which passes the point (a*,c*) and is tangent to the eigenvector for the eigenvalue λ1=-1. Moreover, Wc(a*,c*) is invariant under the coupled system (10) when ϵ=0; that is, cn+1p=h(APDn+1) if (APDn,cnp)∈Wc(a*,c*). One can derive a one-dimensional map of the couple system (10) with ϵ=0 on the center manifold Wc(a*,c*).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 Wc(a*,c*) hold. Then, the unperturbed system (10) undergoes a period-doubling bifurcation at (a*,c*).Numerical simulations for alternans of the coupled system are shown in Figure2. Figure 2(a) gives an example of discordant alternans in APD and cp, where the corresponding Jacobi matrix at the unstable fixed point is J(a*,c*)=[-0.2411,0.0085;-0.0496,-1.6841], with eigenvalues λ1=-0.2413 and λ2=-1.6838, and the eigenvector for λ2 is v=(-0.0059,1), and the slope of v is negative. While Figure 2(b) gives an example of concordant alternans in APD and cp, where the corresponding Jacobi matrix is J(a*,c*)=[-0.5310,-0.0107;-0.2108,-1.6754], with eigenvalues λ1=-0.5290 and λ2=-1.6774, and the eigenvector for λ2 is v=(0.0094,1), and the slope of v is positive. The numerical computations illustrate the above theoretical analysis, and they are performed by the software MATLAB 7.0.Plots of peak calciumcp versus APD for the coupled system (10), with the same function g in Figure 1, and f(Dn)=0.2[1-1/(1+e(Dn-0.02)/0.04)]+0.1e-(Dn-0.06)2/0.8, and  q(Dn)=1-0.5e-Dn/0.02 at the parameters T=0.4,α=2,w=0.4, and ct=2. (a) ϵ=0, and γ=-0.02. (b) ϵ=0.5, under the other conditions same as in (a). (c) ϵ=0 and γ=0.02. (d) ϵ=0.5, under the other conditions same as in (c).
(a)
(b)
(c)
(d)From the viewpoint of cardiology, the long (short) APD corresponds to a small (large) calcium transient for discordant alternans in Figure2(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 andcp [10], which has also been observed experimentally [31, 32]. In the case when (trJ)2-4detJ<0, from which one can derive that detJ>0, the eigenvalues λ1,2 are a pair of conjugate complex numbers
(31)λ1,2=trJ2±i4detJ-(trJ)22,|λ1|=|λ2|=detJ=-(1-w)f′(1-qg′)1-γc*.
Thus, the fixed point (a*,c*) is locally stable if detJ<1, and it is unstable if detJ>1. The fixed point is nonhyperbolic if det=1. Therefore, we obtain the Hopf bifurcation condition as follows:
(32)-(1-w)f′(1-qg′)1-γc*=1,(33)(trJ)2-4detJ=[(1-w)(1-qg′)-γαq′gf(1-γc*)2-f′1-γc*]2-4<0.We now give the nondegeneracy conditions for Hopf bifurcation. Choosect as a bifurcation parameter, let the other parameters be fixed, and denote the solution of (32) by c^t. Rewrite λ1,2 as
(34)λ1,2=detJe±iφ(ct),φ(ct)=arccostrJ2detJ.
Then, the nondegeneracy conditions for Hopf bifurcation are  (iH)
(∂detJ/∂ct)|ct=c^t≠0;(iiH)
eikφ(c^t)≠1, for k=1,2,3,4.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 (a*,c*) when the parameter ct varies in the small neighborhood of c^t.Numerical simulations for Hopf bifurcation are shown in Figures3(a) and 3(b), where the total calcium ct 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 (a*,c*)=(0.2203,3.4375), the Jacobi matrix J(a*,c*)=[-1.0100,0.0038;-3.5759,-0.9333], λ1,2=-0.9717±i0.1109, and detJ=0.9564<1. While in  Figure 3(b), as ct decreases, the fixed point becomes unstable, and an attractive closed invariant curve bifurcates from the fixed point, where (a*,c*)=(0.2206,3.2914), J(a*,c*)=[-1.0159,0.0043;-3.4776,-1.0505], λ1,2=-1.0332±i0.1218, and detJ=1.0823>1. Such computations illustrate well the above theoretical analysis.Hopf bifurcation for the coupled system (10), plots of peak calcium cp versus APD, where T=0.27,w=0.478, and the other parameters are the same as that in Figure 2. (a) ϵ=0, and ct=1.9. (b) ϵ=0, and ct=1.85. (c) ϵ=0.02, under the other conditions same as in (b). (d) ϵ=0.1, under the other conditions same as in (b).
(a)
(b)
(c)
(d)Under noise perturbations, that is,ϵ≠0, we will pay our attention to the effect of the noise on the closed invariant curve. Let S0 denote the invariant curve, and write
(35)Wϵ(S0)=⋂n≥0‍⋃i≥n‍Fi(S0,ΔϵN)¯
for the omega-limit set of S0 under all possible random iterations, where Δϵ=[-ϵ,ϵ] and F denotes the corresponding map for system (10). Note that Wϵ(S0) 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 cp 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 ofx(t), according to Parseval's energy theorem [33]. Let x(t),t=0,1,2,…, and let L-1 be a time series of length L, and then it has a well-defined discrete Fourier transform, X(fk)=F(x(t))=∑t=0L-1‍e-i2πfktx(t), where fk=k/L,k=0,1,…, and L-1 is referred to as Fourier frequency. The power spectrum, defined as P(fk)=|X(fk)|2/L, 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 {cnp:n=0,1,2,…,L-1} be a time series obtained according to subsystem (11); that is, cn+1p=G(cnp)+ϵξn, where ξn is a realization of the uniform noise in [-1,1]. Then the power spectrum for {cnp:n=0,1,2,…,L-1} is
(36)P(fk)=1L|∑n=0L-1‍e-i2πnk/Lcnp|2,k=0,1,…,L-1.
It is clear that P(fk) depends continuously on the noise intensity ϵ as well as the parameters w,α,ct,αg,β, and lc in G. 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 cp 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.(a) Time series of system (11) for different noise intensities ϵ=0,ϵ=0.01, and ϵ=0.1, respectively, where the total calcium ct=1.9 and the other parameters are the same as Figure 1(a). (b) Power spectra corresponding to the time series in (a).
(a)
(b)(a) Time series of system (11) for different noise intensities ϵ=0,ϵ=0.05, and ϵ=0.5, respectively, where the total calcium ct=1.4 and the other parameters are the same as Figure 1(a). (b) Power spectra corresponding to the time series in (a).
(a)
(b)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 (a*,c*), where the orbits of perturbed system can stay in the neighborhood of (a*,c*) for ϵ sufficiently small. Hence, we use the following linearized version of system (10) to approximate its dynamics [34]
(37)(APDn+1cn+1p)=J(a*,c*)(APDncnp)+e2ηn,(APDn+1cn+1p)(APDncnp)n=0,1,…,L-1,
where ηn=ϵξn,e2=(0,1)T, and J(a*,c*) are given in expression (25).Motivated by the work in [34], we obtain the following theorem.Theorem 5.
LetP1,2(fk) denote the power spectra of APDn and cnp in the linearized system (37), respectively; then, one can relate the power spectrum of the noise Pη(fk) to P1,2(fk) as follows:
(38)P1(fk)=Pη(fk)det(B(fk))2J122,P2(fk)=Pη(fk)det(B(fk))2(1+J112-2J11cos(2πfk)).
Moreover, the quotient P1(fk)/P2(fk) is an increasing (decreasing) function with respect to fk in the interval 0≤fk≤0.5 if and only if J11<0 (J11>0).Proof.
DenoteYn=(APDn,cnp)T for simplicity, by applying discrete Fourier transform to both sides of (37); we obtain
(39)ei2πfkX(fk)=J(a*,c*)X(fk)+e2Z(fk),fk=kL,k=0,1,…,L-1,
where X(fk) is the Fourier coefficient with respect to time series {Yn:n=0,1,…,L-1} and Z(fk) is the Fourier coefficient for the uniform noise. Solving (39) with matrix algebra, we have
(40)X(fk)=(ei2πfkI-J(a*,c*))-1e2Z(fk)=B(fk)-1e2Z(fk),
where I is a unit matrix and B(fk)=ei2πfkI-J(a*,c*). Thus, the power spectrum of Yn is
(41)PY(fk)=1L|X(fk)|2=1L|B(fk)-1e2|2|Z(fk)|2=|B(fk)-1e2|2Pη(fk),
with Pη(fk) the power spectrum of the noise η.
Sincee2=(0,1)T, then it is easy to compute that
(42)B(fk)-1e2=1det(B(fk))(J12ei2πfk-J11).
Hence,
(43)PY(fk)=Pη(fk)det(B(fk))2|(J12ei2πfk-J11)|2:=(P1(fk)P2(fk)),
where
(44)P1(fk)=Pη(fk)det(B(fk))2J122,P2(fk)=Pη(fk)det(B(fk))2|ei2πfk-J11|2=Pη(fk)det(B(fk))2(1+J112-2J11cos(2πfk)).
Moreover, the quotient between the two power spectra is(45)P1(fk)P2(fk)=J1221+J112-2J11cos(2πfk),
which is an increasing (decreasing) function with respect to fk in the interval 0≤fk≤0.5 if and only if J11<0 (or J11>0). The proof is finished.Remark 6.
Since the successive noise termsηn are independent and identically distributed, then the noise is white and has equal parts of all frequencies [34]. Note that ηn~U(-ϵ,ϵ), so the variance of ηn is D(ηn)=ϵ2/3, and its expected power spectrum Pη(fk)=D(η)=ϵ2/3.Remark 7.
The quotientP1(fk)/P2(fk) provides a quantitative characteristic to estimate the comparison between the power spectra of APDn and cnp. For instance, if J11<0, then the quotient between the power spectra of APDn and cnp 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, J(a*,c*)=[-1.0100,0.0038;-3.5759,-0.9333] with J11=-1.0100<0. Hence, the quotient between the power spectra of APDn and cnp 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.Power spectra of the coupled system (10), for the two corresponding time series obtained in Figure 3(a). (a) ϵ=0, (b) ϵ=0.005, (c) ϵ=0.05, and (d) ϵ=0.5. To facilitate comparisons, the spectra are normalized to unit variance.
(a)
(b)
(c)
(d)Power spectra of the coupled system (10) for the two corresponding time series obtained in Figure 3(b). (a) ϵ=0, (b) ϵ=0.02, (c) ϵ=0.1, and (d) ϵ=0.5. To facilitate comparisons, the spectra are normalized to unit variance.
(a)
(b)
(c)
(d)


## 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 cp, where period-doubling bifurcation and inverse period-doubling can occur. For the coupled system describing the interplay between APD and peak calcium cp, the three distinct modes of instability, concordant alternans and discordant alternans, as well as quasiperiodic oscillations of APD and cp, 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 [36–38]. 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 cp (see Figure 7). Since time series of APDs and cp 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.



---

*Source: 687472-2013-10-07.xml*