In this study, we analyze a nonlinear map model of intracellular calcium (Ca) and voltage in cardiac cells. In this model, Ca release from the sarcoplasmic reticulum (SR) occurs at spatially distributed dyadic junctions that are diffusively coupled. At these junctions, release occurs with a probability that depends on key variables such as the SR load and the diastolic interval. Using this model, we explore how nonlinearity and stochasticity determine the spatial distribution of Ca release events within a cardiac cell. In particular, we identify a novel synchronization transition, which occurs at rapid pacing rates, in which the global Ca transient transitions from a period 2 response to a period 1 response. In the global period 2 response dyadic junctions fire in unison, on average, on alternate beats, while in the period 1 regime, Ca release at individual dyads is highly irregular. A close examination of the spatial distribution of Ca reveals that in the period 1 regime, the system coarsens into spatially out-of-phase regions with a length scale much smaller than the system size, but larger than the spacing between dyads. We have also explored in detail the coupling to membrane voltage. We study first the case of positive coupling, where a large Ca transient promotes a long action potential duration (APD). Here, the coupling to voltage synchronizes Ca release so that the system exhibits a robust period 2 response that is independent of initial conditions. On the other hand, in the case of negative coupling, where a large Ca transient tends to shorten the APD, we find a multitude of metastable states which consist of complex spatially discordant alternans patterns. Using an analogy to equilibrium statistical mechanics, we show that the spatial patterns observed can be explained by a mapping to the Potts model, with an additional term that accounts for a global coupling of spin states. Using this analogy, we argue that Ca cycling in cardiac cells exhibits complex spatiotemporal patterns that emerge via first or second order phase transitions. These results show that voltage and Ca can interact in order to induce complex subcellular responses, which can potentially lead to heart rhythm disorders.

The heart can exhibit a variety of heart rhythm abnormalities that are responsible for cardiac arrhythmias. One of them is a beat-to-beat alternation in the action potential duration (APD), which is referred to as alternans. Often, its origin lies in a dysregulation of intracellular calcium cycling, which is mediated by stochastic signaling between ion channels in thousands of sub-micron scale domains. We show that these signaling units exhibit complex spatiotemporal dynamics that depends on the interaction between voltage and calcium. We argue that these subcellular patterns may promote arrhythmias by inducing spatial heterogeneity in heart tissue.

## I. INTRODUCTION

The contraction of a cardiac cell is tightly controlled by the voltage across the cell membrane. This process is mediated by Ca ions which flow across the membrane and activate a variety of intracellular signaling cascades.^{1,2} The key players in this signaling process are the L-type Ca channel (LCC), which is a voltage sensitive Ca channel located at the cell membrane, and the Ryanodine receptor (RyR), which controls the flow of Ca from the major intracellular Ca store, the sarcoplasmic reticulum (SR). In heart cells, these channels are situated within thousands of dyadic junctions which are roughly pill box sized regions of height ∼10 nm and diameter ∼100 nm. Within these junctions, there are roughly 10–100RyR channels which form a tight cluster and which face 1–5 LCC channels [Fig. 1(a)]. RyR channels have the important property that they transition from a closed to open state in a Ca dependent manner. Thus, a single opening (or a few openings) of an LCC channel in the dyad is sufficient to induce a domino-like opening of the nearby RyR cluster. This process leads to a release of Ca into the cell which then diffuses and initiates downstream signaling processes. Once the Ca is released into the cell, then powerful transporters pump the Ca back into the SR. This process allows electrical activity on the cell membrane to activate internal cellular processes at each beat of the cardiac cycle.

An important feature of the signaling between LCC and RyR is that it is highly nonlinear.^{3–5} This is because RyRs transition to the open state with a rate that is a nonlinear function of the local Ca concentration. This nonlinearity originates at the single channel level, where multiple Ca binding processes are required to induce a transition to the RyR open state. Thus, a high degree of molecular cooperativity is required to induce channel openings. This nonlinearity is amplified further since RyRs form tight clusters so that if one RyR channel opens, then the Ca flux from the SR, which is four orders of magnitude larger than in the cell, will raise the Ca concentration and stimulate further channel openings.^{4} This domino-like process is highly cooperative and makes the response of the RyR cluster a nonlinear function of the local concentration variables. A second important feature of the signaling process at the scale of the dyad is that it is stochastic. This is because channel numbers within the dyadic junction are small and so that fluctuations of just a few channels can have a substantial effect on the activation rate of the cluster. Therefore, Ca signaling between LCC and RyR channels is inherently a stochastic process. While there has been some work in the literature describing the stochastic features of this system, it is still not completely understood how stochasticity and nonlinearity interact to induce different cell wide phenomena.

In a previous study, Alvarez-Lacalle *et al*.^{6} showed that a population of diffusively coupled dyads, when driven periodically, can exhibit a phase transition to a global alternating phase, where most dyads fire in unison only every other beat. When this occurs, the aggregate behavior of the population exhibits an alternating beat-to-beat response referred to as Ca alternans. This transition to alternans in the global signal is important since an alternating sequence of Ca release can induce a corresponding alternating sequence of the action potential duration (APD) in cardiac tissue. When this occurs, the APD in different parts of tissue can change substantially from beat to beat. This heterogeneity in the voltage activity in cardiac tissue makes the heart more prone to abnormal electrical excitations, which can lead to cardiac arrhythmias. In fact, numerous experimental and clinical studies have shown that these APD alternans are correlated with the risk for ventricular fibrillation.^{7,8} Interestingly, the work of Alvarez-Lacalle *et al.*^{6} suggests that this synchronization transition to alternans is a phase transition in the Ising universality class. The relevance of the Ising transition was first pointed out by Restrepo and Karma,^{9} who showed that during alternans, the local signaling dynamics exhibits an Ising symmetry. To explain this symmetry, we first note that a single stochastic dyadic junction, which is driven periodically, can release Ca due to a spark, or remain refractory and do not release. Thus, during alternans, the local SR load will alternate between two values, which we will denote here by $x1$ and $x2$. During pacing, a local junction can exhibit the two alternating sequences

which differ only by a shift of one beat and are therefore dynamically equivalent. It is this dynamical equivalence of the two types of alternating sequences that endows the system with Ising symmetry. Now, for sufficiently strong coupling, there can be a global phase transition such that the majority of dyadic junctions exhibit one of these sequences. When this occurs, most dyads in the system will alternate with the same phase and exhibit a whole cell period 2 response, in much the same way that a spin system undergoes a ferromagnetic transition to an ordered state. This transition exhibits all the features that are observed near the critical point of second order phase transition, such as a diverging correlation length and critical slowing down. Recent studies^{10} show further that these critical properties are insensitive to the underlying mechanism for calcium alternans and can likely be applied to a broad range of detailed computational models used to describe Ca dynamics in cardiac cells.^{11–15}

