Coupled iterated map models of action potential dynamics in a one-dimensional cable of cardiac cells

Low-dimensional iterated map models have been widely used to study action potential dynamics in isolated cardiac cells. Coupled iterated map models have also been widely used to investigate action potential propagation dynamics in one-dimensional (1D) coupled cardiac cells, however, these models are usually empirical and not carefully validated. In this study, we first developed two coupled iterated map models which are the standard forms of diffusively coupled maps and overcome the limitations of the previous models. We then determined the coupling strength and space constant by quantitatively comparing the 1D action potential duration profile from the coupled cardiac cell model described by differential equations with that of the coupled iterated map models. To further validate the coupled iterated map models, we compared the stability conditions of the spatially uniform state of the coupled iterated maps and those of the 1D ionic model and showed that the coupled iterated map model could well recapitulate the stability conditions, i.e. the spatially uniform state is stable unless the state is chaotic. Finally, we combined conduction into the developed coupled iterated map model to study the effects of coupling strength on wave stabilities and showed that the diffusive coupling between cardiac cells tends to suppress instabilities during reentry in a 1D ring and the onset of discordant alternans in a periodically paced 1D cable.


Introduction
The cardiac action potential is determined by many different types of membrane currents and their couplings with intracellular calcium cycling [1,2]. Mathematical models of action potential are usually high-dimensional and nonlinear, and the newly developed ones are becoming more and more complicated [3]- [12]. This makes it extremely difficult to investigate the mechanisms of action potential dynamics, even in a single isolated cell. On the other hand, low-dimensional iterated maps, which are easy to deal with both analytically and numerically, have been widely used in studies of cardiac action potential excitation and wave propagation dynamics [13]- [34], and have provided great insights for the cardiac action potential dynamics. The earliest map model was developed by Nolasco and Dahlen [13] in their 1968 paper entitled 'A graphic method for the study of alternation in cardiac action potentials'. Subsequent studies by Guevara et al [15] and Chialvo et al [19,20] have used low-dimensional maps to successfully describe alternans and other complex dynamics from their experimental observations, while Vinet et al [18] and Lewis and Guevara [17] have used low-dimensional maps to explain the dynamical mechanisms of their ionic model results. These pioneer studies have led to a widely used map function which is called the action potential duration (APD) restitution curve, defined as the present APD as a function of that preceding diastolic interval (DI), and the following map equation: a n+1 = f (d n ) = f (mT − a n ). (1) In equation (1), a n+1 is the APD at the (n + 1)th beat, d n is the DI at the nth beat, and mT = a n + d n is the stimulation period. m is chosen so that (m − 1)T < APD n + DI min < mT in which DI min is the minimum DI for which an action potential can be elicited by the stimulus. m = 1 corresponds to 1 : 1 capture and m = 2 corresponds to 2 : 1 capture, etc. By iterating equation (1) with a proper APD restitution function f, it has been shown that one is able to obtain very similar bifurcation sequences as observed in experiments [23,35], demonstrating that iterated maps (such as equation (1)) can describe the APD dynamics in real cardiac cells well. Models of coupled iterated maps were used to study reentry dynamics in one-dimensional (1D) rings or cables with the coupling between cells achieved by conduction velocity (CV) but not coupling in APD [21,22,24,33,36,37]. The coupling in APD or repolarization was first introduced by Chialvo [38] and later by Vinet [25] to study reentry dynamics. Same or similar forms were used in other studies [26]- [28], [30], [39]- [40]. The coupling in APD introduced by Vinet [25] for a 1D ring and later adopted by Fox et al [27] for a 1D cable has the following form: where d n (i) = T − a n (i). α is the number of neighbors coupled with cell i, and w j is a normalized weight function which was empirically set as a Gaussian function [25]: where µ is a parameter measuring the spatial scale of coupling, and is referred to as the space constant of coupling in this study. w j is normalized, i.e. α j=−α w j = 1. In a recent study, Echebarria and Karma [41] derived from a simplified action potential model an equation similar to equation (2) for a propagating wave with velocity c in a continuous medium, i.e., where A c is the APD at the bifurcation point leading to APD alternans [41]. If the CV is assumed to be infinite (or equivalently all the cells are activated simultaneously), then equation (4) becomes which is the same as equation (2) but in a continuous system with an infinite spatial dimension. Under certain limited conditions, equation (5) can be approximated by the following coupled map model which was derived by Echebarria and Karma [41] and has also been modified by others to study the dynamics of discordant alternans and conduction block in 1D cables [30,42]. Although several coupled iterated map models [25]- [28], [30], [39]- [42], such as equations (2) and (6), were used to study action potential propagation dynamics in 1D rings and cables, whether they are appropriate for or how accurate they can describe the action potential dynamics remains unclear. There exist several limitations in these models as detailed below: 1. Both equations (2) and (5) require normalization of the coupling strength, i.e. α j=−α w j = 1 or ∞ −∞ w(y) dy = 1. If this normalization is not satisfied, equation (2) or (5) may become invalid. For example, for a homogeneous 1D cable, the coupled map model should degenerate into the single cell map equation (1) when the spatially uniform state 4 is stable in space (i.e. no spatial instability occurs). In other words, if α j=−α w j = ξ = 1 or ∞ −∞ w(y) dy = ξ = 1, the map resulting from equation (2) for a homogeneous state is a n+1 = ξ f (d n ) = ξ f (T − a n ), which is invalid. However, it is not physiologically clear why such a normalization is needed. 2. Even if the normalization is satisfied, the physiological interpretation of these two equations is still problematic. For example, in the case of a 1D cable of diffusively coupled identical cells stimulated uniformly over space, the diffusive currents between cells are zero so that the action potential of one cell is independent of the other cells. But both equations (2) and (5) imply that even in such a homogeneous case, the APD of a cell is the sum of weighted contributions from all the other cells, which is contrary to what happens physiologically. 3. The coupled map model equation (6) was derived under certain limited conditions (e.g. close to the onset of alternans), thus its application to general conditions remains unclear.
One limitation of this model is that the coupling in the original equation (equation (5)) is long range or global while the coupling is short range in equation (6). How this change in coupling affects the dynamics remains unknown in this system. 4. In the 1D cable of ionic model, no-flux boundary condition is usually used, i.e.
where V is the membrane voltage. When the 1D cable is modeled with the coupled iterated maps, it is not clear how this boundary condition should be represented and this issue has not been addressed in previous studies.
In this study, we modified the coupled iterated map model (equation (2)) into a standard diffusively coupled map model for the 1D cable of coupled cardiac cells, similar to a coupled map system called a coupled map lattice (CML) which was first phenomenologically formulated by Kaneko [43,44] and then widely used to study the spatiotemporal dynamics of nonlinear systems [44]- [51]. We also developed another form of a diffusively coupled map model. We developed a method to determine the coupling strength and the space constant. In this method, the coupling strength and the space constant are determined by minimizing the distance between the APD spatial profile from the heterogeneous 1D ionic model and that of the heterogeneous 1D coupled map model. We also used this method to determine the boundary conditions for the 1D coupled map model. To study whether coupled map models can recapitulate the stability conditions of the spatially uniform states, we analytically calculated the spectrum of the maximum Lyapunov exponents of the coupled map model and compared them with the numerically obtained Lyapunov spectra of the ionic model, and showed that they agreed well with each other. Finally, we incorporated conduction into the coupled map model to study how the coupling strength affects the onset of instability for reentry in a 1D ring and discordant alternans in a periodically paced 1D cable.

