New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Turing instability mediated by voltage and calcium diffusion in paced cardiac cells

Edited by Charles F. Stevens, The Salk Institute for Biological Studies, La Jolla, CA, and approved February 8, 2006 (received for review December 22, 2005)
Abstract
In cardiac cells, the coupling between the voltage across the cell membrane (V _{m}) and the release of calcium (Ca) from intracellular stores is a crucial ingredient of heart function. Under abnormal conditions and/or rapid pacing, both the action potential duration and the peak Ca concentration in the cell can exhibit well known perioddoubling oscillations referred to as “alternans,” which have been linked to sudden cardiac death. Fast diffusion of V _{m} keeps action potential duration alternans spatially synchronized over the ≈150μmlength scale of a cell, but slow diffusion of Ca ions allows Ca alternans within a cell to become spatially asynchronous, as observed in some experiments. This finding raises the question: When are Ca alternans spatially inphase or outofphase on subcellular length scales? This question is investigated by using a spatially distributed model of Ca cycling coupled to V _{m}. Our main finding is the existence of a Turingtype symmetry breaking instability mediated by V _{m} and Ca diffusion that causes Ca alternans to become spontaneously outofphase at opposite ends of a cardiac cell. Pattern formation is governed by the interplay of shortrange activation of Ca alternans, because of a dynamical instability of Ca cycling, and longrange inhibition of Ca alternans by V _{m} alternans through Casensitive membrane ionic currents. These results provide a striking example of a Turing instability in a biological context where the morphogens can be clearly identified, as well as a potential link between dynamical instability on subcellular scales and lifethreatening cardiac disorders.
When a cardiac cell is stimulated by a change in the membrane, V _{m}, calcium (Ca) is released from intracellular stores and signals cell contraction (1, 2). This release of Ca occurs within several thousand submicrometerscale junctions distributed throughout the volume of the cell. Within these junctions localized clusters of Casensitive ion channels called Ryanodine receptors, which gate the flow of Ca across intracellular stores, open in response to Ca entry via Ltype Ca channels in the surface cell membrane (1). Signaling within these junctions occurs through a Cainduced–Carelease process (2, 3), which endows the system with rich spatiotemporal properties typical of excitable media.
Rapid pacing rates, or abnormal physiological conditions, are known to produce complex dynamics of both V _{m} and Ca at the cell and tissue level. The most common nontrivial dynamics is the phenomenon of alternans, where the action potential duration (APD) and the Ca transient alternate from one beat to the next (4, 5). Alternans is clinically important because it has been shown to correlate with the onset of lifethreatening arrhythmias and sudden cardiac death (6, 7). Historically, alternans has been attributed to a perioddoubling instability that depends on the electrical restitution properties of cardiac cells (4, 8, 9), which itself is governed by the kinetics of membrane ion currents. More recent studies, however, have demonstrated that alternans also can result from an instability of intracellular Ca cycling (5, 10–12).
Recently, Kockskamper and Blatter (13) used twodimensional (2D) linescan confocal microscopy of cat atrial cells to image the spatial distribution of Ca release when the cell was paced during alternans. They observed that the release of Ca could be spatially outofphase, where one half of the cell alternated with opposite phase with the other half. Outofphase subcellular Ca alternans also has been observed by Diaz et al. (14) in rat ventricular cells. These experimental results demonstrate that Ca release within cells can form complex spatiotemporal patterns. However, what governs the formation and dynamics of these patterns remains unknown.
To address this question, we explore the spatiotemporal dynamics of Ca and V _{m} alternans using numerical simulations of a physiologically based ionic model. A distinguishing feature of this model, in comparison with previous models (10, 15, 16), is that it captures the spatially distributed nature of the Ca cycling system, along with a description of the membrane ion currents that couple to this system. Hence, this model can be used to investigate pattern formation of Ca signaling on subcellular scales. The simulations reveal that spatially distant regions of Ca alternans can either synchronize or spontaneously desynchronize depending on the bidirectional coupling between the intracellular Ca and V _{m}. This bidirectional coupling is mediated by Casensitive membrane ionic currents that gate V _{m}. Because Ca mediates the contraction of the cardiac cell, these asynchronous subcellular patterns will influence contraction on the tissue scale and thus may potentially affect the function of the heart.
To shed light on this pattern formation process, we reduce the multivariable ionic model to coupled reaction–diffusion equations that govern the spatiotemporal evolution of the amplitudes of Ca and V _{m} alternans. These amplitude equations describe the spatiotemporal dynamics of Ca and APD alternans at the subcellular scale. An analysis of these equations reveals that the spatial desynchronization of Ca alternans results from a symmetrybreaking Turinglike instability (17–19) mediated by V _{m} and Ca diffusion. The amplitudes of alternans are shown to play the role of diffusing morphogens, with Ca corresponding to the shortrange activator and V _{m} to the longrange inhibitor. The crucial requirement for a Turing instability, in which the inhibitor diffuses significantly faster than the activator, is ensured by the fact that V _{m} diffuses across the cell membrane ≈10^{4} times faster than intracellular Ca in the myoplasm.
Physiologically Detailed Model
Cell Architecture.
The architecture of a cardiac ventricular cell is shown schematically in Fig. 1 A. The cell membrane forms an array of deep invaginations into the cell (ttubules), distributed along planes (zlines) transverse to the long axis of the cell (20). The ttubules distribute dyadic junctions, where Ca signaling takes place, uniformly over the cell volume. The volume in between zlines is referred to as the sarcomere and contains the sarcoplasmic reticulum (SR). The SR is a diffuse tubular network that sequesters a high concentration of Ca. A typical ventricular myocyte has a length of 150 μm and has ≈75–100 sarcomeres.
We take advantage of the observed structural periodicity in the longitudinal direction to treat each sarcomere as an equivalent subcellular unit. This method lets us represent the cell as a 1D chain of identical subunits, where each subunit is assumed to have an identical distribution of dyadic junctions, and to enclose a roughly equivalent volume of the SR network. We assume further that Ca diffusion is sufficiently fast on the scale of one sarcomere, so that the concentration of Ca can be accurately described by the local spatial average. Note that atrial cells lack ttubules, and consequently Ca release propagates as a wave from the cell membrane toward the central axis of the cell. We expect, however, the present model to also be applicable to these cells as long as this propagation time is shorter than the pacing period, in which case the cell can still be divided longitudinally into equivalent subunits.
V_{m} Dynamics.
The effective diffusion coefficient of V _{m} (D _{V} ≈ 2.5 × 10^{−4} cm^{2}/ms) (2) is large. Therefore, V _{m} equilibrates rapidly over the length scale of a cell on a time scale of ≈0.5 ms that is comparable with the rise time of the action potential (1 ms) but much smaller than the APD. Hence, we can assume that each sarcomere senses the same V _{m}, which is dictated by the ionic currents averaged over the whole cell. The membrane ionic currents, which govern the dynamics of V _{m} across the membrane surrounding a single sarcomere (here the kth subunit), are shown schematically in Fig. 1 B. Here, I _{Na} ^{k} is the fast inward Na^{+} current, I _{K} ^{k} is the total outward K current, and I _{Ca} ^{k} and I _{NaCa} ^{k} are the Cadependent Ltype and NaCa exchange currents, respectively. The ionic currents are formulated by using the canine ventricular cell model of Fox et al. (15), which is based on Luo–Rudy currents (16). The spatially uniform V _{m} across the cell membrane satisfies the standard equation of chargeconservation
where C _{m} is the cell membrane capacitance, and N is the number of sarcomere units. Because the V _{m}dependent Na and K currents do not couple to the intracellular Ca transient, they are spatially uniform along the cell. Therefore, the sums over I _{Na} ^{k} and I _{K} ^{k} in Eq. 1 can be replaced by the whole cell current. In contrast, Cadependent currents are local and vary between sarcomeres. I _{stim} is a stimulus current applied to the whole cell simultaneously, consistent with the fact that V _{m} is assumed to be spatially uniform. A detailed formulation of these ionic currents is in Supporting Text, which is published as supporting information on the PNAS web site.
Ca Cycling Dynamics.
The Ca cycling dynamics within each sarcomere is described by using the model of Shiferaw et al. (10). The essential elements of the local Ca cycling machinery are shown in Fig. 1 B. Ca entry into the sarcomere occurs via Ltype Ca channels (I _{Ca} ^{k}), which triggers a Cainduced–Ca release flux (J _{rel} ^{k}) from Ryanodine receptor channels bound to the SR membrane. The amount of Ca released into the sarcomere depends sensitively, in a nonlinear fashion, on the concentration in the SR (the SRload c _{j} ^{k}). Once the intracellular Ca concentration (c _{i} ^{k}) rises, powerful uptake (J _{up} ^{k}) pumps on the SR membrane are activated and pump Ca back into the SR. The Ca that enters the cell via I _{Ca} ^{k} is extruded out of the cell via the electrogenic NaCa exchange current (I _{NaCa} ^{k}), so that at steady state, Ca entry balances Ca extrusion. The Ca systems in neighboring sarcomeres are coupled by allowing diffusion between the Ca in myoplasm (c _{i} ^{k}) and SR (c _{j} ^{k}). Ca diffusion (D _{Ca} = 150–300 μm^{2}/s) (2) is several orders of magnitude slower than V _{m} diffusion driven by voltage gradients along the cell. Full details of the Ca cycling model are described in Supporting Text.
Bidirectional Coupling.
The coupling of V _{m} to the local Ca transient within a sarcomere is determined by the experimentally established phenomenon of graded SR Ca release (1, 21), where the amount of SR Ca released increases in proportion to the local Ca entry (I _{Ca} ^{k}). The availability of Ltype Ca channels at a given beat depends on the previous diastolic interval, which is the time spent at full repolarization before the next action potential upstroke. Now, a larger diastolic interval gives more time for recovery of Ltype Ca channels at the resting membrane potential. Thus, in our sarcomere model, graded release requires that the peak of the local Ca transient increases in response to an increase of diastolic interval in the previous beat, as illustrated in Fig. 2 A. This coupling will be referred to hereafter as “positive V _{m} → Ca coupling.” In principle, it is also possible to have the opposite effect at slow pacing rate, also illustrated in Fig. 2 A (“negative V _{m} → Ca coupling”). However, we did not investigate this case in the present study.
The local Ca released into a sarcomere also couples to the global APD by its effect on the local membrane ionic currents. As discussed in ref. 22, there are two cases that need to be considered. The first, referred to as “positive Ca → V _{m} coupling,” illustrated in Fig. 2 B, corresponds to the case when a large local Ca release tends to prolong the APD. The second, referred to as “negative Ca → V _{m} coupling” (Fig. 2 B), corresponds to the case when a large local Ca release shortens the APD. The sign of the coupling is dictated by the relative contributions of the local Ltype Ca current (I _{Ca} ^{k}) and the NaCa exchange current (I _{NaCa} ^{k}) to the APD. A larger Ca release tends to inactivate the local I _{Ca} ^{k} more rapidly through Cainduced inactivation, which tends to shorten the APD. However, a large local Ca release also increases I _{NaCa} ^{k}, which tends to prolong APD. For the ionic model used here, the sign of the coupling was controlled by adjusting the degree of Cainduced inactivation of the I _{Ca} ^{k}, as described in Supporting Text.
Alternans.
Steadystate alternans of the Ca transient within a given sarcomere can be induced in the model by increasing the steepness, at high SR loads, of the function relating SR Ca release to SR Ca load, as described theoretically (10, 11) and experimentally (12). In experiments under physiologically normal conditions (5), and in the present ionic model, this bifurcation to alternans occurs at rapid pacing rates (T < T _{c}) when the cell is overloaded with Ca. For low pacing rates T > T _{c}, the steadystate Ca transient is periodic from beattobeat.
Alternans also can be induced by a purely V _{m}dependent mechanism, independently of Ca cycling. This mechanism depends on the kinetics of membrane ion currents, which regulate V _{m}, and are typically described using the restitution relationship between APD and diastolic interval (4, 8). However, in our study we found that alternans due to the latter mechanism always yielded spatially synchronized Ca alternans. The reason for this result is that when the spatially uniform APD alternates, it simply drives local Ca alternans in phase over the extent of the cell by means of the graded coupling described in the previous section. Hence, in this study, we have focused exclusively on alternans due to unstable Ca cycling, which we found leads to nontrivial spatial patterns.
The bidirectional coupling between the APD and the local Ca release determines the relative phase of APD and Ca transient alternans at the subcellular level (23–25). For alternans due to unstable Ca cycling, positive Ca → V _{m} coupling produces steadystate electromechanically concordant alternans, where a large–small–large local Ca transient corresponds to long–short–long APD, whereas negative Ca → V _{m} leads to steadystate electromechanically discordant alternans, where a small–large–small peak Ca transient corresponds to a long–short–long APD. Both of these modes of alternans have been observed experimentally under various physiological conditions (23–25).
Numerical Results
We have simulated a set of N = 75 coupled sarcomeres with a spacing of Δ = 2 μm between Zplanes and a cell length L = NΔ = 150 μm. To characterize the spatiotemporal evolution of Ca alternans, we define the local amplitude of Ca alternans as Δc(x, n) = (−1)^{n}(c _{n+1} ^{k} − c _{n} ^{k}), where c _{n} ^{k} is the peak Ca transient in sarcomere k at the nth paced beat, and where x = kΔ denotes the position along the cell. The factor (−1)^{n} is included so that the phase of alternans does not change sign at every beat.
The parameters of the isolated sarcomere model are chosen such that, at a pacing rate T < T _{c}, alternans due to unstable Ca cycling develops and saturates at a steadystate value. The initial Ca concentration in the SR was chosen to be ≈10 μM less than its steadystate value, so that fully developed alternans emerged only after ≈100 paced beats, i.e., the SR had to load with Ca before it became dynamically unstable. Thus, for the first few beats, alternans amplitude of Ca and APD was always zero for all sarcomeres in the cell. To mimic realistic conditions, the initial conditions of Ca concentration in the myoplasm and the SR were taken from a Gaussian distribution with a standard deviation <1% of their average value.
Our main finding is that the final steadystate pattern of Ca alternans depends sensitively on the bidirectional coupling between Ca and V _{m}. When model parameters were chosen such that the Ca → V _{m} coupling was positive, steadystate Ca alternans were always inphase over the extent of the cell. In Fig. 3 A, we plot the amplitude of Ca alternans in a cell where the degree of Cainduced inactivation is adjusted so that the Ca → V _{m} coupling was positive. The spatial distribution of the amplitude of Ca alternans at the nth beat, Δc(x, n), started at a very small value across the cell and grew as the cell was paced at T = 0.35 s to a fully developed spatially synchronized pattern. To assess the nature of the coupling for these model parameters, we paced an isolated sarcomere at the same rate and found that alternans were indeed electromechanically concordant. In contrast, when the Ca → V _{m} coupling was negative, spatially outofphase alternans developed from an initially small amplitude spatial distribution, as shown in Fig. 3 B. In this case, the final steadystate pattern was such that the two halves of the cell alternated out of phase, and Δc(x, n) changed sign across a node in the center of the cell. Remarkably, a qualitatively similar pattern was observed in the experiments of Kockskamper and Blatter (13). To assess the nature of the coupling for these parameters, we paced again an isolated sarcomere and found that Ca and APD alternans were electromechanically discordant.
To further analyze the patternformation mechanism that produces the asynchronous state for negative Ca → V
_{m} coupling, we carried out numerically a linear stability analysis of the spatially homogeneous state with no alternans. We computed the amplification rate Ω(q) of small sinusoidal perturbations of wavelength λ = 2π/q in a large system of N = 400 coupled sarcomeres. For each mode the amplification rate was obtained from the slope of a linearlog plot of the average magnitude of alternans, defined as
Morphogenetic Roles of Ca and V_{m} Alternans
A basic understanding of this pattern forming instability can be gained by considering the similarity between the amplitudes of Ca and APD alternans and diffusing morphogens in a Turinglike instability. It is useful to discuss this analogy qualitatively to set the stage for amplitude equations in the next section, which provide the formal mathematical framework to make this analogy precise.
Activator.
The fact that Ca diffuses much slower than V _{m} makes the amplitude of Ca alternans the natural candidate for shortrange activator. This role can only be fulfilled, however, if the alternans within an isolated sarcomere are due to an intrinsic instability of Ca cycling rather than V _{m} dynamics. In this case, Ca transient alternans will grow locally as the cell is paced and induce APD alternans. In Fig. 5 A, we show the beattobeat Ca transient and V _{m} during alternans for an isolated sarcomere that exhibits discordant electromechanical alternans. Next, consider a small fluctuation that leads to an increase in amplitude of Ca alternans, as indicated by the red trace. The subsequent feedback of Ca on V _{m} will tend to increase the amplitude of APD alternans, as long as the subsequent changes in APD do not feedback and cancel the effect of the small fluctuation on the next beat. This condition will be made more mathematically precise below but is essentially equivalent to saying that the feedback of Ca on itself should be stronger than the feedback of V _{m} on itself. If this condition is met, then the amplitude of Ca alternans can be viewed as the activator of APD alternans.
Inhibitor.
The extremely fast diffusion of V _{m} in the cell membrane makes the amplitude of APD alternans the natural candidate for the longrange inhibitor. To study the effect of APD alternans on Ca transient alternans, consider a sarcomere in which alternans are electromechanically discordant, as illustrated by the black lines in Fig. 5 B. The mathematical condition for this dynamical mode is that ∂A_{n} /∂c_{n} < 0, where A_{n} and c_{n} denote the APD and peak Ca transient at beat n. This condition, which is equivalent to negative Ca → V _{m} coupling, ensures that an increase in peak Ca transient c_{n} tends to decrease the APD A_{n} , and thus a large Ca transient will correspond to a small APD and vice versa. A more rigorous mathematical derivation of this condition using a nonlinear map reduction is given in Supporting Textand Fig. 6, which is published as supporting information on the PNAS web site. Now consider a fluctuation that increases the amplitude of APD alternans, as illustrated by the red V _{m} trace in Fig. 5 B. Such a fluctuation will cause the Ca transient peak in the next beat to be smaller if ∂c_{n} /∂A _{n−1} < 0. This condition is equivalent to the requirement for positive V _{m} → Ca coupling as illustrated in Fig. 2 A and is precisely the coupling dictated by the phenomenon of graded release incorporated in our detailed ionic model. The later decrease in the peak Ca transient is equivalent to a decrease in Ca alternans amplitude because the larger Ca peak (red line) is smaller than it would have been in the absence of the fluctuation (black line). Combining the above two conditions, we can conclude that for APD alternans to inhibit Ca transient alternans the condition (∂A_{n} /∂c_{n} )(∂c_{n} /∂A _{n−1}) > 0 must be satisfied. This condition also can be satisfied when alternans are electromechanically concordant if the sign of both partial derivatives are reversed, although this case was not explored here.
After identification of the amplitudes of Ca and APD alternans with diffusing morphogens, understanding the formation of the asynchronous state follows straightforwardly from the analogy with a classic Turing instability mechanism. Namely, a fluctuation that causes the amplitude of Ca alternans to increase in a small region of the cell generates APD alternans that in turn inhibit the growth of Ca alternans far from that localized fluctuation (longrange inhibition). Thus, Ca alternans will not synchronize spatially under these conditions and will develop nontrivial spatiotemporal patterns.
Amplitude Equations
A more quantitative analytical understanding of the Turing instability can be obtained by deriving equations that govern the spatiotemporal evolution of the amplitudes of Ca and APD alternans (26). The derivation, which is detailed in Supporting Text, exploits the fact that the local beattobeat dynamics in the detailed ionic model can be described by an iterative map, and that the alternans amplitudes are small and evolve slowly (over many beats) close to the alternans bifurcation.
The starting point of the derivation is the 2D iterative map A_{n} = F _{1}(A _{n−1}, c _{n}), and c _{n} = F _{2}(A _{n−1}, c _{n−1}), which relates the APD and peak Ca transient at one beat to the next. Such a map was recently shown to describe well the beattobeat dynamics of V _{m} and Ca in a situation where the APD and Ca transient are assumed to be spatially uniform (22). Here, we allow these quantities to vary spatially and label each sarcomere by its position x along the cell. We couple the maps spatially using nonlocal spatial kernels that describe the cumulative effect of V _{m} and Ca diffusion in the interval of one beat. For V _{m}, this coupling amounts to simply taking A_{n} equal to the spatial average of A_{n} = F _{1}(A _{n−1}, c_{n} (x)) owing to the fact that V _{m} diffusion is essentially instantaneous. For Ca, we take c_{n} (x) to be the weighted average of F _{2}(A _{n−1}, c _{n−1}(x′)) with a weight that decreases smoothly to zero with increasing x − x′ over a length scale ξ ∼ (D _{Ca} × T)^{1/2}, where D _{Ca} is the average diffusion constant of free Ca inside the cell, and T is the pacing period. Next, we make the substitutions A_{n} = A* + (−1)^{n} a(n) and c_{n} (x) = c* + (−1)^{n} c(x, n) in the spatially coupled maps, where A* and c* are the fixed point of the map, and a(n) and c(x, n) are the amplitudes of APD and Ca alternans, respectively. We use the fact that the amplitudes of alternans vary slowly from beat to beat close to the alternans bifurcation to treat the beat number n, or equivalently the time t = nT as a continuous variable, which allows us to write A _{n+2} − A_{n} ≈ 2T(−1)^{n}∂a/∂t, and similarly c _{n+2} − c_{n} ≈ 2T(−1)^{n}∂c/∂t. In addition, we assume that the amplitude of Ca alternans varies along the cell on a scale larger than the spacing between sarcomeres, such that x can be treated as a continuous variable. With these simplifications, we get where c̄ = (1/L)∫_{0} ^{L} c(x, t)dx denotes the spatial average of the Ca alternans amplitude over the whole cell. For clarity, we have omitted nonlinear terms that are not essential to characterize the linear regime of the instability. The coefficients of the amplitude equations, which are derived in detail in Supporting Text, can be expressed in terms of the matrix elements of the Jacobian matrix of the twovariable map defined by J _{11} = ∂F _{1}/∂A_{n} , J _{12} = ∂F _{1}/∂c_{n} , J _{21} = ∂F _{2}/∂A_{n} , and J _{22} = ∂F _{2}/∂c_{n} , evaluated at the fixed point of the map A* and c*. The effective diffusion coefficient of Ca alternans amplitude is given by D = J _{22} ^{2}ξ^{2}/2T.
It is important to emphasize that the reaction–diffusion Eqs. 2 and 3 have a richer structure than the standard Turing model (17, 19) because of the fact that the diffusing morphogens measure the amplitude of perioddoubling oscillations of two signals (Ca and V _{m}), rather than the concentrations of two chemical species. Consequently, the local rate of generation of the activator depends both on its local value and on its spatial average over the cell, as reflected by the appearance of terms proportional to c and c̄ in Eq. 3 , respectively. In contrast, there is no analog of the global coupling term c̄ in the equation for the activator in the limit of the standard Turing model where the inhibitor diffuses quasiinstantaneously on the scale of the system. It is interesting to note that, beyond the Turing analogy, the breakdown of synchronization leading to the formation of phase domains also has been found to occur generically in oscillatory media in the presence of global coupling, as observed in certain surface catalytic reactions (27).
The conditions for a Turing instability can now be readily obtained from a stability analysis of the amplitude equations. Taking a and c proportional to exp(iqx + Ω(q)t), the condition for a nontrivial solution of Eqs. 2
and
3
to exist yields an eigenvalue equation for Ω(q) that is quadratic in q. The most unstable branch is easily found to be
Physical Interpretation.
To interpret the above conditions for a finite wavelength instability, we express the inequality conditions in terms of the Jacobian of the beattobeat map. The first condition, δ > 0, is equivalent to J _{22} < −1, which corresponds physically to the condition that alternans are due to an instability in Ca cycling. The second condition, δ > Ω(0), is shown in Supporting Text to be equivalent to the two conditions J _{22} < J _{11} and J _{12} J _{21} < 0. Because J _{11} and J _{22} give a measure of the degree of instability of the V _{m} and Ca cycling subsystems, the first condition is satisfied if the Ca subsystem is effectively more unstable than the V _{m} subsystem, which makes the condition for Ca to be the activator more precise, as described in the last section. The second condition, J _{12} J _{21} < 0, is in turn equivalent to the condition for V _{m} to be the inhibitor, (∂A_{n} /∂c_{n} )(∂c_{n} /∂A _{n−1}) > 0, derived in the last section by using qualitative arguments. To see this condition, note that J _{21} = ∂c_{n} /∂A _{n−1}, and J _{12} = (∂c_{n} /∂c _{n−1})(∂A _{n}/∂c_{n} ) and that ∂c_{n} /∂c _{n−1} = J _{22} is negative for alternans to occur.
Analysis of the leading eigenvector of the Jacobian reveals that the relative phase of Ca and APD alternans is determined precisely by the conditions above for the finite wavelength instability. Note that condition J _{12} J _{21} < 0 can be satisfied in two ways: J _{12} < 0 and J _{21} > 0 or J _{12} > 0 and J _{21} < 0. If J _{21} < 0, the condition for a finite wavelength instability is identical to that for electromechanically discordant alternans. This result is precisely the case investigated here where APD alternans inhibit Ca alternans, consistent with the qualitative arguments of the last section. If J _{21} > 0, the instability condition requires that alternans are electromechanically concordant. The latter case was not investigated here.
Comparison with Ionic Model.
The physical interpretation of the amplitude equations is consistent with our numerical simulation results. The simulations revealed that the asynchronous state developed after many beats only when we enhanced the effect of Cainduced inactivation. For the same model parameters, an isolated sarcomere paced at the same rate exhibited electromechanically discordant alternans. This dynamical mode emerged when alternans are due to a steep relationship between SR Ca release and SR load and when the Ca → V _{m} coupling is negative, as illustrated in Fig. 2 B. Because, furthermore, the V _{m} → Ca coupling is positive in the ionic model by virtue of the graded release relationship between Ca release from the SR and the Ltype Ca current, which triggers this release, the asynchronous state developed when (V _{m} → Ca) × (Ca → V _{m}) < 0, which is equivalent to the condition J _{12} J _{21} < 0 for a Turing instability obtained from the amplitude equations.
Discussion
The main finding of this work is that heterogeneity of Ca cycling on the subcellular scale can arise spontaneously from a Turingtype instability mediated by Ca and V _{m} diffusion. In the present context, symmetry breaking is produced by local selfenhancement of Ca alternans and longrange inhibition of these alternans by APD alternans. Amplitudes of Ca and APD alternans are predicted to play a directly analogous role to activator and inhibitor species in the context of reaction–diffusion models in which Turing patterns have been studied (18). The crucial requirement that the inhibitor diffuses faster than the activator is fulfilled here because V _{m} diffusion is several orders of magnitude faster than Ca diffusion.
The conditions for this instability can be summarized as follows. (i) Alternans at the subcellular level must be due to an instability of Ca cycling. (ii) Bidirectional coupling derived from the product of the unidirectional couplings between Ca and V _{m} must satisfy the requirement that (V _{m} → Ca) × (Ca → V _{m}) < 0. The latter condition also determines the relative phase of Ca and V _{m} alternans at the single sarcomere level, i.e., whether alternans are electromechanically concordant or discordant in the absence of spatial coupling. Under normal physiological conditions, the graded release relationship between Ltype Ca current and Ca release from the SR requires that the V _{m} → Ca coupling is positive. Thus, in this case, the condition for instability is determined by the sign of the Ca → V _{m} coupling, which is itself governed by the effects of Ltype Ca current and NaCa exchange current on the APD.
It is important to emphasize that the existence of the Turingtype instability elucidated here does not depend on a specific mechanism for the autocatalytic amplification of Ca alternans. Therefore, even though we have assumed here that Ca alternans are due to a steep relationship between SR release and SR load, other mechanisms of instability also could potentially give rise to a Turing instability, as long as the sign of the coupling between Ca on V _{m} is negative.
Our results provide a framework to interpret experimental observations of outofphase patterns of subcellular Ca alternans (13). In particular, they suggest that such patterns may result from a symmetrybreaking instability, rather than from intrinsic spatial heterogeneities of the Ca cycling machinery inside the cell. However, other mechanisms to form asynchronous patterns of Ca alternans are possible. Specifically, spontaneous Ca release in one region of the cell can desynchronize a spatially uniform state of inphase Ca alternans, in an analogous way that a premature electrical stimulus can initiate outofphase electrical alternans in cardiac tissue (27). Therefore, to distinguish between different potential mechanisms, it is essential to experimentally determine the bidirectional coupling between Ca and V _{m}. The possibility of modulating this coupling is suggested by experiments in cat ventricular cells by Rubenstein and Lipsius (25), who observed a transition from concordant to discordant electromechanical alternans by reducing the temperature of the cell from 35°C to 32°C. Because theory predicts (22) that the relative phases of Ca and V _{m} alternans is determined by the sign of this coupling, these experiments show that this sign can be reversed by a temperature change. It therefore may be possible to validate the existence of a Turing instability experimentally by observing how the degree of synchronicity of subcellular patterns of alternans change with temperature and pacing rate.
The Turing instability elucidated here causes a spontaneous breakdown of synchrony of Ca alternans at a subcellular scale. This breakdown will generally have an influence on the dynamics of membrane voltage via Ca sensitive membrane ionic currents and hence on repolarization alternans at a single cell level. Therefore, this instability may have important implications for understanding the onset of lifethreatening cardiac arrhythmias insofar as repolarization alternans at the cellular level suffice to induce wave instabilities and ventricular fibrillation at a tissue scale, as shown in both experimental and modeling studies (7, 8, 23). However, the link between subcellular and tissue scale instabilities remains to be explored. This exploration may be possible by extending the present study to networks of coupled cells.
The instability mechanism proposed by Turing (17) to explain morphogenesis based on the interaction of diffusing morphogens has been demonstrated experimentally in chemical systems that mimic reaction–diffusion models (28–30). Turing instability mechanisms also have been invoked to explain various experimentally observed patterns in biological systems (31, 32). However, the complexity of these systems makes it generally difficult to identify the morphogens that lead to the observed patterns. From this standpoint, the present work highlights the existence of a rare, and somewhat unexpected, example of a naturally occurring Turinglike instability in a biological context where the diffusing morphogens can be clearly identified. An important distinguishing feature of this instability is that one morphogen is a chemical signal (Ca), whereas the other is an electrical signal (V _{m}). Another distinguishing feature is that both morphogens measure the amplitudes of underlying period doubling oscillations of these signals rather than simply their magnitudes. This feature leads to a more complex interaction of the activator and the inhibitor, reflected in the reaction–diffusion model, whose manifestations in a fully nonlinear regime of this instability are yet to be explored.
Acknowledgments
We thank Daisuke Sato for checking some of the simulation results and Jim Weiss for valuable suggestions at the initial stages of this work. This work was supported by National Institutes of Health/National Heart, Lung, and Blood Institute Grants P50 HL52319 and P01 HL078931.
Footnotes
 ^{‡}To whom correspondence should be addressed. Email: yshiferaw{at}mednet.ucla.edu