In this study, we introduce a two dimensional (2D) coupled map system to explore the dynamics of Ca signaling within dyadic junctions. In this model, the local Ca concentration variables evolve from beat-to-beat using a nonlinear map describing the release and uptake of Ca between the cytosol and the SR. A key feature of this model is that release occurs with a probability that is a nonlinear function of the local concentration and voltage variables. Using this model, we extend the work of Alvarez-Lacalle *et al*.^{6} to analyze the dynamics of the system in the strongly nonlinear regime, where the Ca cycling system can exhibit more complex dynamics. Also, we explore in detail the important role played by the bi-directional coupling between membrane voltage and subcellular Ca. In the case where Ca and voltage are uncoupled, we identify a novel synchronization-desynchronization transition, which occurs when the underlying dynamics exhibits an approximate period-3 or higher response. We show further that when this transition occurs, the system coarsens into synchronized patches of a characteristic length scale that is much larger than the dyadic junction spacing but smaller than the system size. We analyze this length scale and show that it depends only on the extent of diffusive coupling and the nonlinearity of the underlying dynamics. In the case, when voltage is coupled to Ca, we consider two distinct cases. The first case is the scenario of positive coupling, where a large Ca transient promotes a longer APD.^{16} In this case, we find that the coupling between voltage and Ca makes the system more unstable to alternans but tends to eliminate higher order periodicities. Furthermore, voltage coupling tends to spatially synchronize out-of-phase alternans by driving domain walls (or nodal lines) out of the system. This result explains the prevalence of experimental findings that show that cardiac myocytes are prone to alternans at rapid pacing rates but rarely exhibit higher order periodicities.^{17–19} The second scenario that we analyze is the case of negative coupling, where a large Ca transient as the effect of shortening the APD. In this case, we show that the system can exhibit a wide range of complex metastable spatial patterns, which are dependent on initial conditions. We argue further that these complex spatially patterns can be explained using an analogy to the well known Potts model in equilibrium statistical mechanics.^{20} Interestingly, we show that various transitions in our model can be mapped to discontinuous first order phase transitions, in sharp contrast to the transition to alternans, which is a continuous second order phase transition. These results demonstrate that stochastic nonlinear signaling units exhibit rich spatiotemporal behavior that can be understood using the tools of statistical mechanics and which may be relevant to the emergence of heart rhythm disorders.

## II. METHODS

### A. Nonlinear stochastic maps describing Ca signaling

Ca cycling in the cardiac cell is due to the stochastic dynamics of thousands of synapse like junctions, where Ca signaling takes place. The basic architecture of local Ca signaling is illustrated in Fig. 1(a). Here, Ca is released at the dyadic junction (dashed rectangle), where LCC channels are in close proximity to an RyR cluster. Ca released from the junction then diffuses into the cytosol and is then pumped back into the SR. In this study, we will refer to the region surrounding the dyadic junction as a Ca release unit (CRU). To model their spatiotemporal dynamics, we will simplify the system to a two dimensional (2D) lattice of CRUs. To describe the dynamics of each CRU, we will apply a nonlinear map that relates Ca concentration variables from one beat to the next. To describe the state of a CRU, we will denote $xij(n)$ to be the total amount of Ca ions in the SR at site $i,j$ and at the beginning of beat *n* [see Fig. 1(b)], also $cij(n)$ will represent the total Ca in the cell interior before Ca is released, and $cijp(n)$ is the peak of the total Ca. Also, the voltage dynamics of the cell will be characterized by the action potential duration (APD) and the diastolic interval (DI), which are denoted at a given beat *n* as $An$ and $Dn$, respectively. As a starting point, we will first consider the case where there are no diffusive interactions between CRUs. In this case, the beat to beat evolution of the SR load and diastolic Ca at site $i,j$ is described by a mapping

where $Rij(n)$ is the total Ca released at beat *n* and $Uij(n)$ is the total Ca pumped back into the SR. Typically, the release of Ca occurs on a time scale of tens of milliseconds, while the uptake of Ca back into the SR is slower and requires several hundred milliseconds. Here, we assume that the fluxes into and out of the cell are small, so that the total Ca is conserved and can be conveniently fixed at

This constraint is a good approximation close to the periodic fixed point, where Ca entry must balance Ca extrusion so that total Ca is constant from one beat to the next. Here, we point out that when the periodic fixed point becomes unstable to alternans, or higher order periodicities, then this constraint may not hold since Ca entry and extrusion will in general not be equal. However, if the flux of Ca across the cell membrane is much smaller than that across the SR, then the nonlinear dynamics of the system will be driven primarily by the internal Ca cycling. In this case, the model reduction due to Eq. (5) should preserve the main dynamical features of the system. For a more complete analysis, see Qu *et al*.,^{15} who have analyzed a deterministic nonlinear map model without the constraint given by Eq. (5).

The release of Ca into the cell is typically taken to be a function of the SR load before release and also on the amount of LCC current. This is because Ca release from an RyR cluster is initiated by LCC channel openings so that the larger the LCC current then the more Ca released into the cell. However, LCC channels inactivate and must recover from inactivation in order to open on the next beat. Therefore, the LCC is sensitive to the previous DI, which is denoted as $Dn\u22121$. Thus, Ca release at beat *n* is well approximated by a function

To model the release function *R* in more detail, we note that release is stochastic, so that there is a probability that the local RyR cluster will fire or not. To incorporate these physiological observations, we will take the probability of release to be a product

To model the voltage dependence of release, we follow the approach of Qu *et al*.^{21} and use a simple functional form

where $\tau Ca$ is the time scale of recovery of the LCC channel and *A* is an adjustable constant that depends on the current kinetics. To model the SR load dependence, we rely on the previous experimental data showing that the amount of Ca released in the cell increases substantially at high SR loads.^{22–25} To model this nonlinear dependence, we use a functional form

where $\nu $ is the Hill coefficient that controls the strength of the nonlinearity and $x\u2217$ is the threshold SR load. To proceed, we write the release at each beat as

where *r* is a constant and $\eta ij(n)$ is a random variable generated at each beat that satisfies

Thus, local release is modeled as a discrete process with probability that is a nonlinear function of the DI and the SR load.

To complete the map, we note that the total uptake flux $Uij(n)$ is the amount of Ca current that the sarcoplasmic reticulum Ca^{2+}-ATPase (SERCA) pumps back into the SR during one pacing interval *T*. The strength of the uptake pump increases with the Ca concentration in the cell so that uptake can be approximated as a function