Ionic model
We used the following partial differential equation to simulate the action potential excitation in the 1D cable with either no-flux or periodic boundaries: where V(mV) is the membrane voltage, C m = 1 µF cm −2 is the membrane capacitance, and D is the diffusion constant which was set to 0.001 cm 2 ms −1 as control, unless otherwise stated. I ion (µA cm −2 ) is the total ionic current density from the Beeler-Reuter (BR) model [4]: is the slow inward current density in which E s = −82.3 − 13.0287 ln[Ca] i . I x 1 =ī x 1 x 1 is the time-dependent outward potassium current density; and I K 1 is the time-independent outward potassium current density, which is a voltage-dependent function only. The intracellular calcium concentration [Ca] i (mM) is described by m, h, j, d, f and x 1 are the corresponding gating variables satisfying differential equations of the type: where y ∞ = α y /(α y + β y ) and τ y = 1/(α y + β y ). α y (ms −1 ) and β y (ms −1 ) are the corresponding rate constants which are also functions solely of voltage. Detailed forms of I K 1 ,ī x 1 , α y and β y can be found in the original paper by Beeler and Reuter [4]. I sti is the stimulation current density which was applied uniformly in the entire cable with a current strength 30 µA cm −2 and duration 2 ms.

Changing APD and APD heterogeneity in the cable
APD is defined as the time duration of V > −72 mV of the action potential. To change APD of the BR model, we modified the slow inward current equation to I si =Ḡ si d f (V − E s ) and the time-dependent potassium current to I x 1 =Ḡ x 1ī x 1 x 1 . In control,Ḡ si = 0.09 andḠ x 1 = 1. We increase (decrease)Ḡ si or decrease (increase)Ḡ x 1 to increase (decrease) APD. APD heterogeneities in the cable were generated by changingḠ x1 in space as where x 0 is the center of the cable.

