%Electrical Stability of Paced Cardiac Tissue Instability and Spatiotemporal Dynamics of Alternans in Paced Cardiac Tissue
Blas Echebarria and Alain Karma
We derive an equation that governs the spatiotemporal dynamics of small amplitude alternans in paced cardiac tissue. We show that a pattern-forming linear instability leads to the spontaneous formation of stationary or traveling waves whose nodes divide the tissue into regions %separate spatial domains with opposite phase of oscillation of action potential duration. This instability is important because it creates dynamically an heterogeneous electrical substrate for inducing fibrillation if the tissue size exceeds a fraction of the pattern wavelength. We compute this wavelength analytically as a function of three basic length scales characterizing dispersion and inter-cellular electrical coupling.
alternans.tex
\documentstyle[twocolumn,aps,epsfig,floats]{revtex}  

\begin{document} 

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}

\draft  

\twocolumn[\hsize\textwidth\columnwidth\hsize\csname %	
@twocolumnfalse\endcsname

\title{
%Electrical Stability of Paced Cardiac Tissue
Instability and Spatiotemporal Dynamics of Alternans in Paced Cardiac Tissue   
}
\author{Blas Echebarria and Alain Karma} 
\address{Department of Physics and Center for 
Interdisciplinary Research on 
Complex Systems,\\
Northeastern University, Boston, MA 02115}

\date{November 23, 2001}

\maketitle

\begin{abstract}
We derive an equation that 
governs the spatiotemporal dynamics of small amplitude
alternans in paced cardiac tissue.
We show that a pattern-forming linear instability leads to the spontaneous
formation of stationary or traveling waves  
whose nodes divide the tissue into regions 
%separate spatial domains 
with opposite phase
of oscillation of action potential duration.
This instability is important because it creates 
dynamically an heterogeneous electrical substrate 
for inducing fibrillation if the tissue size exceeds a fraction of the
pattern wavelength. We compute this wavelength analytically 
as a function of three basic length scales characterizing
dispersion and inter-cellular electrical coupling.   
\end{abstract}

 \pacs{PACS numbers: 87.19.Hh, 05.45.-a, 05.45.Gg, 89.75.-k}
]

It is well established that the duration of
cardiac excitation can oscillate
from beat to beat at sufficiently short pacing interval \cite{Min}.
Pioneering studies by Nolasco and Dahlen \cite{NolDal}
and Guevara {\it et al.} \cite{Gueetal} have demonstrated that
the generic sequence LSLS... of long and short action
potential duration ($APD$), known as alternans, 
is a direct consequence of the restitution relationship
\begin{equation}
APD^{n+1}=f(DI^n)\label{res}
\end{equation}
between the $APD$ generated by the $n^{th}+1$ 
stimulus, $APD^{n+1}$, and the diastolic time interval
$DI^n$ during which the
tissue recovers its resting properties after
the end of the previous ($n^{th}$) action potential.
If we denote the interval   
between the $n^{th}$ and $n^{th}+1$ stimulus by $T^n$,
we must have $DI^n= T^n-APD^n$. 
Then, for a fixed period: $T^n=\tau$ 
for all $n$, Eq. \ref{res}  
yields the map $APD^n=f(\tau-APD^{n-1})$ that has
a period doubling instability if the 
slope $f'$ of the restitution curve evaluated at its
fixed point exceeds unity, which generically occurs as
the period is decreased.  

Over the last decade, 
the study of alternans \cite{FraSim,Couetal,Kar94,Vin,Pasetal,Foxetal,Quetal,Watetal},
and their control \cite{control}, has become a main
focus of research because of the potentially   
crucial link of this dynamical
instability with cardiac fibrillation \cite{Kar}.  
However, there is presently no simple analytical
understanding of how the bifurcation to alternans is manifested
\emph{spatiotemporally} in paced cardiac tissue.
Analytical progress to date is limited
to the one-dimensional circulation of electrical impulse in
a ring of tissue \cite{Couetal,Kar94,Vin}.

In this letter, we derive an
equation that governs the 
spatiotemporal dynamics of alternans close to the
onset of instability. This enables us to 
obtain a quantitative analytical understanding 
of the formation of recently observed
complex patterns of $APD$ oscillations that can promote
fibrillation \cite{Pasetal,Foxetal,Quetal,Watetal}.
A crucial feature of these patterns is that  
the $APD$ oscillates with opposite phases in two (or more) spatially extended 
regions of tissue, i.e. with a sequence LSLS...
in one region and SLSL... in the other.
These ``discordant alternans'' have been observed experimentally
in both two-dimensional \cite{Pasetal}
and linear strands \cite{Foxetal} of cardiac tissue,  
as well as in ionic model simulations \cite{Foxetal,Quetal,Watetal}. 
Moreover, they have been shown to lead to the formation
of conduction blocks \cite{Foxetal} as well as to the onset of spiral
wave formation and fibrillation \cite{Pasetal}. 
We show here that discordant alternans 
result from a \emph{pattern-forming} linear instability 
that has interesting similarities with classic  
instabilities leading to the spontaneous formation of
spatially periodic patterns in nature 
(such as Rayleigh-B\'enard convection, Taylor-Couette flow,
etc \cite{CroHoh}), but also presents some unique features.  