where the peak of the Ca transient is just the diastolic concentration plus the amount released, i.e., $cijp(n)=cij(n)+Rij(n)$. To simplify further, we see that the peak Ca is just

where

To model the total uptake flux, we will use a simple linear relation

where $u(T)$ gives the dependence of the total uptake flux at beat *n* on the pacing period *T*. Since *U* is the total amount of Ca ions pumped into the SR, then we expect that *U* must increase with period *T*, since this will give more time for the SERCA channels to pump Ca. Here, we will model this effect using a simple linear relationship so that $u(T)=a\u22c5T$, where *a* is a constant. Also, we note that the uptake flux should go to zero when all the Ca is in the SR. Thus, we impose the additional condition that $U=0$ for $xij\u2032\u22651$. Using this formulation, we can now write

which combined with Eq. (14) gives the beat-to-beat mapping of the SR load.

### B. Diffusive coupling between dyadic junctions

In this study, a cardiac cell will be modeled as a 2D lattice of $N=m\xd7m$ CRUs with periodic boundary conditions. To account for diffusive coupling between nearest neighbors, we will include a spatial coupling between lattice sites. For generality, we will introduce a number *l* that is the number of nearest neighbors that Ca can diffuse during a given beat. Now, since Ca release is much faster than the uptake then we can assume that local averaging occurs during the uptake process, so we can replace

where $\u27e8x\u27e9l$ denotes averaging over *l* nearest neighbors in our 2D lattice. In this study, we will average the SR load by computing the local average

where *l* can be interpreted as roughly the distance Ca diffuses during the duration of uptake. Therefore, the beat-to-beat evolution of our SR load obeys the stochastic coupled map

Thus, the SR load on the next beat $xij(n+1)$ is the average release over *l* nearest neighbors, combined with an uptake amount that is determined by this average.

### C. Coupling voltage and Ca

To incorporate voltage dynamics, we rely on the classic restitution function where the APD depends on the diastolic interval on the previous beat $Dn$.^{26,27} However, Ca release can also modulate the APD via two mechanisms.^{16,28} The first is due to the sodium-calcium exchanger (NCX), which is a transporter that pumps Ca ions out of the cell in exchange for the entry of sodium (Na) ions. The effect of this current is to prolong the APD when Ca is released into the cell. The second mode of coupling between Ca release and APD is due to Ca-induced inactivation of LCC. Here, Ca release tends to decrease the whole LCC current and therefore has the effect of shortening the APD. Thus, a given Ca transient can prolong or shorten the APD depending on which effect is dominant in that cell. These cases have been observed experimentally under different physiological conditions and cell types.^{17,29} Here, we follow Shiferaw and Karma^{16} and refer to prolongation (shortening) of the APD as the positive (negative) coupling. To incorporate this coupling between the voltage and Ca, we will use a functional form

where $cnp$ is the peak of the average Ca concentration in the cell at beat *n*. This is just

where *N* is the total number of CRUs in our system. Here, we have made the approximation that $cij(n)\u226aRij(n)$ since the amount of Ca released is typically much larger than the concentration levels before release. We will also consider an explicit model for the APD as a sum of a purely voltage dependent term and a Ca dependent term. So, we will use

where $FV(D)$ represents only the voltage contribution and $\gamma $ represents the coupling between the local Ca release and the APD. In this case, we have that $\gamma >0$ denotes positive coupling and $\gamma <0$ denotes negative coupling. Following Qu *et al*.,^{21} we will model the voltage dependence using a functional form

which models the shape of the restitution curve typically measured in cardiac cells.^{30}

## III. RESULTS

### A. The deterministic limit

An important feature of our coupled map system is that the release is a stochastic variable governed by the random variable $\eta ij$. Here, we will consider the limit where diffusion is fast and the spatial average $\u27e8x\u27e9l$ is over the full lattice so that $l=m$. In this case, at the end of each beat, all CRUs have the same SR load which we denote as $xn=\u27e8xij(n)\u27e9m$. Therefore, we can write

where $\u27e8\eta ij(n)\u27e9m$ is simply the fraction of units in our 2D lattice which fires at that beat. To compute this fraction, we note that the number of dyads which release at a given beat is the number of successes in *N* trails with probability $P=P(Dn\u22121,xn)$. This obeys a binomial distribution $B(P,N)$, where *N* is the number of trials and *P* is the probability of success. In the limit of large *N*, this distribution is approximately Gaussian with average *NP* and variance $\sigma =NP(1\u2212P)$. Therefore, for large *N*, we have that

where $\xi n$ is a Gaussian distributed random variable that satisfies $\u27e8\xi k\xi l\u27e9=\delta k,l$. We can now take the large *N* limit, and the system is described by the deterministic map

Hereafter, this 2D map will be referred to as the deterministic limit of the stochastic coupled map system.

### B. The periodic action potential clamp

We explore now the spatiotemporal dynamics of our lattice in the case where there is no feedback from Ca to voltage, so that $\gamma =0$. In addition, we will consider the case where there is also no coupling between voltage on Ca so that $PCa(D)=1$. Experimentally, this corresponds to the case where the system is driven by a periodic action potential clamp in which the system dynamics is dictated entirely by the Ca cycling system. In this case, the SR load is described by the 1D map

where

Here, we will focus on the strongly nonlinear regime, where the probability of firing $Psr$ is a sharp sigmoid function of the SR load. All model parameters in the simulation are shown in Table I. In Fig. 2(a), we plot the steady state SR load $xn$ for the last 100 beats after pacing to steady state (1000 beats). Here, we observe that the deterministic map possesses a period doubling bifurcation cascade to chaos. This follows from the nonlinearity due to $Psr(xn)$ which yields a nonmonotonic mapping between $xn+1$ and $xn$, so that the system dynamics is qualitatively similar to the standard logistics map.

Parameter . | Description . | Value . |
---|---|---|

$x\u2217$ | Threshold SR concentration | 0.7 |

$\nu $ | Hill coefficient for SR load dependence | $40$ |

$r$ | Release coefficient | $0.7$ |

$a$ | Uptake strength | $10\u22123(ms\u22121)$ |

Parameter . | Description . | Value . |
---|---|---|

$x\u2217$ | Threshold SR concentration | 0.7 |

$\nu $ | Hill coefficient for SR load dependence | $40$ |

$r$ | Release coefficient | $0.7$ |

$a$ | Uptake strength | $10\u22123(ms\u22121)$ |

To explore the dynamics of our stochastic map system, we iterate a system of size $m=200$ with periodic boundary conditions. At each beat, we measure the average SR load of the system defined as