Numerical methods
We used the forward Euler scheme to integrate equation (7), i.e.
. (11) x = 0.015 cm and t = 0.02 ms were used in this study.

Calculation of Lyapunov exponent spectrum in the homogeneous cable
The calculation of Lyapunov exponent was carried out using a method similar to that in the previous study [52]. Assuming a uniform state [V(t), y(t)], a perturbation is then given to this state, i.e. V (x, t) = V (t) + δV (t)e ikx and y(x, t) = y(t) + δy(t)e ikx . Note that y represents the gating variables as described above. Inserting this perturbation into equations (7)-(9), one obtains the linearized equation as: is the Jacobian of the system equations (7)-(9). The maximum Lyapunov exponent λ is calculated as:

Coupled iterated map model I
Rewrite α j=−α w j = 1 into the following form: then equation (2) becomes: Equation (15) is a standard form of diffusively coupled iterated map system, similar to the widely used CML [43]- [51]. The difference between equation (15) and the generally used CML is that the coupling constant in CML is arbitrary, while in equation (15) the coupling strength w j is not arbitrary but normalized. To overcome this limitation, we modified equation (15) to include an arbitrary coupling strength ε, i.e.
and we use the standard Gaussian distribution for w j , i.e.
Equation (16) is now a standard CML with α neighbors of coupling. σ is the space constant which is simply related to µ in equation (3) as σ 2 = 1/2µ. Although equation (16) is a modification of equation (2) or (5), the physiological interpretation is clear: the action potential of one cell is affected by the difference between its own action potential and those of other cells, and normalization of the coupling strength is no longer required.

7
Since y in equation (5) can be represented by y = jdy = j x in which x (which was set to 0.015 cm in this paper) is assumed to be the length of a cardiac cell, then w(y) in equation (5) becomes Comparing this with equation (17), one obtains, Therefore, equation (18) is the theoretically predicted relation between σ 2 and APD and diffusion constant D based on Echebarria and Karma [41]. In the following sections, we will use numerical simulations to determine the coupling strength ε and space constant σ for two types of boundary conditions-periodic and no-flux boundary conditions.
We chose the cable length as an odd number in this case, i.e. L = (2M + 1) with M being a positive integer, so that α = (L − 1)/2 = M is an integer.