We consider  
a one-dimensional (1-d) homogeneous cable 
of length $L$ paced at period $\tau$ from one end ($x=0$). 
Close to the onset of instability, we can expand
the $APD$ and the period in the form
\begin{eqnarray}
APD^n(x)&\approx &APD_c\,+\,a(x,t)\,e^{i\pi n}, \label{aa}\\
T^n(x)&\approx &\tau_c\,-\delta \tau\,+\,b(x,t)\,e^{i\pi n},\label{bb}
\end{eqnarray}
where $APD_c$ and $\tau_c$ are the $APD$ and the period  
evaluated at the bifurcation point of the map 
($f'=1$), $x$ measures the position along the cable, and
$\delta\tau\equiv \tau_c-\tau\ll \tau_c$. 
In this range of period, $a$ and $b$
vary slowly from beat to beat, which allows us  
to treat  the time, $t\equiv n\tau$, as a continuous variable;
the fast beat-to-beat oscillations
are contained in the exponential factor 
$e^{i\pi n}$.  

A relation between $a$ and $b$ can first be
derived by noting that $T^n(x)$ is
the difference of arrival time of two subsequent
action potentials at $x$ \cite{Watetal}, or
\be
T^n(x)= \tau + \int_0^x \frac{d\,x'}{c(DI^n(x'))} -
\int_0^x \frac{d\,x'}{c(DI^{n-1}(x'))}, \label{disp}
\ee
where the first (second) integral on the right-hand-side  
is the time required for the leading front of the
action potential to
travel from the paced end to $x$ at the
$n^{th}$ ($n^{th}-1$) stimulus; concomitantly,
$c(DI)$ is the standard
dispersion curve that relates the propagation
speed of this front with the
local diastolic interval.
The dispersion curve is typically steeply
increasing at small $DI$ and flat at large $DI$. 
Substituting Eqs.
\ref{aa}-\ref{bb} in Eq. \ref{disp} with 
$DI^n(x)=T^n(x)-APD^n(x)$, and expanding to linear
order in $a$ and $b$, we obtain 
$b(x)\simeq \int_0^x a(x') dx'/\Lambda$ 
where we have defined $\Lambda\equiv  c^2/(2\,c')$, 
with $c$ and $c'\equiv dc/dDI$ evaluated 
at the bifurcation, and
assumed that $\Lambda$ is much larger
than the scale over which $a$ varies.

Next, in order to derive an evolution equation for 
the amplitude $a(x,t)$, we first neglect
the influence of the electrical coupling between cells on the
$APD$. This allows us to assume that the restitution 
relationship (\ref{res}), and hence the second iteration of the map
\begin{equation}
APD^{n+2}=f(T^{n+1}-f(T^n-APD^n)),\label{map2}
\end{equation}
continues to be valid even when
the $APD$ is non-spatially uniform; we shall soon see why it is crucial
to relax this assumption. 
We substitute Eqs. \ref{aa}-\ref{bb} in Eq. \ref{map2},
and expand the right-hand-side keeping only the dominant
linear and weakly nonlinear terms. Furthermore, we use the
aforementioned fact that, close to onset, $a$ varies slowly from beat 
to beat, and we therefore expand the left-hand-side as
$APD^{n+2}\simeq APD^n + 2\, \partial a/\partial n \, e^{i\pi n}$,
where $\partial a/\partial n=\tau\partial a/\partial t$.
Finally, we use the integral relation between  
$a$ and $b$ derived earlier. After equating both sides
of Eq. \ref{map2}, we obtain
\be
\tau \partial_t a =
\sigma a -g a^3 - \int_0^x \frac{dx'}{\Lambda}\,a(x')  ,\label{aeq}\\
\ee
where $\sigma\equiv f''(\tau - \tau_c)/2$, $g\equiv f''^2/4-f'''/6$, 
and all derivatives are evaluated at the bifurcation point.

In order to test this evolution equation,
we simulate the standard cable equation  
\be
\partial_t V = D\, \partial_x^2 V 
-\left(I_{\rm ion}+I_{ext}\right)/C_m,
\label{cable}
\ee
for the membrane current
$I_{\rm ion}$ given by the Noble model \cite{Nob} with
time in units of millisecond (ms),
$D=2.5 \times 10^{-4}$ cm$^2$/ms, $C_m=12 \,\mu$F/cm$^2$,
dx=0.01  cm, dt=0.05 ms, 
and $I_{ext}$  
modeling a sequence of stimuli applied   
at $x=0$ at the pacing interval $\tau$.  
To determine the parameters of the amplitude equation,
we compute the restitution and dispersion curves
by pacing Eq. \ref{cable} in a short cable and by using
two subsequent stimuli spaced by different
intervals to vary $DI$; we also use $V=-40$ mV as
threshold of the transmembrane voltage to define the $APD$.
We impose zero gradient boundary conditions on $V$ and
$a$ at the two ends of the cable in all 
our simulations.

Fig. \ref{comp}(a) shows that Eq. \ref{aeq}
produces discordant alternans
consistent with the picture that  
restitution and dispersion suffice to produce this state
\cite{Foxetal,Quetal,Watetal}. 
However, the magnitude
of the spatial gradient of $a$ at the node  
increases with time, and can be shown to diverge in a finite time.  
This yields an unphysical spatial discontinuity of $APD$,
which also occurs if the system of coupled maps 
(Eqs. \ref{res} and \ref{disp}), from which the amplitude equation
is derived, is solved numerically. Note that the 
analogous system of coupled maps for 
a 1-d pulse circulating in a ring produces a smooth $APD$
modulation starting from a smooth initial condition \cite{Couetal},
which highlights the much stronger role of dispersion
at producing heterogeneity of $APD$ during pacing than circulation. 

\begin{figure}[t]
\centerline{
\epsfxsize=5.5cm\epsfbox{fig1.eps}
}
\vspace{0.4 cm}
\caption{Amplitude $a$ of $APD$ oscillation vs $x$ for  
Noble parameters and $\tau=258$ ms.
(a): profiles obtained with Eq. \protect\ref{aeq}  
at different times (dotted lines) and
final stationary profile (solid line). (b): stationary
profiles obtained with Eq. \protect\ref{cable} (solid line)
and Eq. \protect\ref{aeqg} (dashed line). Nodes ($a=0$) separate
tissue regions with $\pi$ out of phase oscillations.
}
\label{comp}
\end{figure}

This discontinuity, which is absent in the simulation
of the cable equation in Fig. \ref{comp}(b), can be cured
by adding spatial derivative terms 
to the amplitude equation. Since the 
cable is paced at one end, the underlying basic state (i.e. traveling pulses)
is not invariant under parity symmetry.
Hence, in addition to $\partial_x^2 a$,   
a term proportional to $\partial_x a$ must  
generally be included, which yields the final
form of our amplitude equation
\be 
 \tau \partial_t a = 
\sigma a -g a^3  
-\int_0^x \frac{dx'}{\Lambda}\,a(x')
- w \, \partial_x a+\xi^2 \partial^2_x a.\label{aeqg}
\label{aeqr}
\ee
To calculate the new lengthscales $w$ and $\xi$, we must
determine how the electrical coupling between cells
modifies the restitution
relationship (Eq. \ref{res}). For complex electrophysiological models
such as Noble, this generally needs to be done numerically using 
a procedure that will be discussed elsewhere.
For a simple two-variable
ionic model, which we study below, the analytical expressions
\begin{eqnarray}
w&=&2D/c,\label{w2v} \\
\xi &=&(D\times  APD_c)^{1/2}, \label{xi2v}
\end{eqnarray}
can be derived by interpreting Eq. \ref{cable} as a diffusion equation
with a source $I_{ion}$, and expressing
$V$ as space-time integral of $GI_{ion}$,
where $G$ is the standard Green's function of the diffusion equation.
This integral can then be evaluated analytically because the action potential  
shape is simply triangular for this model, and
used to calculate $w$ and $\xi$. Eq. \ref{xi2v} has 
the simple physical interpretation that $V$
diffuses a length $\sim \xi$ in the time interval of one 
$APD$. Therefore, the repolarization of a given cell
is influenced by other cells within a length $\sim \xi$ of cable.
In addition, repolarization of this cell
is influenced unequally by its left and right
neighboring cells because these cells are activated at different
times by the propagating wavefront. Clearly, this asymmetry must 
vanish in the limit $c\rightarrow \infty$ where all
cells are activated simultaneously consistent with Eq. \ref{w2v}.
Fig. \ref{comp}(b) 
shows that our regularized amplitude
equation (\ref{aeqg}) now produces a smoothly varying  
stationary profile of $a$ that agrees well with the 
simulation of the cable-Noble equation, where $a$  
is obtained from the $APD$ using Eq. \ref{aa}.



\begin{table}
\caption{Values of various lengths in cm with
$\lambda_{theor}$ (Eq. \protect\ref{lsta} 
or Eq. \protect\ref{ltra}) and $\lambda_{sim}$ (from simulations
of Eq. \protect\ref{cable}).}
\begin{tabular}{l||llllll}
 Model & $\Lambda$ & w & $\xi$ & $\lambda_{theor}/4$ & $\lambda_{sim}/4$ & $L_{min}$ 
 \\  \hline \
Noble &  49.1  & 0.045  & 0.18  & 2.33 & 2.6 & 2.75 \\
Two-variable &  3.55  &  0.031  &  0.235 & 1.33 & 1.1 &  1.15  
\label{lengths}
\end{tabular}
\end{table}

The genesis of discordant alternans can be understood  
by computing the linear stability spectrum
of the spatially homogeneous
state ($a=0$). We have calculated this 
spectrum both numerically for different $L$, 
and analytically for the large $L$ limit. The main result
is that the wave pattern can emerge from the amplification of
\emph{either} a unique finite wavelength 
mode, which yields a stationary pattern,  
\emph{or} from a discrete set of  
complex modes that approach a continuum in the 
limit $L\rightarrow \infty$, and yields a traveling pattern. 
There is indeed experimental evidence for both stationary
\cite{Pasetal} and traveling \cite{Foxetal} waves.

We treat here explicitly the large $L$ limit
since it provides the basis to understand  
finite-$L$ patterns. In this limit, we can analyze
stability by differentiating Eq. \ref{aeqg} with respect
to $x$ and letting $a(x,t) \sim e^{ikx+\Omega t}$,
with both $\Omega\equiv \Omega_r+i\Omega_i$  
and $k\equiv k_r+ik_i$ complex,
which yields at once the eigenvalue equation
\begin{equation}
\Omega=\sigma-\xi^2k^2-i\left[wk-1/(\Lambda k) \right].\label{spec}
\end{equation}
When dispersion is weak, 
a unique real mode with $k_i=0$ 
grows faster than the other complex modes 
and yields a stationary pattern. Its wavelength $\lambda=2\pi/k_r$,
is determined by the condition $\Omega_i=0$,
which yields  
\be 
\lambda = 2\pi (w\Lambda)^{1/2},\label{lsta}
\ee
in good agreement with the wavelength observed
in simulations of the cable-Noble equation
(Table \ref{lengths}). Note that $\cos k_rx$ is
an exact eigenvector of Eq. \ref{aeqg} linearized around
$a=0$ that satisfies  
$\partial_x a=0$ at the two cable ends when $L$ is
an integer multiple of $\lambda/2$.  
The threshold of instability occurs
when $\Omega_r=0$, or for a period $\tau_{th}$ 
defined by $\sigma_{th}=f''(\tau_{th}-\tau_c)/2=\xi^2/(w\Lambda)$.

In the opposite limit where dispersion is strong,  
complex modes that grow exponentially at 
large $x$ are the most unstable. It is simple to
deduce from Eq. \ref{spec} that each $k$-mode
travels towards the pacing end of the cable, but
a wave-packet constructed from a linear superposition of
these modes has a group velocity that makes the packet
move away from the pacing end. This is the signature of a 
convective instability \cite{CroHoh} where perturbations 
are transported as they grow, similarly to
Taylor-Couette vortices developing in an axial flow \cite{Babetal}.
In such a situation, patterns are only transient
unless the group velocity vanishes, 
or $\partial \Omega_i/\partial k_r=0$, and hence
they grow at a fixed point in space.
Moreover, the fastest growing wavelength that dominates at large time must
correspond to a maximum of $\Omega_r$, which yields the additional
condition $\partial \Omega_r/\partial k_r=0$. These two conditions
are equivalent to requiring that $d\Omega/dk=0$. From this, 
we deduce that the threshold of absolute instability occurs
when $\sigma_{th}=(3/2)(\xi/2\Lambda)^{2/3}$, with a pattern of wavelength
\begin{equation}
\lambda=(4\pi/\sqrt{3})(2\xi^2 \Lambda)^{1/3},\label{ltra}
\end{equation}
which travels with phase velocity $\Omega_i\lambda/(2\pi)$
where $\Omega_i=(3\sqrt{3}/2)(\xi/2\Lambda)^{2/3}$. 
It can also be deduced that
traveling waves are favored over
a stationary pattern when $\Lambda=2c^2/c' \ll \xi^4/w^3$,
and hence for large enough dispersion ($c'$),
and vice versa in the opposite limit.

Nonlinearities are mainly 
responsible for saturating the amplitude of the linear
modes, but both the wavelength and evolution of the 
pattern can generally vary with distance from onset.
There is a close analogy between the amplitude equation
(\ref{aeqg}) and the real Ginzburg-Landau equation that has been
extensively studied in the context of phase transitions and
front propagation \cite{CroHoh}. The dynamics is richer here because the 
integral term originating from dispersion 
causes a non-local interaction of
the fronts separating two out-of-phase oscillating
regions with the pacing end of the cable.

\begin{figure}[t]
\centerline{
\epsfxsize=6.6cm\epsfbox{fig2.eps}
}
%\vspace{0.1 cm}
\caption{Space-time plots of $a$  
obtained by simulations of Eq. \protect\ref{aeqg} for
parameters of the two-variable model, showing   
absolutely unstable (left and $\tau=295$ ms) and  
convectively unstable (right and $\tau=298$ ms ) wave patterns. The crosses denote
the positions of the nodes ($a=0$).
}
\label{examples}
\end{figure}


As a final test of our theory, we study
a two-variable model where all the parameters
of the amplitude equation can be calculated analytically
including $w$ and $\xi$ given by Eqs. \ref{w2v}-\ref{xi2v}, and  
for which our theory predicts traveling waves.
In this model, $I_{\rm ion}/C_m$ is the
sum of a slow outward current, $I_o/C_m=\tau_0^{-1}(S+(1-S)V/V_c)$,
and a fast inward current, $I_i/C_m=-\tau_a^{-1} h \,S$, where inactivation 
of the latter is controlled by  
\be
\partial_t h=(1-S-h)/(\tau^-(1-S)+S\tau^+)
\ee
%\begin{eqnarray}  
%& &I_{ion}/C_m=\tau_a^{-1} h \,S-\tau_0^{-1}(S+(1-S)V/V_c),\\
%& &\partial_t h=(1-S-h)/(\tau^-(1-S)+S\tau^+),.
%\end{eqnarray}
In addition, $V$ is dimensionless
and $S\equiv (1+\tanh((V-V_c)/\epsilon))/2$. We choose
$V_c=0.1$, $\tau_0=150$ ms, $\tau_a=6$ ms, $\tau^-=60$ ms, 
$\tau^+=12$ ms, $\epsilon=0.005$, and we simulate Eq. \ref{cable} with
the same $D$ as before, dx=0.01 cm, dt=0.02 ms, and $V=0.1$ as threshold  
for the $APD$.  

\begin{figure}[t]
\centerline{
\epsfxsize=9cm\epsfbox{fig3.eps}
}
\vspace{0.2 cm}
\caption{Stability diagram of
two-variable cable model with domains of no-alternans
(open circles), concordant alternans
(filled squares), and discordant alternans
(filled diamonds); conduction blocks form
at smaller $\tau$ not shown here. Boundaries between the 
same domains obtained by simulations
of the amplitude equation (\protect\ref{aeqg})
are shown by solid lines. The dashed line denotes the bifurcation period
for alternans predicted by the map of Eq. \ref{res}.
On the right panel we show profiles 
of $a$ vs $x$ at times $t$ (solid), $t+2\tau$ (dashed) and $t+4\tau$ 
(dotted lines), for three different cable lengths.}
\label{phasediag}
\end{figure}

As predicted, we observe traveling waves in this model
with a wavelength that agrees well with Eq. \ref{ltra}. 
Fig. \ref{examples} illustrates
sustained and transient wave patterns above (left) and
below (right) the onset of absolute stability, respectively. 
Furthermore, Fig. \ref{phasediag} shows that
the results of simulations of the cable and amplitude
equations are in good agreement over a wide
range of $L$ and $\tau$, including the
fact that the onset of instability occurs at a shorter
period in a larger cable (and, hence, for a slope of the
restitution curve larger than unity).
Fig. \ref{phasediag} is qualitatively similar
for stationary waves, but this 
overstabilization is smaller because dispersion is weaker.  

The $APD$ oscillations  
at a fixed $x$ are periodic when the
wave is stationary and quasiperiodic when it travels. 
For both cases, we find here that 
$\lambda$ is independent of $L$, in contrast to
the oscillations produced by a pulse circulating
in a ring where it is known that
$\lambda \simeq 2L/(2i+1)$ for weak dispersion  
\cite{Couetal} (with $i$ integer and $L=$ ring perimeter). 
As will be discussed elsewhere,
our theory applied to the ring shows that  
the bifurcation to alternans is finite dimensional
with $i=0$ being the most unstable mode, in agreement with
the fact that it is this mode that is generically
selected in experiments \cite{FraSim} and
ionic model simulations \cite{Couetal,Kar94,Vin}.  
In addition, it
shows that the 
gradient term ($-w\partial_x a$) can lead to
quasiperiodicity even in the absence 
of dispersion ($c'=0$).

Our results demonstrate that the formation
of discordant alternans is crucially affected by the 
effect of electrical coupling (diffusion) on repolarization,
in addition to restitution and dispersion. 
Dispersion is responsible for the formation of nodes and
spatial gradients of $APD$ that steepen with time. 
Diffusion, in turn, tends to spread the $APD$ spatially,
and also induces a drift of the pattern away from the pacing site
that is induced by the more subtle gradient term 
($-w\partial_x a$) in the amplitude equation. When
dispersion is sufficiently weak, drift balances dispersion
and produce a stationary pattern. In the opposite limit,
the tendency for dispersion to form steep gradients
of $APD$ is balanced by the spreading effect of diffusion. 
Nodes then travel, cyclically 
disappearing (appearing) at the pacing (opposite)
end of the cable. 

In conclusion, we have derived a simple  
evolution equation that describes the universal spatiotemporal
dynamics of small amplitude alternans in paced cardiac tissue. 
Moreover, we have shown that discordant 
wave patterns \cite{Pasetal,Foxetal,Quetal,Watetal}, 
which are linked to fibrillation \cite{Pasetal}, 
result from a finite wavelength
linear instability. Hence, their formation
requires a minimum tissue size 
$L_{min}\sim \lambda/4$, required for at least
one node to form. The value of $L_{min}$ that we
measure in simulations of reaction-diffusion
models are actually close to $\lambda/4$
with $\lambda$ predicted by Eq. \ref{lsta}
and Eq. \ref{ltra}, respectively 
(Table \ref{lengths}).  
%Several predictions
%made here should be experimentally testable given
%the expected universal character of Eq. \ref{aeqg}
%close to the onset of instability.
This equation can be extended
to higher dimensions as well as to include 
heterogeneities. Our preliminary results in a 2-d
homogeneous tissue are very similar to 1-d.
Finally, this equation is readily applicable to model  
a non-constant pacing interval and provides
a theoretical basis to study the control of 
alternans in spatially extended tissue.

This research was supported in part by NIH 
SCOR in Sudden Cardiac Death P50-HL52319.  
 
\vskip -0.5 cm

\begin{references}

%\bibitem{Wi98}
%A.T. Winfree, Chaos {\bf 8}, 1 (1998). Focus issue on fibrillation in normal 
%ventricular myocardium [Chaos {\bf 8} (1998)].

\bibitem{Min} G. R. Mines, J. Physiol. (Lond.) {\bf 46}, 349 (1913).

\bibitem{NolDal} J. B. Nolasco and R. W. Dahlen, J. App. Physiol, {\bf 25},
191 (1968).

\bibitem{Gueetal}
M. R. Guevara {\it et al.},
%, G. Ward, A. Shrier, and L. Glass
in {\it Computers in Cardiology}, IEEE Comp. Soc., pp. 167 (1984). 

\bibitem{FraSim}
L. H. Frame and M. B. Simson, Circulation {\bf 78},     
1277 (1988).                                                         

\bibitem{Couetal}
M. Courtemanche, L. Glass and J. P. Keener, 
Phys. Rev. Lett. {\bf 70}, 2182 (1993).

\bibitem{Kar94} A. Karma, Chaos {\bf 4}, 461 (1994).

%\bibitem{Gla} L. Glass, Phys. Today {\bf 49}, 40 (1996).

\bibitem{Vin}
A. Vinet, Annals Biomed. Eng. {\bf 28}, 704 (2000).

\bibitem{Pasetal}
J. M. Pastore {\it et al.},
%, S.~D.~Girouard, K.~R.~Laurita, F.~G.~Akar and D.~S.~Rosenbaum, 
Circulation {\bf 99}, 1385 (1999).

\bibitem{Foxetal}
J. J. Fox {\it et al.},
%,M. L. Riccio, F. Hua, E. Bodenschatz, and R. F. Gilmour,
Spatiotemporal Transition to Conduction Block in Canine Ventricle,
Circ. Res., in press.

\bibitem{Quetal}
Z. Qu {\it et al.}, 
%, A.~Garfinkel, P-S.~Chen, J.~W.~Weiss, 
Circ. Res. {\bf 102},
1664 (2000).

\bibitem{Watetal}
M. A. Watanabe {\it et al.},
%, F.~H.~Fenton, S.~J.~Evans, H.~M.~Hastings and A.~Karma, 
J. Cardiovasc. Electrophys. {\bf 12}, 196 (2001).  

\bibitem{control} K. Hall {\it et al.} 
Phys. Rev. Lett. {\bf 78}, 4518 (1997); D. J. Gauthier and J. E. S.
Scoolar, Phys. Rev. Lett. {\bf 79}, 4938 (1997); 
D. J. Christini {\it et al.}, Proc. Nat. Acad. Sci. {\bf 98}, 5827 (2001).

%\bibitem{RoJa94}
%D.~S.~Rosenbaum, L.~E.~Jackson, J.~M.~Smith, H.~ Garan, J.~N.~Ruskin, 
%R.~J.~Cohen, N. Engl. J. Med. {\bf 330}, 235 (1994);
%N.~A.~M.~Estes, G.~Micaud, D.~P.~Zipes, N.~El-Sherif, F.~J.~Venditti, 
%D.~S.~Rosenbaum, P.~Albrecht, P.~J.~Wang, R.~J.~Cohen, Am. J. Cardiol. 
%{\bf 80}, 1314 (1997).

\bibitem{Kar} See: A. Karma, Proc. Natl. Acad. Sci 
USA {\bf 97}, 5687 (2000); D. S. Rosenbaum, J. Cardiovasc. Electrophysiol.
{\bf 12}, 207 (2001); and earlier references therein.

\bibitem{CroHoh} M. C. Cross and P. C. Hohenberg,                       
Rev. Mod. Phys. {\bf 65}, 851 (1993). 

\bibitem{Nob}
D. Noble, J. Physiol {\bf 160}, 317 (1962).

%\bibitem{ToPr98}
%S.M. Tobias, M.R.E. Proctor, and E. Knobloch, Physics D, {\bf 113}, 43 
%(1998). 

%\bibitem{LaGi96}
%K. R. Laurita {\it et al.},
%S.~D.~Girouard, D.~S.~Rosenbaum, 
%Circ Res. {\bf 79},
%493 (1996).

\bibitem{Babetal} K. L. Babcock, G. Ahlers, and D. S. Cannell,
Phys. Rev. Lett. {\bf 67}, 3388 (1991).
    
%\bibitem{Fib} 
%M. L. Riccio, M. L. Koller,
%and R. F. Gilmour, Circ. Res. {\bf 84}, 955 (1999);
%A. Garfinkel {\it et al.}, Proc. Natl. Acad. Sci 
%USA {\bf 97}, 6061 (2000).

\end{references}



\end{document}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



 
A generalization of our results to a two-dimensional tissue paced at one 
point is straightforward, provided the wavefronts originated at the pacing
point remain circular. This will be the case as long as the dimensions of 
the system are large compared with the width of the pulse. In this case,
Eq. \ref{eq.pde.diff} can be generalized to:
\be
T_o \frac{\partial a}{\partial t}=\sigma a -\frac{1}{\Lambda} \in\tau^r
a dr' -ga^3 - w\partial_r a + \xi^2 \nabla^2 a,
\label{eq.2D}
\ee
and the motion of the node becomes nontrivial, 
because of curvature effects. For the Noble model, for instance, the 
node is stationary (Fig. \ref{fig.2Dnodes}a) or 
traveling (Fig. \ref{fig.2Dnodes}b), depending on 
the size of the system.
The minimum length along the diagonal for the node to appear is the same
as in the one-dimensional case, but now the drift term has a contribution 
coming from the curvature of the interface, that can oppose or favor its 
motion. When the system is large enough so that the node would be formed 
in the region with positive curvature (close to 
the pacing point), this extra effect is enough to make it travel.   

\begin{figure}[t]
\centerline{
\epsfxsize=7.5cm\epsfbox{figuras/2Dnodes.eps}
}
\caption{Simulations of Eq. \ref{eq.2D} for the parameters of
the Noble model, with  $\tau=258$ and a) $L=5$cm, b) $L=7$cm. The lines 
in b) are drawn every ten beats.
\label{fig.2Dnodes}}
\end{figure}

where $\sigma\equiv f_{(2)}(\tau - T_c)/2$ 
measures the distance from the onset of the  
bifurcation, $g\equiv f_{(2)}^2/4-f_{(3)}/6$ determines
whether the bifurcation is supercritical ($g>0$) 
or subcritical ($g<0$), and $f_{(i)}$ denotes
the $i^{th}$ derivative of $f$ at $D=D_c$. If $g<0$,
the quintic term $\chi a^5$, with 
$\chi=-(f_{(5)}-15f_{(2)}f_{(4)}+5f_{(2)}^2f_{(3)}+10f_{(3)}^2)/240$,
also needs to be included for $a$ to stay finite. 

The dispersion curve is typically steeply
increasing at small $DI^n$ and flat at large $DI^n$. 
Substituting Eqs.
\ref{aa}-\ref{bb} in Eq. \ref{disp} with 
$DI^n(x)$ defined by Eq. \ref{Ddef}, and expanding to linear
order in $a$ and $b$, we obtain 
$b(x)=\in\tau^x (a(x')-b(x')) dx'/\Lambda$, or equivalently
\be
b(x)=e^{-x/\Lambda} \in\tau^x \frac{dx'}{\Lambda}\,e^{x'/\Lambda }a(x') ,
\label{beq}
\ee
where we have defined  
$\Lambda\equiv  c^2/(2\,c')$ 
with both $c$ and $c'\equiv dc/dDI$ evaluated at $DI=DI_c$.

For the latter, we obtain
$\sigma\simeq  8.85 \times 10^{-3}(T_c - \tau)$ 
and $g=2.34\times 10^{-5}$; the bifurcation is supercritical
($g>0$) and we can therefore set $\chi=0$. For Noble,
we obtain  
$\sigma\simeq 5.7\times 10^{-3}(T_c-\tau)$; the bifurcation
is weakly subcritical with $g\simeq -8\times 10^{-6}$ and 
$\chi \simeq 1.37\times 10^{-8}$.  

This equation provides a simple
picture of the formation of discordant alternans.
For simplicity, assume that $a$ 
is initially uniform along the cable, $a(x,0)=a_0$, where
$a_0>0$ is chosen to make the right-hand-side of Eq. \ref{aeq} vanish
when $b=0$. Solving Eq. \ref{beq} for this case, we easily obtain that
$b(x,0)=a_0(1-e^{-x/\Lambda})$, which shows that dispersion 
sets up a spatial gradient of period oscillation with
$b$ increasing along the cable. Furthermore, Eq. \ref{aeq}
implies that, initially, $\partial_t a=-b<0$, such that this gradient
will cause $a$ to decrease fastest at the end opposite to
the pacing site. If $a$ crosses zero at this end, a
{\it node} will form separating two regions of the cable 
oscillating with opposite phases, with $a>0$ ($a<0$) on the 
side of the node proximal (distal) to the pacing site. 

The numerical integration of the system of  
spatially coupled maps,
 $APD^{n+1}(x)=f(DI^n(x))$ and
Eqs. \ref{disp}-\ref{Ddef}, yields the same discontinuity, which is therefore
not a consequence of the approximations
made in deriving the amplitude equation. Interestingly,
the integration of the analogous coupled maps for the non-paced
case where the action potential circulates in a ring of tissue 
is known not to produce  
spatial discontinuities of $APD$ from a smooth initial condition; dispersion
only produces a slow spatial modulation of $APD$ with a wavelength
that is about twice the ring length. Therefore, the
plots of Fig. \ref{comp}(a) show that dispersion is much
more effective at creating strong spatial heterogeneities of $APD$  
during pacing then during quasi-one-dimensional 
circulation. 

The same discontinuity
 profile obtained from the simulations
of the cable equation is smooth (Fig. \ref{comp}(b)),
consistent with the fact that the diffusive electrical coupling 
term ($D\partial_x^2V$)
prevents $V$, and hence the $APD$, to become discontinuous
in space.

For higher-order ionic models, 
the form of $I_{ion}$ during repolarization is 
considerably more
complicated such that $\xi$ and $w$ 
generally need to be 
calculated numerically. These legnthscales, however,
have similar physical origins as in the two-variable
model, and hence comparable magnitude as shown in 
Table \ref{lengths}. Furthermore, it is important  
to emphasize that for typical physiological values 
of $D$, $A_c$, $c$,
and of the dispersion slope $c'$ at onset of
alternans, the three key lengthscales
$w$, $\xi$, and $\Lambda$ generally span several
decades with $w\ll \xi \ll \Lambda$.

%Furthermore, we discuss the implications of our finding
%for cardiac arrhythmogenesis and experiments.
%It produces a spatially 
%periodic variation of the amplitude of
%alternans, which can be stationary or travel,
%and the nodes of this wave pattern separate regions of tissue where the APD
%oscillate out of phase.
%This implies that the discordant alternans substrate that  
%has been shown to produce fibrillation under rapid pacing will
%spontaneously form if the tissue exceeds some minimum size 
%$L_{min}\sim \lambda/4$, required for at least
%one node to form. In this respect, an essential outcome of
%our study is the identification of basic lengthscales that
%determine $\lambda$, and that depend on
%restitution of conduction velocity and on how
%the electrical coupling between cells influence the APD.
%Another interesting result is that the onset
%of this instability depends on tissue size,
%and occurs counter-intuitively at a smaller period 
%in a larger tissue, a prediction that 
%is validated by our ionic model simulations.

%We will show that it is crucial to keep this term in order to
%obtain the right lengthscale of the system, unlike previous work
%where only the symmetrical part was considered \cite{Vin}. 

cond-mat0111552.extracted_bib
\bibitem{Wi98}
%A.T. Winfree, Chaos {\bf 8}, 1 (1998). Focus issue on fibrillation in normal 
%ventricular myocardium [Chaos {\bf 8} (1998)].


\bibitem{Min} G. R. Mines, J. Physiol. (Lond.) {\bf 46}, 349 (1913).


\bibitem{NolDal} J. B. Nolasco and R. W. Dahlen, J. App. Physiol, {\bf 25},
191 (1968).


\bibitem{Gueetal}
M. R. Guevara {\it et al.},
%, G. Ward, A. Shrier, and L. Glass
in {\it Computers in Cardiology}, IEEE Comp. Soc., pp. 167 (1984). 


\bibitem{FraSim}
L. H. Frame and M. B. Simson, Circulation {\bf 78},     
1277 (1988).                                                         


\bibitem{Couetal}
M. Courtemanche, L. Glass and J. P. Keener, 
Phys. Rev. Lett. {\bf 70}, 2182 (1993).


\bibitem{Kar94} A. Karma, Chaos {\bf 4}, 461 (1994).

%
\bibitem{Gla} L. Glass, Phys. Today {\bf 49}, 40 (1996).


\bibitem{Vin}
A. Vinet, Annals Biomed. Eng. {\bf 28}, 704 (2000).


\bibitem{Pasetal}
J. M. Pastore {\it et al.},
%, S.~D.~Girouard, K.~R.~Laurita, F.~G.~Akar and D.~S.~Rosenbaum, 
Circulation {\bf 99}, 1385 (1999).


\bibitem{Foxetal}
J. J. Fox {\it et al.},
%,M. L. Riccio, F. Hua, E. Bodenschatz, and R. F. Gilmour,
Spatiotemporal Transition to Conduction Block in Canine Ventricle,
Circ. Res., in press.


\bibitem{Quetal}
Z. Qu {\it et al.}, 
%, A.~Garfinkel, P-S.~Chen, J.~W.~Weiss, 
Circ. Res. {\bf 102},
1664 (2000).


\bibitem{Watetal}
M. A. Watanabe {\it et al.},
%, F.~H.~Fenton, S.~J.~Evans, H.~M.~Hastings and A.~Karma, 
J. Cardiovasc. Electrophys. {\bf 12}, 196 (2001).  


\bibitem{control} K. Hall {\it et al.} 
Phys. Rev. Lett. {\bf 78}, 4518 (1997); D. J. Gauthier and J. E. S.
Scoolar, Phys. Rev. Lett. {\bf 79}, 4938 (1997); 
D. J. Christini {\it et al.}, Proc. Nat. Acad. Sci. {\bf 98}, 5827 (2001).

%
\bibitem{RoJa94}
%D.~S.~Rosenbaum, L.~E.~Jackson, J.~M.~Smith, H.~ Garan, J.~N.~Ruskin, 
%R.~J.~Cohen, N. Engl. J. Med. {\bf 330}, 235 (1994);
%N.~A.~M.~Estes, G.~Micaud, D.~P.~Zipes, N.~El-Sherif, F.~J.~Venditti, 
%D.~S.~Rosenbaum, P.~Albrecht, P.~J.~Wang, R.~J.~Cohen, Am. J. Cardiol. 
%{\bf 80}, 1314 (1997).


\bibitem{Kar} See: A. Karma, Proc. Natl. Acad. Sci 
USA {\bf 97}, 5687 (2000); D. S. Rosenbaum, J. Cardiovasc. Electrophysiol.
{\bf 12}, 207 (2001); and earlier references therein.


\bibitem{CroHoh} M. C. Cross and P. C. Hohenberg,                       
Rev. Mod. Phys. {\bf 65}, 851 (1993). 


\bibitem{Nob}
D. Noble, J. Physiol {\bf 160}, 317 (1962).

%
\bibitem{ToPr98}
%S.M. Tobias, M.R.E. Proctor, and E. Knobloch, Physics D, {\bf 113}, 43 
%(1998). 

%
\bibitem{LaGi96}
%K. R. Laurita {\it et al.},
%S.~D.~Girouard, D.~S.~Rosenbaum, 
%Circ Res. {\bf 79},
%493 (1996).


\bibitem{Babetal} K. L. Babcock, G. Ahlers, and D. S. Cannell,
Phys. Rev. Lett. {\bf 67}, 3388 (1991).
    
%
\bibitem{Fib} 
%M. L. Riccio, M. L. Koller,
%and R. F. Gilmour, Circ. Res. {\bf 84}, 955 (1999);
%A. Garfinkel {\it et al.}, Proc. Natl. Acad. Sci 
%USA {\bf 97}, 6061 (2000).