Here, we will first consider spatially uniform initial conditions so that the lattice is initialized at $xij(0)=1$. In a later section, we will explore in detail the system dynamics with non-uniform initial conditions. To plot the bifurcation diagram of the system, we plot $xn$ for the last 100 beats after pacing for 1000 beats, so that the system has reached steady state, at a range of pacing period *T*. In Figs. 2(b) and 2(c), we show the case where the coupling length is $l=2$ and $l=8$, respectively. In the case $l=2$, we see that for $T>T1\u2248600ms$, the steady state $xn$ has a narrow distribution around a well defined average so that the system can be said to converge a period 1 (P1) steady state. At faster pacing rates, so that $T<T1$, then $xn$ fluctuates around two well defined average values, which we will refer to as a period 2 (P2) response. Finally, for $T<T2\u2248350ms$, $xn$ once again returns to an effective P1 state with narrow fluctuations around a well defined average. For stronger nearest neighbor coupling $(l=8)$, the steady state dynamics is even richer and can transition between approximate P1, P2, P3, and P4 regimes. To summarize these results, we measure the steady state values of $xn$ for a range of pacing periods *T* and diffusion lengths *l*. In Fig. 3, we plot the phase diagram of the system identifying the periodicity of the steady state SR load. On top, we show the bifurcation diagram of the deterministic system in order to compare to the stochastic lattice simulation. Here, we see that the stochastic system eliminates the chaotic response observed in the deterministic limit but can still exhibit a broad range of behaviors including period 3 and period 4, depending upon the coupling strength *l*.

#### 1. The spatial distribution of the SR load

In order to understand the phase diagram of the coupled map system, it is necessary to explore the spatiotemporal dynamics of the system. In Figs. 4(a)–4(d), we plot on the left column the spatial distribution of the SR load $xij$ at steady state. On the right, we plot the probability density distribution $\rho (x)$, so that $\rho (x)\Delta x$ gives the probability that $xij$ falls in an interval $[x,x+\Delta x]$. In Fig. 4(a), we consider the case $l=2$ and with $T=650ms>T1$, which corresponds to the P1 state shown in the bifurcation diagram in Fig. 2(b). On the right panel, we plot the density distribution $\rho (x)$ showing that it has a jagged shape with an overlap that is mostly peaked at a well defined average. This jagged distribution can be explained from the fact that we are iterating a discrete beat-to-beat map. Thus, the SR load $xij$ at the end of each beat depends on a discrete number of previous release (or no release) events. In Fig. 4(b), we pace the system into the P2 regime $(T=450)$ and find that the SR load alternates from one beat to the next, and it is effectively spatially uniform with fluctuations around two well defined values. On the right panel, we show the density $\rho (x)$ during steady state, showing that the distribution of the SR loads separates into two distinct peaks. Finally, in Fig. 4(c), we chose $T=275ms<T2$, where the system exhibits a global P1 behavior. Here, we find that the SR load loses its strong spatial correlation observed for $T>T2$ and appears to coarsen into a spatially disordered state with large fluctuations on scales much smaller than the system size. On the right panel, we show the density distribution which shows that $xij$ is broadly distributed. In particular, $xij$ is distributed in the range $[0.45\u22120.75]$, while for $T>T2$, the range of loads is much sharper and in the range [0.7–0.8]. Thus, our numerical simulations reveal that the spatiotemporal behavior of the P1 state above $T1$ is distinct from the state below $T2$. In the case where the spatial coupling has been increased $(l=8)$, we find that the steady state SR load can exhibit even richer spatial patterns. In Fig. 4(d), we show steady state plots of the SR load when *T* is picked so that the system exhibits a global P3 pattern. Here, we see that the SR load is effectively spatially uniform, and the steady state $xij$ is distributed around 3 distinct peaks.

#### 2. Dependence on initial conditions

The spatiotemporal dynamics observed in Figs. 4(a)–4(d) was computed by setting the initial conditions to be spatially uniform $[xij(0)=1]$. Here, we will explore the case where $xij(0)$ is chosen from a uniform distribution in the interval $[0,1]$. In Fig. 5(a), we plot the spatial distribution of the SR load after 1000 iterations using the same model parameters as that used to compute Fig. 4(b), in which the spatially uniform initial condition evolves to a homogeneous P2 pattern. In this case, we find that starting from random initial conditions, the system typically evolves to a heterogeneous spatial pattern in which different regions alternate out-of-phase [Fig. 5(a)]. These out-of-phase regions are demarcated by domains walls (nodal lines) which can form complex spatiotemporal patterns. These domain walls proceed to fluctuate and evolve in the system and can last for thousands of beats. In Fig. 5(b), we show the SR load distribution on alternate beats along the dashed line shown in Fig. 5(a). Indeed, we see that the SR load alternates out of phase across a nodal point, which corresponds to the intersection of the domain wall. Here, we note that while the spatial distribution of $xij$ exhibits rich patterns dictated by the domain wall dynamics, the SR load distribution fluctuates around the two peaks shown in Fig. 4(b). In the case of higher periodicities, we find that random initial conditions robustly evolve to a spatially disordered state. In Fig. 5(c), we show the steady state distribution of the SR load when the system is paced to steady state at $T=270ms$, so that spatially uniform initial conditions evolve to a global P3 pattern. In this case, we find that the system evolves to a complex disordered pattern, where the system coarsens into regions that are out-of-phase or desynchronized. In Fig. 5(d), we show the spatial profile along the dashed line on 3 alternate beats, showing that the steady state spatial distribution is spatially desynchronized.

#### 3. Spatial correlations

A main result of this study is the observation that the transition from P2 to P1, which occurs when $T<T2$ [Fig. 2(b)], yields a spatially disordered state that is distinct from the P1 dynamics that occurs when $T>T1$. To characterize the system more fully, we have also computed the spatial correlations between different points of our lattice. Specifically, we have computed the Pearson correlation between the SR load at sites $xij$ and $xkl$. This quantity is defined as