Author contributions: Y.S. and A.K. designed research, performed research, and wrote the paper.

Conflict of interest statement: No conflicts declared.

This paper was submitted directly (Track II) to the PNAS office.
 Abbreviations:
 APD,
 action potential duration;
 SR,
 sarcoplasmic reticulum.
Abbreviations:
 © 2006 by The National Academy of Sciences of the USA
References
 ↵

↵
 Bers D. M.

↵
 Fabiato A.

↵
 Nolasco J. B. ,
 Dahlen R. W.
 ↵

↵
 Garfinkel A. ,
 Kim Y. H. ,
 Voroshilovsky O. ,
 Qu Z. ,
 Kil J. R. ,
 Lee M. H. ,
 Karagueuzian H. S. ,
 Weiss J. N. ,
 Chen P. S.

↵
 Pastore J. M. ,
 Girouard S. D. ,
 Laurita K. R. ,
 Akar F. G. ,
 Rosenbaum D. S.
 ↵
 ↵
 ↵

↵
 Eisner D. A. ,
 Choi H. S. ,
 Diaz M. E. ,
 O’Neill S. C. ,
 Trafford A. W.

↵
 Diaz M. E. ,
 O’Neill S. C. ,
 Eisner D. A.

↵
 Kockskamper J. ,
 Blatter L. A.

↵
 Diaz M. E. ,
 Eisner D. A. ,
 O’Neill S. C.

↵
 Fox J. J. ,
 McHarg J. L. ,
 Gilmour R. F., Jr

↵
 Luo C. H. ,
 Rudy Y.

↵
 Turing A. M.
 ↵
 ↵
 ↵

↵
 Wier W. G. ,
 Egan T. M. ,
 LopezLopez J. R. ,
 Balke C. W.
 ↵

↵
 Walker M. L. ,
 Rosenbaum D. S.

↵
 Euler D. E.

↵
 Rubenstein D. S. ,
 Lipsius S. L.
 ↵
 ↵
 ↵
 ↵

↵
 Lengyel I. ,
 Epstein I. R.

↵
 Garfinkel A. ,
 Tintut Y. ,
 Petrasek D. ,
 Bostrom K. ,
 Demer L. L.

↵
 Kondo S. ,
 Asai R.