Numerical estimation of the coupling strength and space constant.
To estimate the coupling strength ε and space constant σ in equation (16), we compared the results from the ionic model (equation (7)) with those of the map model (equation (16)). The simulations were carried out as follows (the pacing interval T = 1000 ms): (1) Stimulate every cell simultaneously in a 1D heterogeneous ionic model with no cell coupling (uncoupled array in figure 1(A), or equivalently set D = 0 in equation (7)) to determine an APD profile a 0 = [a 0 (1), a 0 (2), . . . , a 0 (L)] in the 1D space (open squares in figure 1(B)). L is the number of cells. The heterogeneity was generated by changingḠ x 1 according to equation (10). (2) Stimulate every cell simultaneously in the same 1D cable but with diffusive coupling (coupled array in figure 1 (circles in figure 1(B)). (3) Assign the same APD spatial profile a 0 = [a 0 (1), a 0 (2), . . . , a 0 (L)] obtained in step (1) to the iterated map model as initial condition and then iterate equation (16) one step for a given set of ε and σ to result in a spatial APD profile . Note that at the pacing interval we used, iterating equation (16) more steps maintains the same APD profile. (4) Calculate the distance δ(ε, σ ) between the APD profile of the iterated map and that of the ionic model in the presence of coupling, i.e.
We say the APD spatial profile of the coupled map model is best fit to the APD spatial profile of the ionic model when δ reaches a minimum. Examples of APD distribution from an uncoupled array (squares), from the same array but diffusively coupled (circles), and from the coupled iterated maps (triangles). The APD distribution was generated using equation (10) with γ = 0.001. (C) Contour plot of distance function δ(ε, σ ) in the σ -ε parameter space.
ε-σ space, showing that δ reaches a minimum for a certain set of ε and σ . At this optimal case, the APD profile of the map model agrees very well with that from the ionic model (compare triangles with circles in figure 1(B)). Thus, we choose the coupling strength ε and space constant σ at this optimal case for the map model equation (16).
According to the theoretical prediction (equation (18)) based on Echebarria and Karma [41], σ 2 is linearly related to the product of diffusion constant D and APD. In figure 2(A), we show σ 2 versus D obtained by the method mentioned above for the original BR model. The results can be fit to σ 2 = 4.2 × 10 5 D. The APD of the original BR model is about 308 ms. Using equation (18), we obtain which differs substantially from our simulation results, i.e. the slope of the theoretical prediction is 6.5 times of that from the simulation of the ionic model. In figure 2(B), we show σ 2 versus APD for D = 0.001 cm 2 ms −1 by altering APD in two ways: (i) changing the maximum conductance of the slow inward currentḠ si (red squares); (ii) changing that of the timedependent potassium currentḠ x 1 (black squares). Interestingly, in the former case, σ 2 increased first but then saturated as APD increased, while in the latter case, it increased almost linearly with APD in the parameter range we studied. We can linearly fit the data to σ 2 = 1.2a (green line) for the latter case or σ 2 = 2.5a (blue line) for the increase phase of the former one. However, the prediction from equation (18) is  Again the theoretically predicted slope based on Echebarria and Karma [41] is much larger than our numerical results. Similar results (data not shown) as in figure 2(B) were also obtained when we used phase I of the Luo and Rudy model [5].
To gain more insights into the difference between the theoretical predictions and our numerical results as well as the nonlinear relation between σ 2 and APD, we show the voltage gradients at cell 50 (indicated by the red arrow in figure 1(B)) for three different APDs by changingḠ si (red lines) and by changingḠ x 1 (black lines) in figures 2(C)-(E). The voltage gradient was larger (figure 2(C)) when we shortened APD by reducingḠ si than by increasinḡ G x 1 from control, but it became smaller (figure 2(E)) when we increased APD by increasinḡ G si than by decreasingḠ x 1 from control. A larger voltage gradient may effectively increase the distance of voltage diffusion and thus result in a larger space constant σ for the coupled iterated maps. This may explain the difference shown in figure 2(B). Also note that the time span of the voltage gradient is not equal to the duration of the action potential, which may explain why the slopes shown in figure 2(B) are much smaller than the theoretical predictions of equation (18).
Besides the effects of action potential morphology and duration on the space constant σ , cable length may also affect the coupling strength and space constant. Figures 3(A) and (B) show the coupling strength ε and space constant σ versus the cable length L, respectively. ε increases while σ decreases as L decreases, and both can be fit by a single exponential function.

Stability of spatially uniform states.
Diffusive coupling may change the stability of the spatially uniform states [28], [52]- [55]. In a previous study [56], we showed that in a 1D cable of BR model with uniform stimulation, the uniform state could lose stability only in the parameter regime at which chaotic dynamics occurs in the single cell. How this property is maintained in the coupled map model needs to be studied. An important advantage of the coupled iterated map model is that the stability conditions of the uniform states can be analytically derived. When a small perturbation was given to the uniform trajectory a n , i.e., a n (i) = a n + δa n (i), equation (16) can be linearized as: where f = d f (d n ) dd n is the slope of APD restitution curve at d n . If we set δa n ( j) = δ n e i2πjk/L , then from equation (21) we have: where k = 0, . . . , L − 1. By the definition of the maximum Lyapunov exponent, λ = lim n→∞ 1 n ln δ n δ 1 , one obtains the Lyapunov exponent for the kth mode as: where λ 0 is the maximum Lyapunov exponent of the periodic or chaotic state of an isolated cell.
In equation (23), if  ε versus σ for k = 1 (blue), 2 (cyan) and 3 (red) from the relation ε = 1/(2 α j=1 w j sin 2 (π jk/L)) (from equation (25)). (B) ε versus k for σ = 20. (C) The maximum Lyapunov exponent versus wavenumber k from the simulation of the ionic model for the control BR model with stimulation period T = 95 ms (red squares) and from equation (22) for ε = 0.86 and σ = 13.9 (blue line). The wavenumber k is related to k in equation (22) as k = 2π k/L x. then λ k < λ 0 , which indicates that no new instability can be induced by the spatial coupling if the isolated cell dynamics is stable (includes either steady states or periodic states). In other words, if the uniform state is always stable. Figure 4(A) plots the critical ε versus σ for k = 1, 2 and 3, while figure 4(B) plots the critical ε versus k for σ = 20. The limit of ε for large k is 2, i.e., as long as ε < 2, no new instability can be induced by the spatial coupling and the uniform state is always stable. Our numerical estimation of ε = 0.85 indicates that no new instability can occur to the uniform states in the coupled iterated map model equation (16), which agree with our simulation results in the 1D cable of the BR model. But in the chaotic regime, the uniform state can lose stability even though ε < 2. In figure 4(C), we show the maximum Lyapunov exponents (symbols) versus the wavenumber k from the ionic model obtained using the method mentioned in section 2 for a stimulation period T = 95 ms. The maximum Lyapunov exponent for the single cell is λ 0 = 0.28 (λ 0 = λT in which λ was calculated from equation (13)). Inserting λ 0 , ε = 0.85 and σ = 13.9 into equation (23), we obtained the same Lyapunov exponent spectrum (line in figure 4(C)), except for large k . This shows that the coupled iterated map model equation (16) with proper coupling strength and space constant can reproduce the stability conditions of the ionic models well.
To have a direct comparison on this property, we numerically simulated the ionic model and the coupled map model, a cable with 250 cells. Since in the original BR model, the APD restitution curve is monotonic, complex dynamics (periodic states with period higher than 2 and chaos) occurs as pacing interval T decreases only after 2 : 1 conduction occurs [34,56]. In this case, when 2 : 1 conduction occurs, the coupled iterated map is no longer suitable for simulating the losing of stability of the uniform state since 2 : 1 block may occur at one location but not the other. But when the APD restitution curve is nonmonotonic, complex dynamics can occur before 2 : 1 block occurs [34,56]. Therefore, we modified the BR model as previously [57] to result in a nonmonotonic APD restitution curve which is shown in figure 5(A). Figure 5(B) shows the bifurcation diagram from an isolated cell of the BR model, in which APD versus T is plotted. We then transformed the APD restitution data file shown in figure 5(A) into a piece-wise linear function (the increment in DI was 0.05 ms) and then use this function and equation (1) to generate the bifurcation diagram (figure 5(C)) for the iterated map model. Unfortunately, we were not able to obtain the same bifurcation diagram from the two types of models, probably due to the memory and other effects that exist in the ionic model. However, the bifurcation diagrams agree qualitatively. (A better quantitative agreement can be achieved if we are allowed to change APD restitution curve, see another bifurcation diagram in a previous study [34].) To quantify the stability of the uniform state in the 1D cable, we defined a standard deviation of the APD as: δ n = L i=1 [a n (i) −ā n ] 2 /L, where L = 250 is the length of the cable. For the ionic model, a small (1% of the magnitude of each variable) random perturbation to the variables of equations (7)-(9) was given at time zero. For the coupled iterated map model, a small random perturbation was given to the DI of each cell in equation (16) at the beginning of the iteration. Based on our stability analysis, the uniform state in the coupled iterated map model equation (16) is always stable when ε < 2 as long as the dynamics of the single cell is stable. However, the uniform state can lose its stability when ε > 2 even though the dynamics of the single cell is stable. The optimal coupling strength is estimated around 0.8 for the BR model, which satisfies ε < 2. This condition is also satisfied for other action potential models we studied, such as the Luo and Rudy model [5]. However, since ε and σ may be strongly model dependent, further studies are needed to verify whether equation (16) is an appropriate model for action potential dynamics in 1D cable systems.

No-flux boundary conditions.
In the ionic model, the no-flux boundary conditions are: How these boundary conditions are translated into the coupled iterated map model is not clear. The iterated map model equation (2) was originally used with periodic boundary conditions to study reentrant dynamics in 1D rings [25,28], but it was also applied to 1D cables with truncated boundary conditions [27,42], i.e. terms with indices larger than L or smaller than 1 are dropped in equation (2). If we are to find the optimal coupling strength ε and space constant σ for the coupled iterated map model equation (16) with the truncated boundary conditions, large error occurs not only at the two ends of the cable but also at the middle of the cable (compare the triangles with the circles in figure 6(A)). Note that the diffusive coupling in the ionic model corresponds to the space constant σ in the coupled iterated map model, therefore, the no-flux boundary conditions should also correspond to the same spatial scale in the coupled iterated map model. We then applied the following boundary conditions for the coupled iterated map model: a n (1 − k) = a n (k), a n (L + k) = a n (L − k + 1), k = 1, . . . , α, where L is the length of the 1D cable and α the neighbors of coupling and we set α = L. By applying the boundary conditions equation (26), the optimal fit agrees very well with the ionic model (compare the triangles with the circles in figure 6(B)). The optimal ε and σ are almost the same as in the periodic boundary case. The dependence of ε and σ on L is also exponential (figures 6(C) and (D)).
As to the stability of the spatially uniform state, the maximum Lyapunov exponent can be derived according to Heagy et al [54,58] as Time APD n DI n APD n +1 Equation (27) is the same as equation (23) but L is substituted by 2L. The same is true for the ionic model. In other words, if we insert the proper set of ε and σ into equation (27), we obtain the same spectra as in figure 4(C) for both the ionic model and the coupled map model.

Coupled iterated map model II
Although the form of equation (16) has been widely used in studying nonlinear dynamics, it was phenomenologically introduced by Kaneko [43,44] instead of deriving directly from spatially extended systems. The terms used for diffusive coupling in equation (16) are predicted based on the restitution relations, i.e. the DI from the previous beat, instead of the present APD itself. In the cardiac systems, as illustrated in figure 7(A), the adjustment of APD due to diffusive coupling is a consequence of the voltage difference at the present beat, i.e. the diffusive coupling correction for the (n + 1)th beat of the map should come from the coupling of the (n + 1)th beat. By considering this, we modify the coupled iterated map model equation (16) into the following coupled iterated map model: The difference between equation (28) and (16) is that f [d n ( j)] in the diffusive coupling term was substituted by a n+1 ( j). We used the same method to determine the optimal ε and σ . Contrary to the coupled map model equation (16), no single set of optimal ε and σ can be determined. Instead, ε and σ have a power-law relation ( figure 7(B)) at which the distance δ reaches almost the same minimum value ( figure 7(C)). In other words, in this model, for a coupling strength ε, there is always a corresponding space constant σ that the APD spatial profile of the coupled iterated map model can well fit to that of the ionic model.
Similarly, we can analytically derive the stability conditions for the spatially uniform states. The linearized coupled maps are: Inserting δa n ( j) = δ n e i2πjk/L into equation (29), one obtains, which can be written as  From equation (31), one obtains the Lyapunov exponent spectrum as This is different from equation (23) since the logarithm term in equation (32) is always positive and therefore λ k < λ 0 holds for any ε. Interestingly, inserting λ 0 = 0.28 into equation (31), we obtain the same Lyapunov exponent spectrum as the ionic model except for large wavenumber k ( figure 8). This indicates that the coupled iterated map model equation (28) can capture the action potential dynamics in the ionic model well, i.e. the spatially uniform state is always stable unless in the chaotic regime [56].

Wave propagation in 1D tissue
Following the previous studies [22,25,27,28,36,39,42], the CML model developed in this study can be directly extended to study the reentry dynamics in a 1D ring or wave propagation dynamics in a 1D cable.

An integral-delay model for reentry in a 1D ring
The integral-delay model for reentry in a 1D ring was first developed by Courtemanche et al [22], which relates DI (d), CV (θ ), and APD as:  (32) with ε = 2 and σ = 44. However, any ε and σ combination that falls on the relation shown in figure 7(B) can also fit the data well.
where L is the length of the ring. θ (d) is the CV restitution function, and d(x + L) and a(x + L) are the DI and APD at the location (x + L), respectively. Using our CML model equation (16) and changing it for a continuous medium, we have, where w(y) = e −y 2 /2(σ x) 2 / √ 2π σ x is a continuous form of equation (17). We also assumed L to be large so that L −L w(y) dy ≈ 1. Inserting equation (34) into (33), we obtain the revised integral-delay model for reentry in a 1D ring as: Linearizing equation (35) at the steady state d * , i.e., inserting d(x) = d * + δe λx into equation (35), and following the same analysis by Comtois and Vinet [28], we obtain the following characteristic equation: where f is the slope of the restitution curve. At the bifurcation point, Re(Q) = 0 and Im(Q) = q. Inserting Q = iq into equation (36), we obtain: When ε = 1, f c = e ηq 2 which is just the criterion derived by Comtois and Vinet [28]. Figure 9 plots f c versus η from equation (37) for different coupling strength ε, showing that increasing the space constant σ (σ 2 ∝ η) or the coupling strength ε increases the stability of the reentry, i.e. a larger slope of the APD restitution curve is needed for the instability to occur. In other words, in an isolated cell (σ = 0 or ε = 0) the instability leading to alternans occurs at f c = 1, but in a diffusively coupled 1D ring, the instability occurs at f c > 1.

Coupled iterative map model for wave propagation in a 1D cable
When a 1D cable is periodically paced with a period T 0 , the period at any location x is T 0 plus the difference in conduction time to x of two consecutive beats, i.e. [24,30,33,37] Based on equation (16), we have: and since d n (x) = T n (x) − a n (x), we then have the map equation for the pacing cable as: or in a discretized form: where x is the cell length. Here, we use the no-flux boundary conditions equation (26). Using equation (40) or (41), one can study how coupling strength affects the wave dynamics in a paced cable, such as discordant alternans and conduction block, etc.
To study the effects of coupling strength and space constant on discordant alternans, we numerically simulated equation (41) using the following APD and CV restitution functions that were obtained from previous studies [32,33], i.e. and where d c is the critical DI below which a stimulus fails to elicit an action potential. We set a 0 = 220 ms and θ 0 = 0.5 mm ms −1 in this study. Figure 10(A) plots the APD restitution curve of equation (42) and figure 10(B) shows the transition from period-1 to APD alternans under rapid pacing by iterating equation (1) using the function in equation (42). APD alternans occurs when the pacing period T 0 < 208 ms. Figure 10(C) shows two CV restitution curves of equation (43) with different τ θ , corresponding to normal and ischemic conditions [59]. Figures 10(D)-(F) show APD distributions for two consecutive beats (beat 199 and 200) when τ θ = 10 for the coupling strength ε = 0 (D), ε = 1 and the space constant σ = 10 (E), and ε = 1 and σ = 20 (F). As the coupling strength ε and/or the space constant σ increase, the magnitude of APD alternans decreases and the nodal distance increases, indicating that stronger coupling or larger space constant tends to stabilize the system. Note that only one node exists in figure 10(F) but two in the other two panels. In the case of broadened CV restitution curve (τ θ = 50, figures 10(G)-(I)), the number of nodes in the same cable is much larger. The alternans magnitude reduces largely in this case and even more so at the beginning of the cable, indicating that the stabilizing effect of the coupling is even more prominent. An interesting phenomenon is that the alternans amplitude maintains roughly the same in the cable when ε = 0 (figure 10(G)), but it increases from the beginning to the end of the cable when ε = 0 (figures 10(H) and (I)). This same phenomenon was also observed previously in simulations of ionic models [41,59]. To show how the nodes of discordant alternans move in the cable, we plotted the positions of all the nodes versus the pacing beat number n in figure 10(J)-(L) for τ θ = 50. Since the initial condition was a uniform distribution of APD in the cable, the first node formed at the far end of the cable and then more and more nodes formed. When ε = 0, the formed nodes stabilize eventually (figure 10(J)), however, when ε = 0, the formed nodes never stabilize (figures 10(K) and (L)). They form at one end (i = 800) but terminate at the other end (i = 1). Moving nodes have been observed in ionic models and theoretically explained by Echebarria and Karma [26,41].

Conclusions and discussion
In this study, we have developed two coupled iterated map models for the action potential excitation and propagation dynamics of 1D coupled cardiac cells. We first modified previous coupled iterated map models developed by Vinet [25] and by Echebarria and Karma [41] into a standard coupled map form of diffusive coupling with an arbitrary coupling strength (equation (16)). By comparing the spatial APD profiles from the heterogeneous 1D ionic model and the coupled iterated map model, we showed that an optimal set of coupling strength and space constant could be determined and the coupling strength did not obey the normalization requirement as in the previous models [25,41]. Our results show that the theoretically derived equation by Echebarria and Karma [41] from a simple action potential model substantially overestimated the space constant of coupling. In addition, the coupled iterated map model can capture the stability conditions of the spatially uniform states well, i.e., in both the coupled iterated map model and the 1D ionic model the uniform state is always stable unless in the chaotic regime. We also developed another coupled iterated map model in which the coupling strength and the space constant are not independent but have a power-law relationship. Since both the models can fit the APD spatial profile and capture the stability conditions from the ionic model well, we cannot judge whether one is better than the other. However, there is a major difference between the two models, i.e., in the former model, when the coupling strength is greater than 2, spatial instabilities can occur even though the cellular dynamics is stable, while no such instabilities can occur for any coupling strength in the latter one. This needs to be verified by further studies using different action potential models. We also extended the coupled iterated map model to study the dynamics of reentry in a 1D ring and action potential propagation in a 1D cable. Our study showed that the diffusive coupling suppressed instabilities caused by a steep APD restitution curve. In other words, in an isolated cell, alternans occurs when the slope of the APD restitution curve exceeds one at the equilibrium point, but in tissue, a larger slope is needed to cause APD alternans. We also showed that diffusive cell coupling could cause the nodes of discordant alternans to propagate, which agrees with the results by Echebarria and Karma [26,41].
Since cardiac action potential is regulated by complex biophysical and physiological processes, it is difficult to analyze the mechanisms of the action potential dynamics, even in an isolated cell. Therefore, to develop a valid low-dimensional representation of the complex cardiac system is of great importance, and our present study has made one step forward based on the efforts of many other investigators. Further studies are still needed to improve the coupled map models, such as to couple the maps describing the calcium cycling dynamics [29,31,34,60] into the coupled iterated map models to describe the spatiotemporal dynamics of excitation-contraction coupling.