where we consider iterations starting at beat $K1$ to beat $K2$. Here, we consider $K1=500$ and $K2=1000$ in order to capture the correlations at steady state. Similarly, the average is computed at steady state so that $xij=[1/(K2\u2212K1)]\u2211n=K1K2xij(n)$. This correlation function is computed in our lattice for 5 different points, and the correlation length $\xi $ is determined by fitting an exponential function $f(r)\u221dexp\u2061(\u2212r/\xi ),$ where $r=(i\u2212k)2+(j\u2212l)2$. In Fig. 6(a), we compute $\xi $ as a function of *T* for the case $l=2$. Indeed, we find that the correlation length increases substantially near the border to the P2 phase, which indicates that the system undergoes a synchronization transition near $T1$. This result is consistent with previous findings of Alvarez-Lacalle *et al*.,^{6} who show that the P1 to P2 transition occurs via a synchronization transition that is likely in the Ising universality class. Additionally, here we find that $\xi $ once again substantially decreases as $T<T2$, which suggests that the system once again undergoes a phase transition from a synchronized to desynchronized state. This desynchronized state is characterized by large fluctuations in the local SR load, which coarsen according to a length scale $\xi $. This state is distinct from the $T>T1$ state, where the SR load is sharply peaked around a well defined average. In Figs. 6(b) and 6(c), we plot $\xi $ for increasing values of the spatial coupling parameter. For the case $l=4$, we notice that $\xi \u223c14$ for $T>T1$, while for $T<T2$, we find that $\xi \u223c23$ just below the transition. Similarly, when $l=8$, we find that $\xi $ increases substantially near the P2 and P3 transitions.

### C. The free running AP: The case of positive coupling

In this section, we will explore the dynamics of the system when the voltage is free running and coupled to the Ca cycling system. In this case, the deterministic limit is given by the 2 variable map given by Eqs. (27)–(29). In this section, we will consider the case $\gamma >0$, which is the case of positive coupling, where a large Ca transient promotes a longer APD. This coupling is believed to be the typical relationship observed in cardiac myocytes across a range of species.^{17,30,31} All model parameters describing the voltage map are given in Table II. In Fig. 7(a), we plot the bifurcation diagram of the deterministic system, where the Ca cycling system simulated in Fig. 2 is coupled to the voltage model with a coupling parameter $\gamma =50ms$. Here, we observe that positive coupling expands the parameter range of the P2 and P4 regimes, while eliminating higher order periodicities. In Fig. 7(b), we show the phase diagram of the stochastic system where the initial conditions are chosen to be spatially uniform $[xij(0)=1]$. Here, we see that the stochastic system shows a broad range of P2 dynamics, with a small region of P4 dynamics which can be observed for strong spatial coupling $l\u22658$. It should be noted here that the P4 regime is sensitive to fluctuations near the phase boundaries and the 4 distinct states are only intermittently observed. In Fig. 7(c), we repeat the same computations using a larger Hill coefficient $\nu =50$ in order to induce higher order period doubling bifurcations in the deterministic dynamics. In Fig. 7(d), we compute the phase diagram and again observe that the P2 regime is expanded, with a small region of approximate P4 dynamics observed at high spatial coupling. We have repeated these simulations using random initial conditions and find that the steady state phase diagram remains unchanged. Also, visualization of the SR load reveals that at steady state, the system is effectively spatially uniform and domain walls are not observed. These findings show that the coupling between voltage and Ca tends to expand the P2 region of the phase diagram and makes the steady state P2 response independent of initial conditions. We note here that this result is consistent with electrophysiological studies of Ca and voltage alternans, which show that the dominant response pattern at high pacing rates is P2 and that higher order periodicities are rarely observed.^{24}

Parameter . | Description . | Value . |
---|---|---|

$A$ | LCC kinetics constant | 2 |

$\tau Ca$ | LCC channel recovery time | $100ms$ |

$\tau $ | Restitution decay time | $100ms$ |

$Ao$ | Restitution parameter | $250ms$ |

$A1$ | Restitution parameter | $2$ |

$\gamma $ (positive) | Coupling parameter | $50ms$ |

$\gamma $ (negative) | Coupling parameter | $\u221255ms$ |

Parameter . | Description . | Value . |
---|---|---|

$A$ | LCC kinetics constant | 2 |

$\tau Ca$ | LCC channel recovery time | $100ms$ |

$\tau $ | Restitution decay time | $100ms$ |

$Ao$ | Restitution parameter | $250ms$ |

$A1$ | Restitution parameter | $2$ |

$\gamma $ (positive) | Coupling parameter | $50ms$ |

$\gamma $ (negative) | Coupling parameter | $\u221255ms$ |

#### 1. Dependence on initial conditions

An important feature of the positive coupling case is that the steady state spatial distribution of the SR load is independent of initial conditions. This is in sharp contrast to the uncoupled Ca map model, which exhibited long lasting domain walls [Fig. 5(a)]. In Figs. 8(a)–8(d), we plot the spatiotemporal distribution of the SR load for a system paced at a period of $T=500ms$, in which the steady state of the system is in the P2 regime. Here, we consider initial conditions where our $200\xd7200$ system is heterogeneous with a domain wall through the center of the system. To construct these initial conditions, we have simply set the SR load on the left and right sides of the lattice at $xij=0.65$ and $xij=0.80$ on the right side. Upon pacing, we find that the domain wall gradually drifts until after 1200 beats the system is dominated by one alternans phase. In Fig. 8(e), we show the global average SR load $xn$ as a function of beat number. Here, we see that the global alternans amplitude gradually increases as the domain wall is pushed out of the system. We also show that as $\gamma $ is reduced, the rate at which the domain walls drift out of the system becomes slower, which demonstrates clearly that the bi-directional coupling to voltage is responsible for the gradual synchronization of the system. These results demonstrate that the global voltage coupling tends to eliminate domain walls in order to favor the dominant phase in the system. This dominant phase is essentially determined by the initial conditions, which randomly picks a phase during the first few iterations of the system. Therefore, when $xij$ is taken from a random distribution, domain walls will typically form and will in time drift out of the system in order to favor the global average that was established randomly due to the initial conditions.

### D. The free running AP: The case of negative coupling

In this section, we analyze the case of negative coupling, where a large Ca transient tends to decrease the APD. This scenario is not commonly observed in electrophysiological studies but has been noted in the literature in some animal types.^{29} To model this case, we use the model parameters in Tables I and II and set $\gamma =\u221255ms$. To achieve higher order periodicities, we have also increased the instability of Ca cycling by increasing the Hill coefficient to $\nu =50$. In Fig. 9(a), we show the deterministic limit of the map that displays a typical period doubling cascade to chaos. In Fig. 9(b), we show the bifurcation diagram for the stochastic map model with spatial coupling $l=2$. Here, we have considered uniform initial conditions with $xij(0)=1$. In this case, we observe that the steady state global average $xn$ is not unique and can attain a range of values. In effect, the system is multistable and does not evolve to a unique steady state. To uncover the source of this multistability, we have visualized the SR load on the lattice for 3 independent simulation runs at $T=400ms$ (red dashed line) and with an initial APD of 250 ms. We find that these different simulation runs, with identical initial conditions, evolve to distinct steady state patterns. It should be pointed out here that the spatial patterns formed are also sensitive to the initial value of the APD. If the initial APD is chosen such that there are large initial beat-to-beat changes in APD, then these complex subcellular patterns are washed out by the voltage feedback, and the system will tend to evolve toward a spatially homogenous state. However, if the APD does not alternate during the initial beats, then the observed spatial discordant patterns are formed. In Fig. 9(c), we repeat the same simulations using initial conditions so that $xij(0)$ is taken from a uniform distribution in the interval [0, 1], and the initial APD is taken to be 250 ms. On the right hand side, we show the steady state spatial distribution of the SR load at *T*=520 ms, for 3 independent simulation runs with identical initial conditions. In this case, we observe again that the steady state pattern is not unique and varies between simulation runs. Thus, in the presence of negative coupling, subcellular Ca can display complex spatiotemporal patterns which are dependent on the initial conditions.

To uncover the underlying mechanism for these patterns, we study the time evolution of domain walls that we introduce into the system by hand. In Fig. 10(a), we consider the spatiotemporal dynamics with initial conditions such that a domain wall joins two sides of our 2D system. In this case, the yellow region is started with $xij(0)=0.8$, and the purple region is $xij(0)=0.6$. For these initial conditions, we observe that the bifurcation diagram [Fig. 10(b)] displays only very small amplitude alternans in the global average SR load $xn$. Here, we see that the domain wall evolves so that effectively half of the system exhibits alternans in the opposite phase, with two effectively straight domain walls joining two sides of the system. Once the steady state pattern forms, we find that the global alternans amplitude is small due to the almost exact cancellation of the two regions of spatially discordant alternans. In Fig. 10(b), we repeat the same simulation but this time initialize the system with the square pattern shown. In this case, we find that the steady state bifurcation diagram displays a small amplitude alternans that is much larger than in Fig. 10(a). Visualization of the subcellular structure shows that the initial square pattern evolves into a circular region out-of-phase with its surroundings. This region then proceeds to grow until it stops at a fixed radius. A close examination of the steady state pattern reveals that the circular region has a smaller area than the remainder of the 2D system. This explains why the bifurcation diagram [Fig. 10(d)] settled at a small amplitude of alternans at steady state. These results show that the time evolution of patterns is highly sensitive to initial conditions. In particular, here we find that a closed domain wall evolves into a circular region at steady state, while a domain wall joining two sides of the system drifts to the center so as to lead to cancellation in the global alternans amplitude.

## IV. DISCUSSION

### A. Analogy to equilibrium statistical mechanics

In a previous study, Alvarez-Lacalle *et al*.^{6} studied the onset of the alternans (P2 response) and showed that it occurred via a synchronization transition with critical exponents consistent with the Ising universality class.^{32} A detailed numerical analysis of the dynamics near the P2 transition revealed important features consistent with a second order phase transition, such as a diverging correlation length and critical slowing down. These findings reveal that important features of the steady state spatial patterns of the coupled map system can be understood using an analogy to equilibrium statistical mechanics. To extend this approach to the nonlinear regime explored in this study, let us first consider the dynamics of the system in the P3 regime. In this regime, the local SR load will acquire, on average, a sequence of 3 different values which we will denote here by $x1,x2$, and $x3$. For a given set of parameters, we can have three distinct beat-to-beat sequences that we can label as

These sequences differ by a shift of one beat and are therefore dynamically equivalent. Furthermore, diffusive coupling between neighboring units will favor alignment of these sequences, while fluctuations of the SR load will tend to prevent synchronization by randomly introducing a shift in the sequence. In the case of P2 dynamics, the alignment of alternating sequences [Eqs. (1) and (2)] can be viewed as an order-disorder transition, where thermal fluctuations take the role of the beat-to-beat stochasticity, and where decreasing temperature is equivalent to increasing the diffusive coupling between sites. Similarly, in the case of P3 dynamics, we expect that there is a critical coupling strength after which sequences will tend to align in order to form a global P3 response. This expectation is confirmed by our numerical simulations [Fig. 6(c)] which show that indeed the Pearson correlation length increases rapidly near the P1 to P3 transition. To gain a deeper understanding of this transition, it is useful to look for an analogous system in equilibrium statistical mechanics. To proceed in this direction, we can consider a mapping to a discrete system where a site *i* in our lattice can attain 3 possible values represented by a “spin state” $\sigma i$. These spin states will correspond to each dynamically equivalent sequence given by Eqs. (35)–(37). For this 3 state system, we can consider the equilibrium behavior described by the Hamiltonian

where $\delta (\sigma i,\sigma j)=1$, if $\sigma i=\sigma j$, and zero otherwise, and where $\u27e8ij\u27e9$ refers to a summation over nearest neighbors. Here, we will consider an interaction energy $J>0$, so that the system is Ferromagnetic and favors alignment of spin states by introducing an energy penalty for two nearest neighbors to be in different states. Therefore, as the temperature of the system is decreased, we expect a transition to an ordered phase, in direct analogy to a synchronization transition of the 3 distinct beat-to-beat sequences at a critical diffusive coupling strength. To proceed with this analogy, we first note that the Hamiltonian given by Eq. (38) is well known in the literature and is referred to as the Potts model.^{20} In general, the Potts model can be formulated in terms of a discrete number of spin states, denoted as *q*, so that the specific model given by Eq. (38) is the *q* = 3 Potts model, while *q* = 2 gives the standard Ising model. The mean field theory of this model has been studied in detail and yields key insights into the nature of the phase transition in this system.^{20} In particular, the model exhibits a finite temperature phase transition to an ordered phase for all *q *> 1. On a square lattice, the critical temperature for *q* = 2 is $Tc(2)=2J/k$, where *k* is the Boltzmann constant, and $Tc(q)=(2J/k)[(q\u22122)/(q\u22121)log\u2061(q\u22121)]$ for *q *> 2. Evaluation of this expression reveals that the transition temperature occurs at lower temperatures as the number of spin components *q* is increased. This result can be understood intuitively from the fact that the configurational entropy of the system increases with *q* due to the increased multiplicity of states. Thus, a lower temperature is required to induce the phase transition from the disordered to the ordered state. In the context of our stochastic coupled map, the tendency of neighboring sites to align is controlled by the diffusive length scale *l*, so that increasing *l* is analogous to decreasing temperature in the statistical mechanics picture. This result explains why a larger *l* is required to observe the P3 phase, at a fixed pacing rate, than is required to observe the P2 phase. In effect, a larger diffusive coupling is required to overcome the increased entropy of the 3 state system. Now, in the case where *l* is fixed and where the pacing period is decreased, then the P2 to P1 transition occurs because of the increased complexity of the underlying dynamics. In the statistical mechanics picture, decreasing the pacing period is equivalent to increasing the number of spin states *q* at a fixed temperature. In this case, there will be a transition from the ordered to disordered phase since the critical temperature decreases with increasing *q*. This analogy suggests that the P2 to P1 transition that we observe at fast rates is effectively an entropic transition that is driven by the increased complexity of the underlying dynamics.

Here, we mention that when the P2 to P1 transition occurs, then the system coarsens into synchronized regions of length scale $\xi $ which is much smaller than the system size, but larger than the lattice spacing. When these patterns form, we find that the local SR load can vary substantially from beat-to-beat, since the dynamics of a single unit on the lattice becomes highly irregular. In the deterministic limit, the dynamics in this regime likely exhibits a higher order periodicity consistent with the period doubling cascade to chaos. This behavior is in sharp contrast to the P1 to P2 transition at slow pacing rates, where the P1 state above the transition point is characterized by a narrow distribution of the SR load.

To carry the analogy to statistical mechanics, further, we note that the mean field approximation to the Potts model predicts that the order-disorder transition occurs via a discontinuous first order phase transition for *q *> 2.^{20} This is in sharp contrast to the second order continuous phase transition from P1 to P2 at slow pacing rates. This result is consistent with our numerical findings which indicate that the onset of the P3 state is discontinuous and sensitive to initial conditions. In effect, the system is bistable and transitions between the ordered and disorder states occur in an abrupt and discontinuous manner. To date, the role of first order phase transitions in the cardiac system has not been studied. Our findings indicate that these transitions can occur when the underlying dynamics for local Ca signaling is unstable to P3 dynamics or higher. The implications of this finding to the onset of cardiac arrhythmias deserves further investigation. In particular, we point out that a first order phase transition implies bistability. Thus, a P3 response can coexist with a P1 response, and discontinuous transitions can occur between these two dynamical states. In cardiac tissue, this feature is likely highly arrhythmogenic since different parts of tissue can attain different dynamical states and therefore induce abrupt spatial variations of APD, which can form a substrate to induce wave break or triggered activity.

The analogy to statistical mechanics also sheds light on the spatiotemporal dynamics of dyadic junctions in the presence of voltage coupling. Our numerical simulations reveal that turning on a positive voltage coupling $(\gamma >0)$ favors an approximate P2 response (alternans) for a wide range of parameters. Hence, we will first consider the spatiotemporal dynamics of the P2 regime in the presence of voltage coupling. In this case, we can assign explicit spin states $\sigma i=+1$ and $\sigma i=\u22121$ which will corespond to the two possible alternans phase [see Eq. (1)]. Here, we claim that the steady state of the coupled map system can be understood by analogy to an equilibrium statistical system with Hamiltonian

where

and

Here, we will take *h *> 0 so that a global alignment of spin states will be favored by $Hglobal$. This global coupling originates from the fact that the APD alternans phase is determined by the summation of Ca release over all dyads in the system. The summed effect on the APD then feeds back on the local Ca alternans, via the coupling mediated by the term $PV(D)$, in order to drive the local response in phase with the global signal. Thus, $Hglobal$ captures the key interactions between local Ca and global voltage. Note here that this model is distinct from the standard Ising model in an external magnetic field, in which the global alternans phase will tend to align with the fixed direction of the external field. Here, the system can evolve toward either a global spin up or spin down state depending on initial conditions. More precisely, initial conditions will pick a dominant phase, which will tilt the Hamiltonian and proceed to drive more spin variables in that direction. This result is consistent with our numerical simulations, which showed that domain walls are driven out of the system after the Ising symmetry is broken by the choice of initial conditions. Here, we point out that the Hamiltonian in Eq. (40) is symmetric under spin reversal, so that the ground state of the system consists of all spin states to be either $\sigma i=+1$ or $\sigma i=\u22121$. This symmetry has important implications in the tissue setting. In particular, given an array of thousands of electrically coupled cells, we expect that regions that are separated by distances longer than the electrotonic length are effectively independent. Thus, these regions will evolve toward the local dominant phase which is picked by the initial conditions in that region. Therefore, these regions will evolve toward a global phase $\sigma i=\xb11$ with equal probability, so that spatially discordant alternans will form over a length scale consistent with the electrotonic length. This mechanism has been pointed out previously by Sato *et al*.^{33} as a possible mechanism for spatially discordant alternans in cardiac tissue. Here, we point out that from an energetic point of view, these spatially discordant alternans patterns are a direct consequence of the Ising symmetry of the coupled voltage-Ca system.

In the case of negative coupling, the Hamiltonian describing the steady state of our coupled map lattice is given by Eq. (41) but with the requirement that $h<0$. The rationale for this choice is that when $\gamma <0$, then APD alternans, which are driven by the summation of the local Ca response, will alternate out-of-phase with the summed Ca signal. Therefore, the resulting APD alternans will drive local junctions with a phase that is opposite to the phase that induced the global APD alternans phase. This negative feedback between voltage and Ca is captured by the Hamiltonian given by Eq. (41) since $Hglobal$ introduces an energy penalty to form a global alignment of spins. Thus, the global coupling term will make it more energetically favorable for the system to form domain walls. On the other hand, there is also an energy penalty to form domain walls due to the contribution from $Hising$. Thus, the minimum energy configuration of the system will consist of the shortest domain wall such that $Hglobal=0$. This argument explains Fig. 10(a), where a jagged domain wall evolved toward two straight lines which divided the square lattice into two equal areas of opposite phase so that the alternans amplitude at steady state was zero. To proceed with this analogy, we can predict further that on a rectangular lattice, the domain wall will orient itself so as to divide the system into half, while joining the closest distance between two sides of the system. In Fig. 11, we demonstrate this effect by considering a rectangular geometry, with aspect ratio 1:6, in which case the domain walls, which are initially oriented along the long axis, always reorient in order to join the minimum distance between two boundaries. Now, if a small circular domain wall is formed due to initial conditions, then that droplet radius will increase in order to minimize $Hglobal$. However, this reduction in $Hglobal$ due to the size increase comes at the expense of the growth of the domain wall boundary and an increase in $Hising$ in the form of surface tension. In effect, the growth of the droplet will stop when the surface tension energy balances the total energy due to $Hglobal$. However, when this balance occurs, then the system will exhibit a dominant alternans phase such that the total $Hglobal$ will be equal to the total surface tension in the system. This argument explains Fig. 10(c), where the square initial condition always evolved to a finite steady state alternans pattern with a circular domain wall. This result reveals that the case of negative coupling can evolve to a multitude of metastable states, depending on the initial conditions of the system. This explains why we found a hierarchy of complex shapes for both random and spatially uniform initial conditions.

### B. Comparison to the previous work

In a previous study, Shiferaw and Karma^{16} have analyzed the spatiotemporal dynamics of Ca when voltage and Ca are bidirectionally coupled. In that model, each sarcomere in the cell was described by a deterministic ionic model and diffusively coupled to its nearest neighbor sarcomere. In that study, they focused on the negative coupling regime and identified a Turing-like instability, which induced the formation of subcellular nodes. By considering the dynamics of small amplitude voltage and Ca alternans, they cast the pattern formation process as an activator-inhibitor partial differential equation (PDE), with Ca alternans playing the role of a short range activator and APD alternans the role of a long range inhibitor. Here, we point out that this activator-inhibitor picture is consistent with our Hamiltonian given by Eq. (39). In that case, the low energy configurations are precisely those with a single domain wall that reduces the energy penalty due to $Hglobal$. Thus, the equilibrium configurations generated by our Hamiltonian should be consistent with the patterns that form using the PDE describing the motion of small amplitude alternans. However, a key difference between these approaches is that the PDE picture does not include local stochasticity and therefore does not capture the phenomenology associated with phase transitions. We also note here that the patterns shown in Figs. 9(b) and 9(c) are formed with both spatially uniform initial conditions and random initial conditions, respectively. In this case, the system is initialized in the fully nonlinear regime, where the local alternans amplitude needed is not small. Here, the negative coupling between Ca and voltage forces different regions in the cell to alternate out-of-phase via a mechanism analogous to the activator-inhibitor interaction that drives the Turing instability. In the original work of Shiferaw and Karma, the spatial patterns developed from a zero alternans state, where small noise perturbations were amplified by the pattern forming instability. Hence, the steady state patterns attained in that study were obtained using different initial conditions. Using our mapping model, we can also explore the development of spatially discordant patterns from the zero amplitude state. To observe the growth of patterns from this initial state, it is necessary to first pace the system to steady state in the P1 regime, and then to change the pacing rate to the P2 regime so that alternans develops gradually from a spatially homogenous P1 state. Indeed, using this pacing protocol, we find that spatially discordant alternans develop in much the same way as that predicted in the Shiferaw-Karma study. Finally, we mention here that the emergence of multiple steady state patterns, in the negative coupling scenario, has been observed by Zhao^{34} in cable simulations. In that study, it was found that a rapidly paced cable of cells could evolve to a multitude of steady state patterns depending upon initial conditions and pacing history. This result is similar to the multiple steady states shown in Fig. 9. In both cases, it is the negative coupling between Ca and voltage that drives the formation of spatially discordant regions in a manner that is sensitive to the pacing history and initial conditions of the system.

In a later study, Restrepo and Karma^{9} applied a physiologically detailed three dimensional (3D) model to analyze the movement of domain walls. They showed that in the regime of positive coupling, nodes were expelled from the system, while in the case of negative coupling, nodes drifted to the center. These results are consistent with the findings here, which reproduced essentially the same phenomenology for the different coupling signs. Thus, the main spatiotemporal processes observed in detailed computational models are captured using our stochastic coupled map model and the corresponding statistical mechanics analogy.

In a relevant study, Qu *et al*.^{15} applied a detailed stochastic model of subcellular Ca cycling to show that the Ca transient can exhibit higher order periodicities during rapid pacing. As a function of increasing pacing rate, they observe a transition from P2 to P4, which is followed by a clear transition to P1 and P3 [see Fig. 7(e) in Qu *et al*.^{15}]. This result is consistent with the current findings in the case where $l=8$ shown in Fig. 2(c). Thus, the transition to P1 at fast pacing rates, observed in their study, is likely the entropic transition identified in this work. This result confirms the universality of these transitions that should also occur in physiologically detailed models of subcellular Ca. Furthermore, their finding suggests that in the detailed computational model, and perhaps in real myocytes, the spatial coupling between release units is much larger than nearest neighbor. This is likely because of the rapid diffusion of Ca in the SR which may serve to strongly couple release units. In the future, it will be interesting to explore further how Ca diffusion in both the SR and cytosol can influence the spatial organization of Ca release in cardiac myocytes.

## V. SUMMARY AND CONCLUSIONS

In this study, we have explored the dynamics of a stochastic coupled map system that mimics key features of Ca signaling in cardiac cells. The key feature of the model is that Ca release in the SR is modeled as discrete release events which occur with a probability that is a nonlinear function of the SR load. The combination of nonlinearity and stochasticity leads to complex spatiotemporal dynamics. The main new finding is that we identify a novel transition, at fast pacing rates, from a global P2 response to P1 response. We argue further that this transition can be understood using an analogy to the phase transitions described by the *q* state Potts model. This perspective indicates that this transition is driven by the increase in entropy due to the increasing complexity of the underlying dynamics at fast pacing rates. Our mapping to the Potts model also predicts that various transitions that occur on our coupled map lattice can be interpreted as first order phase transitions. This observation may have important implications in the onset of cardiac arrhythmias since these transitions are discontinuous with large changes in system parameters in response to very small parameter changes. This high sensitivity to parameters should lead to pronounced heterogeneities in cardiac tissue and are likely highly arrhythmogenic. Finally, we have also investigated the coupling between voltage and Ca to reveal rich spatiotemporal dynamics depending on the sign of the coupling. Our main finding is that positive coupling tends to synchronize subcellular dynamics via the gradual elimination of domain walls. In the tissue setting, this result implies that subcellular alternans will tend to be locally synchronized and can only be desynchronized on a scale that is larger than the electrotonic length. On the other hand, negative coupling drives the formation of subcellular domain walls, and we expect in this case that APD alternans amplitude on the tissue scale will be small. However, in both cases, it is clear that the detailed implications of the complex spatiotemporal patterns that emerge requires further study.

## ACKNOWLEDGMENTS

This work was supported by the National Heart, Lung, and Blood Institute (Grant No. RO1HL119095) (Y.S.). Part of this work has been completed at the Kavli Institute of Theoretical Physics (KITP). This research was supported in part by the National Science Foundation (NSF) (Grant No. NSF PHY-1748958), by the National Institutes of Health (NIH) (Grant No. R25GM067110), and by the Gordon and Betty Moore Foundation (Grant No. 2919.01). E.A.-L. acknowledges the financial support from Fundació La Marató de TV3 and from the Spanish Ministerio de Economía y Competitividad (MINECO) under Grant Nos. SAF2014-58286-C2-2-R and SAF2017-88019-C3-2-R with the support from European FEDER funds.